Current operators in Bethe Ansatz and Generalized Hydrodynamics: An exact quantum/classical correspondence

Generalized Hydrodynamics is a recent theory that describes large scale transport properties of one dimensional integrable models. It is built on the (typically infinitely many) local conservation laws present in these systems, and leads to a generalized Euler type hydrodynamic equation. Despite the successes of the theory, one of its cornerstones, namely a conjectured expression for the currents of the conserved charges in local equilibrium has not yet been proven for interacting lattice models. Here we fill this gap, and compute an exact result for the mean values of current operators in Bethe Ansatz solvable systems, valid in arbitrary finite volume. Our exact formula has a simple semi-classical interpretation: the currents can be computed by summing over the charge eigenvalues carried by the individual bare particles, multiplied with an effective velocity describing their propagation in the presence of the other particles. Remarkably, the semi-classical formula remains exact in the interacting quantum theory, for any finite number of particles and also in the thermodynamic limit. Our proof is built on a form factor expansion and it is applicable to a large class of quantum integrable models.


I. INTRODUCTION
The description of the collective motion in many body quantum systems is one of the most challenging problems in theoretical physics. There are different possible levels for a theoretical treatment, ranging from the microscopic laws to various effective theories describing mesoscopic or macroscopic physics. For large enough systems one expects that classical behaviour will emerge, at least for certain observables. It is thus important to understand how and under what circumstances the various classical theories can be derived from an underlying quantum mechanical motion [1].
One such classical theory is hydrodynamics: it is known that many quantum systems admit some kind of hydrodynamic description on the mesoscopic and/or macroscopic scales [1,2].
Examples include Bose-Einstein condensates [3] or the quark gluon plasma [4]. Superfluidity is a famous exotic phenomenon, where frictionless flow is realized due to the constraints for the decay of excitations into lower energy modes. Superfluidity has been observed not only for liquid Helium, but also in ultracold bosonic and fermionic gases [5].
An other class of systems with exotic hydrodynamic behaviour is comprised by the one dimensional integrable models. In these models there exist independent conservation laws that constrain the dynamical processes, the number of which grows at least linearly with the volume. As an effect, these models do not thermalize to standard statistical physical ensembles. Instead, the emerging long time steady states can be described by a Generalized Gibbs Ensemble (GGE) that involves all higher conserved charges of the model [6,7]. The conservation laws prevent the decay of quasi-particle excitations, and this leads to dissipationless and factorized scattering, which was already demonstrated by experiments [8]. Dissipationless propagation of the collective modes leads to the emergence of ballistic transport and non-zero Drude weights (DC conductivity) [9].
Generalized Hydrodynamics (GHD) is a recent theory describing large scale non-equilibrium behaviour in integrable models [10,11] (see also [12][13][14][15]). The theory is built on the local continuity equations following from the conservation laws, which lead to a generalized Euler-type equation describing the ballistic transport. The GHD provides exact results for the Drude weights [16,17]. Diffusive corrections to the ballistic transport were also considered in [18,19], including an exact computation of the diffusion coefficients. Remarkably, the predictions of the GHD have already been confirmed in a concrete experimental setup [20].
Despite the successes of the GHD, one of the cornerstones of the theory has not yet been proven. The works [10,11] conjectured an expression for the expectation values of the current operators in local equilibrium, which is central to the derivation of the main equations of motion. Regarding integrable Quantum Field Theories a proof was provided in [10,21], whereas for the spin current of the XXZ model it was proven in [22]. Nevertheless, for arbitrary current operators in interacting lattice models or non-relativistic gases (the models most relevant to experiments) it was completely missing up to now. It is the goal of this paper to provide a proof of the conjecture, valid for a wide class of Bethe Ansatz solvable models.
The problem of the current mean values is also interesting from a purely theoretical perspective, without the immediate application to GHD. A large body of literature has already addressed equilibrium correlation functions in Bethe Ansatz solvable models [23][24][25][26][27][28][29][30][31][32], with particular interest devoted to the asymptotics of two-point functions (for a review see [33]), and to equilibrium mean values of arbitrary short range operators of the Heisen-berg spin chains (for a review see [34]). In contrast, the current operators are very specific short range objects, and as we will show, their finite volume mean values take a remarkably simple form. This has not yet been noticed in the Bethe Ansatz literature, and we believe that it deserves a study on its own right.
In the following subsection we describe the conjecture of [10,11] for the current mean values and explain its role in the GHD, while omitting many technical details. This brief introduction motivates our finite volume investigations.
After that, the remainder of this paper is composed as follows: In Section II we specify the problem and present our main new results, with a semi-classical interpretation given in III. Section IV includes our model independent proof, based on a finite volume form factor expansion. Section V includes known generalities about the charge and current operators in integrable lattice models; also, it gives a short summary of the Algebraic Bethe Ansatz. The proof of our form factor expansion for the XXZ and XXX spin chains is presented in Section VI. In Section VII we point out a connection to the theory of factorized correlation functions. Finally, we conclude in VIII.
A. Foundations of the GHD Isolated integrable models equilibrate to steady states described by Generalized Gibbs Ensembles (GGE's). Each GGE can be characterized by a set of parameters (generalized temperatures), or alternatively, by the mean values of all the local and quasi-local conserved charges in the model [7,35,36].
The time scales of equilibration to the GGE are set by the microscopic laws. It follows that in mesoscopic or macroscopic dynamical processes local equilibration happens much sooner than the characteristic times of the transport processes. This leads to the hydrodynamic description: GHD assumes the existence of fluid cells (regions in space much larger than the inter-particle distance, and much smaller than the variation of the physical observables) such that the state of each fluid cell can be described by a local GGE. The parameters of these local GGE's are then space and time dependent.
The Generalized Eigenstate Thermalization Hypothesis (GETH) [37] states that in local equilibrium the local observables only depend on the mean values of the conserved charges, and not on any other particular detail of the states. Thus, in order to describe the dynamical processes, it is enough to establish flow equations for the conserved charges, which will then determine all other physical observables on the hydrodynamic scales.
For simplicity let us consider here a continuum model in the thermodynamic limit. Let the complete set of local and/or quasi-local conserved charges be Q α = ∞ −∞ dx Q α (x), with α being an index or a multi-index, and Q α (x) being the charge density operators. Conservation of the charges implies that there exist current oper-ators J α (x) satisfying the equation of motion in Heisenberg picture: (I.1) GHD concerns the mean values of these relations: A closed set of flow equations can be obtained if the currents are expressed using only local information about the charges. In hydrodynamics this is performed in a derivative expansion: , t), . . . } β=1,2,... .

(I.3)
In the first approximation we neglect the spatial variations, and express the currents in the fluid cells using the mean values Q β (x, t) in that specific fluid cell only. This approximation describes the ballistic part of the transport. The diffusive part of the transport can also be treated by considering the derivatives ∂ x Q β (x, t) [19], but this will not be considered here. The ballistic flow equations will then follow from (I.2), given that one can compute the exact mean values of the currents in the local equilibrium states, as a function of the charges, or using any alternative description of the local GGE's. This is the problem that we consider in this paper.
A large class of integrable models is solvable by the Bethe Ansatz [38]; prominent examples include the Heisenberg spin chains or the 1D Bose gas with pointlike interaction. In these models the local equilibria in the thermodynamic limit can be characterized by the root densities ̺ n (λ) of the interacting quasi-particles, where n stands for a particle type and λ is the so-called rapidity parameter. The root densities can be understood as generalizations of momentum-dependent occupation numbers in free theory. Dissipationless scattering implies that these densities are well defined concepts even in the presence of interactions. The construction of the GGE is equivalent to specifying all the root densities ̺ n (λ), because they carry all information about the local equilibria. This is the ultimate form of the GETH, and it was understood in [36,39] following the earlier works [40][41][42][43][44][45].
In GHD we thus need to specify the Bethe root densities for each fluid cell, therefore they will depend also on the x, t coordinates of the cell. It is a very fruitful idea of the GHD that instead of concentrating on the equations (I.2) for the charges, one should derive flow equations for the rapidity distribution functions. This can be achieved starting from (I.2), by expressing both mean values in the continuity equations using the densities ̺ n (λ) only.
The mean values of the charges can be computed additively. In local equilibrium we have where q α (λ) is the single particle eigenvalue of the charges, and here we assumed only one particle species for simplicity. In GHD (I.4) is assumed to hold for each fluid cell separately. For the currents it was conjectured in [10,11] that the mean values can be computed using a semi-classical expression, namely by integrating over the carried charge multiplied by an effective propagation speed: Here the effective speed v eff (λ) is a generalization of the one-particle group velocity, which also takes into account the interactions between the particles. It has a physical explanation using a semi-classical argument: the oneparticle wave packets suffer time delays due to the scattering on the other particles, and these time delays accumulate along the orbit, eventually modifying the propagation speed. The resulting effective speed v eff (λ) is a collective property of the local GGE, because for each λ it also depends on the particle density ̺(λ ′ ) for all other λ ′ . For a precise definition of v eff (λ) we refer to [10,11]. Eq. (I. 5) has not yet been proven in Bethe Ansatz. It is our goal to fill this gap and to prove (I.5) starting from a rigorous finite volume computation.
Regarding the interpretation of (I.5) it was explained in [46], that the GHD can be simulated by the so-called "flea gas" model, which describes 1D motion of purely classical particles, subject to time delays (displacements) as an effect of interparticle collisions. Thus one observes a complete quantum/classical correspondence on the hydrodynamic scale.
In this work we show that the functional form of the mean values is the same in finite and infinite volume, therefore the quantum/classical correspondence holds with an arbitrary finite number of particles.

A. Elements of integrability
We consider integrable many body quantum models, including both lattice and continuum theories [23]. In this work we limit ourselves to those theories where particle number is conserved, and which can be solved by the traditional ("non-nested") Bethe Ansatz. Main examples are given by the Heisenberg spin chains, the 1D Bose gas, and also certain integrable QFT's [23]. Regarding other types of integrable models we give a few comments in the Conclusions.
In this Section we present formulas pertaining to lattice models. We consider an integrable Hamiltonian H in a finite volume L: (II.1) Here h(x) is the Hamiltonian density. In the most relevant cases h(x) is a two site operator, but we do not necessarily need this restriction. For simplicity we require periodic boundary conditions. In integrable models there exists a family of conserved operators Q α , where α is an index or multi-index. They mutually commute: and the Hamiltonian is a member of the series. We concentrate on the strictly local operators, which are given as where Q α (x) is a short range operator identified as the charge density. It is important that in a finite volume Q α is well defined for volumes larger than the range of Q α (x), and then the density Q α (x) does not depend further on L. Details on the canonical construction of the charges, and a few concrete examples will be given in V.
The Hamiltonian is a member of the series, thus the global operator Q α is conserved: We are interested in non-equilibrium processes and a hydrodynamic description, therefore we investigate the time evolution of the charge contained in a finite section. This leads to a continuity relation in operator form: (II.5) This relation defines the current operator J α (x) associated to Q α (x) under the time evolution of H. The existence of J α (x) follows simply from (II.4) and locality arguments.
The relations (II.3) and (II.5) do not define the Q α (x) and J α (x) operators uniquely; certain subtleties are discussed in Section V. We just put forward that the additive normalization of the current operators is chosen by the physical requirement that where |0 is the vacuum or reference state with no particles.
Our goal is to determine the mean values where |n is an arbitrary excited state of the finite volume Hamiltonian. In the thermodynamic limit these mean values will enter the flow equations (I.2).
In the models in question the exact eigenstates are found using the Bethe Ansatz [23,38]. The states are characterized by a set of lattice momenta {p 1 , . . . , p N } that describes the interacting spin waves.
The un-normalized Bethe wave function can be written for x 1 < x 2 < · · · < x N as [38] Ψ(x 1 , x 2 , . . . , (II. 8) Here each term in the sum represents free wave propagation with a given spatial ordering of the particles. The amplitude S(p j , p k ) = e iδ(pj ,p k ) is a relative phase between terms with different particle ordering, and it can be interpreted as the two-particle scattering amplitude. It depends on the model in question, and can be determined from the two-particle problem. The wave function is two-particle reducible: any multi-particle interaction is explicitly factorized into a succession of two-particle scatterings.
These wave functions describe eigenstates if they are periodic, from which we obtain the Bethe equations: It is useful to introduce the rapidity parametrization p = p(λ), where λ is the additive parameter for the scattering phase: (II.10) The Bethe equations can then be written as In the following the normalized N -particle Bethe states with rapidities {λ} N = {λ 1 , . . . , λ N } will be denoted as |{λ} N . The total energy and lattice momentum can be computed additively: where the single-particle energy e(λ) is a further characteristic function of the model.
Similarly, the eigenvalues of the conserved charges are where q α (λ) are the one-particle eigenvalues.
For later use let us write the Bethe equations in the logarithmic form: Here I j ∈ Z are the momentum quantum numbers, which can be used to parametrize the states.
In our derivation an important role will be played by the so-called Gaudin matrix G, which is defined as where now the I j are regarded as functions of the rapidities. Explicitly we have There are two interpretations of the Gaudin matrix. First, det G describes the density of states in rapidity space. This follows from the fact that in the space of the quantum numbers the states are evenly distributed, and G is defined as the Jacobian of the mapping from λ j to I j . Second, the Gaudin determinant also describes the norm of the Bethe Ansatz wave function in many integrable models (See Section V and [47,48]).

B. Main result
Our main result for the normalized mean values of the current operators is the following: (II. 18) Here the quantities e ′ and q α are N -dimensional vectors with elements and G −1 is the inverse of the Gaudin matrix.
In the simplest case of N = 1 the Gaudin matrix has a single element G 11 = Lp ′ (λ) and (II.18) gives the anticipated classical result where we introduced the bare group velocity v(λ) = e ′ (λ)/p ′ (λ) = ∂e/∂p. A similar semi-classical interpretation can be given also for higher particle numbers. Using the definition (II. 15) and the additive formula (II.12) we find the alternative expression where we defined the quantities In Section III it is explained, that the v eff can be understood in a simple semi-classical picture as effective velocities, describing the propagation of the individual bare particles in the presence of the others. It is remarkable, that the exact result and the functional form of the effective velocity are so simple in the finite volume situation. In the thermodynamic limit the papers [10,11] conjectured the formula (I.5) with the effective speed given by where ε(λ) and P (λ) are the so-called "dressed energy" and "dressed momentum". These are computed as the energy and momentum differences as we add a particle with rapidity λ into a sea of particles. The "dressing" takes into account the backflow of the other particles, which can be computed from the Bethe equations (II.14). The correspondence between (II.23) and our (II.22) is evident: small changes in the dressed momentum and dressed energy can be traced back to small changes in the momentum quantum numbers and the overall finite volume energy, respectively: (II. 24) This implies that (II.23) is indeed the thermodynamic limit of our (II.22). Furthermore, our (II.21) can be seen as the finite volume origin of the thermodynamic formula (I.5).
Eq. (II.18) is exact in those cases when the Bethe wave function is exact; this holds for integrable spin chains or the 1D Bose gas. On the other hand, in integrable QFT (iQFT) in finite volume the Bethe wave function is only an approximation, and in iQFT (II.18) holds up to exponentially small corrections in the volume. Depending on the model, the bare particles of the Bethe Ansatz can form bound states. These are described by the so-called string solutions of the Bethe equations [49]. It is important that our formula (II.18) is exact on the level of the individual Bethe rapidities, even in the presence of strings. An effective description involving the string centers (describing the rapidities of the composite particles) can be given afterwards using well established methods [50].
The main result (II.18) concerns the physical current operators, that describe the flow under time evolution by the physical Hamiltonian. However, it is also useful to consider certain generalized current operators, that describe the flow of a given charge under time evolution generated by some other charge.
Let us therefore consider two local charges Q α and Q β belonging to the same integrable hierarchy, implying that all three operators H, Q α , Q β commute with each other. We define J β α to be the current of the charge Q α under unitary time evolution dictated by Q β : Locality of the charge densities Q α,β (x) and the global relation [Q α , Q β ] = 0 implies that the operator equation (II.25) can always be solved with some short-range J β α (x). For the mean values of these generalized current operators we have the following result: Here q ′ β is an N -element vector with components q ′ β (λ j ), where q β (λ) is the one-particle eigenvalue of the charge Q β and the prime denotes differentiation. The analogy between (II.26) and (II.18) is evident: only the oneparticle eigenvalues of the time evolution operator are replaced. A special case and certain symmetry properties of this general statement are treated in Appendices A-B.
In the following Section we describe a semi-classical interpretation of these results, whereas the full quantum mechanical proof is provided in Section IV.

III. THE SEMI-CLASSICAL INTERPRETATION
Here we present a semi-classical computation which also gives a simple physical interpretation for the main result (II.18). Our arguments are very similar to those presented in [46], with the main difference being that we consider finite systems: we are looking at the motion of N particles on a finite ring of volume L (see Fig. 1). For convenience we consider here continuum models and a strictly point-like interaction.
Instead of solving the time-dependent Schrödinger equation we employ a semi-classical picture, namely we assume that the particles can be assigned well-defined straight orbits as long as they don't interact with each other. In this picture a particle with rapidity λ is represented by a wave packet, which travels with the group velocity v(λ) = de dp = e ′ (λ) p ′ (λ) . (III.1) In a typical situation all speeds are different and particles meet as they travel around the volume. The scattering events are taken into account by using exact quantum mechanical solution of the two body problem. For pure Fourier modes this would mean that for the scattering of particle j on k the wave function has to be multiplied by the phase S(λ j − λ k ); this is also reflected by the Bethe Ansatz wave function (II.8). This momentum dependent phase results in displacements of the center of the wave packets [51][52][53]. Such time delays are also present in classical integrable models, including models supporting solitons [54] or the hard rod gas [55]. After a displacement process the particles continue their path with their own bare speeds, until a further scattering event occurs. The time delays suffered in each of these events accumulate, and this alters the actual propagation speed of the wave packets, leading to the emergence of effective velocities v eff (λ j ). It is important that due to the higher conservation laws the multiparticle scattering events always factorize into a succession of two-particle scatterings [56], which also implies that the particle rapidities never change during time evolution.
In this semi-classical picture the current mean values are evaluated simply as Our goal is to find the emerging effective velocities. We consider long times such that each pair of particles has scattered on each other many times. In this limit the particular order of the individual scattering events (which depends on the initial positions of the particles) becomes irrelevant. The spatial displacement of particle j caused by the scattering on the particle k (j < k) is given by the derivative [51][52][53] where we used the rapidity parametrization and ϕ(u) is defined in (II.17). This formula is valid when particle j overtakes particle k from the left, i.e. when v j > v k . We can formally extend it as The time elapsed between the two scattering events of the same two particles j and k can be expressed as: For asymptotically long times t the accumulated displacement that particle j suffers is given by: This causes the difference between the effective and bare velocities: Putting everything together we obtain the self-consistent relation (III.9) As an effect of the extra sign σ jk in (III.4) it can be written as (III.10) Multiplying with Lp ′ (λ j ) and using the definition of the bare group velocity we get (III.11) On the r.h.s. we can recognize the action of the Gaudin matrix (II.16), we can thus write (III.12) Multiplying with G −1 , and substituting the v eff (λ j ) into (III.2) leads to which is identical to the full quantum mechanical result (II.18). We have thus demonstrated a complete quantum/classical correspondence with a finite number of particles.
In our derivation we assumed that v eff (λ j ) > v eff (λ k ) whenever v(λ j ) > v(λ k ). However natural this requirement may seem, there can be situations where it does not hold [57]. In those cases the quantum result remains valid, but the semi-classical picture can not be applied.

IV. PROOF USING A FORM FACTOR EXPANSION
Here we present a model-independent proof of our main result for the current mean values. Our technique relies on a finite volume form factor expansion theorem. The proof of this expansion theorem can depend on the particular model, but the computations of this Section are quite general.
In the most part we will concentrate only on the physical currents J α . The generalized currents defined in (II.25) will be considered in IV C.
The starting point is to use the definition of the current operators (II.5), and to consider matrix elements of this operator relation. We intend to compute the mean values of the currents, but taking the mean values of (II.5) gives automatically zero on both sides, due to the Bethe states being translationally invariant and eigenstates of H. Instead, it is useful to take the off-diagonal matrix elements between two Bethe states with non-equal total lattice momentum: In integrable models generic off-diagonal finite volume matrix elements of local operators can be expressed as where the function F O ({λ} N |{µ} M ) is the so-called form factor and det G λ and det G µ describe the norm of the Bethe Ansatz wave function, or alternatively the density of Bethe states in rapidity space. They are N × N and M × M Gaudin determinants computed from the sets of rapidities {λ} N and {µ} M . The form factors describe the transition matrix element for the un-normalized Bethe wave functions (VII.5). They are meromorphic functions and they are completely independent from the volume. The volume dependence of the physical matrix elements only comes through the solution of the Bethe equations and the normalization factors. The properties of the form factors have been investigated both in QFT [58] and using Algebraic Bethe Ansatz (ABA) [23]. A derivation of the analytic properties in the Lieb-Liniger model was also given using coordinate Bethe Ansatz in [59]. The statement (IV.2) is well known in the literature dealing with integrable lattice models [23], and for integrable QFT it was first written down in [60].
As opposed to the transition matrix elements the mean values of local operators in Bethe states can not be expressed directly using the infinite volume form factors.
The reason is the appearance of the so-called disconnected terms: on a mathematical level they arise from the kinematical poles of the form factors, whereas their physical interpretation is that they describe processes when a subset of the particles does not interact with the local operator. On the other hand, those processes when some of the particles do interact with the operator are described by certain diagonal limits of the form factor functions.
In order to describe the finite volume mean values, let us define the so-called symmetric evaluation of the diagonal form factors as There is an other often used diagonal limit, called the connected form factors, but they will not be used in this work and for a thorough discussion we refer to [61].
It is useful to define the functions ρ N (λ 1 , . . . , λ N ) as the N × N Gaudin determinants evaluated at the set of rapidities {λ 1 , . . . , λ N }. In the notations we suppress the index N and write simply We remind that the Gaudin determinants describe the norms of Bethe wave functions for eigenstates, ie. for sets of rapidities satisfying the Bethe equations. On the other hand, the functions ρ({λ} N ) are defined for arbitrary sets of rapidities. It is useful to write down the first two ρ({λ} N ) functions. For N = 1 we have simply and its determinant is Our proof for the current mean values will be based on the following expansion theorem.  This theorem was first formulated in [61] for integrable QFT, where the Bethe Ansatz for the finite volume eigenstates is not exact due to the presence of virtual particles. Therefore, in iQFT the theorem holds up to corrections exponentially small in the volume. On the other hand, it is an exact relation in non-relativistic models including the 1D Bose gas and the Heisenberg spin chains. We believe that the theorem is new in the case of the XXZ and XXX spin chains, and it will be proven rigorously in Section VI. Nevertheless, similar theorems had been known for particular cases in the ABA literature [23].
In [62] it was shown that the LeClair-Mussardo (LM) formula (an integral series developed for thermal mean values) can be considered a thermodynamic limit of this expansion theorem, whereas in [63] it was shown that alternatively (IV.8) can be derived from the LM formula using certain analytic continuations. In [64] it was also shown that in iQFT it can be derived directly using the off-diagonal relation (IV.2). Regarding the continuum gas models, the expansion was proven for certain local operators in the Lieb-Liniger models in [65].
We note that there is an alternative expansion theorem using the connected form factors [61], but it will not be used here.
The continuity relations yield a connection between the symmetric form factors of the charge and current operators. Introducing the short-hand notations we get from (IV.1) (IV.10) Our strategy is the following: First we find the symmetric form factors of the charge densities by comparing the formula of the expansion theorem 1 to the known mean values (II.13). Next we use the above relation to find the J α . Finally we use Theorem 1 for the current operators: we sum up the resulting expansion to obtain (II.18). Essentially the same strategy has already been applied in [10,21] directly in the thermodynamic limit, using the LeClair-Mussardo series. The novelty of our approach is that we perform these steps in finite volume, and that we also provide the proof of our expansion theorem for the XXX and XXZ spin chains (see Section VI).

A. The form factors of the charge densities
We consider (IV.8) in the case of the charge density operator Q α (0) and write it as (IV.11) These are algebraic relations that hold for any particle number N and any finite volume L. The symmetric diagonal form factors can be extracted using a recursive procedure: we consider the above algebraic relations for N = 1, 2, . . . , and at each N we compute the N -particle form factor by subtracting the terms with all Q α ({λ + }) with a lower number of particles, obtained earlier.
At N = 1 the relation immediately gives where we used (IV.5) and substituted Q α (0) = 0. At N = 2 we use (IV.5)-(IV.7). Substituting them into (IV.11) and using also (IV.12) we observe the cancellation of some terms, leading eventually to (IV.13) This procedure does in principle yield all higher symmetric form factors. However, the direct recursive subtractions for N ≥ 3 become more involved and it is advantageous to use an alternative method. For the computation of the Gaudin determinants we will apply a graph theoretical matrix-tree theorem, which has already proven to be useful for different problems [21,66,67].
Let us introduce the following definitions. Given a graph Γ a directed graph F is a directed spanning forest of Γ if it satisfies the following: • F includes all vertices of Γ.
• F does not include any circles.
• Each vertex has at most one incoming edge.
The nodes without incoming edges are called roots. Each spanning forest can be decomposed as a union of spanning trees, which are the connected components of the forest. It can be seen from the above definitions, that each spanning tree has exactly one root. (IV.14) In this case the determinant of G can be expressed as a sum over the directed spanning forests of the complete graph with N nodes. For each spanning forest F let R(F ) denote the set of the roots, and let l jk denote the edges of F pointing from node j to k. Then we have For a proof see for example [68]. The Gaudin determinant given by (II.16) satisfies the requirements of this Theorem with D jj = p ′ (λ j )L and K jk = ϕ(λ j − λ k ). The main idea to obtain the form factors from (IV.11) is to consider the formal L → 0 limit of this relation. To do this we need to observe the L → 0 limit of the various Gaudin determinants. Each term in (IV.15) carries a factor of L r where r ≥ 1 is the number of roots for the particular spanning forest F . On the r.h.s. of (IV.11) we have for every non-empty set of {λ − }. Thus the only term on the r.h.s. which survives the L → 0 limit is the one where On the other hand, on the l.h.s. we need to keep the O(L) term in ρ({λ}), which according to the above theorem gives an expansion over all directed spanning trees F ′ : where for each spanning tree F ′ the index r denotes its root vertex, and we used the abbreviation ϕ jk = ϕ(λ j − λ k ). Each directed spanning tree can be obtained uniquely from a non-directed spanning tree by selecting the single root vertex and choosing the directions of the edges accordingly. In our case the function ϕ is symmetric, thus the direction of the edges does not influence the factors of ϕ(λ u − λ v ). Each vertex has to be chosen exactly one time as a root, thus we obtain where the summation runs over the non-directed spanning trees T . Finally

B. Summation for the current operators
Eqs. (IV. 19) and (IV.10) yield the symmetric diagonal form factors of the current operators: (IV.20) Our task is to sum up the expansion for the currents .
Once again it is instructive to consider the first few cases. At N = 1 we have Using J α (0) 0 = 0 and (IV.5) gives immediately as anticipated. At N = 2 there are 3 non-zero terms in the summation: (IV.24) The two-particle symmetric form factor is (IV.25) Substituting this and also (IV.22) and (IV.5) into (IV.24) we can express the mean value in the product form .
(IV.26) We can recognize the inverse of the two-particle Gaudin matrix (IV.6). Thus in this case we have obtained (IV.27) as stated in our main result (II.18).
The summation for N ≥ 3 is considerably more involved. From the structure of the form factors and the Gaudin determinants we can see that for each pair jk of particles there will be various contributions including the factors e ′ (λ j )q α (λ k ), stemming from different terms in (IV.21). Our main result (II.18) states that the sum of these terms will reproduce the jk element of the matrix G −1 . For this inverse matrix we can use the formula where adj(G) is the so-called adjungate matrix of G. The Gaudin determinant is present in the denominator of (IV.21), thus (II.18) holds if in the nominator the sum of the terms with e ′ (λ j )q α (λ k ) reproduce the jk components of adj(G). The proof of this is not trivial, and it is presented in Appendix C.

C. Mean values of the generalized currents
It is straightforward to repeat the previous calculations for the case of the generalized currents J β α defined in (II.25). For the symmetric diagonal form factors we will use the notation J β α ({λ} N ). The local continuity equation then leads to (IV.30) These form factors are needed to sum up the expansion theorem (IV.8). All of the previous computations can be applied by exchanging the function e ′ (λ) with q ′ β (λ). This has a simple interpretation: for the general currents the time evolution is dictated by Q β instead of the physical Hamiltonian. Performing all the previous steps we obtain our result (II.26).

V. CURRENT OPERATORS AND ALGEBRAIC BETHE ANSATZ IN INTEGRABLE LATTICE MODELS
Here we review the canonical construction of the conserved charge and associated current operators; we also sketch the standard method of the Algebraic Bethe Ansatz to find the eigenstates. We will concentrate on lattice models that are obtained from integrable Lax operators. The treatment below is rather general, with concrete examples given later in V B.
Let H = V 1 ⊗ V 2 ⊗ · · · ⊗ V L denote the Hilbert space of the model, with V j ≈ C D with some D = 2, 3, . . . and let V a ≈ CD with someD denote the so-called auxiliary space. In our main examplesD = D, but this is not necessary. Let L(u) denote the so-called Lax operator acting on V j ⊗ V a . Here u ∈ C is the spectral parameter.
The monodromy matrix T (u) acting on V a ⊗ H is defined as The transfer matrix is given by the trace in auxiliary space, which corresponds to enforcing periodic boundary conditions: The local Lax operators satisfy the exchange property where we introduced two auxiliary spaces V a,b and R(u) is the so-called R-matrix acting on V a ⊗ V b . It satisfies the Yang-Baxter relation which is a relation of operators acting on the triple product V 1 ⊗ V 2 ⊗ V 3 of auxiliary spaces.
As an effect of these relations the transfer matrices form a commuting family of operators [23]: The transfer matrix encodes the hierarchy of the conserved charges, which are obtained by expanding t(u) around certain special points. We consider models where the dimensions of the physical and auxiliary spaces are equal and the local Lax operator can be chosen to be identical to the R-matrix: L(u) = R(u). Furthermore, we will consider cases where the R-matrix satisfies the initial condition R(0) = P with P being the permutation operator, such that the transfer matrix satisfies where U is the translation operator by one site. For any α ∈ N, α ≥ 2 we then define These charges will be extensive and they can be written as a sum over the charge density operators [69]: With this definition Q α (x) spans α sites, and for α = 2 we have with H being the physical Hamiltonian and κ being an L-independent factor, which depends on the particular conventions that are used. If the spectral parameter is chosen appropriately, then the canonical Q α defined above are all Hermitian. We consider models with particle number conservation. In these cases there is a reference state |0 with zero particles present. We require that the overall normalization of the transfer matrix satisfies 0|t(u)|0 = 1, (V. 10) leading to Additional multiplicative factors in t(u) would alter the additive normalization of the charges.
Having found the charge densities, the continuity relations (II.5) and (II.25) uniquely define the current and generalized current operators J α and J β α . It follows from (V.9) that for β = 2 we can identify J 2 α = κJ α .
We note that the relations (V.7)-(V.8) do not define Q α (x) unambiguously: for any choice Q α (x) we can take a local operator D(x) and define an alternative density which leads to the same integrated charge. This transformation can be considered as a "gauge freedom" for the definition of the charge densities, and it is discussed in detail in [19]. This "gauge choice" does not alter the charge mean values, but it changes the definition of the current operators as The additional terms do not affect the mean values of Q α (x) and J β α (x). An alternative way of constructing the charges is with the help of the boost operator [70][71][72][73], which is defined on the infinite chain as the formal expression where the density of the charge Q 2 is simply A formal manipulation shows that [72,73] where t ∞ (λ) is the transfer matrix of an infinite chain. It follows that The additional constant parts are not fixed by the formal computations, and need to be adjusted afterwards. Defining a recursion asQ α+1 = i[B, Q α ] and using (V.11) we get the correct normalization It is also possible to compute the current operators J β α (x) using a generalization of the boost operator. A formal application of (II. 25) gives This relation can be used to obtain J β α (x), but depending on the situation the direct application of (II.25) might be more efficient.

A. The Bethe Ansatz solution
Let us now focus on the case of D =D = 2, furthermore we consider models with U (1) symmetry where the fundamental R-matrix is of the form Specific examples will be given later in V B. The monodromy matrix defined in (V.1) is usually written in the block form where the blocks correspond to the degrees of freedom in auxiliary space, and A(u), B(u), C(u), D(u) are operators acting on the spin chain. It follows from the local relation (V.3) that the monodromy matrix satisfies where a, b refer to two different auxiliary spaces. Commutation relations between the A, B, C, D operators can be derived from (V.22); they are listed in Appendix D.
In our models possessing U (1) symmetry there exists a reference state |0 which is annihilated by all C(u) for all u. Typically |0 is chosen as the state with all spins up. Then the Bethe states can be created in Algebraic Bethe Ansatz as Here σ is a constant which is chosen later such that rapidity parametrization becomes exactly the same as in coordinate Bethe Ansatz. These states are eigenstates of the spin chain Hamiltonian if the rapidities satisfy the Bethe equations where we introduced the function , For on-shell states of the physical chain these are the adjoints of the states (V.23) (for detailed proof see [23]), but in the more general case including certain inhomogeneity parameters (leading to non-Hermitian Hamiltonians) they are only dual vectors. For eigenvectors the norm of the Bethe states is whereκ is a model-dependent constant. This statement was proven in [48] based on the singularity properties of general overlaps. Finally, the eigenvalues of the transfer matrix are

B. Main examples
The SU (2)-invariant XXX Heisenberg spin chain is given by the Hamiltonian where σ a j , a = x, y, z are the Pauli matrices acting on site j. This integrable model can be obtained from the R-matrix of the form (V. 20 the shift parameter is σ = i/2 andκ = 1. The canonical definition (V.7) of the charges gives Q 2 = H/2 [23]; the additive normalization is such that (V.11) is satisfied.
The eigenstates of the model organize themselves into SU (2)-multiplets, the Bethe Ansatz gives the highest weight vectors. The wave functions are of the form (II.8) with the rapidity parametrization given by (V.33) The one particle energy eigenvalue is e(λ) = 2q 2 (λ) with The rapidity parameters take values in the whole complex plain. The solutions of the Bethe equations organize themselves into strings, which describe bound states of spin waves [49]. These bound states can be regarded as different particle types in the thermodynamic limit. However, we perform our finite volume analysis on the level of the individual Bethe roots, therefore we do not treat the string solutions separately. In this model the canonical charges Q α defined by (V.7) can be computed explicitly [74]. For the sake of completeness we present the explicit formulas up to Q 4 in Appendix B; more explicit results are found in [74].
Regarding the one-particle eigenvalues of the charges it follows directly from (V.30) that The transfer matrix commutes with the global SU (2) transformations, therefore the Q α with α ≥ 2 are all SU (2)-invariant operators. The global spin operators are additional conserved quantities, and the traditional choice is to add Q 1 = S z into the commuting family. If the vacuum is chosen as the reference state with all spins up, then the one-particle eigenvalues of Q 1 are simply q 1 (λ) = −1.
Explicit real space formulas for the current operators of the XXX model are not available. Based on [74] it seems plausible that closed form results can be computed, but we have not pursued this direction. Nevertheless we computed the first few currents and generalized currents, the results are presented in B.

VI. PROOF OF THE EXPANSION THEOREM -THE XXZ CHAIN
Here we prove the Theorem 1 using standard methods of the Algebraic Bethe Ansatz (ABA) [23]. In fact, our proof can be considered a generalization of the proof given by Korepin for the norms of the Bethe Ansatz wave functions [48]. Similar ideas have been worked out by one of the authors in the work [62], which considered certain local correlators of the continuum 1D Bose gas. Here we restrict ourselves to the case of the XXZ spin chain, with the R-matrix given by conventions (V.20)-(V.37). The case of the XXX chain can be treated similarly.
For the proof it is useful to define re-normalized operators In order to shorten the notations, in this Section we do not use the rapidity shift −iη/2 that appeared earlier in (V.23).
Our aim is to derive the form factor expansion for the normalized mean values where O is any operator of the finite spin chain. Quite interestingly, for this proof we do not require any locality property from the operator. In fact, locality of an operator is not even a well defined concept in a finite chain. Thus we don't impose any restriction on the range of the operator O.
To be precise, let us consider the elementary matrices E ab , a, b = 1, 2 and let E ab (x) stand for operator which acts as E ab on site x = 1, . . . , L and with the identity elsewhere. We perform the proof for an arbitrary product Each operator of the finite chain is a linear combination of these products. Short range operators are obtained by taking traces in some subset of the indices a x , b x .
In order to compute these mean values we need to embed the operators (VI.3) in the Yang-Baxter algebra. This procedure is called the "quantum inverse scattering problem" and was solved in [75][76][77].
In the case of the homogeneous chain we have Our strategy is that we prove the theorem for an arbitrary product where X(µ) may represent one of the four operators A, B, C, D evaluated at spectral parameter µ. Afterwards we take the µ j → 0 limit, and by continuity we obtain the statement for (VI.4).
Our proof is based on the singularity properties of the matrix elements of operators. We follow closely the proof of Korepin for the norms of Bethe states [48]; in fact, our proof can be considered a slight generalization of the methods of Korepin. For an earlier similar proof see [62]. We introduce the following notation for a matrix element of the arbitrary operator O between two states (not necessarily eigenstates) described by the rapidity (VI.6) These matrix elements depend on the rapidities and the variables l(λ) = a(λ)/d(λ), where a(λ) and d(λ) are the vacuum expectation values of the operators A(λ) and D(λ) respectively [23]. For what follows, it is important, that we can consider arbitrary functions for a(λ) and d(λ) and therefore M O N is the function of 4N independent variables. It follows from the commutativity of the B and C operators that M O N is invariant with respect to simultaneous exchanges λ B j ↔ λ B k , l B j ↔ l B k , and similarly for the rapidities on the left hand side.
The matrix elements have apparent singularities as two rapidities from the two sides approach each other. These apparent poles result from the commutation relations between the B and C operators. We will show that the structure of these poles completely determines the mean values. For simplicity we will focus on the singularities as λ C N − λ B N → 0; the other cases follow simply from the permutation symmetry.
Theorem 3. The matrix elements satisfy the following singularity property: where in the third line the matrix element is calculated with the modified vacuum expectation values The proof is rather technical and it is presented in Appendix D.
In the physical case (when the l B,C N variables are not independent) this expression gives a finite result where z(λ) = i∂ λ log l(λ).
The diagonal evaluation of the matrix element is defined as: where the limit is performed for every k = 1, . . . , N . This quantity depends on 3N independent variables {λ} N , {l} N and {z} N . From (VI.7) it follows that the dependence on z N is linear, and the proportionality factor is To calculate the expectation values in the eigenstates of the system, we have to take the rapidities to the solutions of the Bethe equations. The quantity obtained after this does not depend on the l-parameters anymore, since they are given by the Bethe equations: The dependence on z N is still linear: (VI.14) To continue the calculation we need to introduce the form factors F O N and F O N,s (these differ from the previously used form factors only in an overall normalization) where it is understood that both sets of rapidities solve the Bethe equations, i.e. we express the {l C } N , {l B } N using the rapidities. From (VI.7) it is obvious that this form factor satisfies the following recursion relation (VI. 16) By taking the symmetric, diagonal limit of this quantity, one obtains the symmetric form factor Theorem 4. The symmetric form factor of the operator O is equal to its expectation value in the case when every z j is zero: Proof. From (VI.7) it is obvious that the z dependence of the expectation value arises from the rapidity dependence of the l(λ) function. This means that the z independent, irreducible part can be obtained by choosing On the other hand the symmetric form factor is by definition: This completes the proof.
We define the S function for an arbitrary bi-partition of the set {λ} N = {λ + } n ∪ {λ − } N −n in the following way: . This function depends on z N only through the Gaudin determinant, and its dependence can be calculated easily by expanding the determinant with respect to its N th row or column. Doing so one obtains that It is important to notice that if we take all the z-s in the argument of S to zero, we get zero: This follows from the fact that in this case the Gaudin determinant is zero, which follows from Theorem 2.
With the help of the relations listed above we are now in the position to prove the following theorem: Theorem 5. The un-normalized mean value of an arbitrary operator O in any eigenstate of the system can be calculated in the following way: where the summation goes over every bi-partition of the set of rapidities, {λ} N and for simplicity we did not denote the cardinality of the sets {λ ± } separately.
Proof. We use induction in the particle number N . Let us consider both sides as the multi-linear function of the z j variables. Our goal is to show that both sides depend the same way on every z j and that their z independent parts are also equal. We look at the first N for which the matrix element defined by (VI.6) is not zero, and we denote this number by N min . In the case when N < N min both sides of (VI.25) are zero, therefore the equation is satisfied. If N = N min there is only one non-zero term in the summation in the r.h.s., namely when {λ + } = {λ} N . This means that on the r.h.s. there is only F O N,s ({λ} N ), since S({λ + }, ∅, ∅) = 1, therefore it is z independent. The l.h.s. is also independent of z, which follows from (VI.14). But the z independent parts are equal according to the previously proved theorem. Now let us suppose that (VI.25) is satisfied for every N < M , and examine the N = M case! In the r.h.s. only those terms depend on z j where z j ∈ {λ − }. By taking the partial derivative of it with respect to z j , the initial summation is getting modified to a new one, going over all the bi-partitions of the set {λ} N −1 = {λ} N \{λ j } (this follows from (VI.23)). According to the induction assumption this sum gives exactly O N −1 ({λ} N −1 , {z mod } N −1 ), multiplied by sinh(η) k =j f kj f jk . But using (VI.14), this demonstrates that the z dependence of the two sides are equal. To investigate the part independent of the z variables, we have to take all of them to zero. In this case on the r.h.s. only

VII. CONNECTION TO THE THEORY OF FACTORIZED CORRELATION FUNCTIONS
In the case of the XXZ and XXX spin chains our results for the current mean values are directly related to certain objects in the theory of factorized correlation functions. In the following we describe this connection.
Factorization of correlation functions concerns the equilibrium mean values of local operators (for a review, see [34]). Factorization means that any multi-point correlator can be expressed as sums of products of simple building blocks, which are derived from the two-site density matrix in an inhomogeneous spin chain. On a practical level the theory consists of two parts: the algebraic part and the physical part. The algebraic part deals with the factorization on the level of the operators, and is independent of the actual physical situation. On the other hand, the information about the concrete situation is supplied by the physical part of the computation.
The theory was worked out first for the cases of thermal equilibrium in infinite volume and for the ground states in finite volume [78][79][80][81]. In [82] a conjecture was formulated for the physical part in any excited equilibrium state in infinite volume; these formulas have been used to compute the steady state properties after global quenches. Finally, it was argued in [83] that the known results for the finite volume ground states [81] can be extended naturally to all finite volume excited states, and this leads to the proof of the formulas of [82] in the thermodynamic limit.
In the following we summarize the main statements in the case of the XXX chain.
The Bethe vectors are highest weight vectors with respect to the SU (2) symmetry. Let us consider a Bethe state with S z = 0, thus N = L/2. All such states are SU (2) singlets.
The theory of factorized correlations states [81,83] that in the SU (2) singlet states any multi-point correlation function can be expressed using only a single generating function Ψ(x 1 , x 2 ) that depends on the excited state in question. Defining the coefficients it can be shown that all correlators can be expressed as finite combinations of the quantities Ψ n,m . The alge-braic part of the computation expresses the correlators in terms of Ψ n,m , whereas the physical part supplies their actual values, depending on the Bethe state in question. For example the simplest z − z correlators can be expressed as It was shown in [83] that for the Bethe state |{λ} N the generating function reads where q(x 1 ) is a parameter dependent vector of length N with elements q j (x) = q(λ j −x) with q(λ) = 1/(λ 2 +1/4). Furthermore, G is the Gaudin matrix, which now takes the form The factor of 2 is included in (VII.5) to conform with earlier conventions, see [83]. It follows that the Ψ n,m coefficients can be expressed as Ψ n,m = 2q n+2 · G −1 · q m+2 . (VII.8) The shifts in the indices are due to our conventions, namely that the first member of the series of the charges is called Q 2 .
Comparing to (II. 26) we see that the current mean values are This gives a new interpretation for the generalized currents: they are special operators that are represented by a single Ψ n,m .
In the case of non-singlet states the situation is more involved, and the correlators also involve the so-called moments (see [83,84]). Nevertheless, the mean values of SU (2) invariant operators are still described by the Ψ n,m , and this is in accordance with the fact that in the XXX model the charges Q α and their currents are also SU (2)-invariant.
In the case of the XXZ model spin-flip invariant local correlations are described by two generating functions (traditionally denoted as ω(x, y) and ω ′ (x, y), see [78][79][80][81]). It was shown in [83] that for the finite volume excited states the function ω(x, y) has a form analogous to (VII.5), thus the Taylor coefficients of this function describe the generalized current operators in the XXZ model.
At present we don't have an explanation of the observed coincidences between the factorization and the current mean values. Understanding this connection might lead to an independent proof of our results, at least in the XXZ and XXX models.

VIII. CONCLUSIONS AND DISCUSSION
We computed an exact finite volume formula for the current mean values in Bethe Ansatz solvable systems. The main results are eqs. (II.18) and (II.26). We did not treat the direct thermodynamic limit of these results, but their functional form implies that they reproduce the previously conjectured infinite volume formulas (see eq. (II.21) and the discussion there). We have thus supplied a rigorous foundation for the treatment of ballistic transport in Generalized Hydrodynamics.
Perhaps the most interesting finding of this work is that the semi-classical result for the currents remains exact in the interacting quantum theory, even with a finite number of particles. Ultimately this boils down to the two-particle reducibility of the Bethe Ansatz wave function (see eq. (II.8)), which also enables the applicability of the flea gas model to simulate the dynamics [46].
Our rigorous proof is rather general and it relies on a model independent form factor expansion (Theorem 1). However, the expansion itself has to be proven separately and the techniques to be applied might vary. Here we provided a proof for the case of the Heisenberg spin chains. Whereas we did not treat the Lieb-Liniger model (1D Bose gas) here, the expansion theorem can be worked out with essentially the same techniques (see for example [62]), at least for those local charges and currents which have well-defined real space representations [85]. For integrable QFT the expansion was proven earlier in [64].
It would be highly desirable to develop an alternative, more transparent proof for the current mean values. For special cases this can indeed be achieved, for example the spin current of the XXZ model can be computed using form factor perturbation theory [86], using similar ideas to those of [22]. Nevertheless a simple and general proof is not available. In the special case of the XXZ and XXX models the connection to the factorized correlation functions might lead to a more transparent derivation.
In this work we restricted ourselves to local charges on the lattice. However, it is known that there exist quasilocal charges that are essential for the GGE and thus for GHD [35,36]. For quasi-local charges our results hold asymptotically, and the finite volume formulas (II.18)-(II.26) receive exponentially small corrections. These are due to the long-range contributions to the operators with exponentially decreasing amplitudes.
It would be interesting to extend our results to multicomponent models solvable by the nested Bethe Ansatz. In these systems much less is known about local correlations, and the current operators are good candidates to obtain exact analytic results. In turn, this would give a rigorous foundation for the hydrodynamic treatment of these models [87]. Also, it would be interesting to consider models without the U (1) symmetry responsible for particle conservation. The integrable XYZ spin chain is such a model, where the canonical family of conserved charges [73] does not include the S z operator. Nevertheless there might exist simple exact results for the current operators, and thus also for GHD. We plan to return to these questions in future research. Here we investigate a special case of the main formulas (II.18)-(II.26): we consider the energy current J H , which is itself a conserved operator in integrable lattice models In fact, it is proportional to the canonical charge Q 3 . For completeness we re-derive this correspondence, and show that (II.18) reproduces the mean values of Q 3 .
The continuity equation for the energy flow is where h(x) is the Hamiltonian density. It follows that Using the definition of the boost operator (V.14) and the proportionality relation (V.9) we compute the canonical Q 3 as Therefore we can identify The proportionality factor −κ 2 originates simply in our conventions.
The one-particle eigenvalues of Q 3 are obtained from the transfer matrix construction: (V.35) and (V. 34) give For J H (II.18) takes the form We will show that after re-scaling this formula gives the expected mean value of Q 3 . Let us take an N -dimensional vector u whose elements are equal to 1. It follows from the definition of the Gaudin matrix (II.16) that Multiplying with the inverse and using e(λ) = q 2 (λ)/κ = −p ′ (λ)/κ we get With this we have indeed obtained as expected from (A.4). Eq. (A.9) is a special case of a more general symmetry relation. It follows from (V.35) and the symmetry of the Gaudin matrix, that the following mean values are equal: This interesting relation was already observed in GHD in [12], but a more direct explanation is not yet known. We stress that the relation above does not mean that the operators J β α (x) and J α−1 β+1 (x) are equal; they might differ in total derivatives. The concrete form of the current operators depends on the choice of the charge density representing the global charge operator, and the gauge freedom (V.12) leaves room for re-definitions. It would be interesting to see whether there is a gauge choice, which would guarantee that the above relation holds on the level of the operators. In the XXX Heisenberg chain the canonical Q 2 operator (V.7) is Here we used the short-hand notation that σ x is a vector of operators (σ x , σ y , σ z ) acting on site x. Note that Q 2 is one half of the Hamiltonian. Further charges can be computed easily using the boost operator (V.14). After explicit computations we find where the cross denotes the vectorial cross product. Further examples and a closed form result for all Q α can be found in [74] Let us define the generalized current operators J β α according to the definition (II.25), with the charge densities as given above. Then the real space representations for the first few currents are We can observe the equality J 2 2 (x) = −Q 3 (x), in accordance with the previous section. Similarly, we find the interesting relation which is an analogous identity that follows from (A.10); the minus sign is simply a matter of convention.

Appendix C: The summation of the form factor expansion
In order to sum up the expansion for the current mean values we write the r.h.s. of (IV.21) in a slightly different way where J α ({λ}) are the symmetric diagonal form factors of the currents and we defined Here the summation runs over T , the non-directed spanning trees of {λ + }. Using this notation our task is to show that This matrix product can be written out explicitly using the matrix-tree theorem for the Gaudin determinant where F denotes the directed spanning forests of {λ − }. First we show that the diagonal elements of this matrix product give det G. To do this we consider the j = k case in (C.3) and we split the outer summation over l into two parts: to the l = j and l = j cases.
In the case of l = j the sum of the ϕ terms from the Gaudin matrix appears, multiplying the appropriate terms in A. These can be written in the following way: Here in the second line we denoted the terms inside the summation with (. . . ) for brevity. Using this equation and renaming the summation variable l to s in (C.3), we can finally arrive to the following expression for the diagonal elements of the matrix product (C.3): Looking at this expression from the graph theoretical point of view, and using the matrix-tree theorem it can be shown that this is indeed det G: in the first term for every bi-partition the summation over T (together with the p ′ j L factor) gives the contribution of every such directed spanning tree of {λ + } in which λ j is the root. Together with the summation over F , which is just the contribution from every directed spanning forest of {λ − }, the first term contains all such directed spanning forests of {λ} in which λ j is one of the roots. On the other hand in the second term the summation over T gives the contributions from the free (rootless) spanning trees of {λ + }. These are ,,connected" to one of the spanning trees in the spanning forests of {λ − } by the factor ϕ js . Since {λ + } and {λ − } are disjoint sets, this connection cannot create circles, so the result is still a spanning tree. Furthermore, since we sum over s, this connection is realized in every possible way, which means that the second term gives the contribution of every spanning forest of {λ} in which λ j is not among the roots. Altogether the two terms contain the contribution from all of the directed spanning forests of {λ}, which is just det G according to the matrix-tree theorem.
For the off-diagonal elements of the matrix product (C.3) the same steps can be performed. To show that in this case the result is zero, it is convenient to investigate those terms that contain p ′ j , and afterwards those that do not.
The ones that contain p ′ j are In the second line j ∈ R(F ) denotes that here we consider only those spanning trees of {λ − } in which λ j is one of the roots. Because of this we can pull out a p ′ j L factor from it. Similarly to the previous argument it can be shown that this whole expression is zero: in the second line, for every partition the spanning tree in F that originates from λ j is connected to a spanning tree in T by ϕ js , and the result of this connection is still a spanning tree. Since there is a summation over s, all possible connections are included. This means that the second line contains all such spanning trees that include λ j and λ k , beside all the spanning forests of {λ − }. But the first line is exactly the same, so their difference is zero. The terms that do not contain Here in the first line λ j ∈ {λ + } and λ s ∈ {λ − }, and they are connected by ϕ js , while in the second line they are reversed. Since we sum up for every s and every bi-partition (and since λ j cannot be a root in the second line) the difference is zero. This completes the summation of the expansion.

Appendix D: Proof of the singularity properties
Here we present the proof of the singularity property (VI.7), focusing on the case of the XXZ chain. Following the reasoning in Section VI it is sufficient to show that (VI.7) holds for an operator which is the product of the elements of the monodromy matrix: where X(µ) represents one of the four operators: A, B, C, D. In order to prove (VI.7) we are using the commutation relations of the operators A, B, C, D, and the following singularity property of the scalar products (see Subsection IX. 3. of [23]) Here the elements of the sets {λ C } and {λ B } not necessarily satisfy the Bethe equations, and the mod notation means, that the scalar product is calculated with the modified vacuum expectation values a mod and d mod defined in (VI.8).
We stress an important technical detail already at this point: whereas for the Bethe states we used the re-normalized creation operators B(λ) and C(λ), because these are convenient to study the norms and mean values, for the local operator O we used the original operators A(λ), B(λ), C(λ), D(λ), because these are needed for the solution of the inverse problem.
For the sake of completeness we also present here the commutation relations resulting from the RTT equation (V.22): Instead of considering an arbitrary product in (D.1) it is useful to require a specific ordering of the X = A, B, C, D operators. It follows from the commutation relations above that any operator of the form X(µ 1 ) . . . X(µ M ) can be written as a linear combination of operators in which the order of the A, B, C and D operators is fixed in the following way: where the G({µ}) are the coefficients of the individual terms, which will include combinations of the f (u) and g(u) functions. Because of the linearity of the scalar product, all we need to show is that (VI.7) holds for a single term in this summation. We will first consider products of the form O = D(µ 1 )D(µ 2 ) . . . D(µ m )A(µ m+1 ) . . . A(µ M ). The possible presence of additional C and B operators will be treated later.
First let us consider only one A operator:

(D.6)
Here we only used the commutation relation between the operators A and B. Taking the λ C N → λ B N limit and using (D.2) we arrive at the expected expression: The computation goes similarly for the D operator. Now let's consider the case when O = M k=1 A(µ k ). To calculate the matrix element of this operator, first we have to compute the effect of it on an arbitrary state: It is obvious that the result is the linear combination of states with particle number N . The rapidities of the particles are coming from the set {µ} M ∪ {λ} N . The result can be arranged into a sum, depending on the number of rapidities coming from {µ} M . In every term in this summation we have to take into consideration every possible way how that certain amount of rapidities can be substituted from {µ} M . All this means that we can write the result in the following way: where K goes over every such subset of {λ} N ∪ {µ} M that has the following two properties: the number of its elements is N and it has exactly K elements coming from {µ} M . If µ l1 , µ l2 , . . . , µ lK denotes these elements and λ n1 , λ n2 , . . . , λ nK {λ} N those that are not present in the subset, K can be written out explicitly: . . .  (D.10) The G coefficients can be calculated using only the commutation relations. For example in the case K = 0 none of the rapidities are coming from {µ} M and we have: Let's calculate now a general coefficient, where K rapidities are coming from {µ} M ! To do this we have to keep in mind that every time when we exchange the arguments of the operators A(µ) and B(λ) during commutation there will appear a factor g(µ, λ)d(µ)/d(λ), while if we do not exchange them the corresponding factor will be f (λ, µ). Let µ l1 be the first substituted rapidity and λ n1 the one which is replaced! In this case there will be a g(µ l1 , λ n1 )d(µ l1 )/d(λ n1 ) factor coming from the commutation relation of A(µ l1 ) and B(λ n1 ), describing the replacement. The commutation of A(λ n1 ) and the other B operators will give the factor k=1 k =n1 f (λ k , λ n1 ). Finally the effect of A(λ n1 ) on the pseudo vacuum will result in the factor a(λ n1 ). In the case of the next substituted rapidity (µ l2 replacing λ n2 ) the same factors will appear, expect that now it has to be taken into consideration that one of the commutations will be with B(µ l1 ) and not with B(λ n1 ).
The rest of the substitutions go the same way. The remaining A operators with the non-substituted rapidities commute through the B operators without exchanging the arguments. To take into consideration all of the possible cases, we have to sum over every n and l. But since both the A and B operators commute with each other, the answer does not depend on which λ is replaced by which µ. Therefore, to obtain the right result we have to divide by K!. This leads to the following result for a generic coefficient: G n1,n2,...,nK l1,l2,...,lK ({λ} N |{µ} M ) = 1 K! a(λ n1 )g(µ l1 , λ n1 ) Multiplying (D.9) from the left with 0| N k=1 C(λ C k ) and taking the λ C N → λ B N limit, in the r.h.s. only those terms will give a contribution to the pole, in which λ B N is still among the arguments of the B operators. All these terms will get multiplied by i sinh(η) which means that (VI.7) holds for O = M k=1 A(µ k ). This computation goes the same way for the product of D operators (only the order of the arguments of the g and f functions are reversed). But because of the linearity of the scalar product this also proves that (VI.7) is true for an operator of the form O = D(µ 1 )D(µ 2 ) . . . D(µ m )A(µ m+1 ) . . . A(µ M ): the effect of the products of the A operators on an arbitrary state is a linear combination of states computed above. The effect of the product of the D operators on each term in this linear combination is another linear combination, with coefficients that can be similarly calculated. However, if the original term (the one obtained after acting with the A operators) doesn't contain λ B N among the arguments of the B operators then it won't appear there after acting with the D operators. If it is still an argument of a B operator in the original term, than the appropriate f factor is present next to the v.e.v., and it will appear also after acting with the D operators. This means that in every term that contains B(λ B N ) the appropriate f (λ N , ξ) or f (ξ, λ N ) factor will appear next to the v.e.v. a(ξ) or d(ξ). Therefore (VI.7) holds for O = D(µ 1 )D(µ 2 ) . . . D(µ m )A(µ m+1 ) . . . A(µ M ).
To complete the proof, we consider now the case when there are additional C and B operators. Since they are on the left and right hand side, respectively, they can be considered as additional creation operators in the states. Their presence does not alter the singularity property (VI.7) we intend to prove. To see this, let us first consider the most simple case when O = C(µ C )B(µ B ). The singularities of the form factors of this operator follow from (D.2). First we consider the re-normalized operators: (D.14) Substituting the definition (VI.1) and also using the modification rule (VI.8) for the re-normalization on the r.h.s. we get (D. 15) which has just the desired form. It can be argued similarly, that the presence of the C and B operators does not alter the singularity properties when there are D and A operators present.