Emergent hydrodynamics in integrable quantum systems out of equilibrium

Understanding the general principles underlying strongly interacting quantum states out of equilibrium is one of the most important tasks of current theoretical physics. With experiments accessing the intricate dynamics of many-body quantum systems, it is paramount to develop powerful methods that encode the emergent physics. Up to now, the strong dichotomy observed between integrable and non-integrable evolutions made an overarching theory difficult to build, especially for transport phenomena where space-time profiles are drastically different. We present a novel framework for studying transport in integrable systems: hydrodynamics with infinitely-many conservation laws. This bridges the conceptual gap between integrable and non-integrable quantum dynamics, and gives powerful tools for accurate studies of space-time profiles. We apply it to the description of energy transport between heat baths, and provide a full description of the current-carrying non-equilibrium steady state and the transition regions in a family of models including the Lieb-Liniger model of interacting Bose gases, realized in experiments.


I. INTRODUCTION
Many-body quantum systems out of equilibrium give rise to some of the most important challenges of modern physics [1]. They have received a lot of attention recently, with experiments on quantum heat flows [2,3], generalized thermalization [4,5] and light-cone effects [6]. The leading principle underlying non-equilibrium dynamics is that of local transport carried by conserved currents. Deeper understanding can be gained from studying nonequilibrium, current-carrying steady states, especially those emerging from unitary dynamics [7]. This principle gives rise to two seemingly disconnected paradigms for many-body quantum dynamics. On the one hand, taking into account only few conservation laws, emergent hydrodynamics [8][9][10][11][12] offers a powerful description where the physics of fluids dominates [13][14][15][16][17][18]. On the other hand, in integrable systems, the infinite number of conservation laws are known to lead to generalized thermalization [19][20][21] (there are many fundamental works on the subject, see the review [22]), and the presence of quasi-local charges has been shown to influence transport [23,24] (see the review [25]). However, except at criticality [26,27] (see the review [28]), no general manybody emergent dynamics has been proposed in the integrable case; with the available frameworks, these two paradigms seem difficult to bridge. The study of prethermalization or pre-relaxation under small integrability breaking [22,[28][29][30], the elusive quantum KAM theorem [31,32], the development of perturbation theory for non-equilibrium states, and the exact treatment of nonequilibrium steady states and of non-homogeneous quantum dynamics in unitary interacting integrable models remain difficult problems.
In this paper, using the recent advances on generalized thermalization and developing further aspects of integrability, we propose a solution to such problems by deriving a general theory of hydrodynamics with infinitelymany conservation laws. The theory, applicable to a large integrability class, is derived solely from the funda- mental tenet of emerging hydrodynamic: local entropy maximization (often referred to as local thermodynamic equilibrium) [33][34][35][36][37]. Focussing on quantum field theory (QFT) in one space dimension, we then study a family of models that include the paradigmatic Lieb-Liniger model [38] for interacting Bose gases, explicitly realized in experiments [4,5,[39][40][41]. We concentrate on far-fromequilibrium states driven by heat baths in the partitioning protocol [7,26,27,42] (see Fig. 1). We provide currents and full space-time profiles which are in principle experimentally accessible, beyond linear response and for arbitrary interaction strengths. We make contact with the physics of rarefaction waves, and with the concept of quasi-particle underlying integrable dynamics.
Note added: After a first version of this paper appeared as a preprint, similar dynamical equations as those derived here were independently obtained in the integrable XXZ Heisenberg chain by assuming, in addition to local entropy maximization, an underlying kinetic theory [43]. Solutions to these equations of the same type as those considered here were constructed and confirmed by numerical simulations.

II. SETUP
Let two semi-infinite halves (which we refer to as the left and right reservoirs) of a homogeneous, short-range one-dimensional quantum system be independently thermalized, say at temperatures T L and T R . Let this initial state · · · ini be evolved unitarily with the Hamiltonian H representing the full homogeneous system. One may then investigate the steady state that occurs at large times (see e.g. [28]), If the limit exists, it is a maximal-entropy steady state involving, in principle, all (quasi-)local conserved charges of the dynamics H (see (4) below). Generically, the dynamics only admits a single conserved quantity, H itself: this means that, due to diffusive processes, ordinary Gibbs thermalization occurs. However, when conserved charges exist that are odd under time reversal, the steady state may admit nonzero stationary currents. This indicates the presence of ballistic transport, and the emergence of a current-carrying state that is far from equilibrium (breaking time-reversal symmetry). This is the partitioning protocol for building non-equilibrium steady states. See Fig. 1.
The study of such non-equilibrium steady states has received a large amount of attention recently (see [28] and references therein). They form a uniquely interesting set of states: they are simple enough to be theoretically described, yet encode nontrivial aspects of non-equilibrium physics. They naturally occur in the universal region near criticality described by QFT, where ballistic transport emerges thanks to continuous translation invariance; and in integrable systems, where it often arises thanks to the infinite family of conservation laws.
Works [26,27] open to door to their study at strongcoupling critical points with unit dynamical exponent, obtaining in particular the full universal time evolution. The steady state was found to be homogeneous within a light cone, with the energy current being where c is the central charge of the conformal field theory (CFT) (below we set k B = = 1). This result arises from the independent thermalization of emerging left-and right-moving energy carriers (chiral separation). It was numerically verified [44] and agrees with recent heat-flow experiments [2]. It was generalized using hydrodynamic methods to higher-dimensional critical points [13,14,17,18] and to deviations from criticality [15,16,18]. Under conditions that are fulfilled in universal near-critical regions, inequalities that generalize (2) can be derived [28,45] (here with unit Lieb-Robinson velocity [46]) where e L,R and k L,R are, respectively, the energy densities and the pressures in the left and right reservoirs 1 . Many further results exist in free-particle models (see [28] and references therein), where independent thermalization of right-and left-movers still hold. In contrast, however, only conjectures and approximations are available for interacting integrable models [47][48][49]. In addition, a striking dichotomy is observed between integrable situations and hydrodynamic-based results: for instance, conformal hydrodynamics is expected to emerge in strong-coupling CFT [13,14], leading to shock structures, but generically fails in free-particle conformal models [50], where transition regions are smooth. This points to the stark effect of integrability on non-equilibrium quantum dynamics, still insufficiently understood with available techniques.

III. EMERGING HYDRODYNAMICS IN QUANTUM SYSTEMS
Let us recall the basic concepts underlying the hydrodynamic description of many-body quantum systems, and its use in the setup described above (similar concepts exist in many-body classical systems).
Let Q i , i ∈ {1, 2, . . . , N } be local conserved quantities in involution. These are integrals of local densities q i (x, t), and the conservation laws take the form ∂ t q i (x, t) + ∂ x j i (x, t) = 0, where j i are the associated local currents. A Gibbs ensemble is a maximal-entropy ensemble under conditions fixing all averaged local conserved densities. It is described by a density matrix where β i are the associated potentials. For instance, Q 1 is taken as the Hamiltonian, and β 1 is the inverse temperature. We will denote β = (β 1 , β 2 , . . . , β N ) the vector representing this state, and · · · β the averages. Clearly, the Gibbs averages of local densities q i = q i β (these are independent of space and time by homogeneity and stationarity) may be seen as defining a map from states to averages, β → q. This is expected to be a bijection: the set of averages fully determines the set of potentials. Therefore, the current averages j i = j i β are functions of the density averages, These are the equations of state, and are modeldependent. The averages q can be generated by differentiation of the (specific, dimensionless) free energy f β .
Similarly, one can show [28] (see Appendix A) that there exists a function g β that, likewise, generates the currents, A hydrodynamics description of quantum dynamics is expected to emerge at large space-time scales. This has been exploited, in the present setup, in [13][14][15][16][17][18]. The emergence of hydrodynamics is solely based on the assumption of local entropy maximization (or local thermodynamic equilibrium 2 ). Technically, this is the assumption that averages of local quantities O(x, t) tend uniformly enough, at large times, to averages evaluated in local Gibbs ensembles O β(x,t) with space-time dependent potentials β(x, t). Physically, this is a consequence of separation of scales, as follows (see for instance [37]).
Assume that, after some time, physical properties vary only on space-time scales that are much larger than microscopic scales. This may be referred to as the "local relaxation time". From that time on, microscopic processes such as particle collisions or inter-site interactions give rise to fast, local relaxation: the reaching of a (approximate) steady state on space-time scales small compared to variations but large enough for thermodynamics to be applicable. By Boltzmann's phase-space argument, these local steady states are obtained from entropy maximization, and as usual maximization is under the conditions provided by conservation laws (properties of the microscopic dynamics). That is, on each space-time "fluid cell" a Gibbs state is (very nearly) reached. Neighboring Gibbs states are different, but their variations are small. This is local entropy maximization.
Assume local entropy maximization. On each fluid cell, the Gibbs state is initially characterized by the values of the conserved densities at the local-relaxation time. The large-scale dynamics is thereon obtained from conservation laws, as follows. Consider microscopic conservation in integral form, Since averages of densities and currents, after the local relaxation time, take the form q i (x, t) = q i β(x,t) and j i (x, t) = j i β(x,t) uniformly enough, we have where q(x, t) = q β(x,t) and j(x, t) = j β(x,t) . Here, integrals may be taken to cover a macroscopic number of fluid cells: these become macroscopic conservation equations. Macroscopic conservation equations can be re-written in differential form, with differentials representing small variations amongst fluid cells: These are the pure hydrodynamic (Euler-type) equations, representing the slow, large-scale quantum dynamics of conserved densities and currents flowing amongst neighboring cells.
The problem of emergence of hydrodynamics in manybody systems is one of the most important unsolved problem of modern mathematical physics. Although there are few proofs of emergence of hydrodynamics, there is strong evidence for the validity of emerging Euler equations in many situations; see the books [33][34][35][36][37], and the recent paper [51] for a study of emerging Euler equations in classical anharmonic chains.
Combined with the equations of state (5), Euler equations (8) give where J(q) := ∇ q j is an N by N matrix, the Jacobian matrix of the transformation from densities to currents, Equations (9) are the emergent pure hydrodynamic equations in quasi-linear (or characteristic) form [12]. The complete model dependence, including all quantum effects, is encoded, besides the number N of conserved quantities, in the Jacobian J(q). The density averages q, like the potentials β, correspond to a set of state coordinates. One may choose any other set of state coordinates n, with q = F q (n) and j = F j (n). A similar equation is obtained, where J(n) = ∇ n q −1 ∇ n j. Observe that J(n) and J(q) are related to each other by a similarity transformation: J(n) = ∇ n j −1 J(q)| q=F q (n) ∇ n j. Therefore, the spectrum of J(n) is independent of the choice of coordinates, and is a fundamental property of the model. We will denote this spectrum by {v eff i (n), i = 1, . . . , N }. Choosing coordinates n that diagonalize J(n) one obtains These express the vanishing of the convective derivatives, representing the constance of each fluid mode n i (x, t) on fluid cells. The eigenvalues v eff i (n(x, t)) are therefore interpreted as the propagation velocities of these normal modes. The normal modes interact with each other only through the propagation velocities, which is generically a function of all state coordinates.
Let us now apply the above to the solution of the partitioning problem. For clarity of the following discussion, we come back to the q-coordinates (but it is easy to generalize to any coordinates n). Consider the large-scale limit (x, t) → (ax, at), a → ∞. Because (9) is invariant under this scaling, in the limit, if it exists, the solution is also invariant. Thus we may assume self-similar solutions β(x, t) = β(ξ) where ξ = x/t, and (9) becomes an eigenvalue equation, The initial condition is determined by the state at the local relaxation time (at which the fluid-dynamics description starts to be valid). This state is unknown, as it depends on the full quantum dynamics, but its asymptotic at large |x| is identical to that of the original state.
In the large-scale solution, the initial condition t → 0 + is implemented as asymptotic conditions as ξ → ±∞. Therefore it only depends on the asymptotic form of the initial state, and we impose In the present setup, these involve Gibbs states at potentials β R,L : and the steady-state averages are given by The solution to the eigenvalue equation (13) and initial conditions (14) provides the exact large-scale asymptotic form of the full quantum solution, along any ray x = ξt (see Fig. 1). The eigenvalue equation (13) represents the small changes of averages along various rays, due to the exchange of conserved charges amongst fluid cells. The set of eigenvalues of J(q) -the available propagation velocities in the state characterized by the averages qform a finite, discrete set for finite N .
Solutions to (13), (14) are typically composed of regions of constant q-values separated by transition regions [12]. Transition regions may be of two types: either shocks (weak solutions), where q-values display finite jumps, or rarefaction waves, where they form a smooth solution to (13). Rarefaction waves, the most natural type of solution, cannot, generically, cover the full space between two reservoirs. Indeed, (13) specifies that the curve traced by the solution in the q-plane must at all points be tangent to an eigenvector of J(q). Since eigenvectors -and available propagation velocities -form a discrete set, smooth variations of q along the curve imply a unique choice of eigenvector at each point (except possibly at points where eigenvalues cross). Thus, the curve is completely determined by its initial point, and cannot join two arbitrary reservoir values. That is, in ordinary pure hydrodynamics, shocks are often required.

IV. HYDRODYNAMICS WITH INFINITELY MANY CURRENTS
In integrable systems, there are infinitely many local conservation laws. In fact, this space is enlarged to that of "pseudolocal" conservation laws, where the densities q i (x, t) and currents j i (x, t) are supported on extended spacial regions with weight decaying fast enough away from x. This enlargement plays an important role in non-equilibrium quantum dynamics [20,21,[23][24][25]. Under maximal-entropy principles, Gibbs states are then replaced by generalized Gibbs ensembles (GGE) [19,21,22]: formally the limit N → ∞ of the density matrix (4), involving all basis elements in the space of conserved pseudolocal charges. We choose Q 1 = H (the Hamiltonian) and Q 2 = P (the momentum operator).
Under the influence of infinitely many conservation laws, the picture of local entropy maximization is still expected to hold: all physical principles underlying it stay unchanged, the only difference being the use of GGEs instead of Gibbs ensembles. This, along with the emergence of self-similar solutions in the partitioning protocol, are our working hypotheses; see Appendix B for a discussion. The emergence of a generalized type of hydrodynamics was proven in the classical hard-rod problem [37,52], whose relation with the present quantum problem we will study in a future work. The emergence of self-similar solutions was observed numerically in the quantum XXZ chain in [53]. In free-particle quantum models, hydrodynamic ideas and related semi-classical approximations, as well as ray-dependent local entropy maximization, were studied in various works, see [54][55][56][57][58][59].
Looking for a full solution to the infinity of equations (9), (13) and (14), an appropriate choice of state variables is crucial. A powerful way is to recast them into the quasi-particle language underlying the thermodynamic Bethe ansatz (TBA) method for integrable systems [60]. Using this language, we derive the exact GGE equations of state, and the ensuing generalized hydrodynamics equation. We determine the exact normal modes and propagation velocities, and obtain full ray-dependent solutions.

A. GGE equations of state
We assume that the spectrum of stable quasi-particles is composed of a single quasi-particle species of mass m (see Appendix C for a many-particle generalization). The dispersion relation is encoded via a parametrization E(θ), p(θ) of the energy and momentum: in the relativistic case θ is the rapidity, E(θ) := m cosh(θ), p(θ) := m sinh(θ), and in the Galilean case θ is the velocity, E(θ) := mθ 2 /2, p(θ) := mθ. A differential scattering phase ϕ(θ) fully specifies the dynamics of the model [60]. We denote by h i (θ) the one-particle eigenvalue of the conserved charge Q i ; in particular h 1 (θ) = E(θ) and h 2 (θ) = p(θ).
Let us first recall the basic ingredients of TBA. Three related quantities play important roles: the quasi-particle density ρ p (θ), the state density ρ s (θ), and the quasiparticle occupation number n(θ) := ρ p (θ)/ρ s (θ). The functions ρ p (θ) and n(θ) are two different sets of state coordinates; each can be used to fully characterize the GGE. The former specifies all average densities in a simple way: This can in fact be seen as a definition of ρ p (θ). Here and below, integrations are over R.
As a consequence of interaction, quasi-particle and state densities are related to each other. Using the Bethe ansatz, one finds the following constitutive relation [60]: where p (θ) = dp(θ)/dθ. This relation gives rise to a nonlinear relation between the state coordinates ρ p (θ) and n(θ). The transformation from the former to the latter is direct from the above definitions. In the opposite direction, the transformation is effected by where the "dressing" operation h → h dr is defined by the solution to the linear integral equation The potentials β can be recovered: the occupation number is related to the one-particle eigenvalue w(θ) = i β i h i (θ) of the charge i β i Q i in the GGE (4) via the so-called pseudo-energy w (θ) [60,61]: The above ingredients give exact average densities as functions of GGE states. However, they do not provide expressions for average currents as functions of state coordinates, and for equations of states. Hence they are not sufficient in order to develop generalized hydrodynamics.
We solve this problem by obtaining the following expressions: is the dressed one-particle eigenvalue. These expressions emphasize the role of relativistic or Galilean symmetry: the sole difference between GGE averages of charge densities and currents is the integration measure, determined by the dispersion relation.
The first equation in (22) is well known and is a consequence of (17) and (19). In integral-operator notation (with measure dθ/(2π)), the dressing operation is where N is diagonal with kernel 2π n(θ)δ(θ − α), and ϕ has kernel ϕ(θ − α). Therefore, introducing the symmetric operator U = N (1 − ϕN ) −1 and the bilinear form which leads to the first equation of (22). The second equation in (22) is new. It may be proven, in the relativistic case, using relativistic crossing symmetry, and then obtained by the non-relativistic limit in the Galilean case. In the relativistic case, crossing symmetry says that local currents j i , in the cross-channel, are local densities q i ; therefore the current expression in (22) is obtained from that of the density under an appropriate exchange of energy and momentum. Let C be the crossing transformation (x, t) → (it, −ix), implemented on rapidities by θ → iπ/2 − θ. Note that it squares to the identity, C 2 = 1. Let us denote by q[h] and j[h] the density and current operators, respectively, associated to a one-particle eigenvalue h(θ). Then the statement that local currents j i , in the cross-channel, are local densities . Let us also denote by O w the average of observables O in the state characterized by w(θ). (22), we obtain that for j i = j[h i ]. An alternative proof, using the machinery of integrable systems, is presented in Appendix D.
Expressions (22) have interesting consequences. First, (22), the average current may also be written in terms of a current spectral density ρ c (θ): which takes the forms Here v eff (θ) is the effective velocity, defined by The effective velocity depends on the state via the occupation number entering the dressing operation, and brings out the quasi-particle interpretation of the current expression: since ρ c (θ) = v eff (θ)ρ p (θ), quasi-particles are seen as moving at effective velocities v eff (θ), influenced by the state in which they move.
Second, one may extract explicit GGE equations of state from expressions (22). The equations of states are necessary and sufficient relations between densities and currents, guaranteeing the existence of n(θ) such that both relations in (22) hold for all h i (θ). Assume that q i and j i are averages in a state, not necessarily a GGE. In complete generality, both are linear functionals of h(θ), hence we may still write (17) and (25) for some quasiparticle density ρ p (θ) and current spectral density ρ c (θ). GGE equations of states can therefore be written as relations between ρ p (θ) and ρ c (θ), necessary and sufficient for the existence of n(θ) such that (22) hold. One can show that these relations are These relations are independent of the state: they hold in any GGE, in the model described by the differential scattering phase ϕ(θ − α). They characterize the set of doublets of functions (ρ p , ρ c ) describing available GGEs for this integrable model. The proof of (28) is obtained by isolating n(θ) in both (19) and (26), in the forms 2π(N −1 − ϕ)ρ p = p and 2π(N −1 − ϕ)ρ c = E , and equating the resulting expressions. Finally, recalling (26), the left hand side of (28) is v eff (θ). Simple manipulations of (28) then give a linear integral equation for the effective velocity v eff (θ) in terms of quasi-particle densities: where v gr (θ) = E (θ)/p (θ) is the group velocity. In this form, the equations of state of integrable systems are seen as equations specifying an effective velocity of quasiparticles, as a modification of the group velocity.
We note that the effective velocity derived here agrees with that proposed in [62]. This is interesting, as our derivation is based on comparing current spectral density to quasi-particle density, while the concept proposed in [62] is based on stationary-phase arguments 3 .

B. Generalized hydrodynamics
The basic equation of generalized pure hydrodynamics is derived from (8) along with the quasi-particle expressions (17) and (25). The fact that the space of pseudolocal charges is complete [21] suggests that these hold for a complete set of functions h i (θ), and thus (here and below we suppress explicit x, t dependences for lightness of notation): Using the equations of state (28), this is an integrodifferential system for the space-time dependent state characterized by the particle densities ρ p (θ). Alternatively, using the dressed-velocity formulation (26) and (29), Equation (30) may be written as This is the conservation form of generalized hydrodynamics. It is a density-type conservation equation, and identifies ρ p (θ) as a conserved fluid density. The state coordinates ρ p (θ) are, however, not the most convenient. We show that the occupation numbers n(θ) diagonalize the Jacobian J(n) in the quasi-linear form (11): the space-time dependent occupation number n(θ) satisfies the following integro-differential system, the vanishing of the convective derivative of n(θ): Here (27) may be used to express the effective velocity in terms of n(θ). Hence, n(θ) are the normal modes of generalized hydrodynamics, and further, the eigenvalues -the propagation velocities -are exactly the effective velocities v eff (θ). The proof of (32) is as follows. Using the integraloperator relations 2πρ p = Up and 2πρ c = UE , we have which gives (32) using (23). Observe that using (31) and (32), it is simple to show that the state density ρ s (θ), as well as the hole density ρ h (θ) := ρ s (θ) − ρ p (θ), also satisfy the same density-type conservation equation (31). Further, as a consequence, the entropy density [60], also satisfies this conservation equation, ∂ t s(θ) + ∂ x v eff (θ)s(θ) = 0. Conservation of entropy density is a fundamental property of perfect fluids, as no viscosity effects are taken into account.
One can show that the solution for the nonequilibrium, ray-dependent occupation number n(θ) is the discontinuous function where Θ(· · ·) is Heavyside's step function. The position of the discontinuity θ depends on ξ and is selfconsistently determined by v eff (θ ) = ξ; equivalently, it is the zero of the dressed, boosted momentum p ξ (θ) := p(θ − η) where ξ = tanh η (relativistic case) or ξ = η (Galilean case), The GGE occupation numbers n L,R (θ) entering (35) guarantee that the asymptotic conditions on ξ correctly represent the asymptotic baths as per (14). They are obtained using (21) with w = w L,R (θ) the one-particle eigenvalues characterizing the GGE of the left and right asymptotic reservoirs; for instance, with reservoirs at temperatures T L,R , we have w L,R (θ) = T −1 L,R E(θ). Indeed, the solution (35) of the scaled problem holds since v eff (θ) is monotonic and covers the full range of θ (which is [−1, 1] in the relativistic case and R in the Galilean case): therefore there is a unique solution to v eff (θ) = ξ, thus a unique jump; and θ is monotonic with ξ, hence the asymptotic conditions are correctly implemented.
The system of integral equations (22), (20), (35) and (36) can be solved numerically using Mathematica, yielding extremely accurate results. Integral equations in (21) and (20) can be solved iteratively, a procedure that converges fast [60]. The hydrodynamic solution is obtained by first constructing the thermal occupation numbers n L,R (θ) (21). Then, the non-equilibrium occupation number is evaluated by solving the system (35), (36): one first chooses θ = η in order to construct n(θ), and then evaluates p dr ξ (θ). The zero of p dr ξ (θ) is numerically found -we observed that p dr ξ (θ) always has a single zero. The process is repeated until the zero is stable -we observed that this is a convergent procedure. Finally, the nonequilibrium occupation number is used in (22), (20). The solving time increases slowly with the numerical precision demanded, thus this allows arbitrary-precision results.
The solution presented may be interpreted as a single space-covering rarefaction wave, in the sense that it is a solution to the eigenvalue equation (13) where all physical observables q i , j i are continuous and interpolate between the two reservoirs. With relativistic dispersion relation, the solution is smooth within the light cone, beyond which the states are constant and equal to the initial baths' states; while in the Galilean case, the solution is generically smooth on the whole space. In this solution, every normal mode n(θ), seen as a function of ξ for fixed θ, is discontinuous exactly at its propagation velocity. Every normal mode therefore displays a "contact discontinuity" (a discontinuity without entropy production) [12]. Hence, the rarefaction wave may be seen as being composed of infinitely many contact discontinuities. In contrast to the finite-dimensional case, this single rarefaction wave can account for generic reservoirs, and no shock need to develop. This is because in the infinite-dimensional case, the eigenvalues of J(n) form a continuum: all propagation velocities v eff (θ) are available as conserved charges guarantee a large number of stable excitations, providing an additional continuous parameter tuning the smooth state trajectory and guaranteeing its correct asymptotic-reservoir values. Since weak solutions (shocks) are not necessary to connect the asymptotic reservoirs, they do not appear.

V. ANALYSIS AND DISCUSSION
Concentrating on pure thermal transport, we have analyzed the above general system of equations for two related models: the relativistic integrable sinh-Gordon model and its non-relativistic limit [63,64], the (repulsive) Lieb-Liniger model. We have also verified that our hydrodynamic equations reproduce the known results for the case of free particles.

A. The relativistic sinh-Gordon model
One of the simplest integrable relativistic QFT with non-trivial interactions is the sinh-Gordon model. It is defined by the Lagrangian [65,66]: where φ is the sinh-Gordon field and m is the mass of the single particle in the spectrum. The model is integrable and therefore the only non-trivial scattering matrix is that associated to two-particle scattering. It is given by [67][68][69] The parameter B ∈ [0, 2] is the effective coupling constant which is related to the coupling constant β in the Lagrangian by under CFT normalization [70]. The S-matrix is obviously invariant under the transformation B → 2 − B, a symmetry which is also referred to as weak-strong coupling duality, as it corresponds to B(β) → B(β −1 ) in (39). The point B = 1 is known as the self-dual point. At the self-dual point the TBA differential scattering phase is simply Contrary to the Lieb-Liniger model which we will discuss later, the general features of any quantities of interest in the sinh-Gordon model are very similar for any values of the coupling B. For this reason, in this paper we will concentrate our analysis solely on the self-dual point in the understanding that similar results hold for other values of B.
We have evaluated the energy density e := q 1 , energy current j := j 1 and pressure k := j 2 . Typical profiles are shown in Figs. 2, 3. Fig. 2 shows smooth interpolation within the light cone between the asymptotic baths at ξ = −1 and ξ = 1 (the speed of light is set to 1). Fig. 3 shows how, as temperatures rise, the current approaches the plateau (2) predicted by CFT [26,27]. Further, in Fig. 4, the relative deviation of the steady-state current from its bounds (3)   The numerical data have been obtained by solving the integral equations recursively until convergence is reached. Sources of error are the discretization and finite range of θ for numerical integration. Adjusting the number of divisions and the range, we estimate the error to be less than 0.1%.

B. The Lieb-Liniger model
The Lieb-Liniger (LL) model, in the repulsive regime (λ > 0), can be regarded as a non-relativistic limit of the sinh-Gordon model, as shown in [63,64]. The Hamiltonian of the model is given by This is obtained from the Hamiltonian of the sinh-Gordon model by a double limit where c is the speed of light (which was implicit in (37)) 4 . This limit can be performed within the TBA formalism [64], and accordingly, the density and current averages q i , j i are given by (22), with the non-relativistic dispersion relation. There, the occupation number is given by n LL (θ) = 1/(1+e w (θ) ), and the pseudo-energy w (θ) and the dressed one-particle eigenvalues h dr i (θ) are defined in the same manner as in Equations (21) and (20) (where θ = p/m is the velocity), with scattering matrix given by and corresponding differential scattering phase A uniform chemical potential µ is introduced, associated to the conserved charge Q 0 that counts the number of quasi-particles (with h 0 (θ) = 1). The energy current is chosen to be the current associated to of charge H −µQ 0 , Below we present some numerical results for several values of the coupling λ and for m = 1. Current profiles obtained for λ = 3 and various values of µ are displayed in Fig. 5. The main difference between the relativistic and non-relativistic cases is the lack, in the latter, of any sharp light-cone effect. Nevertheless, at low temperatures T L,R µ, Luttinger Liquid physics emerges [71], including an emerging light-cone due to the Fermi velocity. This can be seen in Fig. 5: a plateau forms whose height is again in agreement with the general CFT result (2). The plateau lies between nearly symmetric values ξ/v F ≈ ±1 fixed by the Fermi velocity v F . Thermal occupation numbers n L,R (θ) are very sharply supported between Fermi points ±θ L,R F with θ L,R 2µ/m, and the Fermi velocity, which depends on ξ very weakly, is the effective velocity v eff (θ R F ) associated to the lowest temperature (T R < T L in the present example). In agreement with general CFT results [26,27], a light cone thus builds up (despite the model having no intrinsic maximal velocity), and the full state is in fact homogeneous between the Fermi velocities.
In the LL model the coupling constant may take any values between 0 and ∞ and the limits λ → 0 and λ → ∞ are of particular interest.
For λ → 0 the differential scattering phase (44) becomes heavily peaked around θ = 0. Formally, lim λ→0 ϕ LL (θ) = 2πδ(θ). The resulting TBA equations, with this differential scattering phase, admit no solution for the pseudoenergy if µ > 0, but for µ < 0 they can be solved exactly and reproduce the free Boson solution (for which µ > 0 would make no physical sense). In particular the energy current takes the free Boson form, where α L,R = β L,R ( ξ 2 2 − µ). In Fig. 6 we compare numerical values for λ = 0.05 and µ = −1 to this analytical expression. The agreement is very good, confirming that a free Boson theory is smoothly recovered in this limit. With µ > 0, as λ becomes small the TBA equations gradually breakdown. How this occurs is subtle, and will be discussed in [81].
The qualitative change in behaviour of the TBA solutions as λ → 0 is related to the two distinct regimes observed at small values of λ [72]. Consider the dimensionless coupling γ := 2mλ/q 0 (where we recall that q 0 is the gas density, which may be taken in the initial baths for instance) and the reduced temperature τ := 2mT /q 2 0 . The "decoherent regime", with large phase and density fluctuations, occurs for γ min(τ 2 , √ τ ). In this regime, ideal Bose gas physics is recovered, and we have indeed verified that the inequality is satisfied in the parameter space where good agreement with (46) is observed (small λ, negative µ). On the other hand, the "Gross-Pitaevskii regime" occurs for τ 2 γ 1, a quasi-condensate with large phase fluctuations but suppressed density fluctuations. It is such quasi-condensate physics that strongly affects TBA solutions as λ → 0 with µ > 0.
The other interesting limit is lim λ→∞ ϕ LL (θ) = 0. In this case we can also find an analytical expression for the current: . (47) This corresponds to a free Fermion, in agreement with the expected Tonks-Girardeau physics occuring in the regime γ max(1, √ τ ) [72]. For ξ ≈ 0 and µβ L,R 1 it is easy to show that the integral above gives π 12 (β −2 L − β −2 R ) so that we recover the CFT result for the current with c = 1 (Dirac Fermion). Fig. 7 shows a comparison between numerical values of the current for λ = 50 and the formula above.
Let us now consider the particle current. Naturally, in the LL model, equilibrium states at higher temperatures have lower particle densities. Therefore, although the energy current flows from the left to the right in the present setup (with T L > T R ), the initial particle density imbalance would naively suggest a particle flow from the right (higher density) to the left (lower density). The opposite occurs: Fig. 9 shows that the particle current is positive, hence flows form the left to the right. This   (47) for the same temperatures and chemical potential, whose profile is not dissimilar to the plots shown in Fig. 5. As before, the bold horizontal line is the CFT value π 12 1 − 1 25 . The agreement is extremely good.
means that the fluid flow produced by the temperature difference drags particles with enough force to counteract the particle imbalance and bring particles towards the higher-density bath. The fact that heat carries particles along its motion is a thermoelectric effect. It has been experimentally demonstrated in a quasi-two-dimensional fermionic cold atoms channel [3], and theoretically shown in CFT in dimensions higher than one [17]. It is nontrivial in integrable models, as the large amount of conservation laws allows for independent currents to coexist, and our result gives the first theoretical prediction of this effect in the integrable one-dimensional Bose gas. An additional consequence of the thermoelectric effect is that the particle density q 0 (ξ) shows particle accumulation around v F and depletion around −v F (see Fig. 8). For instance, the start of the dip can be explained by the fact that, in any local spacial region originally in the left reservoir, the first particles to start moving towards the  9. A characteristic profile of the Lieb-Liniger particle current for TL,R µ, λ = 3 and µ = 6. The local maxima/minima are located around ξ = ±vF . The dashed curve is an interpolation.
right are those on the right of the region, escaping and thus depleting it. Since time evolution at fixed position is obtained by scanning Fig. 8 from left to right, this explains the initial dip on the left. This depleting effect continues as long as the outgoing current on the right of the region is higher then the incoming current on its leftthat is, until the region lies in the steady state. However, as time goes on, the effective local temperature decreases, and this tends to increase the particle density. This effect eventually overtakes the depleting effect, accounting for the rebounce to the higher steady-state value. The behavior of the current j 0 in Fig. 9 is then a consequence of the continuity equation ξ∂ ξ q 0 = ∂ ξ j 0 . This is a nonuniversal effect, not present in the density q 1 (ξ) − µq 0 (ξ) controlled by low-energy processes, where the physics of chiral separation dominates and monotonic transition regions occur.

C. General features
The form of the non-equilibrium occupation number indicates that quasi-particles are thermalized according to the initial GGE's, in a way that depends on the rapidity. It connects with the picture, proposed in [26,47], according to which in the steady state (ξ = 0), quasiparticles traveling towards the right (left) are thermalized according to the left (right) reservoir. However, in the present solution, what determines the traveling direction is the effective velocity in the steady state: quasi-particles with positive (negative) dressed velocities, reaching the point x = 0 at large times, will have travelled mostly towards the right (left) (after a relatively small transient).
In the sinh-Gordon model with T L > T R , the effective velocity behaves as in Fig. 10. We observe that it is greater than the bare velocity tanh θ for small or negative rapidities, and smaller for large positive rapidities. This is in agreement with the intuition according to which the quasi-particles are effectively carried by the flow, which transports them towards the right, for small enough rapidities, but slowed down by dominant "friction" effects of thermal fluctuations at large rapidities. A similar effects occur in the Lieb-Liniger model.
The generalized hydrodynamic result differs from previous proposals in interacting integrable models [47][48][49] (while all results agree in noninteracting cases). The original proposal [47] was later shown [45] to break the second inequality in (3), while the second proposal [48], based on similar ideas, gave slight disagreements with numerical simulations. The conjecture [49] which corresponds to taking θ = 0 in our framework, seems to give good agreement with numerical simulations. This may be due to the fact that taking θ = 0, gives very small errors in wide temperature ranges, of the order of 0.5-1% (we have confirmed this numerically).

VI. CONCLUSIONS
In this paper we developed a hydrodynamic theory for infinitely-many conservation laws, and applied it to the study of heat flows in experimentally relevant integrable models. It would be interesting to study further the nonequilibrium physics of the Lieb-Liniger model, including the effects of the Gross-Pitaevskii quasi-condensate and transport between different regimes. The emerging physical picture and solution we have given can be applied to any Bethe-ansatz integrable model, where infinitelymany conservation laws exist and a quasi-particle description is available. This includes quantum chains (see [43]), as continuity of space on which the microscopic theory lies is not needed for emerging hydrodynamics. It also includes relativistic models with non-diagonal scattering such as the sine-Gordon model, where, for instance, our TBA construction may be generalized along the lines of the famous approach of Destri and de Vega [74,75]. Of course, the hydrodynamic ideas do not require a quasiparticle description, and it might be possible to develop generalized hydrodynamics using a variety of techniques from integrability. We note that it is remarkable that independent quasi-particle mode thermalization agrees, in integrable models, with local entropy maximization. The dynamical equations derived can be used to describe more general situations in ultracold gases such as the release from a trap (see e.g. [73]). This new theory and its extensions, including viscosity effects and forcing, should also allow for efficient studies of integrability breaking and related problems in any dimensionality, as well as for exact descriptions of dynamics in smooth trapping potentials [4] at arbitrary coupling strength.
In the first line we used the fact that dx q m (x) is a conserved quantity and thus commutes with the density matrix ρ; in the second we used space translation invariance, in the third integration by parts and the fastenough vanishing of correlation functions; in the fourth current conservation; in the fifth time-translation invariance; in the sixth current conservation, in the seventh space-translation invariance; and in the eight integration by parts. Therefore, and thus showing Equation (6).

Appendix B: Emergence of generalized hydrodynamics
The only principle at the basis of hydrodynamics, and of the derivation we provide, is that of the emergence of local generalized thermalization (local entropy maximization). Technically, this is the assumption that averages of local quantities O(x, t) tend uniformly enough, at large x and t, to averages evaluated in GGEs (infinite-volume, maximal-entropy states, under conditions on infinitelymany conservation laws), with space-time-dependent potentials. This assumption is sufficient to derive the explicit dynamics for all single-point averages of local conserved densities and currents: no ad-hoc kinetic principle is needed.
In the case of infinitely-many conservation laws, one delicate point is the consideration of quasi-local densities and currents, which are involved in generalized thermalization. Such a quantity is not supported on a finite region, but on an extended region, with a weight (as measured by, for instance, the overlap with any other local observable) that decays away from a point. However, since hydrodynamics is concerned with large-scale space-time regions (the fluid cells), it is natural to consider them on the same footing. This is implicitly done in the derivation presented in this paper by assuming a completeness property of conservation laws.
Another delicate point concerns the definition of GGEs itself. In finite systems, such states depend on the boundary conditions imposed, and in general, these boundary conditions may still have an effect in the infinite-volume limit. For instance, walls simply preclude any nonzero potential associated to the momentum operator, as they break translation invariance. Nevertheless, given a set of allowed conserved charges, at large volumes, boundary conditions have little effect on local averages of conserved currents and densities (as they do not affect specific free energies). Further, periodic boundary conditions, at the basis of the TBA formalism, appear to provide the maximal set of conserved charges. It is in fact possible to construct GGEs directly in infinite volumes [21]. We expect local thermalization, and the full set of available conserved charges, to be correctly described by such constructions; and we expect these to agree with the TBA formalism used here.
We finally mention that the classical hard-rod problem, proven to give rise to a form of hydrodynamics [37,52], has strong connections with the integrable systems investigated here, which we will investigate in a future work.
This exactly coincides with (22). Similar arguments give rise to current averages associated to flows i[Q k , q i ] + ∂ x j (k) i = 0 with respect to any local conserved quantity Q k (with odd spin): A full derivation will be given in [81].