Quasilocal charges in integrable lattice systems

We review recent progress in understanding the notion of locality in integrable quantum lattice systems. The central concept are the so-called quasilocal conserved quantities, which go beyond the standard perception of locality. Two systematic procedures to rigorously construct families of quasilocal conserved operators based on quantum transfer matrices are outlined, specializing on anisotropic Heisenberg XXZ spin-1/2 chain. Quasilocal conserved operators stem from two distinct classes of representations of the auxiliary space algebra, comprised of unitary (compact) representations, which can be naturally linked to the fusion algebra and quasiparticle content of the model, and non-unitary (non-compact) representations giving rise to charges, manifestly orthogonal to the unitary ones. Various condensed matter applications in which quasilocal conservation laws play an essential role are presented, with special emphasis on their implications for anomalous transport properties (finite Drude weight) and relaxation to non-thermal steady states in the quantum quench scenario.


Introduction
Local conservation laws are amongst the most important fundamental concepts in theoretical physics. In generic systems these usually comprise of energy, momentum, particle number, etc., and correspond to Noether charges connected to rather obvious physical symmetries. On the other hand, in systems which are exactly solvable, or integrable, the number of conservation laws and the corresponding conserved charges can be much larger and the underlying symmetries sometimes quite hidden. According to a widespread belief, integrability should provide us with a 1-to-1 correspondence between conserved charges and physical degrees of freedom. However, such a definition is only really applicable -or unambiguous -in classical deterministic (Hamiltonian) systems with a finite number of degrees of freedom where it amounts to the historical, Liouville-Arnold integrability.
Interacting quantum systems, where local degrees of freedom (quantum spins, fermions, or bosons) are arranged in a regular 1D lattice, are typically considered integrable in one of the following cases: Firstly, there may exist a canonical (Bogoliubov) transformation which maps the local degrees of freedom to non-interacting quasiparticles. Such is, for example, the situation with quantum transverse field Ising model, or XY spin-1/2 chain [1]. These systems, which are reducible to a single particle picture and are often referred to as quasi-free, shall not be of interest in this article, even though they allow for an illustration of some non-trivial many-body phenomena, such as area laws for entanglement [2]. Secondly, there exist systems exhibiting genuine interparticle interaction whose dynamics is representable in terms of quasi-particles which undergo non-diffractive scattering without particle production. A central feature in such a case is factorizability of an arbitrary multi-particle scattering process in terms of subsequent 2-particle scattering events, mathematically phrased in the form of the celebrated Yang-Baxter (or star-triangle) equation. One of the most remarkable physical consequences of that mechanism is the emergence of a macroscopic number of local integrals of motion (conservation laws). One of these charges, usually the first one in the series, is considered as the Hamiltonian (with local interactions). Here locality means that the densities of these charges act non-trivially only on a finite number of adjacent lattice sites. Integrability in the sense of Yang and Baxter, which is universally understood within the framework of algebraic structures known as quantum groups [3][4][5][6], is perhaps the most general widely acceptable definition of integrability known to date. Besides defining and describing integrability in closed quantum many-body systems in 1D [7], it also covers 2D equilibrium classical statistical systems [8], nonequilibrium classical driven diffusive 1D systems [9], as well as classical Hamiltonian systems [10,11], and since more recently, also integrable nonequilibrium steady states of open quantum interacting systems [12].
The fact that certain integrable many-body systems can already be routinely controlled in a concrete experimental setup also underlies a remarkable degree of structural stability for some of their dynamical properties with respect to model imperfections (perturbations), in spite of the fact that strict integrability technically requires precise (or fine-tuned) cancellations of most of generically allowed processes. This may hint to an existence of a yet undisclosed quantum analogy of KAM (Kolmogorov-Arnold-Moser) scenario [34]. To our opinion this is one of the potentially most exciting problems for future research [35,36].
As discussed above, Yang-Baxter integrability for a lattice system with N sites, guarantees a macroscopic number ∝ N of local conservation laws and the corresponding local currents. By a local conservation law one understands an operator-valued continuity equation, involving a charge and a current density being operators supported on a finite number of, say n N , physical sites. The summation of the local charge density over the whole volume of N sites then defines an extensive local conserved charge of an integrable model. One might wonder whether such local conserved charges represent a complete set, meaning that any extensive conserved operator which scales linearly with N can be represented as a linear combination of these local charges. Some formal completeness results for specific models have been put forward a while ago [37], and one might have been tempted to conclude that local charges (derived from fundamental Yang-Baxter transfer matrix) are all the conserved operators needed to understand local physics. However, certain unconventional phenomena discovered later in studies of paradigmatic examples of interacting integrable systems gave, in spite of a missing formal understanding, quite the opposite indications. Firstly, it has been discovered [38,39] that the spin Drude weight in the integrable anisotropic Heisenberg chains (XXZ model) is finite at finite temperature, despite the fact that contributions of all hitherto known local charges to spin current were zero. In more recent works it has been found [40,41] that a Generalized Gibbs Ensemble (GGE) formed of the same standard set of local conserved charges fails to describe thermalization after a quantum quench in the gapped XXZ model. These results hinted at the existence of additional effectively local conserved charges linearly independent from the strictly local ones. One should note that in studies of infinite quantum (and even classical) lattice systems, extensive observables form a vector space rather than the full algebra, so it is the linear independence and not functional independence that matters.
The first progress along the above lines came, unexpectedly, with the solution of an open XXZ model [42] driven out of equilibrium with effective magnetic (particle) reservoirs at the boundary formulated in terms of Lindblad master equation. The steady state solution in the perturbative (weak-coupling) regime turned out to be tightly related to a novel effectively local (or quasilocal) conservation law which in turn explained the controversial problem of the ballistic conductivity by providing a rigorous nontrivial lower bound on the spin Drude weight and thus confirmed previous results of several numerical studies [43][44][45][46] and bosonization techniques [47,48]. In a subsequent study [49], a connection to certain non-standard solutions to Yang-Baxter equation has been uncovered, permitting a systematic construction of a large set of quasilocal conservation laws directly from commuting transfer matrices associated to complexspin (non-unitarty) representations and yielding a further improved Mazur bound on the Drude weight. Generalizations of the results to periodic boundary conditions were simultaneously obtained in Refs. [50,51]. A distinguished property of these socalled 'non-unitary' quasilocal charges is that they do not exhibit the spin-reversal invariance of the XXZ Hamiltonian and hence may have a nonvanishing overlap with observables which are odd with respect to spin reversal, such as the spin current. Very recently, even more exotic non-unitary quasilocal charges have been discovered where even the particle conservation (U (1)-symmetry) is broken [52]. Similar constructions of quasilocal charges and consequent Drude weight bounds can be performed also in other gapless integrable quantum spin models, for example in spin-1 Fatteev-Zamolodchikov chain [53]. We should remark, however, that it is the compactness of q-deformation rather than masslessness of the elementary excitations which plays the essential role in the construction of current carrying quasilocal charges which break the parity symmetry of the model (e.g., spin reversal). This observation should make it possible to extend these concepts to massive integrable models like the sine-Gordon theory.
In spite of all rather profound implications mentioned above, the family of nonunitary quasilocal conserved operators could not offer the answer to the puzzling findings of Refs. [40,41,54] which cast doubts on the applicability of the concept of a Generalized Gibbs Ensemble which was vividly debated about at the same time. In particular, it became clear that in a generic case the GGE has to be appropriately extended by incorporating quasilocal conservation laws which are viable for the whole range of anisotropies, invariant under spin-reversal transformation (i.e., of even parity), but still distinct from the canonical ones obtained from expanding the fundamental transfer matrix. Such quasilocal charges have been constructed (for the isotropic case) in Ref. [55], invoking transfer matrices built from unitary but non-fundamental spin representations of the auxiliary spin. Soon after, a study [56] confirmed that those charges exactly explain the GGE conundrum.
Outline. The present review article aims at a coherent and pedagogical (i.e. nontechnical) introduction to the notion of quasilocal conserved charges and various physical applications in which they take the center stage. As the focus is primarily to elucidate the main ideas and their interrelations, a reader seeking for a more detailed and rigorous exposition is referred to the cited literature. Sec. 2 consists of a minimal technical background for getting familiar with the main concepts presented in this article. Sec. 3 is devoted to the construction of what we call 'unitary' quasilocal charges, namely conservation laws arising from the unitary representations of an underlying symmetry group. In Sec. 4 a more intricate case of 'non-unitary' quasilocal charges which break the spin reversal (or, in general, some other Z 2 parity) symmetry is presented. Sec. 5 is dedicated to the exposition of several physical applications: Sec. 5.1 discusses rigorous Mazur bounds on the spin Drude weight. Sec. 5.2 makes a link to quantum quenches from spin-reversal symmetric initial states and highlights the duality between the spectra of quasilocal charges and Bethe root distributions which describe bound states in the formalism of the Thermodynamic Bethe Ansatz. Sec. 5.3 illustrates the connection to integrable nonequilibrium steady states of boundary-driven quantum master (Lindblad) equations. In this review, all the concepts are presented explicitly on a concrete example of the XXZ chain and the associated U q (sl(2)) quantum symmetry. We conclude in Sec. 6 where certain possible generalizations to other integrable models and some questions which enter in the broader context are briefly discussed.

Prerequisites
In this section we introduce the framework and technical tools that shall be used in our paper. In the Sec. 2.1 we introduce the concepts of quantum spin systems on the lattice and the corresponding operator (C * ) algebra, and define the notions of locality, extensivity, pseudolocality and quasilocality. In Sec. 2.2 we define the main concepts of Yang-Baxter integrability: R-matrices, Lax matrices, transfer matrices, and fusion hierarchies which allow one to build unitary representations of these objects from the fundamental one. These concepts enable us to reformulate Bethe's original 'coordinate ansatz' [57] in an entirely algebraic language, a technique which is nowadays typically referred to as the quantum inverse scattering method or the algebraic Bethe ansatz [7,58,59].
The point of our review is to show that one can develop a new perspective on nonequilibrium quantum physics by combining the concepts from Yang-Baxter integrability with the notions of pseudo-and quasilocality of extended quantum lattice systems.

Pseudolocal and quasilocal operators over quantum lattices
The main theme of this article are conserved charges of integrable lattice models which comply with a certain weaker version of locality. As such, they extend beyond the orthodox concept of local charges, derived from logarithmic derivatives of the fundamental transfer matrix [7,[58][59][60], and exhibit physical relevance for computing time-averaged values of dynamical response functions.
Since we are only concerned with integrable systems, we can limit our discussion to a one-dimensional lattice Λ = Z, although the concepts of this section can be readily extended to a D-dimensional lattice Λ = Z D . The total Hilbert space, formed by a tensor product of d−dimensional single-site Hilbert spaces, will be denoted by H. The Hilbert space of a lattice subinterval between sites x and x , x ≤ x , will be denoted by H [x,x ] ⊂ H and the corresponding operator subalgebra by A [x,x ] . The entire quasilocal C * operator algebra A is obtained as the limit of a sequence {A [−n,n] ; n = 1, 2, 3 . . .}, closed in the operator norm topology [61]. We shall refer to an observable represented by an operator a ∈ A as local, if it acts nontrivially only on a finite subinterval [x, x ], The smallest such interval is referred to as the support of a, and its length r = x − x + 1, as the order of locality. Denoting by Tr [x,x ] the trace over H [x,x ] , one defines the tracial state ω 0 as and extends it over an entire A by continuity (of ω 0 ). The tracial state can be interpreted as the infinite temperature Gibbs state, satisfying ω 0 (ab) = ω 0 (ba) and having the strongest clustering property, namely being separable: ω 0 (ab) = ω 0 (a)ω 0 (b) for any pair of local observables a, b with disjoint supports. We define the Hilbert-Schmidt (HS) inner product as and denote the corresponding HS norm ‡ by a HS ≡ (a, a). The latter satisfies the standard Cauchy-Schwartz inequality and a mixed inequality in relation to the operator norm • , Equipped with these structures we can define an orthonormal basis of local observables. A choice of an on-site basis such that (σ α x , σ α x ) = δ α,α , induces the HS orthonormal basis of algebra A [x,x ] consisting of elements of the form For example, in case of 2-dimensional local Hilbert space, σ α≥1 x are just the Pauli matrices, while for 3-dimensional local space they are the Gell-Mann matrices, etc. In all cases we choose σ 0 We furthermore define a lattice shift automorphism byŜ y (a [x,x ] ) = a [x+y,x +y] and associate to each element a ∈ A a translationally invariant sum which represents an extensive observable of a translationally invariant infinite quantum spin chain. Note that A is not an element of quasilocal algebra A, but the above sum can ‡ Note that, strictly speaking, (a, b) and a HS become a proper HS product and HS vector norm, respectively, only after one takes the identity operator 1 out of the algebra A. Otherwise they yield the HS product and HS norm of the corresponding 'nearest' traceless observables. In other words, any operator of the form c1, c ∈ C, has 'zero length'.
still be attributed a precise mathematical meaning as a sequence of operators {A (N ) } acting on finite lattices of increasing lengths N . For example, the Hamiltonian of locally interacting translationally invariant models, as well as other strictly local charges, are precisely of such form. In this sense, a local operator a is called a density of an extensive local observable A. The above sequences have the following properties: (i) volume scaling extensivity and (ii) a finite overlap lim N →∞ (b, A (N ) ) = 0 with at least one local operator b (say b = a). In what follows, the upper index N will be left out, since an extensive operator A is always identified with the corresponding sequence. By definition, any operator sequence A, satisfying extensivity (i) given by Eq. (2.7), and the finite overlap criterion (ii), shall be referred to as pseudolocal. This relaxes the constraint on the strict locality of the densities and generalizes the concept in a physically meaningful way. As we shall argue later, pseudolocality of conserved charges is the decisive property responsible for ballistic (or non-ergodic [62][63][64]) scaling of dynamical response functions. Note that if the density a can be written as a sum of mutually orthogonal terms a [1,r] , for which a stronger condition, known as quasilocality [49], a [1,r] HS < Ce −ξr , ξ > 0, (2.9) holds, A is automatically pseudolocal.
Here we have considered lattices with open boundaries. For systems with periodic or twisted boundary conditions, the same concepts can be introduced by making the shift operatorŜ periodic [50].
The definition of pseudolocality and quasilocality can be generalized (see Ref. [65]) to an arbitrary sufficiently strongly clustering state ω (say Gibbs, or generalized Gibbs state, etc.) simply by replacing the HS inner product by with the main conclusion that the set of all pseudolocal observables forms a Hilbert space.

Yang-Baxter relation, quantum transfer matrices, and fusion hierarchies
A distinguished feature of integrable models is an existence of a macroscopic number of conservation laws. They arise as a consequence of an exceptional amount of symmetry which is governed by algebraic structures known as quantum groups [3][4][5][6]. The central element in the story is the so-called quantum R-matrix, an operator acting on a tensor product of a pair of vector spaces, that can be considered as representations V 1 and V 2 of an underlying symmetry algebra, which we here for simplicity assume to be su(2) or its quantum deformation. In addition, R(λ) depends analytically on a spectral parameter λ ∈ C. The cornerstone equation of quantum integrability is obtained by embedding R-matrices into a three-fold tensor product space V 1 ⊗ V 2 ⊗ V 3 , by making use of a suggestive notation R 12 (λ) = R(λ) ⊗ 1, and imposing the requirement where we have omitted the indices of vector spaces on which the operators act trivially. This condition is the celebrated Yang-Baxter equation [8,66,67] (YBE). Physically speaking, YBE expresses equivalence of two distinct sequences of two-particle collisions which, as a consequence, give factorization property of the whole many-particle scattering process [66,68]. What is perhaps even more remarkable is, that such an equivalence automatically generates an infinite number of conserved quantities. The procedure is outlined below. The simplest solution to YBE (2.12) is obtained when the R-matrix acts in two fundamental spin representations V 1/2 ∼ = C 2 , where P is a permutation operator, P |ψ 1 ⊗ |ψ 2 = |ψ 2 ⊗ |ψ 1 . Furthermore, we introduce the Lax operator L(λ) by interpreting one fundamental space of the R-matrix as a local physical spin while the second fundamental space is referred to as an auxiliary space, L 12 (λ) ≡ R 12 (λ) = λ + 2i s 1 · s 2 , or 14) The spin generators fulfil the su(2) algebraic relations, [s + , s − ] = 2s z and [s z , s ± ] = ±s ± , and in terms of the Pauli matrices read s z = 1 2 σ z and s ± = σ ± = 1 2 (σ x ± iσ y ). For clarity of notation, we shall here and below use bold-roman fonts to denote all operators which act nontrivially in auxiliary (non-physical) spaces. From YBE (2.12) it follows that the Lax operator Eq. (2.14) by construction obeys the local fundamental commutation relation (also known as the RLL relation [4,58]) over the auxiliary vector space H a ⊗H a , which can be extended to the entire physical Hilbert space H p ∼ = (C 2 ) ⊗N of the N -spin lattice Here and subsequently we use a compact notation of ⊗ N to denote a 'partial' tensor product, i.e. an operation where the tensor product only affects the physical components, whereas for the auxiliary components ordinary matrix multiplication applies. Finally, by tracing over the auxiliary space of Eq. (2.17) we produce the fundamental transfer matrix T (λ) = Tr a M(λ), (2.18) acting over the spin chain Hilbert space H p . An infinite set of conservation laws is a consequence of commutativity property which follows directly from the definition (2.17) in combination with the YBE (2.12). In fact, by considering higher-dimensional irreducible unitary representations of auxiliary spaces (s > 1/2), one sees that the entire construction also holds for higher-spin transfer operators. These are constructed from Lax operators L s (λ) associated with (2s + 1)dimensional auxiliary spaces H a = V s ∼ = C 2s+1 and satisfy (2.20) 2.2.1. Lax operator for the anisotropic Heisenberg model. In this work we discuss the properties of quasilocal conservation laws in the anisotropic Heisenberg spin-1/2 chain (XXZ model), where, unless otherwise stated, periodic boundary conditions are assumed. Including the anisotropy requires employing a one-parametric deformation of the su(2) symmetry algebra, which formally gives rise to a quantum-deformed (quantized) enveloping algebra U q (sl (2)). The suitable deformation is achieved through the deformation parameter q = exp(η), yielding the Lax operator of the following form (see e.g. [58]) Three regimes are to be distinguished with respect to the anisotropy parameter ∆: • gapped regime, corresponding to anisotropy ∆ = cosh (η) > 1 with η > 0, • gapless regime, corresponding to |∆| < 1, which we shall write as ∆ = cos(η) with q-parameter lying on the unit circle q = exp (iη) for η ∈ (0, π). In this regime, replacement η → −iη and λ → −iλ is needed in (2.22) to restore the notation that is most often used (equivalent to exchanging sin and sinh in Eq. (2.22)), and that is used below.
• isotropic point, ∆ = 1, is obtained from either of the regimes by taking the scaling limit, namely to write the spectral parameter as λ → λη and then take η → 0.

(2.24)
In addition to finite-dimensional unitary representations of U q (sl(2)) algebra, YBE (2.12) in fact admits a much larger class of solutions which pertain to generic complex-spin highest-weight representation V + s , s ∈ C (see e.g. Refs. [69][70][71]). These are of infinite dimension for a generic value of s. For values of deformations corresponding to η = π l/m, with l, m, l < m, being co-prime positive integers -or equivalently, for q being a primitive root of unity -we shall be interested in irreducible finite-dimensional sub-representations V [n + 1] q |n n + 1| , [2s − n] q |n + 1 n| .

(2.25)
Here the state |0 designates the highest-weight vector, alias the 'vacuum', s + s |0 = 0. Highest-weight transfer operators T hw s with s ∈ C are defined according to the same prescription as in Eq. for all distinct spin labels and pairs of spectral parameters λ, µ ∈ C.
The standard set of local charges is generated by an expansion of log T 1 2 (λ) around the so-called shift point, where H (2) ∼ H is the Hamiltonian (2.21). The locality of conserved operators H (k) is manifested in the fact that each H (k) admits an expansion in terms of homogeneous sums of local densities h (k) of order k, i.e.
for any finite length N . Let us now switch the focus to the properties of higher-spin transfer matrices T s and their spectra, which play a vital role in the construction of unitary quasilocal conserved charges. These properties will only be used later in the 'fusion approach' (Sec. 3.2) and for obtaining closed-form results in the quantum quench problem (Sec. 5.2.4).

Quantum Hirota equation.
The quantum Hirota equation [72][73][74][75][76][77], also known as the T -system [78,79], is a bilinear difference equation which takes the form with bar denoting complex conjugation. This relation can be formally understood as the quantized version of Weyl's formula for characters of classical representations [73,80], while physically it represents fusion rules on an underlying algebra in a covariant way.
Higher-spin transfer operators T s represent the canonical solution to the Hirota equation.
In this case, the scalar potentials have to be identified as φ(λ) = T 0 (λ + iη 2 ) and φ(λ) = T 0 (λ − iη 2 ), where T 0 (λ) = (sin (λ)/ sinh (η)) N . There exists some (gauge) freedom in choosing the operators T s , which is the reason for defining their gauge-invariant combinations known as the Y -operators. They are defined through the non-linear transformation where the following compact notation is introduced: in the isotropic case (after applying a scaling limit λ → λη and sending η → 0). We shall write f ± (λ) ≡ f [±1] (λ). The Y -operators obey the Y -system functional relations where the boundary condition Y 0 = 0 is assumed. In this article, Hirota equation appears in two different (but related) contexts: (i) as the fusion relation among higher-spin transfer operators T s which is automatically inherited by their eigenvalues, and (ii) as an analytic closed-form description of certain solutions of equilibrium states which typically arise in the scope of quantum quench applications (cf. Sec. 5.2).
The Hirota equation (2.29), can be understood as a discrete integrable classical system of its own. A central relation in this regard is the Baxter's T Q-equation [81][82][83][84] which represents a discrete second-order difference equation for the fundamental transfer matrix T 1/2 . The operator Q stands for Baxter's Q-operator. We do not derive it here explicitly (see e.g. Ref. [82]), but make use of its spectral representation which will provide the connection to Bethe eigenstates (cf. Eq. (3.34)).
The Q-operator allows us to linearize Eq. (2.29), i.e. enable us to express T s (λ) explicitly as a combination of Q-operators where the scalars are provided by .

(2.34)
Since the T Q-equation (2.32) is of a second order, it admits two (linearly) independent solutions, Q and Q, whose independence requires the Wronksian determinant to be non-degenerate, By virtue of commutativity of T s (λ) and Q(µ), for all s ∈ 1 2 Z + , λ, µ ∈ C, all previously stated identities can be taken at the level of their eigenvalues. To distinguish commuting operators from their eigenvalues, we write the latter with the calligraphic font. Bethe roots λ j are by definition zeros of eigenvalues of Q, i.e. solutions of Q(λ) = 0. Bethe ansatz equations can be obtained algebraically by eliminating Q through the combination of Eq. (2.32) and the Wronskian condition (2.35), yielding an equation for the eigenvalues Similarly, Eq. (2.33) turns out to be useful in studying the large-N limit spectra of the transfer operators T s . We shall exploit this trick later on in Sec. 5.2.

Quasilocal charges from unitary representations
In this section we construct quasilocal charges from half-integer representations of the auxiliary algebra (2.24), extending the standard family of local charges. In the first part we formulate the pseudolocality condition in terms of auxiliary transfer matrices and subsequently demonstrate its equivalence to the inversion identity. Furthermore, the construction allows us to obtain a representation of conserved charges, which is useful for computation of their norms and subsequently performing orthogonalization procedure. Subsequently we present an alternative approach to obtain the inversion identity by resorting to previously discussed Hirota equation. The latter enables us to identify quasilocal charges which pertain to the gapless regime.

Auxilliary transfer matrix approach
Initially, we consider the |∆| ≥ 1 regime of the XXZ model and show that an infinite tower of conserved operators generated from the higher-spin transfer operators T s are indeed quasilocal conserved charges. The sketch of the proof given below is based on establishing the inversion formula derived in Ref. [55], which allows for an alternative representation (or definition) of the charges Eq. (3.1) in a more convenient product form Subsequently we will adopt Eq. (3.3) as a working definition when proving quasilocality property of operators X s (λ). Initially, we shall not rely on the apparatus of integrability but rather employ a direct technique using auxiliary transfer matrices. By doubling the auxiliary space the operator product on the left hand-side of Eq. (3.2) can be represented as where the trace takes place in V s ⊗ V s and L ± s (λ, µ) are composite Lax operators acting over V s ⊗ V s ⊗ C 2 given by with the index set J = {x, y, z, 0}. For later convenience, we have introduced the normalization factor . Schematic depiction of a quasilocal charge X s (λ) for a spin chain composed of N sites: Each vertex represents a copy of an irreducible spin-s Lax operator L s . Each row represents one copy of an auxiliary space V s , carrying their own rapidity variables (λ and µ). Horizontal stacking pertains to tensor multiplication with respect to physical spaces V 1 2 ∼ = C 2 , while vertical stacking should be understood as tensor multiplication with respect to auxiliary spin spaces (for physical components ordinary multiplication applies). The dashed lines denote partial tracing with respect to auxiliary spaces V s . The upper row is acted upon by the derivative operation ∂ µ (magenta), where Leibniz chain rule should be assumed. In addition, a reducible two-component Lax matrix L s (λ, µ) sits on every vertical rung (shown in blue only for the 3rd site). Notice that to generate a quasilocal charge X s (λ) one has to finally set µ = λ.
The central object to establish pseudolocality of the family X s (λ) to be considered is the normalized Hilbert-Schmidt kernel (HSK) Evaluation of expression (3.7) requires the introduction of an auxiliary transfer operator Equipped with this result, the quasilocality condition for X s (λ) is equivalent to demanding that is finite and non-zero. The goal is to obtain K s,s (λ, µ) by calculating the dominating (i.e. the largest in modulus) eigenvalues of auxiliary transfer matrices T s,s and L ±0 s .
Let τ j s (λ, µ) denote the eigenvalues of L +0 s (λ, µ) = L −0 s (λ, µ), while for coinciding parameters we put τ j s (λ) ≡ τ j s (λ, λ) (and similarly L ± s (λ) ≡ L ± s (λ, λ), and N ± s (λ) ≡ N ± s (λ, λ)). In the normalization we use, the dominating eigenvalues τ 0 s (λ) of L ± s (λ) are equal to 1, while the rest of the spectrum is sub-unitary, |τ j s (λ)| < 1 for j = 0. Moreover, by analyzing the spectra of matrices L ± s one can learn that the left/right eigenvector L ±0 |ψ 0 = |ψ 0 , ψ 0 | L ±0 = ψ 0 | (corresponding to the leading eigenvalue), is the spin-singlet state (3.10) The singlet vector |ψ 0 obeys ( S 1 + S 2 ) |ψ 0 = 0 where (following Ref. [55]) auxiliary spins are given by S 1 = ( s ⊗ 1 s ) and S 2 = 1 s ⊗ s and act over V s ⊗ V s . For the remaining These relations imply that the product state The last step to perform in order to show that the kernel from Eq. (3.9) is finite, is to rigorously show that τ s,s (λ, λ, µ, µ) = 1 is indeed the leading eigenvalue. This statement can be conveniently phrased by defining the operator 13) and showing that it is a positive-definite operator on the orthogonal complement of the singlet state |Ψ 0 . The SU (2) symmetry of the isotropic point ∆ = 1 makes the task of demonstrating that the matrix (3.13) represents a contracting map much easier. The scalar component of double Lax operator L +0 s (λ, µ) can be readily expressed in terms of the Casimir from where we conclude that the eigenvalues are while the dominating vector is clearly the spin singlet state |Ψ 0 . A complete proof and further details on this part are presented in Ref. [55] and the Supplementary material attached to it. Note that factorizability of the leading eigenvalue, Eq. (3.12), in fact implies the inversion identity (3.2). Similar inversion formulae have been discussed earlier in the literature [8,85,86]. Quasilocality then follows essentially as a corollary of Eq. (3.12). To finalize the proof it remains to be shown that the kernels K s,s (λ, µ) given by Eq. (3.9) are well-defined and can be evaluated directly by accounting only for the contributions from the leading eigenvalues of auxiliary transfer matrices T s,s (λ, λ , µ, µ ) and L ±0 s (λ, µ). Using arguments based on the first order perturbation theory in combination with factorizability of the leading eigenvalue results in (3.16) 3.1.1. Local operator expansion. An important practical advantage of the present formulation is that X s (λ) can be readily expanded in terms of local operators. This step is of main interest in applications where evaluation of local correlation functions plays the primary role. To see how this works, we consider the resolution of operators X s (λ) with respect to local clusters of r adjacent spins (2.5), by summing over all projections onto the finite sublattices of length Λ, .

(3.17)
Of course the 'limits' have to be understood in the sense as discussed in Sec. 2.1. Here operators d r (λ) represent projections of X s (λ) onto local densities with support size (order) r, where by virtue of Eq. (2.9) the HS norms d r (λ) HS decay exponentially with r. We note that strictly local charges H (k) are, ignoring irrelevant constant prefactors, just the Taylor series coefficients generated by expanding X 1/2 (λ) around λ = 0. Thanks to the factorizability of the leading eigenvalue and the corresponding eigenvector, all k-point amplitudes (σ α [1,k] , X j (λ)) can be efficiently computed by introducing a set of auxiliary vertex operators, one for each α ∈ J . This allows us to write a matrix product representation This formula is exact in the thermodynamic limit (N → ∞, see Eq. (3.17)) while in finite lattices there are corrections which vanish exponentially in N and can be estimated in terms of subleading eigenvalues of T s,s . The boundary vectors in Eq. (3.19) are set as Here we wish to note that, in order to produce a non-vanishing amplitude, the µderivative which is included in the definition of X s (λ) (cf. Eq. (3.3)) must necessarily act on the first site in the matrix product representation of operators X j (λ) (see Eq. (3.4)) due to Eq. (3.11).

Computation of Hilbert-Schmidt kernel.
Quasilocal charges X s (λ) are linearly independent, but not manifestly orthogonal with respect to HS inner product. Below we show how to obtain explicit expressions for kernels K s,s , and subsequently use them to carry out the 'Gram-Schmidt orthogonalization'. For simplicity we restrict our discussion to the isotropic point ∆ = 1, where we find while boundary vectors given in Eq. (3.20) can now be chosen symmetrically and take the form A direct route to evaluate HSK K s,s (λ, µ) as defined in Eq. (3.16) is to rewrite the initial representation (3.9) in terms of the resolvent of the auxiliary transfer matrix (see Ref. [55] for details) which can be rewritten in terms of a geometric series where |Ψ = α∈{x,y,z} |ψ α ⊗ |ψ α . In the above sum, each term Ψ| [T s,s (λ, µ)] k |Ψ actually corresponds to a contribution of an order-k density d k (λ), which is finite since it obeys the quasilocality condition. A key point in this calculation is to recognize that the leading eigenvalues reside in an invariant singlet subspace Noticing that F s,s does not couple |Ψ 0 to the remaining states from V 0 allows to cast Eq. (3.16) expressed as Eq. (3.23) in terms of a solution to a linear system of 2s equations, introducing the restriction of F s,s to subspace V 0 denoted by F where a 2s (λ) = s/(s 2 +λ 2 ) are Cauchy-Lorentz kernels. Kernels a 2s play the central role as quasi-particle scattering phase shifts of the underlying scattering theory, as briefly explained in Sec. 5.2.2.
3.1.3. Orthogonalization procedure. The aim here is to construct mutually orthogonal families of quasilocal operators X s (λ). By considering a generic charge with s > 1 2 we set (3.27) and minimize the inner product by solving the following variational problem: This yields a linear system of 2s − 1 coupled Fredholm integral equations, . Explicit results for the solutions of Eq. (3.30) can be found in Ref. [55].

Fusion hierarchy approach
We have previously highlighted the meaning of the inversion identity Eq. (3.2) and learned about its importance for identifying quasilocal conserved quantities. In this section, we explore a different route and show how to consistently retrieve the inversion formula from Eq. (3.2) by resorting to an algebraic diagonalization of higher-spin operators T s (λ).
In Sec. 2.2.2 we explained how the entire set of canonical T -operators can be simultaneously diagonalized by means of Baxter's Q-operator. Assuming that the large-N behaviour of Eq. (2.33) can be read from the N -dependent scalars ζ 2s,k , the sum is dominated by the highest term at index k = 2s, This manifestly produces the inversion formula (3.2) on the level of operators. We therefore expect that the formula (3.31) also makes sense on the level of typical eigenvalues and can therefore be used to obtain the action of X j (λ) on (Bethe) eigenstates.
In view of Eq. (3.31) we, in addition, conclude that the 'quasilocality domain' can be analytically continued from the real axis to the whole 'physical strip' in the complex plane, We note that the charges X s (λ) are Hermitian for λ ∈ R, but they become non-Hermitian for Im(λ) = 0.
As a consequence of Eq. (3.31), the general version (for arbitrary anisotropy ∆) of the unitary quasilocal charges from Eq. (3.1) admits a useful compact representation in terms of the Q-operator The charges X s (λ) can now be effectively diagonalized using the fact that eigenvalues of the Baxter's Q-operator (denoted by Q(λ)) are q-deformed polynomials with zeros coinciding with the set of Bethe roots {λ j } §, where c is an inessential scalar prefactor. At this point the identification with the spectrum of the model has been made, which shall play a central role in the subsequent discussion of applications in the area of 'quantum quenches'. Further details are presented in Sec. 5.2.

Gapless regime
In this section we generalize the results for the isotropic and gapped cases derived in the previous section to the gapless regime. Without loss of generality we restrict our considerations to the positive side of the critical interval ∆ ∈ (0, 1). For technical reasons we exclude the non-interacting point at ∆ = 0, which due to the exceptional degeneracy requires a special treatment.
In the gapless regime we introduce a three-parametric family of conserved operators An important difference with respect to the family of charges used in the gapped regime is that T -operators now acquire another quantum label, the so-called (string) parity number u ∈ {±1}. The latter merely represents a π/2 displacement of the spectral parameter in the imaginary direction, namely It is important to stress that operators from Eq. (3.35) do not automatically inherit quasilocality from the gapped counterparts. Even though in the present case the § Here we ignore a subtle fact that Baxter's Q-operator becomes singular in the presence of periodic boundary condition and requires to be regularized in some way [82]. In our formulae, Q-s always appear in certain ratios which are always well-behaved. Apart from this, we do not rely on an operatorial construction of Q-operator, but merely use its spectrum which pertains to Bethe string configurations.
structural form of the solution Eq. (2.33) to the Hirota equation remains unaffected, the scalar functions undergo the following modification , k = 0, 1, . . . 2s. (3.37) For the inversion identity to hold, the following condition should be satisfied In stark contrast to the gapped (and isotropic) case, given a root of unity deformation q = exp (iπl/m), only a finite number of (linearly) independent charges with quantum labels (s, u) can satisfy this condition. For instance, for the simple roots of the form η/π = 1/m, there are precisely m − 1 charges with labels (s, +) for s = 1 2 , 1, . . . m−1 2 . On the other hand, at generic roots of unity identifying the complete set of charges becomes more involved [87]. To give a flavour, at η/π = 3/7, we have four independent families of charges corresponding to the set While the total number of quasilocal charges at a given value of η and their associated quantum labels might seem a bit arbitrary at a first glance, it is explained below in Sec. 5.2, that the labels can be matched to the known and well-established quasi-particle thermodynamic content of the model.

Quasilocal charges from non-unitary representations
Here we turn our attention to the construction of quasilocal conserved charges, using non-unitary representations of U q (sl (2)). In the first part we consider the highest-weight representations as elaborated on in Ref. [50] (see also [51]), building on previous results [42,49]. This construction yields conserved operators which break the spin reversal symmetry of the model and which are used for establishing the ballistic transport property of the high-temperature anisotropic Heisenberg model. The second part discusses an analogous construction, this time with semi-cyclic representations which, interestingly, break even the U (1) symmetry of the model, following Ref. [52].

Charges from highest-weight representations
Let us remain in the gapless regime and keep the root of unity parametrization of the anisotropy given as ∆ = cos(η), or q = e iη , with η = πl/m, and l, m ∈ Z + co-prime.
In what follows, the basic building block of our construction is a reparametrized Lax operator Eq. (2.22), where for our convenience (and to comply with Refs. [49,50]) we perform a rescaling by a factor sinh (η)/ sin (λ) and subsequently make a substitution η → −iη (but refraining from substituting λ → −iλ as in Sec. 2.2.1). This results in a trigonometric form of the Lax operator: . (4.1) Considering the m-dimensional highest-weight auxiliary space representation (2.25), the commuting transfer operators are given in accordance with the standard prescription Without further ado, we define the following family of commuting operators by differentiating T hw s (λ) with respect to continuous spin s, Note that in this way the contribution of the magnetization M = x∈Λ σ z x cancels from Z(λ), and hence, by construction, only the operator terms acting non-trivially on two or more sites remain. With the aid of Lax operator components we can, following the logic presented in Sec. 3.1, expand the family of conserved operators Z(λ) in the large-N limit in terms of r-spin clusters, Ŝx (dr(λ)) . (4.5) The amplitudes are now encoded as matrix product expressions while the boundary vectors are given as L| ≡ sin λ sin η 0| L − , |R ≡ sin λ 2η L + |0 (in addition to that, (σ − ⊗ σ + , Z(λ)) = 1). By inspecting the Lax components (cf. Eq. (4.8) below) we learn that all amplitudes which violate the selection rule α 1 = − and α r = + vanish. Another remark that we would like to make is that in any finite-N lattice the expression for the conserved operators Z(λ), as given by Eq. which gets exponentially suppressed with N with respect to HS norm (see Ref. [50]).
Let us briefly comment on the technical part of what steps have been made to arrive at Eq. (4.5). Due to translational invariance, each term in the operator expansion of Eq. (4.3) has been rearranged so that the right-most position in the product of Lax operators always belongs to the differentiated Lax operator, L(λ). The trace in ∂ s T s (λ) is then split into two parts, a sum over states |n = 0 , producing the correction (4.7), and the projection onto the 'vacuum' |0 part which results in Eq. (4.5). Explicit form of the amplitudes given by Eq. The associated auxiliary transfer matrix T is an operator on the reduced auxiliary space lsp{|n |n ⊗ |n ; n = 1, ..., m − 1} of the form | sin (nη) sin ((n + 1)η)| 2 sin (λ) sin (µ) (|n n+1| + |n+1 n|) . This matrix is contracting when parameters λ and µ lie inside the strip (4.11) Quasilocality of conserved operators from Eq. (4.5) is then an immediate consequence of this statement [50]. The reader has to be reminded that we have disregarded the The exact bijective correspondence, used to produce this symmetrized matrix form is |n ⊗ |n ↔ | sin (nη)| |n , n| ⊗ n| ↔ | sin (nη)| −1 n|. correction term (4.7). Using an equivalent procedure to the one described above, it can be shown that the contribution of this term to HSK is exponentially suppressed in the system size [50]. In order to see that, one must examine the action of T on invariant subspaces of V (m) s ⊗ V (m) s which are spanned by elements |n ⊗ |n + k , for different fixed k. Such a decomposition reduces the auxiliary transfer matrix into the block diagonal form. One then proceeds by proving that each block itself is a contracting matrix.
Evaluating Eq. (4.9) amounts to solving the linear equation for the components ψ j = j|ψ of |ψ . The final result is The construction from above can also be applied to the case of twisted boundary conditions. The Hamiltonian then consists of an open boundary part and a two-site term, acting on the first and the last site of the chain introducing a flux parameter φ, such that the φ = 0 case corresponds to the Hamiltonian with periodic boundary conditions. The transfer operator in case of twisted boundary conditions takes the following form, while the conserved charges are generated similarly as in Eq. (4.5), with the prescription (4.3), but using a modified s-derivative, ∂ s → ∂ s + iφ. In this case the HS kernel from Eq. (4.9) remains independent of φ and hence quasilocality is preserved. This concludes the review of highest-weight conserved charges.

Charges from semi-cyclic representations
After having discussed how to obtain quasilocal charges from the highest-weight auxiliary modules, we now turn our attention to another family of representations of U q (sl (2)) at roots of unity -the semi-cyclic representations. To this end we retain the m-dimensional auxiliary spaces, V  Here we have introduced the 'coupling' parameter α, linking the first and the last basis states. ¶ Since the action of ladder operators is periodic only in one direction, such a representation is referred to as semi-cyclic. The algebraic relations (2.23) are still satisfied.
There are other possible alterations of the representation of the algebra generators, all of them resulting in a certain kind of periodicity [5]. In the following we will, for the sake of simplicity, only consider the above example. Since all other semi-cyclic representations generate the same quasilocal charges, up to trivial transformations, this means no loss of generality [52].
As we will see, the coupling of the lowest and highest-weight vectors in V (m) s results in a family of conserved charges which do not conserve the total magnetization M (i.e. they break the U (1) symmetry). Apart from this, the charges considered here only exist for odd dimensions m. While non-conservation of magnetization is obvious from the explicit expressions, non-existence of these charges for even m stems from the mismatch between the canonical U q (sl(2)) relations (2.23) and slightly modified relations which directly imply commutativity of the transfer operators with the Hamiltonian, see ¶ In our notation, the dependence on additional parameter α will not be explicitly written. One should nevertheless bear this dependence in mind.
Refs. [52,88]. The allowed values of anisotropy parameter are: for k, l ∈ N, l < k. (4.17) 4.2.1. Constructing the semi-cyclic charges. The same transfer matrix as in the case of highest-weight representations can be used, but this time we differentiate it with respect to the coupling parameter α, at α = 0 and s = 0. We now put with the only non-trivial component of L being The conserved charges are this time defined as where T sc s (λ) is the semi-cyclic transfer matrix defined with auxiliary space generators (4.16). Once again the formula (4.5) applies, thereby the amplitudes can be expressed in a canonical way with L| ≡ sin λ sin η 0| L − , |R ≡ sin λ sin η L − |0 . The remaining string of Lax components in the LHS of Eq. (4.20) must connect 1| to |m − 1 so the second sum in the expansion (4.5) actually starts at r = m. Because each term of Z(λ) consists of a surplus of exactly m operators σ − over operators σ + , these charges do not conserve magnetization M . A diagrammatic presentation of semi-cyclic Z-charges is shown in Fig. 2 (panel (b)).

Quasilocality
. What remains to be done is to derive the quasilocality property. The latter follows from a slightly modified calculation with respect to the situation which we had previously with the highest-weight charges. A careful inspection shows that the same auxiliary transfer matrix as given by Eq. (4.10) for a highest-weight representation, can be used to express the semi-cyclic HSK as Again, a solution of a simple tridiagonal system (4.12) of equations yields an explicit expression K(λ, µ) = 1 4 ψ m−1 = sin (λ) sin (µ) sin (λ + µ) 2 sin 2 (η) sin (m(λ + µ)) .

(4.22)
To produce ψ m−1 as defined previously in Sec. 4.1, the states |1 and |m − 1 have to be exchanged. To this end we conjugate Eq. (4.21) and recall that T(λ, µ) is symmetric.

Applications
Let us finally focus on various physical applications of quasilocal conserved charges in the domain of non-equilibrium quantum physics. Here both classes considered above, i.e. unitary and non-unitary charges, will be examined. We shall begin with non-unitary Z-charges and show how they directly relate to non-equilibrium states with currents. On the flip side, unitary X-charges will play an instrumental role for understanding equilibration in quantum quenches. But before heading on, we need to clarify an important role of the spin reversal parity symmetry and its breaking.
Spin reversal and CPT symmetry of generic transfer matrices. We wish to elaborate on an important Z 2 symmetry of all finite-dimensional unitary representations of the quantum group U q (sl(2)), and consequently of the XXZ Hamiltonian itself, which is manifestly broken for non-unitary representations. This symmetry breaking has some remarkable physical implications which shall be presented in the following. The Z 2 symmetry under scrutiny is a parity generated by the spin-reversal canonical transformation In fundamental representation the latter amounts to applying the product of σ x , where A can be any observable on the entire Hilbert space H. It is easy to show that all transfer matrices belonging to finite-dimensional irreducible unitary representations are manifestly P -invariant, implying the same property also for the corresponding local and quasilocal charges, For the root of unity deformations q = exp(iπl/m) there exists another class of irreducible representations. These are non-unitary m-dimensional highest-weight representations of U q (sl(2)) discussed previously in Sec. 4. They are distinguished by the property, which can be readily verified, that no similarity transformation x → GxG −1 of the auxiliary space representation of the algebra exists which would generate the spinreversal canonical transformation (5.1). These non-unitary representations (2.25) are labelled by a complex-spin parameter s ∈ C and are henceforth not P -invariant. We note that existence of an invertible G, such that Gs z G −1 = −s z , Gs ± G −1 = s ∓ , is equivalent to a spin-reversal symmetry of the Lax operator (4.1) P L s (λ)P −1 = GL s (λ)G −1 , where P acts nontrivially only on the physical space and G only on the auxiliary space, and consequently implies Eq. (5.3).
The highest-weight transfer matrices for complex spins and the quasilocal charges they generate instead exhibit a weaker symmetry, P T hw s (λ)P −1 = (T hw s (π − λ)) T , P Z(λ)P −1 = (Z(π − λ)) T , s ∈ C. (5.5) As the transposition can be associated with time-reversal operation, while reflection of the spectral parameter λ → π − λ can be thought of as the 'charge conjugation' (after a suitable rotation and a shift of the spectral parameter it would correspond to λ →λ), the relation (5.5) can in fact be interpreted as a CPT symmetry of a generic highest-weight transfer matrix. The fact that complex-spin transfer matrices T hw s (λ) break spin-reversal symmetry can be fruitfully explored for the analysis of ballistic spin transport in anisotropic Heisenberg chains as will be demonstrated in Sec. 5.1.
An equivalent CPT symmetry (5.5) holds also for the semi-cyclic transfer matrices and the corresponding quasilocal charges as discussed in Sec. 4.2.

Mazur bounds on Drude weights
5.1.1. Ballistic linear response. The main motivation for constructing pseudolocal conservation laws originated from the idea of using such objects to estimate the ballistic contribution to transport coefficients, such as Drude weights or, more generally, zero frequency dynamical susceptibilities [63,64]. It is perhaps worth noticing that related indicators of ballistic transport are nowadays directly experimentally accessible [30][31][32][33].
By considering an extensive current J = xŜ x (j) with a local density j, say the spin/particle/energy/etc. current, the Kubo linear response formula for the nondissipative (real) part of the respective conductivity is of the form Under certain mild assumptions on analyticity of local correlation functions, which are discussed in Ref. [89], the order of the limits for D J can in fact be reversed and using time-invariance of the thermal state ω β the Drude weight gets expressed in terms of time-averaged current as A nontrivial value of the Drude weight D J > 0 signals the ballistic (ideal) DC transport and is equivalent (cf. Eq. (2.7)) to the statement that the time-averaged current is a pseudolocal operator with respect to the Gibbs state ω β (see also Ref. [65]). We have thus related pseudolocality of time-averaged observables to ballistic linear response.

Mazur bound.
Computing time-averages of current operators seems a highly nontrivial task in interacting models. One can instead estimate the Drude weight from below using a bound due to Mazur [90] and Suzuki [91] in terms of some conserved Hermitian operator I = I † , [H, I] = 0. We start by writing out the expectation value of a nonnegative operator (J − αI) 2 , where α ∈ R is a free parameter, We used the fact that ω β (JI) = ω β (JI), which is due to the time-invariance of ω β and conservation of I. After optimizing Eq. (5.12) with respect to α, we obtain Dividing by 2N and taking the limit N → ∞, we produce the Mazur bound on the Drude weight, which has first been pointed out in Ref. [38], . (5.14) In summary, a conserved pseudolocal operator I which satisfies ω β (jI) = 0 implies ballistic transport and consequently allows to put a strict lower bound on the Drude weight. For example, by taking a translationally invariant extensive conserved operator I = x S x (q), with density q satisfying ω β (q 2 ) < ∞, one finds D J > 0 if x ω β (jŜ x (q)) = 0, where the last sum always converges due to exponential clustering of Gibbs states in one dimension [92].
In addition, as a consequence of an effective causality on the locally interacting lattice (i.e. Lieb-Robinson bounds [61]) it can be shown that the above Mazur bound holds even when I is not exactly conserved on any finite lattice with open boundaries but the commutator [H, I] contains terms localized near the boundary sites [89].
When dealing with a larger set of pseudolocal conserved operators, say a countable set {I k , k = 1, 2 . . .}, the Mazur bound can be further improved. To see how this works, we study the operator (J − k α k I k ) 2 , which after repeating the above reasoning results in . In this sense, if the above bound gets saturated for all local currents j, it would be meaningful to regard the set of pseudolocal charges {I k } as being complete. It is presently not known if such complete sets of pseudolocal conserved operators can be systematically identified in interacting models.
In previous sections we have defined and discussed certain continuous families (rather than discrete sequences) of pseudolocal charges which were referred to as quasilocal (cf. Eq. (2.9)). They comprise the charges X s (λ) and Z(λ) which are analytic in λ ∈ C and become quasilocal when restricted to suitable domains D ⊂ C. Since all X s (λ) are even under spin-reversal transformation, while the spin current is odd, we immediately conclude that all the charges coming from unitary representations are irrelevant for the Drude weight, namely ω β (jX s (λ)) ≡ 0. For this reason we subsequently consider only the set {Z(λ); λ ∈ D}. Similarly as in the previously considered discrete case, we start by studying the following operator where the integration is over the quasilocality domain D. It is worth stressing that in general Z(λ) are not Hermitian. Nevertheless, the expectation value of B † B is always nonnegative We proceed by defining the overlap coefficients of an extensive observable J along the conserved operators in terms of the holomorphic function ω β (jZ(λ)), (5.19) assuming the limit N → ∞ exists. For infinite temperature β → 0 the existence of the limit and consequently holomorphicity of Z(λ) simply follow from the explicit matrix product operator expression (4.6). The limit in the last term of Eq. (5.18) exists as well, due to pseudolocality of Z(λ), and can be written in terms of a Hermitian kernel Therefore D J should satisfy the inequality The bound is manifestly real due to the hermiticity of the kernel. It is noteworthy that the lower bound (5.26) agrees exactly with the Thermodynamic Bethe Ansatz (TBA) calculation [39,93] at the special (isolated) points of anisotropy at η = π/m, corresponding to q-deformation at simple roots of unity (l = 1). Since TBA calculation for other values of l seems to be highly nontrivial and has not yet been performed, we can only conjecture that the bound (5.26) is in fact saturating the exact value of high-temperature spin Drude weight for a dense set of commensurate anisotropies ∆ = cos (πl/m). Such a conclusion can also be based on the comparison with numerical results of the state-of-the-art density matrix renormalization group (DMRG) methods [45,46] which indicate no significant deviations from the lower bound D K [94]. One obtains similarly good agreement by comparing to exact real-time dynamical simulations with random initial wave-function sampling on smaller systems and perform appropriate finite size scaling analysis [95].  [49]. In comparison we show (in red) the bound optimized for a single charge obtained initially in Ref. [42]. In either case the bound exhibits a pronounced fractal-like (nowhere continuous) dependence on parameter ∆. Note that one can use the concept of operator time averaging to formally describe the steady state of XXZ model pierced with a flux φ and undergoing a small flux quench φ → φ + δ φ , namely starting from a thermal density matrix β , one may show [96] that after-quench current carrying steady state is given by the density operator Furthermore, the concept of time-averaged extensive local operators has been used to implement a useful numerical algorithm to search for unknown quasilocal charges of an arbitrary locally-interacting lattice model [97]. One should simply recall that for any operator O, which is an extensive translational invariant sum of traceless local operators, O is by construction a pseudolocal conserved operator, or it vanishes in a suitable norm if O is ergodic. Taking a maximal linearly independent set of such local extensive operators {O n } up to some maximal order of locality M N , enumerated with n = 1 . . . M, M ∼ (d 2 ) M , one can define a nonnegative definite HS kernel as the matrix K n,n = (Ō n ,Ō n ) = (O n ,Ō n ). The number of independent pseudolocal conserved operatorsŌ n , with effective support not larger than M , can thus be determined as an effective rank of the matrix K with eigenvectors yielding the quasilocal charges expanded in {O n }. Implementation of this method in the case of isotropic XXX model [97] gave the first constructive empirical evidence on existence of unitary quasilocal charges X s [55].

Quantum quenches
Motivated by recent experimental progress in optical lattices [21,[23][24][25][26][27][28][29] and a plethora of numerical simulations of strongly correlated matter in low dimensions, a very popular setup studied over the last decade is the problem of a 'quantum quench' [98][99][100][101][102][103][104][105][106][107]: at initial time, an ideally isolated (closed) system is prepared in an initial state |Ψ , and subsequently, by a sudden change of interactions, let to evolve according to a unitary evolution generated by a post-quench Hamiltonian H. The situation which is particularly appealing from the theoretical viewpoint is when H is integrable. Many aspects regarding quantum quenches, ranging from classical field theories [17], conformal field theories [13,14], disordered systems [19], Luttinger model [16], to integrable lattice systems [15,18,20] are discussed in the reviews of the present volume.

Complete Generalized Gibbs Ensembles.
One of the pivotal questions is to understand the process of equilibration from the microscopic perspective [108,109]. In homogeneous quantum systems with generic interactions the relaxation towards canonical Gibbs ensemble is typically explained in the framework of the eigenstate thermalization hypothesis [108,110,111], which states that eigenstates which are close in energy give approximately the same values of local correlation functions. The situation with integrable interactions is however different as time-evolution is severely constrained due to the existence of a macroscopic number of local (and quasilocal) conserved quantities.
It has been conjectured in Refs. [112,113] that statistical properties of local quantities in many-particle quantum systems which possess an 'extensive number' of conserved local charges I n should comply with predictions of a Generalized Gibbs Ensemble [20,103,105,107,[114][115][116] , given by a formal expansion ρ GGE ∼ exp − n β n I n . (5.29) The 'GGE conjecture' asserts that the ergodic average of an operator A with a finite support can be reproduced by tracing with respect to an appropriate GGE of the form (5.29), with the 'chemical potentials' β m being determined from expectation values of the charges with respect to the initial state. A great body of work has already been devoted to applicability of the GGE in non-interacting models [103,105,106,117,118], and a closely related phenomenon of prethermalization [26,115,[119][120][121].
Explicit verification of the GGE paradigm in a truly interacting quantum integrable models required a bit more effort though. Initial studies focused on Heisenberg XXZ chain and compared predictions of truncated GGEs made of hitherto known local charges against numerical results for the time-evolved local observables [122][123][124]. First exact results have been obtained for the case of the Lieb-Liniger model in Ref. [125] by resorting to the so-called quench action method, developed previously in [126] (cf. [15] for a review). In this approach, a generalized free energy functional is constructed which incorporates the restrictions imposed by the initial condition in the form of an exact overlap coefficient. By employing TBA framework [127][128][129][130], the saddle-point of such a functional yields the sought for steady-state ensemble via coupled non-linear integral equations for a set of variational variables. These thermodynamic variables are, as we shall shortly discuss, a set of analytic functions representing distributions of Bethe strings.
Sometimes, e.g. for certain simple product states, the overlap formulas which enter as an input to quench action method can be evaluated explicitly [131]. Two independent studies [40,41] unambiguously demonstrated that GGEs composed from only the hitherto known local charges fail to recover the exact results (see also Refs. [132,133]). The failure has been related to the fact that strictly local charges Eq. (2.27) do not provide enough information to determine the distributions of the bound states which are present in an initial state [54]. The results of these studies hinted on the presence of additional (sufficiently local) conservation laws in the unitary (or spin-reversal symmetric) sector.

String-charge duality.
Here we explain, following Ref. [87], the connection between the spectra of quasilocal charges X s and distributions of Bethe strings. The latter should be interpreted as thermodynamic particle content of an integrable lattice theory. Hence, the main task shall be to extract the large-N behaviour of eigenvalues of T -operators. A convenient tool to achieve this is to employ the Baxter Q-operator [8,81,82] and exploit the fact that its eigenvalues are given by a (deformed) polynomial with zeros coinciding with Bethe roots. Below we present the main steps by specializing to the gapped regime.
Bethe equations and string hypothesis. To set the stage we need to briefly describe how to characterize the spectra of integrable lattice models in the thermodynamic regime. The elementary building block of an integrable model is the single-particle S-matrix S 1 which for the XXZ model reads .
From a scattering theory point of view, the spectral parameters λ and µ pertain to rapidities of the two quasi-particles involved in a scattering event. For composite objects which consist of j excitations -commonly referred to as the j-strings -a set of fused scattering matrices S j are introduced Scattering among two different types of strings is governed by string-to-string scattering matrices With the aid of scattering matrices, the Bethe Ansatz equations, representing a quantization condition for quasi-particle rapidities λ j in a periodic system, are cast in the form Here M is the number of Bethe roots (related to the magnetization of the eigenstate) and p(λ) encodes the momentum of an elementary excitation on top of a ferromagnetic vacuum state, .
The string hypothesis [128][129][130]134] states that in the large-N limit the Bethe roots (i.e. solutions λ j to Eq. (5.34)) for a typical eigenstate become equidistantly displaced in the imaginary direction in the rapidity complex-plane, Such string formations physically correspond to bound states of magnons. By partitioning the Bethe roots in terms of strings, Bethe equations (5.34) can be rewritten in terms of string centres λ k α ∈ R. Thus, taking their logarithmic form and considering the thermodynamic limit when string centres get smoothly distributed along the real axis, we arrive at the following non-linear coupled integral equations [128,130,134] known as the Bethe-Yang equations for the strings. The integral kernels in Eq. (5.37) are given by the derivatives of scattering phase shifts and the corresponding string-to-string phase shifts a j (λ) = −i∂ λ log S j (λ), a j,k (λ) = −i∂ λ log S j,k (λ), (5.38) in the respective order. One of the advantages of Eq. (5.37) in comparison to the finite-volume counterpart is that we no longer have to deal with a complicated set of quantized quasi-momenta (encoded by Bethe roots λ j ). Instead, now quasi-momenta take values in the continuum which allows us to cast the description in terms of analytic distributions ρ j (λ) which count the number of Bethe strings whose centres occupy an interval [λ, λ + dλ]. Similarly, ρ j (λ) denote the complementary variables, parametrizing distributions of Bethe holes (the positions of string centres which are in principle available, but remain unoccupied).
Thermodynamic spectra. To obtain the spectra of charges X s we make use of representation (3.33). By neglecting the contributions which are subleading in N we have where |{λ j } denote a Bethe eigenstate parametrized by a set of roots {λ j }. Working under the 'string hypothesis' (cf. Sec. 5.2.2), the spectra of quasilocal charges X s , can be readily expressed in terms of densities of string centers ρ j (λ). Specifically, by plugging the expression for the spectrum (cf. Eq. (3.34) in Eq. (5.39)), we arrive at [87] The set of kernels G 2s,k can be expressed using scattering matrices among the strings Let us introduce a discrete d'Alembert operator , whose action on any set of objects f s ≡ f s (λ) (with s = 1 2 Z + ) which are analytic inside the physical strip P η is prescribed by By acting with the d'Alembertian on the kernel functions from Eq. (5.42) we conclude that G j,k (λ) = δ j,k δ(λ), j = 2s ∈ N. (5.44) This result allows us to interpret G j,k as a discrete 2D Green's function of the 'wave operator' . The relation Eq. (5.41) can be readily inverted, enabling to express the entire set of density functions ρ j (λ) in terms of eigenvalues of the charges X s (λ) as [87] ρ 2s (λ) = X s (λ). (5.45) The distributions of holes ρ 2s (λ) can be obtained in a similar fashion [56,87], In the scope of quantum quench applications, a set of densities ρ 2s provides a complete description of local correlation functions (cf. [133,135]). Finally, let us make a brief account on the gapless regime as well. Although the string hypothesis in the |∆| < 1 regime can still be formulated, taking the deformation parameter q = e iη from the unit circle makes the analysis rather cumbersome and technically involved. The string content in the gapless regime for an arbitrary value of anisotropy has been derived in Ref. [129]. Due to limited space we do not attempt to review it here. We nevertheless wish to point out the three principal differences in comparison with the situation in the gapped case: (i) string configurations acquire (beside the string length) an additional parity label u ∈ {±1} (see Sec. 3.3), (ii) the allowed string lengths depend strongly (and discontinuously) on η, and (iii) at root of unity value of q the number of allowed distinct string types is always finite. Moreover, in the spirit of string-charge duality, the number of (dynamical) strings should still be in a bijective correspondence with the number of quasilocal charges, as discussed in Sec. 3.3. To complete our example for η = 3π/7, where the charge content is given by a set (3.39), we provide the corresponding string content: (1, +), (1, −), (3, +), (5, −).
Below we explain a computational scheme to determine the densities of Bethe roots from the eigenvalues of X s (λ). This can be done, in contrast to a more common practice, without ever resorting to the variational approach based on a generalized free energy functional. The manifest locality of quasilocal charges X s (λ) in the spin basis (cf. Eq. (3.3)) greatly simplifies this task and allows us to resort to rather standard techniques.

Evaluation of charges.
In this section we address the problem of computing expectation values of the quasilocal charges X s with respect to a generic + pure state + Strictly speaking, we are implicitly assuming that our reference state is 'local', i.e. is compatible with cluster decomposition principle [136,137]. In this case we are able to express |Ψ in the thermodynamic limit as a single macrostate (a state given by prescribing distributions of Bethe strings).
|Ψ . While performing this task in full generality remains out of reach at the moment, we make a restriction to a class of matrix product states where an efficient implementation is possible. In what follows we essentially recast the results of Refs. [122,124] in the present language.
In order to keep the level of technicality at a minimum, we shall in addition restrict ourselves only to periodic product states where |ψ is a state on the block of N p spins and N p ∈ N is the periodicity of the state. Our aim is to compute Due to the product structure of |Ψ we can make use of standard transfer matrix techniques. The first step is to introduce a boundary partition function 49) which is given by iterating a one-step auxiliary propagator, Subsequently we evaluate We note that partition functions given by Eq. (5.49) are in essence merely the contracted quantum transfer operators X s (λ, µ) from Eq. (3.4) (depicted in Fig. 1) where in the vertical direction we project onto components determined by the reference state |ψ . Such a contraction over one period N p yields the propagator from Eq. (5.50). The construction sketched above can be adapted for general translational invariant matrix product states (see Refs. [87,122,124]).

5.2.4.
Closed-form results. In Sec. 2.2.2 we already mentioned that higher-spin Toperators constitute the canonical solution to Hirota difference equations (alias the T -system). However, Hirota difference equations admit different solutions as well.
Remarkably, there exists a class of initial conditions which relax to equilibrium steady states (specified by a collection of density functions ρ Ψ j ) which can be cast as distinct solutions of the Hirota equations. Below we mention two particular examples, which have been previously studied in the literature, when equilibrium states admit simple representative product states: (i) a spin-singlet dimerized state |D = 1 √ 2 (|↑↓ −|↓↑ ) ⊗N/2 and (ii) Néel state |N = |↑↓ ⊗N/2 . In other words, these two states can be understood as members of a basin of attraction for equilibrium states which assume parametrizations in terms of non-canonical solutions to the functional relation of the T -system.
In the following, we use a small font to explicitly distinguish non-canonical tfunctions and q-functions, t s (λ), q(λ), from the canonical objects, i.e. fused transfer matrices T s (λ) and Baxter Q-operator Q(λ) defined in Sec. 2.2.2. By relaxing the constraint t 0 = ϕ − , the linear auxiliary problem associated to the Hirota equation takes the form and can be explicitly solved as (5.53) In the present case, q-functions can be considered as auxiliary complex-valued functions which contain the information about the equilibrium state at hand, closely related to auxiliary functions which enter in non-linear integral equations in the scope of the Quantum Transfer Matrix (QTM) method [138][139][140].
To keep things simple, we specialize below only to the isotropic point ∆ = 1. For the dimerized state |D the solution is remarkably simple and reads q D (λ) = λ 2 . These results generate the entire tower of t-functions which can, in turn, be mapped to y-functions y Ψ , j ∈ Z + . (5.55) As we have already explained (cf. Eqs.(5.45), (5.46)), the y-functions can be related to expectation values of the charges on the state |D , We remark that in practice one should work in the opposite direction: by computing a few initial values of the charges and employing the string-charge relationship one can explicitly check whether the y-functions fulfil the Y -system hierarchy. It is not clear if a general systematic procedure exists to directly determine which states admit a description in the Y -system format. The analogous expressions for the Néel state (including the expressions for the gapped case) are provided in Ref. [87].
To conclude this section, let us stress that the unitary charges X s from the compact sector cannot be sufficient for characterizing non-equilibrium steady states, i.e. states which exhibit particle currents. In this situation, the quasilocal charges Z(λ) which break the spin-reversal invariance have to be included [96,141]. To the best of our knowledge, it remains presently unknown how the Z-charges act on Bethe eigenstates.

Steady states of boundary-driven chains
Quantum transport is typically studied in the framework of the linear response theory. An alternative way is to adopt an open system perspective. A simple effective setup for that is to use the approach of non-unitary evolution equations which are commonly referred to as quantum master equations. A central concept here is a Markovian (timelocal) evolution (t) = eL t (0), (5.57) which preserves the trace and positivity of density operators at any time. The generator L is of Lindblad form and acts linearly on density matrices aŝ where H encompasses all interactions attributed to the unitary part of the process, and the set {A k } contains the Lindblad 'jump operators' which are used to model dissipative processes. For a comprehensive introduction on the Lindblad equation formalism we refer the reader to Refs. [142,143]. In Refs. [42,[144][145][146][147][148][149] Lindblad equation has been used to 'drive' a quantum manybody system far from equilibrium. Two common scenarios describe the situations where the Lindblad bath operators simulate (a) dephasing noise due to uncontrolled degrees of freedom in the bulk, or (b) particle/magnetic reservoirs with different chemical potentials/magnetizations attached at the system's boundaries.
General instances can be studied by adapting a time-dependent DMRG technique to the Liouville dynamics [150]. On the flip side, certain interesting situations permit an exact analytic description, the most notable example being non-interacting particles experiencing Gaussian noise which can be treated in a unified manner within the formalism of 'third quantization' [151,152]. While deriving exact solutions for the full Liouvillian dynamics of an interacting system remains an open challenge up to date, certain steady state density operators, i.e. fixed points of Liouvillian dynamics, which allow for an efficient matrix product form have been found and investigated (the first non-trivial example being perhaps the situation of noninteracting particles with bulk dephasing noise [153,154]). In some sense, one can understand these as quantum counterparts of their more popular classical cousins known as asymmetric simple exclusion processes [155,156].
For the Heisenberg spin-1/2 chain, a model under scrutiny in this review, driven by incoherent in/out boundary processes: the steady state in the weak-coupling limit has been constructed first in Ref. [42], and later on extended to the non-perturbative regime in Ref. [146]. What is remarkable, and perhaps somewhat surprising as well, is that the density operator of the current carrying steady state found in Ref. [42] is a fully mixed state perturbed with an operator of the non-unitary quasilocal family, namely The only distinction from the conserved operators given by Eq. (4.3) is that instead of taking the trace over the auxiliary space the adequate transfer matrix is now defined as an expectation value in the highest-weight state (vacuum) where the Lax operator is taken from Eq. (4.1). Consequently, the local operator expansion of the open boundary charge Z vac is given with the same formula as before, Eq. (4.5), where the shiftŜ x no longer acts periodically (meaning that the sum over x runs only up to N − r). While the vacuum transfer matrices and the derived quasilocal charges still mutually commute, [T vac s (λ), T vac s (λ )] = 0, and [Z vac (λ), Z vac (λ )] = 0, the manifest absence of translational invariance breaks the conservation property, where the Hamiltonian of the anisotropic Heisenberg chain is now taken with open boundary conditions and b = σ z sin(λ) sin(ηs)−σ 0 cos(λ) cos(ηs) is a boundary operator. The first identity follows straightforwardly from the RLL relation (2.15), while the second one follows from the first one after taking the derivative ∂ s | s=0,λ=π/2 . Note that the second line of Eq.(5.61) has a form of a conservation law, i.e. time-derivative of Z vac in a finite volume equals net surface currents, where the spin density σ z x plays the role of the formal 'current'. In spite of 'almost-conservation' in a finite volume, it has been rigorously shown in Ref. [89], resorting to quasilocality and Lieb-Robinson causality bounds, that Eq. (5.61) yields a conserved quantity in the thermodynamic limit and in effect provides an equivalent set of quasilocal conservation laws as those introduced in Sec. 4.1.
Moreover, it can be shown (see [49,146,147], and [12,157] for a review) that the vacuum transfer matrix generates an exact, non-perturbative steady state density operator via the purification ansatz , if one sets the spectral parameter and identifies the noise strength Γ with a complex auxiliary spin s, in either one of the following two ways λ = π 2 , tan(ηs) = iΓ 2 sin(η) or λ = 0, cot(ηs) = iΓ 2 sin(η) .
These two assignments yield identical steady-state density operator (5.62). In light of the fact that in the canonical σ z x -basis the amplitude operator Ω becomes a strictly upper-triangular matrix [146], the ansatz (5.62) can also be understood as a many-body Cholesky factorization. The ansatz (5.62) in fact exactly solves a much larger set of boundary-driven Lindblad equations, namely taking an arbitrary pair of asymmetric (left/right) noise strengths Γ L,R and adding arbitrary boundary magnetic fields in z−direction h L,R uniquely parametrizing two complex variables s, λ (see Ref. [12]). We note that the notation of this section was adapted for the regime |∆| < 1 where quasilocal Z-charges have an effect and the corresponding transport is ballistic. To find the gapped counterparts one has to make a substitution η → iη, or replace trigonometric functions with the corresponding hyperbolic functions.
It is also perhaps instructive to stress that the vacuum charges Z vac (λ) are manifestly nondiagonalizable objects with a nontrivial Jordan structure. For example, the spectrum of Z vac (π/2) only consists of {0}, hence the operator is nilpotent for any finite volume, but nevertheless generates a highly nontrivial steady state. The approach of generating quasilocal almost-conserved charges as perturbative solutions to boundarydriven Lindblad equations may be useful also in other integrable models (see Ref. [12] for a review) and should perhaps be further explored in future.

Future prospectives
Spin chains. Even though applications of quasilocal conservations laws which we covered in this review have been fully concentrated on the paradigmatic case of the Heisenberg XXZ model, it is natural to expect that analogous quantities exist for a much broader class of integrable models (see e.g. Ref. [53] for a recent application to gapless spin-1 chains). The simplest extensions should involve quantum lattice models associated with the so-called fundamental solutions to Yang-Baxter equation, with underlying symmetry algebras based on Lie algebras of higher rank and their quantizations (deformations). Additionally, supersymmetric cousins (e.g. t − J model, EKS model) shall be of interest in paving the way towards the celebrated Hubbard model [158][159][160]. Note that a novel family of transfer matrices which violate particlehole symmetry and correspond to non-unitary auxiliary representations has recently been proposed for the Hubbard model [161], based on preserving the integrability of the associated boundary driven master equation [162]. A possibility of generating quasilocal conserved quantities remains to be explored.
For all models mentioned above it is well-known that thermodynamic spectra can be partitioned into Bethe root compositions (strings) which pertain to bound states of elementary excitations. In order to ensure that macrostates (e.g. thermal states and their generalizations) are mutually distinguishable, the number of distinct particle types (see Refs. [75][76][77]79,163,164]) has to be matched with the number of independent families of (quasi)local charges.
Another example of an integrable theory which has recently drawn a great deal of attention due to its experimental significance is a Bose gas with δ-like repulsive interactions, known better as the Lieb-Liniger model [165] (Nonlinear Schrödinger equation in the language of second quantization). Yet, in spite of its wide popularity, the second-quantized form of the entire tower of local charges have not been obtained explicitly so far [166,167]. Besides, there also exist certain obstructions which are intimately related to pathological UV divergences as discussed in Refs. [125,168]. In Ref. [169] the authors attempted to overcome the difficulties by 'mildly' relaxing the conventional form of locality. Alternatively, there is an option to employ a suitable integrable regularization (e.g. by introducing a UV cutoff) allowing to treat the lattice counterparts in the thermodynamic limit first, then construct/compute the observables, and take the continuum limit only at the end (see e.g. Ref. [170]). The effectively local, or quasilocal conserved charges could then be derived using the methods presented in this review.
Integrability in AdS/CFT correspondence. One of the hallmarks of theoretical physics of the last two decades is the discovery of the gauge-gravity duality [171,172]. The most prominent example of this is the celebrated N = 4 superconformal Yang-Mills theory which is conjectured to be dual to a certain type of the superstring theory [173]. One of its surprising features is that the Heisenberg spin Hamiltonian arises in the scalar sector as the one-loop approximation of the dilatation operator. The scattering matrix behind the scenes has an exceptional structure and turns out to be tightly related to the famous Fermi-Hubbard model and some other related models of strongly-correlated electrons dubbed as the Hubbard-Shastry model [174]. Constructions and physical applications of quasilocal charges have not yet been explored in this context.
Correlation functions. In this work we have not devoted any attention to the problem of calculating static and dynamic correlation functions of local observables, a task which typically represents an ultimate goal of any successful computational framework. A systematic procedure to encompass a wide range of interacting integrable theories in a universal and robust language still awaits to be developed. In this review, we have only addressed the problem of determining Bethe root distributions which parametrize a (non-thermal) equilibrium state. A mapping between the string densities and local correlators for the gapped regime of the XXZ model has been conjectured in Refs. [133,135,175]. An alternative route is to follow the Quantum Transfer Matrix approach [139,176,177] which was pursued in Ref. [124].

Beyond quasilocality
We have discussed at length the implication of pseudolocality of conserved quantities on several observable physical properties, such as ballistic (ideal) high-temperature transport and equilibration to non-thermal states. However, in some other rudimentary integrable models, a normal, diffusive spin or particle transport has been observed by numerical simulations, e.g. in the gapped Heisenberg model [144,145,150,178,179], or half-filled Fermi-Hubbard model [180][181][182].
Diffusive high-temperature transport in an integrable model can be considered as an indication of the absence of relevant pseudolocal charges, i.e. linearly extensive charges with non-vanishing overlap with a current operator. In the opposite case, the Mazur bound is strictly positive, implying ballistic conductivity. Even then, however, one may obtain other interesting bounds employing conserved operators with different volumescaling properties. For example, if there exists a conservation law Q with quadratic volume scaling of the Hilbert-Schmidt norm Q 2 HS = qN 2 + O(N ), then a rigorous derivation [183], in spirit very similar to the proof of Mazur bound for quantum spin lattice systems [89] but with appropriately balanced limits N → ∞, t → ∞, yields a rigorous lower bound on the diffusion constant Here j is a local current operator and v is the Lieb-Robinson velocity [61,184].
Simple examples of such bounds have been elaborated in Ref. [183] for the XXX and Fermi-Hubbard models, where Q in fact corresponds to a level-1 generator of Yangian symmetry [185]. Systematic exploration of quadratically extensive conserved charges in integrable systems and their applications to diffusive transport and quantum relaxation has not yet been undertaken.

Conclusions
This review is devoted to certain types of effective localities in the context of quantum integrable lattice models termed pseudolocality and quasilocality. The notion of locality indisputably plays a monumental role in the foundations of statistical mechanics, both on the classical and quantum level. We have presented and exemplified the meaning of quasilocal conserved quantities by discussing various applications of nonergodic phenomena in a paradigmatic interacting system, the anisotropic Heisenberg model. Specifically, we have elaborated on the importance of quasilocal charges in the description of generalized (non-thermal) equilibria on one hand, and their vital role in understanding certain anomalous transport characteristics such as divergent hightemperature spin conductivity on the other hand.
A key observation is that statistical ensembles, given by reduced density matrices which emerge in the steady-state limit after a relaxation process starting from any 'physical' initial state, are, due to effective dephasing, only capable of retaining a part of information about the initial condition which is encoded in local and pseudolocal conservation laws. Identification of a complete set of such charges provides us with a complete description of local correlation functions in generalized equilibria.
This naturally brings us to an elusive question which has been posted at the beginning, namely a controversial issue of the proper counting of degrees of freedom in an integrable lattice model. As we have explained, spectra of integrable lattice models in the thermodynamic limit organize in an astounding way and permit to cast our description in terms of stable quasi-particles [68,127,128,134]. This picture is in principle valid for any equilibrium state and even for elementary quasi-particle excitations on top of them [186]. Quasi-particles are labelled by a representation label (auxiliary spin in our example) and a continuous rapidity variable (corresponding to quasi-momentum).
Having this in mind it should not be difficult to understand why higher-spin transfer operators, despite fulfilling a system of functional identities, are nonetheless linearly independent variables. Therefore, the naive proposal of matching the number of degrees of freedom with the number of local Hilbert spaces of the lattice system cannot be correct.
We have furthermore discussed an interesting (although somewhat atypical) situation when the above picture is incomplete and needs to be appropriately extended. This happens when the underlying symmetry algebra becomes enlarged which implies extra degeneracies in the spectrum. Perhaps the simplest example of that occurs in the gapless regime of the XXZ model, governed by U q (sl (2)) quantum symmetry at root of unity deformations where an enriched symmetry led to the discovery of an extra family of quasilocal charges pertaining to non-unitary representations of the quantum group in the auxiliary space. In this review we exposed some of their implications on non-decaying currents and associated anomalous transport properties and presented a rigorous non-trivial bound on the singular contribution to the spin conductivity (Drude weight).
The last type of applications which we presented briefly were integrable spin chains subjected to Markovian dissipative boundaries. The time evolution in such cases is governed by a non-unitary process described by Lindblad master equation and generally leads to a unique steady state which is far from canonical thermal equilibrium. We owe to stress however that such 'integrable instances' which emerge as a consequence of an effective evolution describing an open system can be profoundly different from the conventional non-equilibrium settings in the scope of isolated systems which evolve according to the unitary evolution law and consequently the relevant class of symmetries which become important might be quite different. In addition, we notice that dissipation processes are strictly only well-defined in a finite volume while studies of equilibration in isolated systems typically deal with extended systems.
Aside of several novel theoretical insights which have been outlined in this review, it is worth mentioning that our formulation can also prove advantageous from a practical computational standpoint. An obvious example of that are explicit matrix product representations of quasilocal charges X s (λ) and Z(λ) which do not only admit a unified abstract representation but also enable a direct and efficient computation using methods from the standard statistical mechanics toolbox. In essence, this lifts the Bethe ansatz concepts to operator level right away in the thermodynamic regime, circumventing a long-standing challenge of achieving this by pursuing the programme of algebraic Bethe ansatz, see e.g. [187,188].
In conclusion, apart from a few successful physical applications in the realm of quantum quenches and quantum transport, much of the formal origin and grouptheoretic interpretation is still missing at the moment. A notable example is the question of completeness of the Z-charges and their reconciliation with the spectrum and the quasi-particle content. We hope that this review can provide a source of inspiration for the ongoing investigation of open directions.