Statistical mechanics of integrable quantum spin systems

This script is based on the notes the author prepared to give a set of six lectures at the Les Houches School"Integrability in Atomic and Condensed Matter Physics"in the summer of 2018. The school had its focus on the application of integrability based methods to problems in non-equilibrium statistical mechanics. The lectures were meant to complement this subject with background material on the equilibrium statistical mechanics of quantum spin chains from a vertex model perspective. The author was asked to provide a minimal introduction to quantum spin systems including notions like the reduced density matrix and correlation functions of local observables. He was further asked to explain the graphical language of vertex models and to introduce the concepts of the Trotter decomposition and the quantum transfer matrix. This was basically the contents of the first four lectures presented at the school. In the remaining two lectures these notions were filled with life by deriving an integral representation of the free energy per lattice cite for the Heisenberg-Ising chain (alias XXZ model) using techniques based on non-linear integral equations.


Preface
This script is based on the notes the author prepared to give a set of six lectures at the Les Houches School "Integrability in Atomic and Condensed Matter Physics" in the summer of 2018. The responsibility for the selection of the material is partially with the organisers, Jean-Sébastien Caux, Nikolai Kitanine, Andreas Klümper and Robert Konik. The school had its focus on the application of integrability based methods to problems in non-equilibrium statistical mechanics. My lectures were meant to complement this subject with background material on the equilibrium statistical mechanics of quantum spin chains from a vertex model perspective. I was asked to provide a minimal introduction to quantum spin systems including notions like the reduced density matrix and correlation functions of local observables. I was further asked to explain the graphical language of vertex models and to introduce the concepts of the Trotter decomposition and the quantum transfer matrix. This was basically the contents of the first four lectures presented at the school. In the remaining two lectures I started filling these notions with life by deriving an integral representation of the free energy per lattice cite for the Heisenberg-Ising chain (alias XXZ model) using techniques based on non-linear integral equations.
Up to small corrections the following sections L1-L6 display the lectures almost literally. The only major change is that the example of the XXZ chain has been moved from section L5 to L2. During the school it was not really necessary to introduce the model, since other speakers had explained it before. But for these notes I thought it might be useful to introduce the main example rather early. I also supplemented each lecture with a comment section which contains additional references and material of the type that was discussed informally with the participants.
I am grateful to my colleagues at the University of Wuppertal, Hermann Boos, Michael Karbach and Andreas Klümper, as well as to my long-term collaborators Karol Kozlowski and Junji Suzuki for sharing their considerable insight into the subjects of these lectures. L1 Statistical mechanics of quantum chains 1

.1 Introduction
Spin systems are the simplest conceivable quantum mechanical systems. In nature the spin occurs in first place as an internal degree of freedom of elementary particles. When many particles bind together in a many-body quantum system like a crystalline solid, the spin may also take the role of a discrete quantum number of collective excitations. In certain experiments on such systems, e.g. on Mott insulators at low temperatures or on ultra-cold atomic gases trapped in optical lattices, it is possible to create 'pure spin excitations'. Such systems are well described, in a certain energy range, by (generalised) Hubbard or Heisenberg models [16], which are in the class of models to be considered in these notes.
Spin systems can be used to illustrate the characteristic properties of quantum mechanics, like its probabilistic nature or the superposition and entanglement of states. For this reason quantum-spin systems are ubiquitous in introductory text books on quantum mechanics [18] and familiar to all graduate students in physics.
• A single spin- 1 2 is the simplest possible quantum system. Its Hilbert space is H = C 2 equipped with the Hermitian scalar product.
• Two spins- 1 2 describe the simplest interacting quantum systems with space of states H 2 = C 2 ⊗ C 2 . A two-spin state may be entangled.
• N many interacting spins- 1 2 constitute the simplest many-body quantum systems with space of states H N = (C 2 ) ⊗N . Depending on the interaction, these systems may possess complicated highly entangled ground states and may carry collective excitations of various types.
In the thermodynamic (or 'infinite volume') limit, N → +∞, quantum-spin systems may exhibit critical behaviour. They can be used to study phase transitions and quantum criticality. Apart from the number of constituents N , quantum-spin systems typically depend on several interaction parameters in their Hamiltonians. Considering certain scaling limits, in which these parameters depend on N and N is send to infinity, quantumspin systems may be used to realise quantum field theories on the lattice. Conversely, one may think of quantum-spin systems as of 'fully regularised quantum field theories' [17], meaning that ultra-violet and infra-red regularisations are provided by the fact that the number of spins is finite, and the Hilbert space is regularised due to the fact that spins have a finite number of degrees of freedom.
Arguably, all many-body quantum physics can be phrased in terms of spin systems of sufficiently general type. This should provide enough motivation to thoroughly study their statistical mechanics. It seems to indicate, on the other hand, that the statistical mechanics of quantum-spin systems in general should be too big as a subject for an introductory lecture course. For this reason, after developing part of a general theory, we shall restrict ourselves to integrable quantum-spin systems in these notes. Integrable systems are defined in one spatial dimension (1d) and have a rich algebraic structure underlying, which makes it possible to obtain more or less explicit results, at least for the thermodynamics of some non-trivial quantum spin systems in the infinite-volume limit.
As far as the general theory is concerned we shall introduce and explain what we call the quantum 'transfer matrix approach' to quantum spin systems. In our understanding this approach is a clever and systematic way of attaching an equilibrium statistical operator with any type of local interaction. The equilibrium statistical operator of 1d spin systems will be composed of transfer matrices of a 'classical vertex model' on a two-dimensional lattice. Such a procedure is non-unique. The non-uniqueness may be utilised to optimise the properties of the transfer matrix for various purposes. e.g. for an efficient calculation of the partition function by quantum Monte Carlo methods [44]. We shall use a construction, introduced by A. Klümper in [34], which is optimised for the use with integrable quantum spin models of Yang-Baxter type. As we shall see, this construction cannot only be used to calculate the free energy per lattice site of such models, but seems to be optimised as well for the calculation of their correlation functions [21].

States and operators
In our mathematical set-up we shall consider systems that are slightly more general than spin-1 2 systems in that they have d ≥ 2 degrees of freedom on a 'local Hilbert space' H = C d , equipped with the canonical Hermitian scalar product. We fix a basis in this space.

Operators on local Hilbert space
In order to introduce a space of local observables we start with a set of operators e α β ∈ End C d , α, β = 1, . . . , d, defined by their action on the basis (1.1), Then, for any A ∈ End C d , where (1.2) was used in the last equation. In (1.3) we have also employed the common 'summation convention', implying that Greek indices that occur twice in an expression are summed over from 1 to d. We shall keep this convention throughout these notes. Comparing left and right hand sides of (1.3) and taking onto account that the set {e α } d α=1 is a basis, we conclude that is a basis of End C d . Any basis element e α β will be called an elementary endomorphism. The action of a product of two elementary endomorphisms on the basis (1.1) can be computed by means of (1.3), e α β e γ δ e ϕ = e α β δ γ ϕ e δ = δ α δ δ γ ϕ e β = δ α δ e γ β e ϕ . (1.6) Comparing the left and the right hand side of this equation and using again that is a basis of C d , we obtain the relation providing a set of structure constants for the algebra End C d . From the first equation (1.3) we see that the identity operator I d ∈ End C d is represented by the matrix A β α = δ β α . Hence, by (1.4), (1.8)

Local basis of L-site Hilbert space
A multi-spin system is defined on a regular or irregular lattice consisting of N points x k ∈ R n , k = 1, . . . , N , such that a local Hilbert space H = C d is associated with every point. The Hilbert space of the multi-spin system is then H N = (C d ) ⊗N . Since we will soon focus on large integrable lattice systems, we shall assume that n = 1 and We define the embedding of the basis of elementary endomorphisms (1.5) into the lattice,

Examples
(i) If A ∈ End C d , then A j is a 'single-site' (or ultra-local) operator acting on 'site j'.
Then P x ⊗ y = e β α ⊗ e α β x γ e γ ⊗ y δ e δ = x γ y δ e β α e γ ⊗ e α β e δ = x γ y δ δ β γ δ α δ e α ⊗ e β = y ⊗ x . (1.12) Thus, P induces the transposition of factors in a tensor product. In physical terms, it interchanges the states on two sites. This operator, called the transposition or exchange operator, is an important object in the theory of spin systems and occurs in many places.
Most prominently, perhaps, it occurs as 'exchange interaction' in the Heisenberg Hamiltonian where J > 0 and P −L,−L+1 = P L,−L+1 by definition. This Hamiltonian defines one of the simplest and most generic interacting quantum-spin systems. It is simple in the sense that only neighbouring sites interact and also because P ( The operator P also plays in important role for the implementation of the action of spatial symmetries on quantum-spin systems, since it induces the action of the symmetric group S 2L on H 2L . For j, k, l ∈ {−L + 1, . . . , L} mutually distinct The symmetry group of a spin chain with an even number of sites is the dihedral group D 2L = C 2L × C 2 ⊂ S 2L which is the symmetry group of a regular polygon with 2L edges. Being a product of two cyclic groups it has two generatorŝ the 'shift operator'Û and the 'parity operator'P . Note thatÛ 2L =P 2 = id. Here is a third example for the occurance of P in the theory of quantum-spin systems. It provides a family of rational (or 'Yangian') solutions of the Yang-Baxter equation which is behind the integrability of the GL(d) invariant Hamiltonians (1.13). Define Then, using (1.15c), it is easy to see that if j, k, l ∈ {−L + 1, . . . , L} are mutually distinct.

Interactions
In the following we shall focus on quantum-spin chains with Hilbert space H 2L and with local interactions h ∈ End(C d ) ⊗m , where m ∈ {2, . . . , 2L} will be called the range of the interaction. Setting

Statistical mechanics of quantum-spin systems
After preparation in an experiment any quantum-spin system will be in a state described by a density matrix (a statistical operator) ρ L ∈ End H 2L with properties These properties guarantee that ρ L is diagonalizable and that the spectrum of ρ L is a discrete probability distribution. We may think of ρ L as representing an ensemble. Subsequent experiments then measure ensemble averages In general ρ L is time-dependent and its time dependence (in the Schrödinger picture!) is determined by the von-Neumann equation (1.24) Hence, a stationary density matrix should be a function of the conserved quantities commuting with the Hamiltonian.

Examples
In these notes we restrict ourselves to stationary density matrices. Some important examples are listed below.
(i) Many-body quantum systems cannot be separated from their environment forever. Eventually they relax to the canonical ensemble, However, transients and long relaxation times are possible, particularly for integrable quantum-spin systems, and stationary non-equilibrium ensembles may be realised in driven systems.
(ii) Two more examples are the zero and infinite-temperature limits of ρ L (T ).

Integrable quantum-spin systems
These notes will focus on large integrable quantum-spin systems. The problems we are going to address comprise: (i) A description of density matrices ρ L in a way compatible with the integrable structure, (ii) the calculation of the free energy per lattice site in the thermodynamic limit, (iii) the calculation of two-point functions of local operators in the thermodynamic limit, (1.30)

Comments
It is an interesting problem to describe the class of density matrices that may be generated as a result of the relaxation of an isolated integrable quantum spin system towards equilibrium. Due to the existence of a large number of additional local conserved quantities, the class of such density matrices must be much larger than the one-parametric family (1.25). In [38] the concept of generalised Gibbs ensembles was suggested with an inverse-temperature like Lagrange parameter for every additional conserved quantity. Conceptual difficulties arise from the fact that it is hard to identify 'a complete set of conserved local operators' in an infinite integrable system. We recommend [26] for further reading.
L2 The quantum transfer matrix

Fundamental models
A sufficiently general setting for a statistical mechanics of integrable quantum-spin systems is that of 'fundamental Yang-Baxter integrable models'. Fundamental integrable spin systems are entirely defined in terms of a matrix R(λ, µ) : Assuming differentiability in λ, µ in a vicinity of (0, 0) equations (2.1) imply another property of the R-matrix which is called unitarity: There is a function g : C 2 → C, differentiable in a neighbourhood of (0, 0), g(0, 0) = 1, g(λ, µ) = g(µ, λ), such that We may therefore assume in the following that R is normalised in such a way that 3) The proof of the existence of the function g is left as an exercise to the reader. With any R-matrix satisfying (2.1), (2.3) we associate two transfer matrices and a Hamiltonian where (P R) −L,−L+1 = (P R) L,−L+1 by definition and where h R ∈ C is a constant which may be used to render H Hermitian and to set the energy scale (alternatively one may rescale the spectral parameter). In order to obtain the second equation in (2.5) one has to use (2.3).

Trotter formula
Let and observe that, due to (2.5), It is not difficult to see that where · is the operator norm. Hence, This way the (unnormalised) density matrix of the canonical ensemble is represented as a product of transfer matrices. Equation (2.9) is sometimes called 'the Trotter formula', N 'the Trotter number'.

External fields
Assume there is Θ(α) = e αφ withφ ∈ End(C d ), α ∈ C such that which is called a U (1) symmetry of the R-matrix. Then this allows us to couple an external field to the Hamiltonian without spoiling its integrability.

Quantum transfer matrix
For N even introduce 'vertical spaces'1, . . . , N . By definition is the 'staggered and twisted inhomogeneous monodromy matrix' of the fundamental model. Here (Exercise: Prove it!). These are the 'Yang-Baxter algebra relations' for the (operatorvalued) matrix elements of T a (λ|α) considered as a matrix in 'auxiliary space a'. We shall call the transfer matrix associated with T a (λ|α) the 'quantum transfer matrix' of the fundamental model and denote it by t(λ|α) = tr a T a (λ|α) . (2.18) Using (2.9) and the Yang-Baxter equation (2.1a) we obtain . . , N/2, and α = κ/T . It is instructive (and not too hard) to derive (2.19) by algebraic means (see e.g. [21]). In the following, however, we shall introduce the graphical language of vertex models and use it for an intuitive and easily memorizable proof.

Example: the XXZ chain
The basic example of a fundamental integrable model is the XXZ spin-1 2 chain. Its Rmatrix (the R-matrix of the six-vertex model [2]) can be understood as a 'q-deformation' of the rational R-matrix (1.18) for d = 2. With a rescaling appropriate for our purposes it becomes .
It is a simple exercise to verify that (2.20) describes a one-parameter family of solutions of the Yang-Baxter equation (2.1a) that is regular (2.1b) and unitary (2.3).
In order to generate the Hamiltonian we differentiate where σ α , α = x, y, z, are Pauli matrices and ∆ = ch(η) by definition. Setting h R = 2J sh(η) in our general formula (2.5) we obtain the XXZ Hamiltonian where σ α −L = σ α L by definition. H XXZ is hermitian for all real J and ∆. A closer inspection of its discrete symmetries reveals that we may restrict ourselves to J > 0. Clearly meaning that R(λ, µ) has a U (1)-symmetry generated by Θ(α) = e αφ withφ = σ z /2. Thus, in this case 14)) and κ has the meaning of a Zeeman magnetic field coupling to the individual spins.

Graphical representation of integrability objects
A rather efficient way of dealing with the relation between the various 'integrability objects' introduced above, the R-matrix and the transfer and monodromy matrices, utilises a certain graphical representation [2].
(i) We identify every matrix element of R(λ, µ) with a 'vertex', (ii) Every arrangement of a finite number of directed, crossing lines is then in one-to-one correspondence with a product of R-matrix elements. We identify the connection of lines with summation over indices, e.g. λ ν µ (2.33) Then the U (1) symmetry (2.10) becomes

Comments
The graphical notation is generally useful in tensor calculus, not only in integrable systems. Other examples are Feynman diagrams, tensor networks and matrix product states.
L3 Partition function, density matrix and static correlation functions 3.1 Statistical operator of the (grand) canonical ensemble The graphical notation is appropriate for discussing the various types of (reduced) density matrices considered in the literature. We start with the most fundamental one which is the statistical operator of the grand canonical ensemble. Define Since, by (2.19), lim we call the finite Trotter number approximant to the statistical operator.
ρ N,L has the graphical representation Regarding the graph directly and in a reference frame rotated by π/2, the equality of left and right hand side (an hence the proof of (2.19)) becomes obvious and needs no further explanation. Note that the object in the dashed frame is the staggered, inhomogeneous monodromy matrix introduced above.

Partition function and free energy per lattice site
The thermodynamics of a quantum spin system is determined by its partition function or by its free energy respectively. Let us denote the eigenvalues of t(λ|κ/T ) by Λ n (λ|κ), n = 0, . . . , d N −1, and assume that they are ordered in such a way that If (i) the limits N → ∞ and L → ∞ commute and (ii) |Λ 0 (0|κ)/Λ n (0|κ)| < 1 for all N ∈ 2Z + and all n = 1, . . . , d N − 1, the expression of the free energy per lattice site greatly simplifies in the thermodynamic limit: It is determined by a sequence of non-degenerate eigenvalues Λ 0 of the quantum transfer matrix. We shall call these eigenvalues of largest modulus for fixed N the dominant eigenvalues, the corresponding eigenvectors |κ the dominant eigenvectors. Proving statements such as the commutativity of the Trotter limit and the thermodynamic limit usually falls into the realm of mathematical analysis. Physicists often take a pragmatic attitude toward more sophisticated mathematical questions, assuming everything to be alright until a counterexample appears. This is perhaps why statements (i) and (ii) above have not yet been rigorously justified. Still we have good reasons to believe that they are true for all fundamental integrable models and for all T > 0.
(i) Equation (3.8) was used to study the thermodynamics of many fundamental integrable lattice models, most notably of the Heisenberg XXX and XXZ spin chains [34,43] and of the Hubbard model [30]. These studies have been compared with other independent methods and give the same results within the available numerical accuracy. They reproduce the correct high and low-temperature behaviour and have the correct free Fermion limits in those cases where they exist.
(ii) For the special case of the XXZ chain a rigorous proof was recently provided for high enough finite temperatures [25]. Since the proof uses basically only the locality of the interaction of the model, it is expected to be generalizable at least to all fundamental integrable models.
(iii) The commutativity of limits was proved for a slightly differently defined quantum transfer matrix (which is unfortunately less compatible with the integrable structure imposed by the Yang-Baxter equation) by M. Suzuki in 1985 [44].
We would like to point out that the existence of the first limit in (3.8), defining the free energy per lattice site, was proved long time ago (see e.g. [39]). Thus, the commutativity of the limits would also imply the existence of the limit on the right hand side of the equation.
Remark. So far we did not use the integrability in any essential way. What we used was and the regularity where H (2) is a two-site Hamiltonian and the dots denote terms quadratic in λ and µ.
Simply defining for a given two-site Hamiltonian H (2) R H (2) (λ, µ) = P + (λ − µ)P H (2) (3.11) and using the U (1) symmetry of this object if it exists we obtain a quantum transfer matrix with dominant eigenvalue Λ 0 (0|κ) and (3.8) remains valid. Equation (3.8) can then be used as the starting point for the implementation of a numerical algorithm for the calculation of thermodynamic properties of infinite spin chain systems [41]. The integrability enters the game only when we calculate Λ 0 .

Density matrix of a chain segment (reduced density matrix) and expectation values of local operators
For any A ∈ End H 2L other than the identity operator there exist k, l ∈ {−L + 1, . . . , L}, k ≤ l and X ∈ End C d ⊗(l−k+1) such that A acts non-trivially on sites k and l and A = X k,k+1,...,l . The number of sites in [k, l] will be called the length of X, (X). Here (X) = l − k + 1. These notions still make sense for L → ∞. An operator, whose length stays finite for L → ∞ will be called local. For any local operator of length m, with non-trivial part X on the infinite chain its grand canonical expectation value is . .

= lim
N →∞ . (3.14) Here The reduced density matrix D m (T, κ) is obtained in the Trotter limit N → ∞, and, as can be read off from (3.14), Using the sequence D m (T, κ) m∈N of reduced density matrices we can calculate the expectation value of any local operator on the infinite chain (a procedure which is called an 'inductive limit'). Equation (3.15) is central for the theory of correlation functions of integrable models, in particular of the XXZ chain to be considered in more detail below.

Comments
Coming back to algebraic expressions we see that For any two ultra-local operators x, y ∈ End C d we define X(λ|κ) = tr{xT (0|κ/T )} and Y (λ|κ) = tr{yT (0|κ/T )}. Employing this notation and using (3.16) and (3.17) we obtain the following expression for the two-point correlation function of x and y: In the second line we have expanded t(0|κ/T ) m−1 tr{Y (0|κ)}|κ in an eigenbasis {|κ, n } of t(0|κ/T ). The series on the right hand side is what we call a 'thermal form factor series' [15]. It is a useful tool for analysing the large-distance asymptotic behaviour of the two-point function.

L4 The characterisation of reduced density matrices 4.1 Expectation values of local operators in pure states
Using the regularity relation (2.1b) we obtain the following graphical identities This generalises to If now |Ψ is any eigenstate of t ⊥ (λ) with eigenvalue Λ(λ), then (4.2) implies that These are two graphical representations of the reduced density matrix of the interval [1, m] associated with the eigenstate |Ψ . Remarkably, the expression on the right hand side is of the same form as in (3.14) (a twist can also be included).
Remark. We would like to emphasise the following: (i) We still did not use integrability here. The construction works for non-integrable models as well.
(ii) The above is a graphical version of 'the solution of the inverse problem' [23,32].
(iii) The problem of calculating the reduced density matrix from (3.15) or (4.3) is still largely unsolved, even for integrable models. Most of the results available in the literature refer of a few simple example systems related to the spin-1 2 XXZ chain on which we concentrate in the following.

Further generalisations of the reduced density matrix with the example of the XXZ chain
An idea in the spirit of Baxter [2], which was rather helpful for actually calculating the reduced density matrix of a chain segment, was to generalise the definition by attaching spectral parameters to the vertical lines [27,28] and a twist to one of the states [9]. We shall denote the corresponding un-normalised reduced density matrix by  It turns out that the new parameters ξ j , j = 1, . . . , m, and α = (κ − κ)/T regularise the mathematical expressions for D (N ) . The parameter α acquires a physical meaning in certain scaling limits, e.g. in the conformal limit [6].
(iii) By an 'exponential form', involving a double integral and the annihilation part of the so-called Fermionic basis [10]. * (iv) By the 'JMS theorem' using the creation part of the Fermionic basis [29].
Every single of these charactersations is technically involved and would justify a series of lectures on its own. For this reason we can only provide a brief description, which will necessarily stay somewhat vague, and a short guide to the literature. The first explicit description of the reduced density matrix of a chain segment of length m of the XXZ chain was obtained in [27] and had form of an m-fold multiple integral. It was derived for the ground state of the spin chain in the massive antiferromagnetic regime. The derivation relied on the construction of representations of a deformed vertex-operator algebra.
Subsequently the multiple-integral formula was rederived and generalised by different methods and different authors. An extension to the massless groundstate phase at vanishing magnetic field was obtained in [28]. The derivation relied on the use of functional equations of qKZ-type [19,42] that had been introduced before in order to characterise form factors of integrable massive quantum field theories.
A derivation by Bethe Ansatz [33] made it possible to take into account a finite magnetic field and opened the way to treat the finite temperature and finite length cases in [14,22]. A multiple-integral representation for the most general inhomogeneous and twisted case (4.6b) was eventually obtained in [3]. The non-vanishing density matrix elements are of the form where we have used the notation dm(λ) = dλ 2πi ρ(λ|κ, κ )(1 + a(λ, κ)) , dm(λ) = a(λ, κ)dm(λ) , (4.8) with ε + j the jth plus in the sequence (ε j ) m j=1 , ε − j the jth minus sign in the sequence (ε j ) m j=1 and p the number of plus signs in (ε j ) m j=1 . The definition of the integration contour C depends on which parameter regime is considered. For |∆| < 1, for instance, we may choose a rectangle centered around the origin of the complex plane with sides parallel to the real and imaginary axes and of height |η| − 0 + and large finite width. The functions a and G are solutions of convolution type integral equations with respect to the contour C and integration kernel K α (λ) = e −α cth(λ − η) − e α cth(λ + η) . (4.9) The function a solves the non-linear integral equation a(µ, κ)) , (4.10) whereas G is the solution of the linear integral equation Here and in (4.8) is the ratio of dominant eigenvalues with different twist parameters. The multiple-integral formula (4.7) is compact and memorizable but, as any true multiple integral, not very efficient for the actual computation of the density matrix elements [13]. Remarkably, however, it turned out that the multiple integrals factorise into sums over products of single integrals. This was explictly worked out with examples [3,4,11,12]. It motivated efforts to directly calculate the density matrix elements in factorised form. This is what is the main point behind items (ii)-(iv) above. In [7,8] a reduced form of the qKZ equation (rqKZ) was derived for the ground state of the XXX and XXZ chains and solved in a form corresponding to the factorised integrals. In [1] the equation was generalised to the finite temperature case. The effort to uncover the structure behind the factorisation led to the discovery of the Fermionic basis in [5,10] and eventually to a proof of the factorisation under very general conditions ('JMS theorem' [29]). The latter claims that all density matrix elements and hence all static correlation functions can be expressed in terms of only two basic functions, the function ρ, equation (4.12), and a function ω which for the inhomogeneous finite temperature case was characterised in terms of the functions a and G, equations (4.10), (4.11), in [3]: Factorisation has been used to calculate short-range static correlation functions of the XXZ and XXX chains (see e.g. [37,40] and references listed therein). Here is an example for the XXX chain at T = 0, h = 0 that was first obtained in [12], The Riemann ζ-functions arise from the function ω in the necessary limits. A more generic example of a finite temperature correlation function of the XXZ chain is Factorisation becomes rapidly inefficient for operators X of length (X) 10 or larger as the number of terms in the factorised form of the correlation functions grows very rapidly. In particular, it was so far inefficient for the calculation of the asymptotics of multi-point correlation functions. Yet, for these other methods, based e.g. on effective field theoretical descriptions [36] or on the use of form factor series [31], are available. It is an interesting open question, if factorisation of static correlation functions is a peculiarity of the XXZ quantum spin chain or rather a generic feature of integrable lattice models.

Properties of the generalised density matrix
While the derivation of the multiple-integral formula and the construction of the Fermionic basis are rather technical and certainly exceed the scope of these lecture notes, a characterisation of the reduced density matrix by its properties is comparatively easy.
These properties hold for any fundamental model (with U (1)-symmetry). Further properties follow from special properties of the R-matrix, e.g. group invariance, asymptotics as a function of the spectral parameter, crossing symmetry.
Remark. The left reduction fixes the one-point functions connected with the U (1)symmetry. Setting m = 1 and taking the derivative of (4.19b) with respect to κ at κ = κ we obtain Sending N → ∞, setting ξ 1 = 0 and using (3.8) we conclude that which is consistent with the standard thermodynamic relation.

Comments
Using the graphical notation it is not difficult to explain the reduced qKZ equation which has be used in order to calculate the reduced density matrix of the XXZ chain. We define Here we suppose that there are arbitrarily many horizontal lines and the spectral parameter on the lowest line is u. We further assume that |Ψ is a transfer matrix eigenstate. Then If the R-matrix exhibits crossing symmetry, like in case of the XXZ chain, then the arrow direction of the last transfer matrix on the right hand side can be reversed and the second equation can be interpreted as a discrete version [1] of the reduced qKZ equation [7]. In this section we consider the staggered and twisted inhomogeneous monodromy matrix (2.15), where R is the R-matrix (2.20) of the XXZ spin-1 2 chain. Then the auxiliary space 'a' is two-dimensional and T a can be interpreted as a 2 × 2 matrix with operator-valued entries acting on C 2 ⊗N , The Yang-Baxter algebra relations (2.17) are a set of quadratic relations for these entries. These relations allow one to construct a set of eigenvectors of the quantum transfer matrix generated over a pseudo vacuum |0 which has the properties for some complex functions a(λ), d(λ).
The existence of a pseudo vacuum is a non-trivial requirement. There are representations of the Yang-Baxter algebra which do not have a pseudo vacuum. For the quantum transfer matrix of the XXZ chain, however, a pseudo vacuum does exist. This can be easily inferred from the structure of the R-matrices composing the staggered monodromy matrix (2.15). They take the form we see that where a(λ) = e α/2 For the density matrix of the grand canonical ensemble we set as before ν Then we have the following Theorem. Algebraic Bethe Ansatz [35].
if {λ} is chosen in such a way that the 'Bethe Ansatz equations'

Auxiliary functions
We may assume that α and the ν j are generic. Otherwise we slightly change their values to make them generic. Then all eigenstates and eigenvalues of the quantum transfer matrix t(λ|α) can be labeled by solutions {λ and the auxiliary functions a n (λ) = d(λ)Q n (λ + η) a(λ)Q n (λ − η) . (5.14) By construction these functions have the important property that 1 + a n (λ

Nonlinear integral equations
Equation (5.15) allows us to characterise the auxiliary functions by means of nonlinear integral equations. Let C n be a simple closed contour that encircles {λ (n) j } Mn j=1 and the N/2-fold pole of a n at − h R N T , but no other poles or roots of 1 + a n . Then the following 'monodromy condition' holds, Cn dλ 2πi ∂ λ ln 1 + a n (λ) = M n − N 2 = −s n .
We shall call s n the '(pseudo-) spin of the nth excited state'.
In order to avoid case differentiations we restrict the parameter η from now on to η = −iγ, γ ∈ (0, π/2]. We would, like to point out, however, that a similar analysis is possible in all physically relevant parameter regimes. Define and Assuming that λ ± iγ is outside the contour C n defined above, for all λ ∈ C n and for all µ on an inside C n we obtain the identity Here we have to supply several comments and explanations. Generally, some care is necessary, when we take the logarithm of a meromorphic function and even more if we integrate it up along a contour. The first logarithm under the integral on the left hand side is defined by its principal branch. Then, due to our prerequisites, it defines a holomorphic function of µ inside and on the contour C n for all λ ∈ C n . The logarithmic derivative under the contour is meromorphic with simple poles with residue 1 at the Bethe roots and a simple pole with residue −N/2 at − h R N T . This explains the first equation. The second equation may be understood as fixing the branch in the definition of the logarithm of the function a n (λ). In the third equation we perform a partial integration of the integral on the left and side of the equation. This requires that we define the function ln(1 + a n ) as a holomorphic function, having no jumps of 2πi, as we move along the contour C n . For this purpose we fix any point x n ∈ C n and define a contour C λ xn running from x n to λ in positive direction along C n . Then ln(1 + a n )(λ) = C λ xn dµ ∂ µ ln(1 + a n (µ)) + ln 1 + a n (x n ) , (5.20) where the rightmost logarithm is defined by its principal branch, has the required properties. Equation (5.19) can be interpreted as nonlinear integral equation for the auxiliary function a n . For this purpose we rewrite it in the form ln a n (λ) = − s n ln sh(λ − x n + η) Note that the explicit dependence on x n vanishes for states with s n = 0. Another possibility to simplify the appearance of equation (5.21) occurs if the contour C n can be deformed in such a way that we can send Re x n → −∞. Then ln a n (λ) = 2iγs n − ε

Back to the roots
We just argued that every solution {λ (n) j } Mn j=1 of the Bethe Ansatz equations corresponds to an auxiliary function a n subject to the monodromy condition (5.16) and solving the nonlinear integral equation (5.21). This reasoning can be reversed, performing the following steps.
(i) Exponentiate (5.21). It follows that a n has an N/2-fold pole at λ = − h R N T (located inside C n ).
(ii) The monodromy condition (5.16) then implies that 1 + a n has precisely M n = N/2 − s n roots inside C n .
(iii) Going backwards through partial integration we see that these roots must satisfy the Bethe Ansatz equations (5.12).
Summing up, we have seen that the Bethe Ansatz equations are in one-to-one correspondence to pairs s n , C n , where the C n denote equivalence classes of simple closed contours.

Comments
(i) We did not explain the algebraic Bethe Ansatz in any detail, since this has been done elsewhere in the literature (e.g. [17,35]) and during the Les Houches 2018 summer school on Integrability in Atomic and Condensed Matter Physics.
(ii) For the description of the thermodynamics of the spin chain only the dominant eigenvalue Λ 0 and the corresponding auxiliary function a 0 are needed. These will be identified in the next section.
(iii) There are other methods for deriving non-linear integral equations and representations of the eigenvalues involving their solutions, most importantly the 'method of functional equations'. For further reading we recommend [43]. by considering the special cases ∆ = 0, T → 0, T → +∞ and by performing numerical calculations for small Trotter numbers N . For space-time limitations we restrict ourselves to the consideration of the high-T limit here in which we obtain a particularly clear and simple picture.

Bethe Ansatz at high temperature
Let us have a closer look at the explicit form of the auxiliary functions a n introduced in (5.14). Recalling that h R = −2iJ sin(γ) for the XXZ chain and that we agreed upon restricting γ to the interval (0, π/2] for simplicity, we see that and that this number becomes arbitrarily small for large T . Inserting it into the explicit expression for a n we obtain a n (λ) = e −κ/T sh(λ + iε) sh(λ − iε) is a solution to the Bethe Ansatz equations a n (λ (n) j ) = −1, j = 1, . . . , M . We would like to find 'perturbative solutions' for small ε > 0. One particular such solution is almost obvious. We shall see that it describes the dominant state in the high-temperature limit. This solution is 'generated' by the second factor on the right hand side of (6.2). The latter has an N/2-fold pole at iε and an N/2-fold zero at −iε. If ε is small, pole and zero are very close to each other and a model of this function for small |λ| is Close to the zero and close to the pole there are N/2 directions in which the phase of this function is iπ. They are connected by lines on which f ε is real negative and goes from zero to minus infinity. Since f ε (λ) takes values on the unit circle for λ on the real line, there are N/2 solutions λ j , j = 1, . . . , N/2, of the equation f ε (λ) = −1 on the real axis which all go to zero for ε → 0. This is sketched in figure 1. Thus, setting M n = N/2 and inserting the λ j for λ (n) j into the fourth factor on the right hand side of (6.2), we see that the product of third and fourth factor goes to 1 for ε → 0. Since the first factor goes to 1 as well in the high-T limit, we have obtained a special high-temperature solution.

A special high-temperature solution
In order to formalise this we set . . , N/2. We will look for a high-T solution of the Bethe equations with |x j | < R for some R > 0. Setting λ = x/T , inserting (6.4) into (6.2) and sending T → +∞ the Bethe Ansatz equations turn into 5) or, equivalently, Now p is a polynomial of order N/2 with asymptotics p(x) ∼ 2x N/2 for x → ∞. Thus, there are x 1 , . . . , x N/2 ∈ C such that

The corresponding eigenvalue
The corresponding eigenvalue is Here we have used (6.6) and (6.7) in the last equation.

Full spectrum the in high-temperature limit
Observe that due to the regularity (2.31) of the R-matrix. Clearly |u v| is a one-dimensional projector. Moreover, v|u = = 2 . (6.10) Thus, the spectrum of t ∞ is {2, 0}, where the eigenvalue 0 is 2 N − 1-fold degenerate. This means that the dominant state in the high-temperature limit is non-degenerate and has eigenvalue the Λ(0) = 2. Comparing with (6.8) above we see that the corresponding Bethe roots are λ j = x j /T , j = 1, . . . , N/2, where the x j are the roots of the polynomial p.

Bethe roots of the dominant state in the Trotter limit
It is not difficult to calculate these roots explicitly. For this purpose we have to solve the Bethe Ansatz equations in the high-temperature limit, 11) j = 1, . . . , N/2. Clearly, if x j is a root, then −x j is a root, and if N/2 is odd, then x j = 0 is a root. Taking the logarithm of (6.11) and setting we obtain, for any non-zero root x j , Using once more (6.12) and solving for x j we arrive at where, due to the π-periodicity of the tangent function, we may restrict the range of j to −N/4 + 1 ≤ j ≤ N/4 if N/2 is even or −N/4 + 3/2 ≤ j ≤ N/4 − 1/2 if N/2 is odd. In the Trotter limit, N → +∞, the roots get confined in the interval (2J sin(γ)/π) × [−1, 1] and accumulate at the origin. The outer roots converge to 2J sin(γ) (2j−1)π . This behaviour is illustrated with an example in figure 2.
Exercise: Find the other roots of the equation Answer: For T → +∞ we have N/2 roots close to +η and N/2 roots close to −η, which can be seen by setting λ = z/T ± η and sending T → +∞. with R > 0 large enough. As required, C 0 encloses all Bethe roots of the dominant state in the high-temperature limit, but no other root of the equation a 0 (λ) = −1 and no pole of this function other than the pole at λ = − h R N T . Since s = 0 for the dominant state, as seen above, we conclude that the auxiliary function of the dominant state satisfies the nonlinear integral equation with C 0 according to (6.16).
6.2 Trotter limit and free energy per lattice site 6.2.1 Trotter limit As follows from 6.1.5 the Bethe roots stay confined close to λ = 0 for N → +∞. We therefore obtain the auxiliary function in the Trotter limit by replacing ε = ln 1 + a 0 (λ) + κ 2T ln(1 + a 0 )(µ) mod 2πi = ln 1 + a 0 (λ) + κ 2T = ln Λ 0 (λ) mod 2πi . (6.20) Here C 0 is a modification of the contour C 0 such that C 0 − C 0 is a small positively oriented circle around λ. In the partial integration in the second equation we have used that s = 0, implying that there are no boundary terms. Equation (6.20) determines the eigenvalue in the Trotter limit. Recalling (3.8) we obtain the free energy per lattice site of the XXZ chain in the thermodynamic limit, where a 0 is the solution of the nonlinear integral equation (6.19). For the identification of the dominant state and the corresponding auxiliary function we have considered the high-temperature limit here. This brought us to the conclusion that s 0 = 0 and that a possible contour C 0 is the contour defined in (6.16). There are many good reasons to believe that (6.21) and (6.19) with the same choice of the contour hold for all T > 0.

Comments
(i) As we mentioned in the introduction, the latter claim is supported by numerical studies at finite Trotter number (for a pedagogical review see [24]), by a low-temperature analysis (see e.g. [15]) and by considering the XX chain (this is recommended as an exercise, for some information see [20]). In addition we would like to recommend the work [25], where the case of high but finite temperature was treated with full mathematical rigour.
(ii) The high-temperature analysis presented for the dominant state can be extended to obtain a large class of excited states in the high-temperature limit. We shall only sketch the calculation and leave the details as an exercise. Let us look for a solution {λ j } M j=1 of the Bethe Ansatz equations (5.12) that has the following high-temperature asymptotics: lim T →+∞ λ j = 0 for j = 1, . . . , n , (6.22a) λ j ∼ x j T with |x j | < R for some R > 0 for j = n + 1, . . . , M . (6.22b) Inserting this high-temperature Ansatz into the Bethe Ansatz equation (5.12) and performing the limit T → +∞, we see that the first n equations decouple and become sh(λ j + iγ) sh(λ j − iγ) s+n n k=1 sh(λ j − λ k − iγ) sh(λ j − λ k + iγ) = −1 , (6.23) for j = 1, . . . , n. These equations can be interpreted as a set of so-called higher-level equations for the high-T limit. They resemble the Bethe Ansatz equations of the spin-1 XXZ chain.
Taking the product over all j = 1, . . . , n in (6.23) we obtain the 'momentum quantisation condition' for some ∈ {0, 1, . . . , s + n − 1} that depends on {λ j } n j=1 . Inserting the λ j , j = n + 1, . . . , M , into the Bethe Ansatz equations (5.12), performing the limit T → +∞ and using (6.24), we obtain a set of equations that determine the x j , s+n . (6.25) Depending on N, , s, n this equation may admit a root x j = 0. All other roots are given by where k(j) is integer, if s + n is even, or half-odd integer, if s + n is odd. This means that we may choose the M − n roots x j from a set of N/2 inequivalent values, giving N/2 M −n different solutions. In [25] the high-T limit was worked out on more rigorous grounds, starting from the nonlinear integral equations for the excited states.