The equilibrium landscape of the Heisenberg spin chain

TL;DR

This paper characterizes the full set of local equilibrium states in the integrable Heisenberg spin chain using two frameworks, revealing a hidden structure and clarifying subtle features of equilibrium ensembles.

Contribution

It establishes the equivalence of the Thermodynamic Bethe Ansatz and transfer matrix approaches for describing equilibrium states in the Heisenberg chain, and uncovers a hidden structure in the Y-system.

Findings

Demonstrates equivalence of two frameworks for equilibrium states

Clarifies the breakdown of the canonical Y-system

Reveals hidden structures in equilibrium ensemble parametrization

Abstract

We characterise the equilibrium landscape, the entire manifold of local equilibrium states, of an interacting integrable quantum model. Focusing on the isotropic Heisenberg spin chain, we describe in full generality two complementary frameworks for addressing equilibrium ensembles: the functional integral Thermodynamic Bethe Ansatz approach, and the lattice regularisation transfer matrix approach. We demonstrate the equivalence between the two, and in doing so clarify several subtle features of generic equilibrium states. In particular we explain the breakdown of the canonical Y-system, which reflects a hidden structure in the parametrisation of equilibrium ensembles.

Equations

$$\mathbf{H}=J\sum _{i=1}^L\left(\frac{1}{4}-{\vec{\mathbf{S}}}_i\cdot {\vec{\mathbf{S}}}_{i+1}\right),$$

$$\Big{[}\mathbf{H},\mathbf{S}^{\alpha}_{\rm tot}\Big{]}=0,$$

$$e^{\mathrm{i}k(u_i)L}\prod _{j=1,j\ne i}^{M}S(u_i-u_j)=1.$$

$$k(u)=\frac{1}{\mathrm{i}}\log \left(\frac{u+\frac{\mathrm{i}}{2}}{u-\frac{\mathrm{i}}{2}}\right),$$

$$S(u)=\frac{u-\mathrm{i}}{u+\mathrm{i}},$$

$$E=J\sum_{i=1}^{M}\big{(}1-\cos{k(u_{i})}\big{)}.$$

$$\bigcup _{i=1}^M\{u_i\}⟼\bigcup _{j=1}^{\infty }\bigcup _{i=1}^{M_j}\bigcup _{a=1}^j\left\{u_{j,i}+(j+1-2a)\frac{\mathrm{i}}{2}\right\},$$

$$e^{\mathrm{i}{k}_j(u_{j,i})L}\prod _{j^{\prime }=1}^{\infty }\prod _{i^{\prime }=1}^{M_j}{S}_{j,j^{\prime }}(u_{j,i}-u_{j^{\prime },i^{\prime }})=-1,$$

$$k_j(u)=\frac{1}{\mathrm{i}}\log \left(\frac{u+j\frac{\mathrm{i}}{2}}{u-j\frac{\mathrm{i}}{2}}\right),$$

$$S_{j,\ell }(u)=\prod _{a=1}^j\prod _{b=1}^{\ell }S\left(u+(j-\ell -2a+2b)\frac{\mathrm{i}}{2}\right).$$

$${\rho }_j+{\bar{\rho }}_j=\frac{1}{2\pi }\left|\frac{\mathrm{d}{k}_j}{\mathrm{d}u}\right|-K_{j,\ell }\star {\rho }_{\ell },$$

$$K_{j,\ell }(u)=\frac{1}{2\pi \mathrm{i}}{\partial }_u\log S_{j,\ell }(u),$$

$$F_{j,\ell }\star f_{\ell }\equiv \sum _{\ell =1}^{\infty }{\int }_{-\infty }^{\infty }\mathrm{d}w{F}_{j,\ell }(u-w)f_{\ell }(w),$$

$$f\star g\equiv {\int }_{-\infty }^{\infty }\mathrm{d}wf(u-w)g(w),f\circ g\equiv {\int }_{-\infty }^{\infty }\mathrm{d}wf(w)g(w).$$

$$\mathfrak{s}[{\rho }_j,\overline{{\rho }_j}]=({\rho }_j+{\bar{\rho }}_j)\log ({\rho }_j+{\bar{\rho }}_j)-{\rho }_j\log {\rho }_j-{\bar{\rho }}_j\log {\bar{\rho }}_j,$$

$$(1-R{)}_{j,\ell }\star (1+K{)}_{\ell ,k}=1,$$

$$R_{j,\ell }(u)=I_{j,\ell }s(u),$$

$$s(u)=\frac{1}{2\cosh (\pi u)},$$

$$I_{j,\ell }={\delta }_{j-1,\ell }+{\delta }_{j+1,\ell }.$$

$$(s^{-1}\star f)(u)=f(u+\frac{\mathrm{i}}{2}-\mathrm{i}ϵ)+f(u-\frac{\mathrm{i}}{2}+\mathrm{i}ϵ),$$

$$\log\tau(u;w)\equiv\log\tanh{\big{(}\tfrac{\pi}{2}(u-w)\big{)}},\qquad w\in\mathcal{P},$$

$$\mathcal{P}=\left\{u\in \mathbb{C}:|\mathrm{Im}(u)|\le \frac{1}{2}-ϵ\right\},$$

$$s(u)=\mp \frac{1}{2\pi \mathrm{i}}{\partial }_u\log \tau (u;\pm \frac{\mathrm{i}}{2}).$$

$$\square =s^{-1}\star (1-R)=s^{-1}-I,$$

$$(\square f{)}_j(u)=f_j(u+\frac{\mathrm{i}}{2}-\mathrm{i}ϵ)+f_j(u-\frac{\mathrm{i}}{2}+\mathrm{i}ϵ)-f_{j-1}(u)-f_{j+1}(u),$$

$$G_{j,k}=(1+K{)}_{j,k}\star s.$$

$$K_j(u)=\frac{1}{2\pi \mathrm{i}}{\partial }_u\log S_j(u)=\frac{1}{2\pi \mathrm{i}}\left(\frac{1}{u-j\frac{\mathrm{i}}{2}}-\frac{1}{u+j\frac{\mathrm{i}}{2}}\right),S_j(u)=\frac{u-j\frac{\mathrm{i}}{2}}{u+j\frac{\mathrm{i}}{2}},$$

$$\frac{1}{2\pi }\left|\frac{\mathrm{d}{k}_j}{\mathrm{d}u}\right|=K_j(u),(1-R{)}_{j,\ell }\star K_{\ell }={\delta }_{j,1}s.$$

$$G_{j,k}(u)=\frac{1}{2\pi \mathrm{i}}{\partial }_u\log {\Psi }_{j,k}(u)=\sum _{a=1}^{\min (j,k)}{K}_{j+k+1-2a}(u),{\Psi }_{j,k}(u)=\prod _{a=1}^{\min (j,k)}{S}_{j+k+1-2a}(u),$$

$$K_{j,k}(u)=I_{j,\ell }{G}_{\ell ,k}(u),\log S_{j,k}(u)=I_{j,\ell }\log {\Psi }_{\ell ,k}(u).$$

Full text

The equilibrium landscape of the Heisenberg spin chain

Enej Ilievski1, and Eoin Quinn1,2

1 Institute for Theoretical Physics Amsterdam and Delta Institute for Theoretical Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

2 LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 91405 Orsay cedex, France

e.ilievski@uva.nl + eoin.quinn@u-psud.fr

Abstract

**We characterise the equilibrium landscape, the entire manifold of local equilibrium states, of an interacting integrable quantum model. Focusing on the isotropic Heisenberg spin chain, we describe in full generality two complementary frameworks for addressing equilibrium ensembles: the functional integral Thermodynamic Bethe Ansatz approach, and the lattice regularisation transfer matrix approach. We demonstrate the equivalence between the two, and in doing so clarify several subtle features of generic equilibrium states. In particular we explain the breakdown of the canonical $\mathcal{Y}$-system, which reflects a hidden structure in the parametrisation of equilibrium ensembles. **

Contents

1. 1 Introduction

2. 2 The Heisenberg spin chain

3. 3 Equilibrium ensembles

4. 4 Thermodynamic Bethe Ansatz

5. 5 Integrable lattice regularisation

6. 5.1 The two-dimensional vertex model

7. 6 The mirror system

8. 6.1 The scaling limit

9. 7 Discussion

10. 8 Conclusion

1 Introduction

The equilibration phenomena of quantum many-body systems have become a vigorous research topic for both theoretical and experimental studies of condensed matter systems in recent years. For generic interacting systems a central role has been played by the Eigenstate Thermalisation Hypothesis [1, 2, 3], which offers a unifying framework for characterising ergodic behaviour. The unconventional equilibration exhibited by (nearly) integrable systems has also drawn substantial interest, leading to the notion of the Generalized Gibbs Ensemble (GGE) [4, 5, 6]. The anomalous behaviour of integrable systems is due to the existence of infinitely many local conservation laws, whose essence is to protect quasi-particle excitations against decay [7]. The majority of the literature on equilibration in integrable systems has focused on non-interacting models, where the concept of a generalised Gibbs ensemble is synonymous to prescribing the occupations of single-particle modes [8, 9]. Interacting integrable systems on the other hand exhibit rich spectra of stable excitations which undergo non-trivial completely factorizable scattering [10]. This places interacting integrable systems in a distinguished position, and raises the question whether interactions induce physically discernible features among equilibrium states. The objective of this paper is to establish a framework to address this.

We focus our study on the isotropic Heisenberg spin-$1/2$ chain, a paradigmatic model of exactly solvable quantum many-body dynamics, due to both its simplicity and physical relevance. Our main result is an explicit construction of the entire manifold of equilibrium states, which helps expose a rich structure intrinsically linked to inter-particle interactions. We term this the ‘equilibrium landscape’.

Studies of the thermodynamic properties of exactly solvable models have been traditionally focused on canonical Gibbs equilibrium. Only in recent years has interest shifted towards the Generalized Gibbs ensembles, predominantly discussed in the context of quantum quench dynamics in several physically relevant models such as the anisotropic Heisenberg model and Lieb-Liniger Bose gas [11, 12, 13, 4]. Here prominence was given to simple initial states of potential experimental relevance, while more recent works have considered more generally a special class of ‘integrable’ product states [14, 15, 16, 17]. In the present work we pursue a general and systematic characterisation of the entire space of equilibrium ensembles without appealing to any initial state specific considerations.

Throughout the work we shall employ a range of techniques from the integrability toolbox, combining the algebraic, thermodynamic, and functional Bethe ansatz approaches. We begin by formulating an explicit algebraic construction of the GGE, and then proceed to analyse two complementary routes for evaluating equilibrium partition sums. On the one hand, the celebrated Thermodynamic Bethe Ansatz (TBA) approach [18, 19, 20] casts the partition function as a functional integral, invoking a spectral resolution through coupled interacting quasi-particle modes. A saddle-point of the functional integral yields an infinite set of coupled integral equations encoding the equilibrium state. On the other hand, we achieve a regularisation of a general partition function as a two-dimensional classical vertex model where, similarly as in the Gibbs canonical ensemble [21, 22, 23, 24, 25, 26, 27], the main subject of study is the dominant eigenvalue of a column transfer matrix. To our knowledge, no previous work on the statistical mechanics of exactly solvable models, including a large body of work on solvable classical vertex models, has achieved a similarly comprehensive description of the entire equilibrium manifold of a model.

For the case of canonical Gibbs equilibrium, compatibility of the two approaches to thermodynamics has been already demonstrated previously in [26]. By means of an integrable Trotterisation of the density operator, usually referred to as the Quantum Transfer Matrix, the Gibbs free energy can be expressed as a solution to the non-linear integral equation [24, 28, 25, 27]. In this regard, a special and seemingly non-generic analytic structure of the transfer matrix spectrum turns out to be crucial. In this work we demonstrate that typical macrostates from the equilibrium landscape have a much richer structure which necessitate going beyond the conventional Trotterisation techniques. Through achieving this we establish compatibility with the TBA formalism on the general grounds, and find that this yields a clear and instructive picture of the emergence of the equilibrium landscape.

2 The Heisenberg spin chain

In this work we outline our construction for perhaps the most widely studied interacting integrable quantum system, the one-dimensional isotropic spin-$1/2$ Heisenberg model [29, 19, 30, 20]

[TABLE]

with exchange coupling $J$ (here $J>0$ corresponds to a ferromagnetic ground state) and periodic boundary conditions $\vec{\mathbf{S}}_{L+1}=\vec{\mathbf{S}}_{1}$. The $\vec{\mathbf{S}}\equiv(\mathbf{S}^{\rm x},\mathbf{S}^{\rm y},\mathbf{S}^{\rm z})$ are the local generators of the $\mathfrak{su}(2)$ algebra, $[\mathbf{S}^{\alpha},\mathbf{S}^{\beta}]={\rm i}\,\epsilon_{\alpha\beta\gamma}\mathbf{S}^{\gamma}$. The model possesses a manifest global $\mathfrak{su}(2)$ symmetry,

[TABLE]

with $\mathbf{S}^{\alpha}_{\rm tot}=\sum_{i=1}^{L}\mathbf{S}_{i}^{\alpha}$. In this section we summarise the integrable structure of the model, introducing the concepts and notations which form the foundation of the work.

Spectrum.

The degrees of freedom of the spin chain are magnon excitations, corresponding to spin waves with respect to a reference fully polarised ferromagnetic vacuum. A key feature of integrability is that the magnons undergo non-trivial yet non-diffractive scattering, implying that any interaction process can be reduced to a sequence of two-particle scatterings. With respect to the vacuum state, the quantisation conditions for the magnons are known as the Bethe equations [29]

[TABLE]

Here the magnons are conveniently parametrised by a rapidity variable $u$, through which their momentum is

[TABLE]

and the two-magnon scattering amplitude is given by

[TABLE]

which depends only on the difference of rapidities. Due to the $\mathfrak{su}(2)$ symmetry of the model the eigenstates arrange into degenerate $\mathfrak{su}(2)$ multiplets, and each highest-weight state (with respect to the orientation of the reference ferromagnetic vacuum) corresponds to a set of magnon rapidities $\{u_{i}\}$ satisfying Eqs. (2.3). The corresponding energy eigenvalue has the additive form

[TABLE]

The descendant states in a multiplet are obtained by adding zero momentum magnons, i.e. rapidities with $u=\infty$, for which the scattering amplitude trivialises.

Bound states.

The scattering between magnons induces bound state formation. These correspond to the collections of complex magnon rapidities aligning into ‘string’ patterns in the complex rapidity plane. In the large-$L$ limit the Bethe roots are classified according to the ‘string hypothesis’,

[TABLE]

with all $u_{j,i}\in\mathbb{R}$. Bound states of $j$ magnons are accordingly called $j$-strings, and the set of $j$-strings provide the thermodynamic particle content of the model, i.e. $M=\sum_{j=1}^{\infty}M_{j}$. The string rapidities are subject to the ‘string Bethe equations’,

[TABLE]

valid up to corrections which are suppressed in system size $L$. Here the bare momentum of a $j$-string is

[TABLE]

and the scattering amplitudes between a $j$-string and a $\ell$-string are

[TABLE]

The factor $-1$ on the right-hand side of Eq. (2.8) compensates the self-scattering factor $S_{j,j}(0)=-1$ from the left-hand side.

Macrostates.

In the thermodynamic limit, defined as $L,M\to\infty$ with ratio $M/L$ kept fixed, the string rapidities distribute densely along the real line. Macrostates are characterised by the complete set of the string rapidity densities $\{\rho_{j}(u)\}$, with $L\rho_{j}(u){\rm d}u$ being the number of occupied $j$-string modes in an infinitesimal rapidity interval ${\rm d}u$ around $u$. These obey the log-differential form of Eq. (2.8), the Bethe–Yang equations,

[TABLE]

where $\bar{\rho}_{j}(u)$ denotes the corresponding density of holes, i.e. unoccupied modes, and the scattering kernels

[TABLE]

are the differential scattering phases. Here we use the following short-hand notation for matrix convolutions

[TABLE]

and adopt the summation convention for the repeated indices. In addition, we also use a short-hand notation for scalar integrations as follows

[TABLE]

For each macrostate there is an associated entropy, which is the logarithm of the number corresponding microstates. This is expressed through the entropy density functional [18, 20]

[TABLE]

where \exp\big{(}L\mathfrak{s}\big{[}\rho_{j}(u),\bar{\rho}_{j}(u)\big{]}\,{\rm d}u\big{)} counts the number of ways of distributing $L\rho_{j}(u)\,{\rm d}u$ particles between the L\big{(}\rho_{j}(u)+\bar{\rho}_{j}(u)\big{)}\,{\rm d}u many $j$-string mode numbers on an infinitesimal rapidity interval ${\rm d}u$ centred at $u$.

Kernel identities.

The scattering kernels $K_{j,\ell}$ are differential scattering phase shift which encode interactions at the level of macrostates. They exhibit a rich structure which we will exploit throughout this work. Firstly, the Fredholm operator $(1+K)$ admits a pseudo-inverse $(1-R)$ through

[TABLE]

with $1\equiv\delta_{j,k}\delta(u)$, where the Fredholm resolvent,

[TABLE]

is defined through the $s$-kernel

[TABLE]

and the nearest-neighbour incidence matrix $I$,

[TABLE]

Here $(1-R)$ admits a non-trivial nullspace and is thus not the true inverse of $(1+K)$. In particular, the relation $(1-R)_{j,\ell}\star n_{\ell}=0$, with boundary condition $n_{0}=0$, has a one-parameter solution $n_{j}=h\,j$. Similarly the $s$-kernel admits a pseudo-inverse $s^{-1}$, which is a left-inverse under convolution, i.e. $s^{-1}\star s\star f=f$, and which also possesses a non-trivial nullspace. Explicitly it is,

[TABLE]

where $\epsilon\equiv 0^{+}$ is an essential positive infinitesimal which prescribes the avoidance of the poles of $s(u)$ at $u=\pm\tfrac{{\rm i}}{2}$. The nullspace of $s^{-1}$, generated by functions $\zeta$ obeying $s^{-1}\star\zeta=0$, is a linear span of basis functions $\log\tau(u;w)$ of the form

[TABLE]

where $\mathcal{P}$ is a strip in the complex plane defined as

[TABLE]

commonly referred to as the ‘physical strip’. We highlight the explicit dependence on the infinitesimal regulator $\epsilon$ here, which ensures that the boundaries ${\rm Im}(u)=\tfrac{1}{2}$ are excluded from the strip, as it will prove useful in later sections. The functions $\log\tau(u;w)$ are related back to the $s$-kernel through the identity

[TABLE]

A related object is the discrete d’Alembertian operator

[TABLE]

or explicitly,

[TABLE]

for a set of functions $f_{j}(u)$. The associated Green’s function, obeying $\square\,G=1$, is given by

[TABLE]

The d’Alembertian $\square$ inherits a non-trivial nullspace from both $(1-R)$ and $s^{-1}$. From $(1-R)$ the functions $n_{j}=hj$ obey $\square\,n=0$, while given a set of functions $\zeta_{j}$ in the nullspace of $s^{-1}$, there exist related functions $\nu_{j}=(1+K)_{j,\ell}\star\zeta_{\ell}$ which also satisfy $\square\,\nu=0$. It is further useful to define kernels $K_{j}$ and their associated amplitudes $S_{j}$ as follows

[TABLE]

with the kernels obeying the identities

[TABLE]

These provide convenient explicit expressions for the matrix elements of the Green’s function $G_{j,k}$ and its associated amplitude $\Psi_{j,k}$ as follows

[TABLE]

and in turn for the scattering kernels and their amplitudes through

[TABLE]

Finally, we emphasise that the $\epsilon$ regulator is tied to the pseudo-inverse $s^{-1}$, and in particular does not appear in the string compounds in Eq. (2.7), nor in the general definitions of kernels and amplitudes, e.g. Eqs. (2.10), (2.27). In the following we employ a compact notation for half-unit imaginary shifts not involving a regulator

[TABLE]

Transfer matrices.

The algebraic formulation of integrability is founded upon the Lax representation [31, 32, 33], see also [34, 35, 36, 37] and references therein. The Lax matrices are a family of operators $\mathbf{L}_{k,1}(v,u):\mathcal{V}_{k}\otimes\mathcal{V}_{1}\to\mathcal{V}_{k}\otimes\mathcal{V}_{1}$, where $\mathcal{V}_{k}$ denotes the $(k+1)$-dimensional irreducible unitary representation of $\mathfrak{su}(2)$, of the form

[TABLE]

These provide local building blocks for the transfer matrices, a commuting family of operators acting on the Hilbert space $\mathcal{H}\cong\mathcal{V}^{\otimes L}_{1}$, which are given as traces over path-ordered products of operators $\mathbf{L}_{k,1}$,

[TABLE]

where the trace is taken over the common auxiliary space $\mathcal{V}_{k}$, with $k\in\mathbb{N}$. Trivially, $\mathbf{T}_{0}(v)=v^{L}$. Integrability ensures that the transfer matrices $\mathbf{T}_{k}(v)$ mutually commute,

[TABLE]

for all values of $k,k^{\prime}\in\mathbb{N}$ and $v,v^{\prime}\in\mathbb{C}$.

Fusion hierarchy.

The eigenvalues $T_{k}(v)$ of the transfer matrices $\mathbf{T}_{k}(v)$ are called $T$-functions. They are polynomial and satisfy the Hirota equation111The Hirota equation (2.35) can be understood as the ‘quantum’ counterpart of the fusion identities for characters of the ‘classical’ Lie algebra $\mathfrak{su}(2)$, i.e. fusion identities amongst unitary irreducible representations of $\mathfrak{su}(2)$. Additional details on fusion identities can be found in e.g. [38, 24, 39, 34, 40, 35].

[TABLE]

with initial conditions $T_{-1}\equiv 0$, $T_{0}(v)=v^{L}$, and boundary ‘scalar potentials’,

[TABLE]

The Hirota equation exhibits a gauge freedom corresponding to the overall normalisation of the Lax matrix. There however exist $Y$-functions,

[TABLE]

which are gauge-invariant quantities satisfying the canonical $Y$-system hierarchy [41, 42]

[TABLE]

with initial condition $Y_{0}(v)=0$.

Thermodynamic inversion identity.

In the large-$L$ limit, the entire family of transfer matrices $\mathbf{T}_{k}(v)$ satisfy the useful identity [43]

[TABLE]

which allows for their inversion. This property can equivalently be expressed as the large-$L$ decay of the physical $Y$-functions

[TABLE]

Local charges.

The transfer matrices serve as generating operators for the local charges through their logarithmic derivatives [43, 4],

[TABLE]

These charges are well-defined only on the physical strip, $v\in\mathcal{P}$, specified above in Eq. (2.22). When $v$ approaches the boundary of the strip, the charges acquire a divergent localisation length [43, 44] and thus become singular at the boundaries $\partial\mathcal{P}={\rm Im}(v)=\tfrac{1}{2}$.

The Hamiltonian Eq. (2.1) is given by $\mathbf{H}=\pi J\,\mathbf{X}_{1}(0)$.

String charge-duality.

In the thermodynamic limit, the eigenvalues of the charges $\mathbf{X}_{k}(v)$ are expressible as a linear functional of the rapidity densities [7]

[TABLE]

where $G$, given in Eq. (2.26), is the Green’s function of the d’Alembertian $\square$. The inverted relation, promoted to the level of operators,

[TABLE]

bears the name ‘string-charge duality’. We highlight that this defines the mode operators $\boldsymbol{\rho}_{j}$ in a manner which is independent of the of the orientation of the reference ferromagnetic vacuum.

Even though the positive infinitesimal $\epsilon$ does not appear in the definition of the charges $\mathbf{X}_{k}(v)$, the mode operators $\boldsymbol{\rho}_{j}$ inherit the $\epsilon$-prescription through the left-inverse $s^{-1}$ which enters in the d’Alembertian $\square$, Eq. (2.25). The important consequence of the regulator $\epsilon$ is that the boundaries $\partial\mathcal{P}$ at ${\rm Im}(u)=\tfrac{1}{2}$ are avoided, ensuring a finite localisation length of the $\boldsymbol{\rho}_{j}$, i.e. as $\epsilon$ is strictly positive the localisation lengths are strictly finite, cf. [43, 44]. Thus $\epsilon$ admits a physical interpretation as a regulator which governs the notion of locality in the large-$L$ limit.

3 Equilibrium ensembles

The purpose of this article is to characterise the equilibrium landscape of an interacting integrable model, specifically the Heisenberg spin chain. Equilibrium ensembles emerge dynamically in the long-time limit of unitary evolution from generic initial states in thermodynamically large systems. These can be viewed equivalently in either the canonical sense as density matrices governing the expectation values of all local observables, or in the microcanonical sense as unbiased collections of eigenstates sharing the same values of all local charge densities. There are two key mechanisms underlying local equilibration: (i) decoherence, causing the dynamical phases between individual eigenstates to average out during the relaxation process, and (ii) the eigenstate thermalisation hypothesis which supposes the local equivalence of distinct eigenstates from the same microcanonical shell [1, 2, 3, 45]. Local correlation function are thus expressible as functionals of the quasi-particle densities characterising a macrostate, see e.g. [46, 47, 48].

For the Heisenberg spin chain, a microcanonical shell corresponds to the set of microstates associated with a given macrostate, parametrised by the full set of occupied mode distributions $\rho_{j}(u)$. In this context, the statement of the eigenstate thermalisation hypothesis is substantiated by the Yang–Yang form of the entropy density, Eq. (2.15), see e.g. [49, 6, 50]. The corresponding (unnormalised) density matrix is [51]

[TABLE]

where $\mu_{j}(u)$ provide a general set of chemical potentials for the mode operators $\boldsymbol{\rho}_{j}$. The term $\vec{h}\cdot\vec{\mathbf{S}}_{\rm tot}$ incorporates a general Cartan charge of the global $\mathfrak{su}(2)$ symmetry, which serves to specify the polarisation direction $\vec{h}=(h_{\rm x},h_{\rm y},h_{\rm z})$ with respect to which the Bethe magnons are defined, i.e. with respect to a ferromagnetic vacuum oriented in the $\vec{h}$ direction. For consistency, the chemical potentials must not diverge with $j$, i.e.

[TABLE]

The above functional parametrisation of the density matrix differs from the more commonly used definition involving a formal infinite sum over a discrete basis of local conservation laws (see e.g. [52, 4, 6, 5]) or a ‘truncated’ GGE [53, 54, 55]. We argue however that gauge-invariant formulation (3.1) not only clearly conveys the physical picture underlying the GGE concept, it moreover provides a natural and convenient starting point for analysis of the ensemble as developed in the following sections.

4 Thermodynamic Bethe Ansatz

In this section we revisit the formalism of Thermodynamic Bethe Ansatz, a functional integral formulation of thermodynamics [18, 19, 20]. Partition sums are cast in the basis of Bethe eigenstates, which in the thermodynamic limit translates to functional integration in the space of macrostates. In particular, the (generalized) free energy density,

[TABLE]

is reformulated as a functional integral over the string rapidity distributions,

[TABLE]

with $h=|\vec{h}|$, and the entropy functional $\mathfrak{s}$ is specified by Eq. (2.15). Identifying the saddle point of the equilibrium partition sum through $\delta f/\delta\rho_{j}(u)=0$, subject to the constraints (2.11), yields the celebrated TBA equations,

[TABLE]

expressed in terms of the ‘thermodynamic $\mathscr{Y}$-functions’,

[TABLE]

The equilibrium free energy density then becomes

[TABLE]

which can be equivalently expressed in the form

[TABLE]

obtained by inserting Eq. (2.16) in Eq. (4.5), and making use of the identities Eq. (2.28), the TBA equations Eq. (4.3), and that $hj$ belongs to the nullspace of $(1-R)$.

The TBA equations provide the link between the chemical potentials $\mu_{j}$ and the thermodynamic functions of general equilibrium states. To further elucidate the underlying structure, we next bring the TBA equations to a local form by convolving with the left-inverse $(1-R)$ of the Fredholm operator $(1+K)$, resulting in

[TABLE]

with $\mathscr{Y}_{0}\equiv 0$, and source terms

[TABLE]

Here information about $h$, which has dropped from the source term as it belongs to the nullspace of $(1-R)$, instead appears implicitly through the large-$j$ asymptotics

[TABLE]

which must be supplemented to Eqs. (4.7) in order to unambiguously fix a solution. The information stored in $\mu_{j}$ is preserved by $(1-R)$, as condition Eq. (3.2) forbids a nullspace contribution, and so Eq. (4.8) can be readily inverted

[TABLE]

The TBA equations are integral equations defined on the real rapidity axis. We now analytically continue them to the complex rapidity plane, obtaining an equivalent set of functional relations called the ‘thermodynamic $\mathscr{Y}$-system’ [51]. This is achieved by convolving both sides of Eqs. (4.7) with the pseudo-inverse $s^{-1}$ and subsequently exponentiating, resulting in

[TABLE]

with

[TABLE]

This provides a natural decomposition of the source terms,

[TABLE]

where $\zeta_{j}$ satisfy

[TABLE]

Owing to Eq. (2.21), we adopt the following generic parametrisation

[TABLE]

in terms of parameters $\alpha_{j,i}\in\mathbb{R}$, and complex-conjugate $w_{j,i}$ and $\bar{w}_{j,i}$ located in $\mathcal{P}$. Although the sum involves a finite number of terms, infinite convergent sequences of simple poles and zeros ($\alpha_{j,i}=\pm 1$) which condense along certain contours are also permissible and can be understood as limits of Padé approximants for complex functions with branch points. We adopt the convention that all branch cuts in $\mathcal{P}$ extend vertically away from the real axis.

The decomposition Eq. (4.13) can also be expressed at the level of the chemical potentials

[TABLE]

where here $\lambda=\square\,\mu$, and $\nu_{j}$ encodes the nullspace of $\mu_{j}$ inherited from $s^{-1}$ through

[TABLE]

Indeed the thermodynamic $\mathscr{Y}$-system can be obtained directly from the TBA equations Eq. (4.3) by applying the d’Alembertian $\square$ and exponentiating.

To this point the analysis appears completely formal. It however reveals a structure in the space of equilibrium states. The information from the nullspace of $s^{-1}$, Eq. (2.21), which appears to have dropped out from Eq. (4.11), enters instead implicitly via analytic properties of functions $\mathscr{Y}_{j}$. Specifically, $w_{j,i}$ and $\bar{w}_{j,i}$ are nothing but branch points of $\mathscr{Y}_{j}$ of degree $\alpha_{j,i}$. To establish this, we now re-obtain Eq. (4.7) from Eq. (4.11). In order to undo the complex shifts and transform Eq. (4.11) to the real axis, we introduce adapted $\mathscr{Y}$-functions,

[TABLE]

which obey ${\mathscr{Y}}^{{\rm adp}}_{j}(u+\tfrac{{\rm i}}{2}-{\rm i}\epsilon){\mathscr{Y}}^{{\rm adp}}_{j}(u-\tfrac{{\rm i}}{2}+{\rm i}\epsilon)={\mathscr{Y}}_{j}(u+\tfrac{{\rm i}}{2}-{\rm i}\epsilon){\mathscr{Y}}_{j}(u-\tfrac{{\rm i}}{2}+{\rm i}\epsilon)$ due to $\tau^{+}\tau^{-}=1$, and are (by construction) analytic on $\mathcal{P}$ with constant large-$u$ asymptotics. Substituting these into Eq. (4.11), taking the logarithm and convolving with $s$ we readily recover Eq. (4.7). The analytic properties of the thermodynamic $\mathscr{Y}$-functions in $\mathcal{P}$ are therefore completely determined by the data $w_{j,i}$ and $\bar{w}_{j,i}$ and $\alpha_{j,i}$. If all $\alpha_{j,i}\in\mathbb{Z}$ then the $\mathscr{Y}$-functions are meromorphic on $\mathcal{P}$, while more generally they possess non-integer branch points.

5 Integrable lattice regularisation

In this section we derive an explicit lattice regularisation of a generic equilibrium ensemble. This is achieved this by re-expressing the density matrix as a product of transfer matrices, thereby casting it as a two-dimensional classical vertex model on a cylinder.

The mode operators $\boldsymbol{\rho}_{j}$ are related to the transfer matrices via string-charge duality, Eq. (2.43). The density matrix Eq. (3.1) can thus be presented in the form

[TABLE]

where the charges $\mathbf{X}_{j}$ are the logarithmic derivatives of the transfer matrices, Eq. (2.41). To proceed it is necessary to be careful about the potential nullspace of $\mu_{j}$ under action of the d’Alembertian $\square$. As described in Section 2, this nullspace is spanned by functions $n_{j}$ in the nullspace of $(1-R)$, along with functions $\zeta_{j}$ in the nullspace of $s^{-1}$. The condition Eq. (3.2) forbids a contribution $n_{j}$, and so the chemical potential decomposes as

[TABLE]

where $\lambda=\square\,\mu$, and $\nu_{j}=(1+K)_{j,\ell}\star\zeta_{\ell}$ encode the nullspace inherited from $s^{-1}$, $\square\,\nu=0$. This essentially repeats the derivation of Eq. (4.16) in Section 4. In the following we will adopt the same generic parametrisation of $\zeta_{j}$ as in Eq. (4.15).

The decomposition of $\mu_{j}$ induces a factorisation of the density matrix Eq. (3.1) into three commuting factors

[TABLE]

given respectively by

[TABLE]

where the first factor is obtained through (G\star\lambda)\circ\big{(}\square\,\mathbf{X}\big{)}=\big{(}\square\,G\star\lambda\big{)}\circ\mathbf{X}=\lambda\circ\mathbf{X}. For convenience we will refer to $\boldsymbol{\varrho}_{\lambda}$ as encoding the ‘node data’, $\boldsymbol{\varrho}_{\nu}$ as encoding the ‘analytic data’, and $\boldsymbol{\varrho}_{\vec{h}}$ as encoding the Cartan charge. We now proceed to analyse each factor in turn.

The node data.

The first factor of the ensemble encodes the node data $\lambda_{k}$,

[TABLE]

To begin we regularise the $\lambda_{k}$, by invoking a large-rapidity cut-off $\Lambda^{\infty}$, and casting $\lambda_{k}(u)$ as a discrete sum on the remaining interval. First it is useful to introduce, for each $k$ separately, two disjoint domains $\{-\Lambda^{\infty},\Lambda^{\infty}\}=\mathcal{I}^{(+)}_{k}\cup\mathcal{I}^{(-)}_{k}$, such that $\lambda_{k}(v)\geq 0$ on $\mathcal{I}^{(+)}_{k}$ and $\lambda_{k}(v)<0$ on $\mathcal{I}^{(-)}_{k}$. That is, we decompose the $\lambda_{k}$ as

[TABLE]

where $\lambda^{(\pm)}_{k}(v)=0$ on $\mathcal{I}^{(\mp)}_{k}$. Each part can then be regularised via a discrete sum of Dirac $\delta$-distributions,

[TABLE]

where $n^{(\pm)}_{k}\in\mathbb{N}$ control the resolution of the discretisation, the parameters $x^{(\pm)}_{k,i}\in\mathbb{R}$ are chosen uniformly from the distributions $\lambda^{(\pm)}_{k}(v)$, and the overall normalisation is fixed by

[TABLE]

This readily translates to a further factorisation $\boldsymbol{\varrho}_{\lambda}=\boldsymbol{\varrho}^{(+)}_{\lambda}\cdot\boldsymbol{\varrho}^{(-)}_{\lambda}$, with

[TABLE]

Referring back to the definition of the charges Eq. (2.41), each exponential factor is written as

[TABLE]

Using the inversion identity Eq. (2.39), which is exact in the large-$L$ limit, this becomes

[TABLE]

Now, for each factor in Eq. (5.9) we couple the finite difference to the corresponding discretisation parameters $n^{(\pm)}_{k}$ through the identifications $v\to x^{(\pm)}_{k,i}$ and $\Delta v\to{\rm i}\,\xi^{(\pm)}_{k}$, with parameters

[TABLE]

resulting in

[TABLE]

Let us highlight that while Eq. (5.13) may in a sense be viewed as an ‘integrable Trotterisation’, the Suzuki–Trotter formula [56, 57] has not been invoked. Instead, we coupled the resolution of the discretisation of $\lambda_{k}(u)$ in Eq. (5.7) with the finite-difference approximation $\Delta v$ of the derivative in the definition of the charges.

The analytic data.

We next analyse the factor

[TABLE]

where $\nu=(1+K)\star\zeta$, with $\zeta_{j}$ given generically by Eq. (4.15). Here we substitute $\boldsymbol{\rho}=\square\,\mathbf{X}$, bringing it to the form

[TABLE]

Then exploiting the null property $\square\,\nu=0$, the exponent becomes recast as a sum of two contour integrals

[TABLE]

where contours $\mathcal{C}_{+}$ and $\mathcal{C}_{-}$ encircle the upper and lower half of the physics strip $\mathcal{P}$ in the counter-clockwise direction, respectively. Here we have shifted the integration contours using

[TABLE]

along with the asymptotic behaviour $\lim_{|u|\to\infty}\nu_{j}(u)=0$. Introducing the functions

[TABLE]

and integrating by parts, we obtain

[TABLE]

Specifically, the functions $\gamma_{j}$ are given through

[TABLE]

using identities (2.23) and (2.26).

Now we evaluate the above contour integrals. Recalling the explicit expression for $G$ from Eq. (2.29), and noting that the poles of kernels $K_{j}(u)$ are of the form, $2\pi{\rm i}\,{\rm Res}_{u=0}K_{j}(u\pm\tfrac{{\rm i}}{2}k)=\pm\delta_{jk}$, we observe that the poles of $\gamma^{\pm}_{j}(u)$ with residues $\pm\alpha_{j,i}$ are located at the null components shifted by integer multiples of $\tfrac{{\rm i}}{2}$. Given that all the null components lie within the physical strip $\mathcal{P}$, the only contribution which does not shift the poles outside the contours then comes from the term $K_{1}(u)$ which appears only in $G_{j,j}(u)$. In effect, the analytic data gets shifted by $\pm\tfrac{{\rm i}}{2}$, giving a pole in $\mathcal{C}_{-}$ ($\mathcal{C}_{+}$) for each null component in $\mathcal{C}_{+}$ ($\mathcal{C}_{-}$). Thus Eq. (5.19) simplifies to

[TABLE]

upon using the inversion identity Eq. (2.39).

Finally, Eq. (5.21) must be further regularised in order to convert to a product of transfer matrices. To achieve this we note the following property: any term $\alpha\log[\tau(u;w)\tau(u;\bar{w})]$ with $0<\alpha<1$ and $w,\bar{w}\in\mathcal{P}$ can be systematically approximated by a sum

[TABLE]

where the sets of zeros $\{\omega^{\rm z}_{a},\bar{\omega}^{\rm z}_{a}\}$ and poles $\{\omega^{\rm p}_{b},\bar{\omega}^{\rm p}_{b}\}$ are obtained as the zeros and poles lying within the physical strip $\mathcal{P}$ of a Padé approximation of $(u-w)^{\alpha}(u-\bar{w})^{\alpha}$ about $w_{0}={\rm Re}(w)={\rm Re}(\bar{w})$. Applying this to each contribution of $\zeta_{j}$ individually we obtain a regularised form

[TABLE]

where both the distinction between distinct zero modes and dependence on the degree of the Padé approximation are left implicit. In this way we obtain the factor $\boldsymbol{\varrho}_{\nu}$ of the ensemble as a product of transfer matrices

[TABLE]

The Cartan charge.

The final factor of the ensemble is simply given by

[TABLE]

There is no trace over an auxiliary space associated with this factor.

5.1 The two-dimensional vertex model

Now recombining the three factors the ensemble takes the compact form

[TABLE]

where we have collected the impurity parameters as follows

[TABLE]

along with their complex conjugates $\{\bar{\theta}_{k,i}\}$.

The regularised equilibrium ensemble is naturally cast as a two-dimensional classical vertex model wrapped around a cylinder of circumference $L$ and height $2N$, with

[TABLE]

as illustrated in Figure 2. The ‘Boltzmann weights’ are formally identified with the matrix elements of Lax matrices, through the Lax representation of transfer matrices provided in Eq. (2.33). By virtue of the involution property of Eq. (2.34), all stacking configurations of the row transfer matrices are equivalent. Let us emphasise that the resulting vertex model is not the common six-vertex model, but can be regarded instead as a fused variant thereof. Indeed, in general there is no upper bound on the ‘number of vertices’, as a generic equilibrium state requires arbitrary large spin representation labels. The contribution $\boldsymbol{\varrho}_{\vec{h}}$ appears as an additional factor. It plays the role of a boundary twist when the cylinder is wrapped to a torus (i.e. when the ensemble is traced over).

6 The mirror system

The vertex-model regularisation of a generic equilibrium ensemble achieved in the previous section allows for a description of equilibrium states in manner which is complementary to the TBA analysis of Section 4. Here the free energy density is

[TABLE]

which can be viewed as an inhomogeneous vertical iteration of row transfer matrices $\mathbf{T}_{k}$. Alternatively this can be re-expressed as a horizontal iteration of an inhomogeneous column transfer matrix, as illustrated in Figure 3. As the two-dimensional vertex model is homogeneous in the horizontal (i.e physical) direction, the iteration of the column transfer matrix is homogeneous, with the consequence that its dominant eigenvalue yields the free energy density. In this section we analyse this process in detail. We then pay particular attention to the re-emergence of the thermodynamic $\mathscr{Y}$-system, Eq. (4.11), in the large-$N$ limit.

Given a two-dimensional vertex model, we introduce the corresponding ‘mirror system’ as the combination of the auxiliary spaces: $\mathcal{H}_{\rm M}\cong\bigotimes_{k}\mathcal{V}^{\otimes 2N_{k}}_{k}$. The fundamental column transfer matrix 222The canonical version of the column transfer matrix appears under several different names in the literature, including the ‘virtual-space’ [21], ‘crossing’ [23], ‘column-to-column’ [58], ‘quantum’ [27] and ‘Matsubara’ [59] transfer matrix. acts on the mirror system as follows

[TABLE]

with spectral parameter $u$, and the trace taken over the common fundamental space $\mathcal{V}_{1}$ of all the $\mathbf{L}_{k,1}$ inherited from a lattice site of the original spin chain. That is, in switching to the column transfer matrix the role of physical and auxiliary spaces is interchanged. We emphasise that unlike the row (physical) transfer matrices defined in Eq. (2.33), the column transfer matrix is inhomogeneous in terms of both impurities and representation labels, and has a boundary twist encoded by the factor $\mathbf{D}(\vec{h})$ which also acts on $\mathcal{V}_{1}$. For notational convenience we proceed by suppressing explicit dependence of $\boldsymbol{\mathcal{T}}_{1}(u)$ on impurities $\theta_{k,i}$, twist $h$ and dimension $N$.

The fundamental column transfer matrix $\boldsymbol{\mathcal{T}}_{1}(u)$ allows for computation of the free energy density through its dominant (i.e. largest) eigenvalue according to the prescription

[TABLE]

where $\widehat{\mathcal{T}}_{1}(u)$ denotes the dominant eigenvalue. Care must be taken when interchanging the thermodynamic $L\to\infty$ and the scaling $N\to\infty$ limits (cf. [21, 22]), and in our analysis we justify this step by establishing full consistency, as summarised in Figure 1.

One way to proceed to determine $\widehat{\mathcal{T}}_{1}(0)$ is to obtain the Bethe equations which diagonalise $\boldsymbol{\mathcal{T}}_{1}(u)$ as a route towards obtaining their spectrum $\mathcal{T}_{1}(u)$, and thereby its dominant eigenvalue. Here we instead take a different route and employ the framework of functional Bethe ansatz [60, 24, 39, 35], which conveniently allows for the mirror-system Bethe equations to be completely bypassed.

With aid of the fusion rule for the Lax matrices, we proceed by embedding $\boldsymbol{\mathcal{T}}_{1}(u)$ into the infinite fusion hierarchy of commuting transfer matrices $\boldsymbol{\mathcal{T}}_{j}(u)$. To this end, we extend the family of Lax matrices $\mathbf{L}_{k,1}(v,u)$ to $\mathbf{L}_{k,j}(v,u):\mathcal{V}_{k}\otimes\mathcal{V}_{j}\to\mathcal{V}_{k}\otimes\mathcal{V}_{j}$. These are readily obtained through fusion [38, 34, 61, 62]

[TABLE]

where each of $j$ copies $\mathbf{L}^{(a)}_{k,1}$ acts on $\mathcal{V}_{k}\otimes\mathcal{V}^{(a)}_{1}$, $\boldsymbol{\mathcal{P}}_{j}:\bigotimes_{a=1}^{j}\mathcal{V}^{(a)}_{1}\mapsto\mathcal{V}_{j}$ is the totally-symmetric projection operator,

[TABLE]

with $\vec{\mathbf{S}}=\sum_{a=1}^{j}\vec{\mathbf{S}}^{(a)}$, and $\mathcal{N}_{k,j}(v,u)=\prod_{a=1}^{j-k}(v-u+(k-j-1+2a)\tfrac{{\rm i}}{2})$ is the common polynomial produced by fusion. Correspondingly, the higher-spin column transfer matrices are given explicitly as

[TABLE]

where here $\mathbf{D}(\vec{h})$ acts on $\mathcal{V}_{j}$.

In particular, analogously to the physical $T$-functions Eq. (2.35), the hierarchy of column $\mathcal{T}$-functions $\mathcal{T}_{j}(u)$ satisfy the Hirota equation

[TABLE]

with initial conditions $\mathcal{T}_{-1}\equiv 0$, $\mathcal{T}_{0}=1$, and the scalar potentials

[TABLE]

with $\Psi_{j,k}$ given by Eq. (2.29). By construction, $\mathcal{T}_{j}(u)$ are meromorphic functions of rapidity variable $u$. The boundary twist manifests itself as non-trivial large-$u$ asymptotics

[TABLE]

parametrised by $h=|\vec{h}|$ through the $\mathfrak{su}(2)$ characters $\chi_{j}(h)$ obeying the ‘classical limit’ of the Hirota equation

[TABLE]

The associated ‘gauge’-invariant $\mathcal{Y}$-functions

[TABLE]

are again meromorphic functions, obeying the canonical $\mathcal{Y}$-system relations [42]

[TABLE]

with large-$u$ asymptotics

[TABLE]

The above functional relations, Eqs. (6.7), (6.11) and (6.12), are valid for the full spectrum of $\boldsymbol{\mathcal{T}}_{j}(u)$ at finite $N$. To specify a particular eigenvalue, it is necessary to identify their analytic data in the physical strip $\mathcal{P}$. This is made transparent by transforming the functional relations to integral equations on the real rapidity axis. To perform this step, it is useful to introduce adapted $\mathcal{T}$-functions,

[TABLE]

obtained by factoring out their (possible) zeros $\mathlcal{t}^{z}_{j,i}$ in $\mathcal{P}$, so that $\log\mathcal{T}^{\rm adp}_{j}(u)$ are analytic on $\mathcal{P}$ with constant large-$u$ asymptotics. The absence of poles of $\mathcal{T}_{j}(u)$ in $\mathcal{P}$ can be seen from the denominator in Eq. (6.6). Now manipulating Eq. (6.11), by taking the logarithm and convolving with $s$, we obtain

[TABLE]

We similarly transform the $\mathcal{Y}$-system relations, Eq. (6.12), to the real axis. The procedure is analogous to that described in Section 4, specialised to meromorphic data only. Introducing adapted $\mathcal{Y}$-functions

[TABLE]

by again multiplying out all (possible) zeros $\mathlcal{y}^{\rm z}_{j,a}$ and poles $\mathlcal{y}^{\rm p}_{j,b}$ on the physical the strip $\mathcal{P}$, so that $\log\mathcal{Y}^{\rm adp}_{j}(u)$ are analytic on $\mathcal{P}$ with constant large-$u$ asymptotics. Substituting this into Eq. (6.12), taking the logarithm and convolving with $s$, we end up with the coupled integral equations

[TABLE]

with source terms

[TABLE]

The spectrum of the $\mathcal{T}$-functions is thus given by Eq. (6.15), where the $\mathcal{Y}$-functions obey Eq. (6.17), and eigenvalues are specified by the locations of the analytic data $\mathlcal{t}^{\rm z}_{j,a}$, $\mathlcal{y}^{\rm z}_{j,a}$, $\mathlcal{y}^{\rm p}_{j,a}$. To determine this state-dependent data, one has in general to resort to the Bethe equations. The $\mathcal{Y}$-functions however inherit certain zeros and poles directly from the impurity data through the scalar potentials in Eq. (6.11). To identify these we need the following properties which can verified from the definition of $\Psi_{j,k}$ in Eq. (2.29): if $w$ belongs to the lower-half of the strip $\mathcal{P}$ then $\Psi_{j,k}(u-w)$ has a zero at $u=w+\tfrac{{\rm i}}{2}$ if and only if $j=k$, and has no poles in $\mathcal{P}$, while if $w$ belongs to the upper-half of the strip $\mathcal{P}$ then $\Psi_{j,k}(u-w)$ has a pole at $u=w-\tfrac{{\rm i}}{2}$ if and only if $j=k$, and has no zeros in $\mathcal{P}$. Then from Eq. (6.8) one deduces that zeros of the $\mathcal{Y}$-functions are located at

[TABLE]

along with their complex conjugates, and we remind that $\xi^{(-)}_{k}<0$. Correspondingly, the poles of the $\mathcal{Y}$-functions appear at

[TABLE]

along with their complex conjugates, with $\xi^{(+)}_{k}>0$.

The eigenstate for which these are the only zeros and poles of the $\mathcal{Y}$-functions inside $\mathcal{P}$ is a distinguished state of the mirror system. In particular, the corresponding $\mathcal{T}$-functions $\mathcal{T}_{j}(u)$ are analytic and free of zeros in the strip $\mathcal{P}$ for this state, a further consequence of Eq. (6.11). As any zero of a $\mathcal{T}$-function in $\mathcal{P}$ yields a negative definite contribution to Eq. (6.15), i.e.

[TABLE]

it is natural to anticipate that this state gives the dominant eigenvalue determining the free energy density through Eq. (6.3). In the following we establish this assertion by taking the large-$N$ limit, and demonstrating that it reproduces precisely the TBA description of Section 4.

6.1 The scaling limit

In the preceding analysis we considered the mirror system at finite $N$. At this level, there is no strict distinction between the impurities encoding the node and analytic data from Eq. (5.27). Now we take the large-$N$ scaling limit in which this distinction becomes manifest. We focus on the distinguished state identified above for which all the $\mathcal{T}$-functions are analytic and free of zeros in the physical strip $\mathcal{P}$.

The task is to inspect the large-$N$ scaling limit for both Eq. (6.15) and Eq. (6.17). For Eq. (6.15) the non-trivial $N$-dependent term is \log\big{(}\varphi_{j}\bar{\varphi}_{j}\big{)}. For clarity, we consider below the contributions coming from the node and analytic data separately, that is we split

[TABLE]

The node data $\lambda_{j}$ involves a sum over pairs of impurities

[TABLE]

yielding

[TABLE]

Individual contributions, which can be deduced with aid of Eq. (2.29), are of the form

[TABLE]

In the large-$N$ limit, the net contribution of Eq. (6.24) yields

[TABLE]

upon removing the large-rapidity cut-off $\Lambda^{\infty}$ and converting sums to convolution-type integrals. In the process, the impurities re-condense in accordance with Eq. (5.7).

The contribution from the analytic data stemming from zeros $\{\omega^{\rm z}_{k,i}+{\rm i}\epsilon-\tfrac{{\rm i}}{2},\bar{\omega}^{\rm z}_{k,i}-{\rm i}\epsilon+\tfrac{{\rm i}}{2}\}$ and poles $\{\omega^{\rm p}_{k,i}+{\rm i}\epsilon+\tfrac{{\rm i}}{2},\bar{\omega}^{\rm p}_{k,i}-{\rm i}\epsilon-\tfrac{{\rm i}}{2}\}$ inside $\mathcal{P}$, reads

[TABLE]

Employing the identity $\log\Psi^{\pm}_{j,k}=\pm(1+K)_{j,k}\star\log\tau$ which follows from Eqs. (2.23) and (2.26), this becomes

[TABLE]

with $\zeta^{\rm reg}_{j}$ given by Eq. (5.23). Here the infinitesimal regulator $\epsilon$ can be safely dropped. As the degrees of the Padé approximations tend to infinity with $N$, we then obtain

[TABLE]

Thus combining the contributions from the node and analytic data we recover the TBA chemical potentials (4.16),

[TABLE]

Next we take the large-$N$ scaling limit of Eq. (6.17). Here we must examine the source terms $\mathlcal{d}_{j}$ given by Eq. (6.18). The node data consists of pairs of zeros and poles, cf. Eq. (6.19),

[TABLE]

each contributing to $\mathlcal{d}_{j}$ a term

[TABLE]

Since there are $n_{k}$ terms, in the $n_{k}\to\infty$ limit we are left with

[TABLE]

employing Eq. (2.23). The full contribution to $\mathlcal{d}_{j}$ in the large-$N$ scaling limit then retrieves the node data, that is $s\star\lambda_{j}$. The contribution from the regularised analytic data, which is for finite $N$ straightforwardly given by $\zeta^{\rm reg}_{j}$, and in the large-$N$ limit converges to $\zeta_{j}$. In conclusion,

[TABLE]

are precisely the source terms (4.15) of the local form of the TBA equations Eq. (4.7).

Free energy.

Having taken the continuum scaling limit, we can now determine the eigenvalues of $\mathcal{T}$-functions $\mathcal{T}_{j}(u)$ encoding the distinguished state. In particular, employing Eq. (6.30) in Eq. (6.15) we obtain

[TABLE]

Finally, substituting this into Eq. (6.3) as the dominant eigenvalue we re-obtain precisely the expression for the free energy given by Eq. (4.6). As the large-$N$ scaling limit of the canonical $\mathcal{Y}$-functions $\mathcal{Y}_{j}$ satisfy the local TBA equations Eq. (4.7) due to Eq. (6.34), we have fully recovered the TBA description of Section 4.

In conclusion,

for a general equilibrium state of the Heisenberg spin-$1/2$ chain, the thermodynamic $\mathscr{Y}$-functions coincide with the large-$N$ scaling limit of the distinguished canonical $\mathcal{Y}$-functions associated with the dominant eigenvalue of the mirror system.

7 Discussion

In the preceding sections we considered generic equilibrium states and developed the dual approaches of TBA and lattice regularisation, as summarised in Figure 1. The central role in establishing their equivalence is played by the mirror system, which provides a compact regularisation of a given equilibrium state. The mirror system exhibits a ‘universal’ canonical structure, seen here in Eqs. (6.7)-(6.13), which is deep-rooted in the fusion properties of the underlying quantum symmetry algebra. One of the key findings of this work is however that in thermodynamic/large-$N$ limit the canonical structure is superseded by an emergent equilibrium landscape, encoded in non-trivial node terms $\lambda_{j}$ and generically non- meromorphic thermodynamic $\mathscr{Y}$-functions involving branch points inside $\mathcal{P}$. There are several discernible aspects which deserve a brief discussion.

Breakdown of the canonical $\mathcal{Y}$-system.

The thermodynamic $\mathscr{Y}$-system generically contains state-dependent node terms $\lambda_{j}$, see Eq. (4.11). This non-canonical form, which is seemingly conflict with the expected universality of the canonical $\mathcal{Y}$-system [42], has been previously observed the context of quantum quenches where several explicit examples have been identified [7, 63, 16]. Here we clarify this aspect.

The $\mathcal{Y}$-functions at finite $N$ obey the canonical $\mathcal{Y}$-system, Eq. (6.12), with $n^{(\pm)}_{j}$ zeros and poles stemming from the regularised $\lambda_{j}$ residing are finite distance $\xi^{(\pm)}_{j}=|\Lambda^{(\pm)}_{j}|/2\pi n^{(\pm)}_{j}$ from the boundary of $\mathcal{P}$. A subtlety however arises in the thermodynamic scaling limit $N\to\infty$, where $\xi^{\pm}_{j}$ invariably become smaller than the positive infinitesimal regulator $\epsilon$ which controls the width of the physical strip $\mathcal{P}$ as in Eq. (2.22). 333We remind that in distinction to $N$, $\epsilon$ is not an adjustable parameter of the regularisation scheme. In this event, a subset of the analytic data leaves the nullspace of $s^{-1}$ by escaping from the integration contours $\mathcal{C}_{\pm}$ and finally collapsing onto the boundary of $\mathcal{P}$ at ${\rm Im}(u)=\tfrac{1}{2}$. This piece of information emerges through the re-condensed $\lambda_{j}$ as in Eq. (6.34), rendering the thermodynamic $\mathscr{Y}$-system of the non-canonical form. Indeed, a simple example of a non-canonical $\mathscr{Y}$-system is provided by the canonical Gibbs equilibrium state, specified by $\lambda_{j}(u)=\pi J\,\beta\,\delta_{j,1}\delta(u)$ and $\zeta_{j}(u)=0$. The corresponding TBA source term thus reads $d_{j}=s\star\lambda_{j}=\pi J\,\beta\,\delta_{j,1}s(u)$, in agreement with [26].

Non-meromorphic $\mathscr{Y}$-functions.

As observed in the TBA analysis, in general the thermodynamic $\mathscr{Y}$-functions are non-meromorphic complex functions. This is to be contrasted with the canonical $\mathcal{Y}$-functions of the mirror system describing regularised macrostates which are manifestly meromorphic on $\mathcal{P}$. The non-meromorphic structure appears in the large-$N$ scaling limit, when the analytic data from the interior of $\mathcal{P}$ form macroscopic condensates. A degree-$N$ Padé approximation of a generic non-meromorphic $\mathscr{Y}$-function in terms of a canonical (i.e. meromoprhic) $\mathcal{Y}$-function provides a particular algorithm of information compression.

Non-canonical asymptotics.

A third subtlety concerns the large-$u$ asymptotics of the thermodynamic $\mathscr{Y}$-functions. At finite $N$ the canonical $\mathcal{Y}$-functions possess a one-parameter family of large-$u$ asymptotics parametrised by $h$, cf. Eq. (6.13), expressed through the solution to the classical limit of the Hirota equation Eq. (6.10). In the large-$N$ limit however, upon removing the regulator $\Lambda^{\infty}\to\infty$, these asymptotics may decouple, with the $\mathscr{Y}$-functions acquiring non-canonical large-$u$ asymptotics provided the large-$u$ asymptotics of the chemical potentials $\mu_{j}(u)$ are non-trivial. A particular instance of this is an infinite-parametric family of ‘dispersionless’ states, that is the class of states characterised by constant $\mu_{j}$.

8 Conclusion

In this work we have undertaken a comprehensive study of the equilibrium landscape of the isotropic Heisenberg spin-$1/2$ chain. We have developed a robust and unified framework which encompasses both the Thermodynamic Bethe Ansatz and the two-dimensional vertex-model regularisation approaches to thermodynamics. In particular we have explained the emergence of a splitting of the chemical potentials $\mu_{j}$ into two contributions: the node data $\lambda_{j}$ which determine the thermodynamic $\mathscr{Y}$-system, and the analytic data $\zeta_{j}$ which encode the branch points of the $\mathscr{Y}$-functions in the physical strip $\mathcal{P}$. These characterise equilibrium states in distinct ways, endowing the equilibrium landscape with a structure that is deserving of further exploration.

There are several novel features of the work worth highlighting. Firstly, we express the generic density matrix Eq. (3.1) in a form which reflects the underlying $\mathfrak{su}(2)$ symmetry of the model, which clarifies how the polarisation of the mode operators $\boldsymbol{\rho}_{j}$ is set. In our TBA analysis, we demonstrate on general grounds that the equilibrium $\mathscr{Y}$-functions are generically non-meromorphic in the physical strip $\mathcal{P}$. Our lattice regularisation of a generic ensemble treats the node and analytic data separately. For the node data we invoke a discretisation of the $\lambda_{j}(u)$, and achieve a variant of ‘Trotterisation’ without actually appealing to the Suzuki-Trotter formula. For the analytic data we develop a contour integration procedure, and introduce Padé approximants to regularise generic branch points of $\mathscr{Y}$-functions. In Section 6 we reconnect the TBA and vertex model approaches by identifying a distinguished eigenvalue of the mirror system, and use it to demonstrate complete equivalence of the descriptions. The interconnected nature of our analysis is summarised in Figure 1.

A technical point emerging from our study is an explanation of the breakdown of the canonical $\mathcal{Y}$-system, i.e. the emergence of the thermodynamic $\mathscr{Y}$-system. Put simply, this is due to competition between the infinitesimal regulator $\epsilon\equiv 0^{+}$ which is tied to locality of the $\boldsymbol{\rho}_{j}$, and the regularisation parameter $N$ which is finite for any vertex model approximation of the ensemble and diverges in the scaling limit. At finite $N$ all the poles and zeros of the $\mathcal{Y}_{j}$ lie on the physical strip $\mathcal{P}$ of Eq. (2.22) and the corresponding $\mathcal{Y}$-system is canonical, Eq. (6.12). The poles and zeros resulting from the node data $\lambda_{j}$ however approach the boundary of the strip as $1/N$, so that in the large-$N$ limit they escape $\mathcal{P}$, resulting in the appearance of node terms in the thermodynamic $\mathscr{Y}$-system, Eq. (4.11).

Throughout the work we have adopted a universal language, which should facilitate extensions to other models solvable through the Bethe ansatz framework. The simplest generalisation of the model considered here is its uniaxial anisotropic deformation [19]. The extension to the ‘hyperbolic’ (easy-axis) regime, with quantum deformation parameter $q\in\mathbb{R}$, is straightforward and readily follows from an analytic continuation of the local degrees of freedom (i.e. the scattering matrix) of the isotropic model considered here. In distinction, in the ‘trigonometric’ (easy-plane) regime with the deformation parameter $q=e^{{\rm i}\gamma}$ at the roots-of-unity $\gamma=m/\ell\in\mathbb{R}$, the total number of local degrees of freedom depends quite intricately on the values of co-prime integers $m$ and $\ell$, see [64, 20, 40]. To implement our approach, the compete set of unitary local charges is required, already identified previously in [7]. The closure of the functional hierarchy and the string-charge duality requires here an extra family of non-unitary charges [65, 44, 66].

Going forward, an important open question is whether there exist physically discernible features associated with the structures identified here. Addressing this will require analysing the correlation functions of local observables [59, 67, 68, 69, 47]. We anticipate that it will be interesting to examine the splitting between node and analytic data in this context, and hope that the framework we provide offers a solid foundation upon which to proceed.

Acknowledgements.

We thank O. Gamayun for valuable remarks on the manuscript. E. I. acknowledges support by VENI grant number 680-47-454 by the Netherlands Organisation for Scientific Research (NWO). E. Q. acknowledges support from ANR IDTODQG project grant ANR-16-CE91-0009 of the French Agence Nationale de la Recherche, the Foundation for Fundamental Research on Matter (FOM), the Netherlands Organization for Scientific Research (NWO), and the European Research Council under ERC Advanced grant 743032 DYNAMINT.