A note on generalized hydrodynamics: inhomogeneous fields and other concepts

Generalized hydrodynamics (GHD) was proposed recently as a formulation of hydrodynamics for integrable systems, taking into account infinitely-many conservation laws. In this note we further develop the theory in various directions. By extending GHD to all commuting flows of the integrable model, we provide a full description of how to take into account weakly varying force fields, temperature fields and other inhomogeneous external fields within GHD. We expect this can be used, for instance, to characterize the non-equilibrium dynamics of one-dimensional Bose gases in trap potentials. We further show how the equations of state at the core of GHD follow from the continuity relation for entropy, and we show how to recover Euler-like equations and discuss possible viscosity terms.


I. INTRODUCTION
The emergence of hydrodynamics in many-body extended systems is based on the phenomenon of local entropy maximization (often referred to as local thermodynamic equilibrium) [1][2][3][4][5]. This is the phenomenon according to which, at large times, the system decomposes into slowly varying local "fluid cells" where homogeneous Gibbs states exist. At leading order in a derivative expansion, the ensuing dynamics on the Gibbs potentials is completely fixed by the local conservation laws -this is often referred to as "pure hydrodynamics", as viscosity terms are absent. This is a powerful description, replacing the full many-body evolution, either quantum or classical, by differential equations for the few (or at least fewer) relevant local state parameters. It allows for the precise description of large-scale structures and the unearthing of exact results, and its universal applicability has been demonstrated in various situations and models [6][7][8][9]. In particular, it provides striking results in the context of quantum transport far from equilibrium [10][11][12][13][14][15][16][17] (see also the review [18]).
Recently [19], see also the related work [20], the hydrodynamic idea was extended to many-body integrable systems, where infinitely-many conservation laws are present. In this context, entropy maximization is conjectured to generate states in the infinite-dimensional variety of so-called generalized Gibbs ensembles 1 (GGEs) [22,23], which therefore are used to characterize fluid cells. In [19], it was shown, in general diagonal-scattering integrable models of quantum field theory including the Lieb-Liniger and sinh-Gordon models, that the infinite system of conservation laws -for the infinite number of GGE potentials -can be recast into a system of hydrodynamic equations for quantities characterizing occupations and densities of quasi-particle states. In [20], the same equations were obtained in integrable quantum Heisenberg chains (the derivation making use of an additional assumption about the underlying dynamics). Interest-ingly, as will be studied in a coming work, these equations appear to give a universal description of quantum and classical quasi-particle elastic scattering; they widely generalize, for instance, hydrodynamic equations proven to emerge in the classical hard-rods model [5,24]. In the same context, the effect of a localized defect on the non-equilibrium transport in quantum chains was also analyzed in [25,26]. In fact even in free models, the hydrodynamic idea, as a semi-classical approximation, has found many applications [27][28][29][30][31][32].
The purpose of this letter is to extend this "generalized hydrodynamics" (GHD) theory further, within the quantum framework. We start by reviewing the main results of GHD in Section II. In Section III, we show that the GGE equations of state, at the core of GHD, are consequences of hydrodynamic entropy conservation. In Section IV we show how to represent the dynamics associated to all conservation laws, not just the Hamiltonian. Using this, in Section V we derive GHD equations in the presence of external inhomogeneous fields, including force fields. Finally, in Section VI we connect with aspects of ordinary fluid dynamics, including a derivation of Euler equations and a proposition for possible viscosity terms. We provide additional details in appendices, and especially in Appendix A we discuss how, in general free models, weak space-time variations of local densities and currents at large times guarantee the emergence of local GGEs, hence of GHD.
We emphasize that the force-field equations obtained can serve as a powerful tool in describing the late time non-equilibrium dynamics of a one-dimensional Bose gas confined in a weakly-varying trap potential. We also note that, recently, an alternative method to incorporate the inhomogeneity introduced by an external potential in a one-dimensional conformal field theory was proposed in [33].

II. REVIEW OF GHD
In this section we recall some of the basic concepts developed in [19] and [20], concentrating on the approach taken in the former, which puts emphasis on hydrody-namics ideas. The basic objects in the hydrodynamic theory of many-body extended systems are the local conserved densities q i (x, t) and local currents j i (x, t). These are quantum operators satisfying, under unitary dynamics, the continuity relations, or conservation laws, ∂ t q i (x, t) + ∂ x j i (x, t) = 0 (1) as a consequence of the total charge Q i := dx q i (x, t) being conserved ∂ t Q i = 0. The set of such local conservation laws is a characteristics of the many-body system. In integrable systems, this set is infinite, and the charges Q i relevant to the problem span the space of pseudolocal conserved charges [34]. In particular, entropy maximization of local subsystems under constraints of these conservation laws, as occurs under unitary dynamics, gives rise to GGEs, formally described by density matrices of the form 2 exp [− i β i Q i ]. It was in fact shown rigorously [35] that, in homogeneous systems, the existence of the long time limit implies that the stationary state is a GGE, where the completeness of the space of pseudolocal charges plays an important role. We will denote averages in such GGEs as · · · β (with β = (β i ) i ), and, for lightness of notation, The problem of pure generalized hydrodynamics, as formulated in [19], is a direct generalization of usual pure hydrodynamics (without viscosity): it is the continuity problem applied to local cells where independent entropy maximization has occurred. That is, one assumes β = β(x, t), and writes where q i = q i β(x,t) and j i = j i β(x,t) . A convenient way of fixing the hydrodynamic problem for a given model is to provide the equations of state: relations connecting averages of currents to averages of densities. The thermodynamic Bethe ansatz (TBA) formulation of GGE averages offers a powerful way of obtaining these equations of state. In this formulation, the most natural objects are the quasi-particles. Quasi-particles are parametrized by their internal quantum numbers a (parametrizing the spectrum of the model) and a continuous "rapidity" parameter θ. In this letter we concentrate on Galilean and relativistically invariant models, wherefore θ will be identified with the velocity (Galilean) or the rapidity (relativistic). We will use the combined parameter The fundamental object that complete the full specification of the model is the differential scattering ϕ(θ, θ ′ ), describing the scattering between particles. By relativistic or Galilean invariance, it depends on the rapidities or velocities only through their differences θ − θ ′ . In this paper we keep the discussion general and do not specify any particular model (any particular choice of particle spectrum and differential scattering phase), except when stated otherwise. A conserved charge Q i is characterized, in terms of quasi-particles, by its one-particle eigenvalue h i (θ). It will be convenient to consider the linear space of pseudolocal conserved charges as a function space spanned by the h i s: we will denote Q[h] the conserved charge (a linear functional of h) associated to one-particle eigenvalue h(θ), and likewise q[h] and j[h] for the density and current. The density and current operators are also linear functionals of h 3 . Therefore, in any state (it does not even need to be homogeneous or stationary), the averages q[h] and j[h] are linear functionals of h, and we may consider the kernels ρ p (θ) (a "quasi-particle density") and ρ c (θ) (a "current spectral density") where here and below dθ = a dθ. These kernels are characteristics of the state. One may conveniently introduce the effective velocity v eff (θ) which relates them: The GGE equations of state, which is the requirement, obtained from the TBA quasi-particle picture, of the existence of β such that both q for all h, are the following integral relations for these kernels [19] (here and below prime symbols ( ′ ) represent rapidity derivatives ∂/∂θ): where E(θ) and p(θ) are the energy and momentum of a particle of type a at velocity or rapidity θ. These relations are independent of the state itself, they characterize the family of GGE states for a given model. In terms, instead, of the doublet ρ p (θ) and v eff (θ), the GGE equations of state can be represented as [19] v eff (θ) = v gr (θ) + dα with the group velocity v gr (θ) := E ′ (θ)/p ′ (θ) (that is, (7) and (8) are equivalent under (6)). In this form, the equations of state of integrable systems are seen as equations specifying an effective velocity of quasi-particles, as a modification of the group velocity that depends on both the model and the state.
GGE equations of states mean that ρ p (θ) completely determine the state, as both q[h] and j[h] may be evaluated once it is known. Hence the function ρ p (θ) is a state variable. Other state variables exist. A particularly useful one is the occupation number n(θ) (taking values in [0, 1]). Given n(θ), consider the symmetric bilinear form 4 where the dressing operation is defined by solving Charge densities and currents are expressed in terms of n(θ) as [19] q The nonlinear relation between the state variables ρ p (θ) and n(θ) is 2πρ p (θ) = n(θ)(p ′ ) dr (θ), and the effective velocity takes the simple form (see [19] for more details). Finally, as a consequence of completeness of the set of functions h(θ), the GHD equations (3) can be expressed in various forms, using either state variables: The first form is immediate, and the second form can be derived from the first using the equations of state. The second form, involving occupation numbers, is particularly useful to solve initial-domain-wall problems (again see [19] for details).
Showing that GGE equations of state do indeed emerge in nontrivial integrable models is of course extremely difficult, and will not be discussed here. See however appendix A for a general discussion of how the GGE equations of state may emerge in free-particle models.

III. GGE EQUATIONS OF STATE FROM HYDRODYNAMIC ENTROPY CONSERVATION
It was noted in [19,20] that the density of available states ρ s (θ) = ρ p (θ)/n(θ), which takes the form 2πρ s (θ) := p ′ (θ)+ dα ϕ(θ, α)ρ p (α) = (p ′ ) dr (θ), (15) and the density of holes, defined as ρ h (θ) := ρ s (θ) − ρ p (θ), also satisfy the continuity equation (13) (that is, the equation holds with the replacements ρ p (θ) → ρ s (θ) and ρ p (θ) → ρ h (θ)). Further, the fact that (13) holds for the densities of particles, states and holes with the same effective velocity implies that the Yang-Yang entropy density also follows the same continuity equation. The entropy density is (16) Its integral dθ s(θ) gives the specific entropy of the fluid cell at position x, t (that is, the specific von Neumann entropy of the local GGE). It is found in [19] that The statement (17) provides an interesting physical interpretation of GGE equations of states. Indeed, in ordinary pure hydrodynamics (with finitely-many conservation laws), any function of the state identified as an entropy must obey a similar, natural conservation law; the exact form of the entropy is therefore related to the fluid equations of state 5 . One may then postulate that local conservation of entropy s(θ) is a basic principle, in some way equivalent to the GGE equations of state.
Assume the following.
Let us replace, in the continuity equation for ρ s (θ) , the constitutive relation (15). We obtain (19) Using the continuity equation for ρ p (α), we then find Therefore the expression in the square brackets on the right-hand side of (20) must be independent of x. Since this holds for any x-dependent ρ p (α), it must be independent of it. Using the condition that the limit ρ p (α) → 0 of the effective velocity gives the group velocity, we finally find (8) as claimed.
In relation to the above result, it has recently been pointed out that entropy conservation can be seen as the conservation of an effective Noether current associated to a certain symmetry emerging at late times [38,39]. It would be illuminating to understand if similar concepts can be applied to the specific fluid entropy s(x, t) = dθ s(θ) (note that if we integrate (17) over θ, we obtain a conservation law for the specific entropy s(x, t)). This might be the case, as entropy conservation is a dynamical symmetry, and emerges only when the GHD description becomes sensible. In the context of classical many-body systems, for models that follow trajectories consistent with quasi-static processes in thermodynamics, a symmetry whose conserved charge is the entropy was found recently in [39]. Applying this finding to the present situation might shed some light on the role of entropy in non-equilibrium dynamics.

IV. EQUATIONS OF STATES AND GHD ON COMMUTING FLOWS
In integrable systems, one may consider flows generated not only by the Hamiltonian, but also by any other conserved quantity Q k ; this will be useful when studying the effect of force fields in the next section. The goal of this section is to report on the main equations that generalize GHD to such commuting flows. Since the conserved charges Q k are linear functionals of the one-particle eigenvalues h k , we will also use the notation Let us denote by t k the associated "time", ∂ t k O := i[Q k , O] (with t 1 = t the ordinary time, under Hamiltonian evolution Q 1 = H). By involution, all flows commute, wherefore conserved quantities Q i are also conserved with respect to all t k evolutions. There are associated currents j k,i : which are bilinear functionals of h k and h i , denoted by j k,i = j[h k , h i ] (we will also use the notations j k,i and j[h k , h i ] for averages in GGEs) 6 . Generalized hydrodynamics may also be applied to all these flows. Below we assume that the set of local conserved charges in involution, from which GGEs are formed, with respect to time t k , is the same as that with respect to the original time t -this is usually case in integrable systems. Thus local entropy maximization is described by the same set of GGE states. By commutativity of the flows, under local entropy maximization, local GGE potentials are welldefined functions simultaneously of all time variables, β = β(x, {t}), and we have Note that the currents j k,i are fixed by conservation, (21), only up to the addition of a constant times the identity operator. We fix this gauge freedom, implicitly, by providing explicit expressions for these currents in GGE states below. Bilinearity implies, in general states, the existence of the kernel ρ c (γ, θ) (by abuse of notation, we use the same symbol ρ c as in (5) but with two rapidity arguments in order to represent this new kernel) such that The GGE equations of state encompass relations between this kernel and ρ p (θ), generalizing (7) in a natural manner. This can be obtained following the derivation of [19] and using the results of [19,App D]: This is the most general form of the equations of state, as integrating over γ against E(γ) reproduces the GGE equations of state for the usual time evolution. Likewise, one may define a γ-dependent group velocity v gr (γ, θ) := ∂ θ δ(θ − γ)/p ′ (θ) and a γ-dependent effective velocity and the equations of state (24) Using the bilinear form (9), results of [19, App D] also enable us to express the density and current associated to a conserved charge Q k , in any GGE state parametrized by the occupation number n(θ), as follows 7 : Note that integrating ρ c (γ, θ) (resp. v eff (γ, θ)) against h(γ) gives a current spectral density ρ c [h](θ) (resp. the effective velocity v eff [h](θ)) corresponding to a flow produced by Q[h], and we have with the usual effective velocity being The generalized hydrodynamic problem (22) including all commuting flows of a given integrable model can then be recast as follows. Consider times t h generated by with equations of state (29). Following the derivation of [19], the GHD equations for arbitrary flows can also be written in terms of occupation number variables n(θ), Finally, commuting-flow continuity equations hold for state and hole densities, as well as for the density of the entropy: 7 The second of Equations (27) was explicitly obtained in [19, App D] (see Eq. (D13)), but only for h(θ) with certain propertiescorresponding, in the sinh-Gordon model, to time evolution with respect to local charges. It is natural, however, to assume that the same form holds for any quasi-local charge, and it is under this assumption that (27) is written in this general form.

V. EVOLUTION IN INHOMOGENEOUS FIELDS
It is natural and physically meaningful to consider how external potentials, temperature fields, or inhomogeneous fields associated to other conserved quantities modify the GHD equations (13), (14).
Recall the basic hydrodynamic assumption that averages of local densities and currents, in an inhomogeneous state, may be approximated by averages in a homogeneous, entropy-maximized state, with inhomogeneous potentials. This approximation leads to Euler-type hydrodynamic equations. These equations are first-order differential equations for hydrodynamic variables, which may be taken as the densities, or as the inhomogeneous potentials themselves. The hydrodynamic assumption is expected to be a good approximation when variations of densities and currents occur on large scales, so that locally, the state looks homogeneous. The Euler-type hydrodynamic equations are in fact the leading terms in a derivative expansion; higher derivative terms would give rise to viscosity and other effects. It is within this picture that we may consider external inhomogeneous fields affecting the time evolution. We assume that these external fields also only display variations on large scales, so that locally the evolution looks homogeneous. The Euler-type hydrodynamic equations obtained will therefore again be leading-order terms in a derivative expansions, neglecting any term containing derivatives of hydrodynamic variables or potentials with a total order of 2 or more. The equations are simply obtained by deriving leading-order evolution equations at the operator level, and then using the hydrodynamic approximation of local entropy maximization.
Inhomogeneous fields, of course, break the integrability of the dynamics. However, since locally the dynamics still looks like a homogeneous evolution with respect to an integrable Hamiltonian, the local GGE approximation stays valid under time evolution. This is made clear below, as we show that the first-derivative-order approximation of the dynamics leads to a consistent equation for local GGE states (the initial state does not know about the dynamics, and thus can be chosen within the space of local GGE states). This is certainly not surprising. In usual hydrodynamics (such as for water waves), local fluid cells are approximated by Galilean boosts of equilibrium states, thus involve the momentum operator. In an inhomogeneous field, this description still holds and Euler and Navier-Stokes equations with force terms give a good description. This is so even though inhomogeneity breaks translation invariance (thus the momentum operator is not a conserved quantity of the dynamics). The derivation below is simply a generalization of this fact to infinitely-many conserved charges. Naturally, the higher-order derivatives terms neglected would modify this picture, and may be expected to lead to integrability breaking effects at large times. But this is beyond the scope of this work.
To start with, let us briefly recall a typical case in relativistic one-dimensional quantum field theory with U (1) symmetry: coupling the particle current J µ (x, t) to an external electric field, is the electric potential 8 . Here in order to fix the notation, we assume the particle current is associated to some conserved charge Q 0 (that is, J 0 (x, t) = q 0 (x, t) and J 1 (x, t) = j 0 (x, t))), and we take Q 1 = H to be the total energy without external field. The external field deforms the evolution Hamiltonian in a familiar fashion: Accordingly the hydrodynamic conservation equations become [42] (keeping the (x, t)-dependence implicit), at the first order in a derivative expansion, where T µν is the energy-momentum tensor, F 01 = −F 10 = ∂ x V 0 and F 00 = F 11 = 0, and averages are taken in local fluid cells. Alternatively this can be written as We now generalize this, as well as more complicated external fields, to GHD. In order to have a clearer general framework, we divide the external field, in general, into two types. We first understand an external force field, arising from a potential V 0 (x), as a field coupled to a conserved density q 0 (x) = q[h 0 ](x) which has the property the associated conserved charge Q 0 = dx q 0 (x) commutes with all conserved densities q i (x): This is a sensible definition of an external potential V 0 , as it implies that physical quantities in GGEs only depend on potential differences.
and as a consequence of (39), evolution of local densities with respect to H + V 0 Q 0 , and averages of local densities with respect to density matrices of the form e − i βiQi−V0Q0 , are independent of V 0 . Note that thanks to (39), all currents associated to the Q 0 evolution must vanish 9 , j 0,i = 0. Therefore, using (27), the one-particle eigenvalue h 0 (θ) must be independent of the rapidity, h 0 (θ, a) = h 0 (a) (that is, h ′ 0 (θ) = 0). As an example, in the Lieb-Liniger model (a Galilean model with one particle type only) one may choose Q 0 to be the number operator, which counts the number of quasi-particles, h 0 (θ) = 1. In a model with an internal charge a ∈ {+1, −1}, such as the (relativistic) sine-Gordon mode, one may take Q 0 to be the total charge, with h 0 (θ, a) = a.
We are thus interested, in a first instance, in deriving a force-field, pure hydrodynamic equation describing the time derivative of local conserved densities under the time evolution with respect to the force-field Hamiltonian, We show in Appendix B that the infinite set of force-field hydrodynamic equations, under the assumption both of local entropy maximization and of weak spacial variations of the potential V 0 (x), are We see that the force term, proportional to the space derivative ∂ x V 0 of the potential, involves the charge current associated to the time evolution with respect to Q i (see (22)). Specializing to the energy Q 1 = H (choosing i = 1), we observe that the force term controlling the continuity equation for the energy density is proportional to the usual particle current j 1,0 = j 0 , as is intuitively clear and in agreement with (38). Equation (41) is to be seen as the leading part of a derivative expansion, where neglected terms are higher space derivatives in the potential and in conserved densities and currents. In a second instance, we consider more general external fields, associated to general conserved densities. These are perturbations of the type dx k V k (x)q k (x): For instance, as q 1 (x) is the energy density (according to our convention), the term dx V 1 (x)q 1 (x) may be understood as a perturbation by an inhomogeneous temperature field, with x-dependent temperature (V 1 (x)) −1 (this interpretation being valid under the hydrodynamic assumption, with weak variations). It is useful to introduce the one-particle potential which is the one-particle eigenvalue function of the operator k V k (x)Q k (in this notation, W (x) is, implicitly, a function of θ). Using q k (x) = q[h k ](x), the perturbation is written in a somewhat more general way in terms of any W (x): We show in Appendix B that the infinite set of hydrodynamic equations in inhomogeneous fields, again under the assumption both of local entropy maximization and of weak spacial variations of the potentials V k (x) (weak spacial variations of W (x)), are These generalize (41), which is recovered by choosing V k (x) = 0 for all k ≥ 1 and using j 0,i = 0.
Equations (41), (45) and (46) are derived without invoking integrability, which is only included in the fact that there are infinitely-many of these equations. They are valid under the following assumptions. First, there is the usual hydrodynamic assumption that local averages of densities and currents are well approximated by averages within local GGEs. Second, all higher-derivative terms occurring from the quantum dynamics with respect to H field are neglected. These are terms composed of products of the first or higher derivative of the potentials V k (x), times local fields and their derivatives, with, in total, two or more space derivatives. As long as the potentials are varying in a smooth enough fashion, such higher-derivative terms are indeed negligible. We recall that the assumption of local GGEs (which comes from that of local entropy maximization) gives rise to a continuity equation, which is a first-derivative equation. It therefore only gives the leading first-derivative terms in a derivative expansion of the full hydrodynamics (thus neglects higher derivatives of hydrodynamic variables such as viscosity terms). Assuming that variations of the potentials are of the order of the variations of the hydrodynamic variables, it is thus consistent to neglect higher derivative terms as above, and to keep the total number of derivatives, of hydrodynamic variables and potentials, to a maximum of 1. Of course, such derivative expansions, common in hydrodynamic problems, are not controlled approximations, and it is difficult to evaluate the corrections.
Further, we show in Appendix B that (45), (46) can be recast, in the quasi-particle basis, into the following equivalent equations for the occupation number n(θ) and for the densities (here keep implicit the x and t dependencies), which holds for ρ = ρ s , ρ p and ρ h . Recall that E is the function of θ giving the one-particle energy (the Hamiltonian one-particle eigenvalue). Here the effective acceleration is The space-dependent force function F (θ) is the derivative of the total energy E(θ) + W (θ) with respect to space; since E(θ) is independent of space, this is The effective velocity v eff [E + W ](θ) depends on x both through the one-particle potential W (θ) = W (x; θ), which modifies the local energy to E(θ) + W (x; θ), and thus modifies the local group velocity; and through the (x, t)-dependent occupation number n(θ), or particle density ρ p (θ), which determines it (see (29) and (31)). Likewise, the effective acceleration depends on x both through the potentials and through the dressing operation.
Equations (47) and (48) invoke integrability in the use of the TBA formalism, and of the completeness of the space of pseudolocal conserved charges. They are otherwise both direct consequences of (45) (or equivalently (46)), without further approximation. They use the quasi-particle expressions of the local GGE densities and currents that are involved in (45), (46).
We see that the effects of the potential W (x) (or equivalently V k (x)) are twofold. First, there is a modification of the effective velocity to , which takes into account the local potential W at the position x. Second, there is an extra term involving θ derivatives, which takes into account the acceleration due to spacial variations of W around the position x. We note that since v eff [E + W ](θ) only involves θ-derivatives W (x) ′ of the one-particle potential, and since h ′ 0 = 0, it is clear that the force-field potential V 0 (x) does not affect the effective velocity. A force field only leads to an acceleration, without modifying the local effective velocity. Other external fields such as temperature fields, however, do modify the local effective velocity.
Consider a pure force field in a Galilean model with a single-particle spectrum (such as the Lieb-Liniger model). In this case, we have θ = θ, h 0 (θ) = 1 and p(θ) = mθ. Then, the effective acceleration a eff (θ) simplifies to the usual acceleration, independently of θ, a eff (θ) = −∂ x V 0 /m (Galilean, single-particle spectrum, pure force field).
(51) Equation (47) (or equivalently (48)) represents evolution in the presence of space-dependent external fields; it is valid in the limit of weak variations of both the hydrodynamic quantities and of the potentials themselves. As it is a pure-hydrodynamic equation, it does not take into account any viscosity effects, which give rise to terms with higher derivatives of the hydrodynamic variables, or, similarly, any effect related to the presence of nonzero higher derivatives of the potentials. In a pure force field, V k≥1 = 0, the effective velocity is not affected, and if the force field is constant, ∂ 2 x V 0 (x) = 0, the effective acceleration does not depend on space. In this case, one may argue that, as usual, at large times variations of hydrodynamic variables become smaller, and thus pure hydrodynamics provides a good description 10 . Otherwise, spacial variations of potentials are present in the pure hydrodynamic equations, and as they do not change with times, they will fix a minimum spacial-variation scale for the hydrodynamic variables. Thus, in this case, the pure hydrodynamic equations cannot become more accurate at large times, and we must understand (47) as being valid for a finite period of time, whose extent depends on the size of spacial variations of the potentials. During this time, discrepancies between the predictions of (47) and the actual evolution, due to neglected terms whose amplitude does not decay, accumulate. Beyond this time, one might expect the integrability-breaking effects of the presence of space-varying potentials to become important.
Let us now investigate stationary solutions to the force-field equations (47). Consider a fluid state which is, at every position x, the Gibbs state associated to H + k V k (x)Q k at the temperature β −1 (independent of x). This is the local density approximation of the finite-temperature, inhomogeneous state e −βH field . We show that the one-parameter family of such local-Gibbs states, parametrized by the temperature β −1 , is indeed a stationary solution to (47).
Clearly ǫ(θ) satisfies the same equation (47) as does n(θ). Note that ∂ x ǫ(θ) = (∂ x w) dr (θ), and that, using the fact that ϕ(θ, α) depends on the rapidities through their difference θ−α only, ∂ θ ǫ(θ) = (∂ θ w) dr (θ) (we recall that the superscript dr indicates dressed quantities as per (10)). Using these statements and setting ∂ t n(θ) = 0, one finds that in terms of the local-GGE one-particle eigenvalue w(θ), a stationary solution satisfies the equation (we also used (49), (15) and (12)). It is simple to see that is a solution to this equation for any β (recall that E = E(θ) depends on θ = (θ, a) but not on x, and that W = W (x) = W (x)(θ) depends on both x and θ). This is the local density approximation of the state e −βH field . The above statement is very natural physically. Assuming that the spatial variations of the potential occur only on large distance scales, we do not expect these inhomogeneities to lead to localization (the latter would of course break the hydrodynamic assumptions). Yet, we expect inhomogeneities to lead to integrability-breaking effects. Therefore, at very large times, after integrabilitybreaking effects have arisen, we expect the stationarystate density matrix to be of the thermal form e −βH field for some β. In such a state, variations of all densities and currents occur on large distance scales, hence this is well approximated by a local density approximationthe "fluid form" of this state. Further, since the forcefield hydrodynamic equations should approximate well the dynamics when variations are on large scales, we expect this approximation to be a stationary solution to these equations, as indeed shown above. That is, although the force-field equations do not contain all integrability breaking effects and might not by themselves lead to thermalization, the thermalized state should be stationary with respect to it. This is much like the fact that the ideal-gas distribution is invariant under the freeparticle evolution, although it may only arise, physically, as a consequence of the small interactions between the particles of the gas.
We have not established uniqueness of this stationary solution to the force-field equations -in particular, it is simple to see that in the case of free-particle models, any function f (E+W ) is a solution. One may wonder if, similarly, there are additional stationary solutions in interacting integrable models, and if these make physical sense. One may also wonder what, if any, stationary solution is actually reached at long times from solving the pure hydrodynamic equations (47) without higher-derivative terms. If it is not the local-Gibbs state above, then this might correspond to a "pre-thermalization" plateau, which appears before the integrability-breaking effects of the inhomogeneous potential become important. We leave these questions for future works.

VI. EULER AND NAVIER-STOKES EQUATIONS
An important ingredient in conventional hydrodynamics is what is often referred to as the Euler equation: this is a continuity equation relating the fluid velocity v to the internal pressure P and the fluid's mass density ρ fl : It is a simple consequence of conservation of the mass density and mass current, and expresses the variation of the fluid's velocity as a convection term and a term due to pressure variations. In generalized hydrodynamics, such equations also arise in a natural fashion. It is obvious from the symmetry of the bilinear form (9) that, in any GGE state, the current associated to the conserved quantity with oneparticle eigenvalue h(θ) = p ′ (θ) is equal to the density associated to h(θ) = E ′ (θ): For instance, in Galilean invariant systems, p ′ (θ) = m a is the mass of the particle, and E ′ (θ) = p(θ) is its momentum, and this is equality between mass current and momentum density. In relativistic system, p ′ (θ) = E(θ) and E ′ (θ) = p(θ), so this is instead equality between energy current and momentum density (which amounts to the fact that the energy-momentum tensor is symmetric). Let us then define the fluid velocity v as follows: This is the velocity for the mass current (Galilean) or energy current (relativistic). The quantity v depends on x and t (but is of course independent of θ). Conservation We may then define the fluid mass density and pressure as and we recover (56), using E ′ (θ) = p(θ). The interpretation of the above identification is particularly clear with Galilean invariance. In this case ρ fl is indeed the physical mass density, and the second equation in (60) is the correct relation between the momentum current j[p] and the pressure P: it identifies the momentum current as the internal pressure plus v times the current associated to the displacement of the fluid cell ρ fl v. In the relativistic case, ρ fl is the energy density, and P has a similar interpretation. Notice that the physical pressure P is not equal to the generalized specific free energy (free energy per unit volume) f = dp(θ)/(2π) log(1 + e −ǫ(θ) ) (where the pseudoenergy is defined in (52)); this is of course natural in states that are not thermal Gibbs states. As such, unlike the case in conventional hydrodynamics, in the Galilean case the continuity equation for the energy ∂ t q[E] + ∂ x j[E] = 0 is no longer expressible in terms of the fluid velocity and thermodynamic variables.
It is also straightforward to generalize the above to the forced equation (41). Repeating the above derivation with the conservation laws ∂ t q[p ′ ] + ∂ x j[p ′ ] + j[p ′ , h 0 ] = 0 and ∂ t q[E ′ ] + ∂ x j[E ′ ] + j[E ′ , h 0 ] = 0, and using the following identities (see (11) and (27)): and (63) where q 0 is the charge density and j 0 is the charge current (the densities and current of the charge Q 0 associated to the force term). In the Galilean case with a singleparticle spectrum (such as the Lieb-Liniger model), we have ρ fl = mq 0 and thus we find the usual forced Euler equation, (64) So far we have considered only ideal fluids, that do not account for viscosity. An accurate consideration of viscosity terms corresponding to the underlying many-body model requires an analysis of how the unitary dynamics approaches pure hydrodynamics. However, one may consider a simple, possible correction to (13) that could account for the presence of viscosity effects. Let us exemplify in the Galilean case with a single-particle spectrum. From standard hydrodynamic arguments, the Navier-Stokes equation in one-dimensional non-relativistic systems reads where ζ is the (mass-normalized) bulk viscosity (note that we do not have the kinematic viscosity as there occurs no shear flow in one dimension). A continuity equation for ρ p (θ) that gives the above Navier-Stokes equation is This might or might not correspond to any underlying quantum model, but in any case it could provide a way of regularizing the GHD equations for numerical purposes 11 . It would be interesting to analyze further such viscosity corrections.

VII. CONCLUSIONS
In this letter, we further developed the generalized hydrodynamics (GHD) theory first proposed in [19]. We showed that the GGE equations of state, at the basis of GHD, follow from a principle of hydrodynamic conservation of entropy. We provide in Appendix A arguments for the emergence of GHD in general free-particle relativistic models under smoothness assumptions, which we expect could be extended to interacting models using the form factor program. Then, we generalized to flows generated by arbitrary conserved charges, and employed this in order to establish the conservation equations (45), (46) and continuity equations (47), (48) within an external field, be it a force field, a temperature field or any other field associated to conserved quantities of the model. We expect that these equations should effectively capture the large-scale (long-wavelength) dynamics of a Lieb-Liniger model in an external potential, such as a harmonic potential (see e.g. [43]). This, we believe, is particularly interesting: indeed, despite a lack of full justification, conventional hydrodynamics has been exploited in analyzing the quench dynamics of one-dimenional bose gases in a trap potentials [44][45][46], and we believe our equations might lead to more accurate results. In particular, the consideration of all conservation laws in the forcefield GHD might give rise to a more accurate theoretical description of the notable "quantum Newton's cradle" experiment [47]. All equations hitherto derived within GHD are, however, for ideal fluids: no dissipation effect has been taken account of. For a precise treatment one has to add viscosity terms. We proposed one possibility from considering the Navier-Stokes equation, but we expect a more in-depth study will be necessary in order to clarify this aspect.

VIII. ACKNOWLEDGEMENT
TY and BD are indebted to B. Bertini, M. Fagotti and especially J. Dubail for a careful reading of the manuscript and for discussions. BD thanks F. Essler, as well as the group of Statistical Physics at the Scuola Internazionale Superiore di Studi Avanzati (SISSA), Trieste, Italy, for discussions. TY and BD are grateful to J.-S. Caux for sharing ideas during the conference "Entanglement and non-equilibrium physics of pure and disordered systems", International Centre for Theoretical Physics (ICTP), Trieste, Italy. BD thanks H. Spohn for discussions during the conference "Non-equilibrium dynamics in classical and quantum systems: from quenches to slow relaxations", Pont-à-Mousson, France. TY acknowledges support from the Takenaka Scholarship Foundation for a scholarship. BD acknowledges support from SISSA, ICTP, the University of Bologna and the University of Lorraine where parts of this work were done.
Appendix A: Emergence of GGE equations of state in free-particle models The problem of showing the emergence of hydrodynamics in many-body systems is notoriously difficult, see [40,41] for recent progress. This is particularly true because usual hydrodynamics requires, as per its principles, strong interactions, by their nature hard to treat analytically. The interaction should provide the mixing necessary in order for all degrees of freedom that do not follow a conservation law to thermalize; thus minimizing, locally, the free energy under the conditions of all local conservation laws, and rendering applicable, locally, the thermal equations of state. As explained in [19], the sole assumption at the basis of GHD is, likewise, the emergence, in a uniform enough fashion and at large enough times, of the GGE equations of state at every point in space-time. GHD follows from this, simply by combining it with the conservation equations (1) of the model's unitary dynamics. In this respect, GHD offers a unique opportunity in that it accounts for infinitely-many conservation laws: as a consequence much less interaction effects, or mixing, is required for the emergence of the GGE equations of state. This is particularly evident in "quadratic models", or models whose asymptotic particles do not interact. In such models, GGE equations of states should still emerge, although the interaction between fundamental degrees of freedom is quadratic and amenable to exact treatment. Thus, in these models, we may analyze with much more depth these fundamental principles, making use of the large simplification afforded by the triviality of the scattering matrix.
An important question is therefore what basic properties either of the initial state or of the large-time evolution guarantee that the GGE equations of state emerge in free-particle models. Although hydrodynamic ideas have been used successfully in the past in such models [27][28][29][30][31][32], to our knowledge, no general assessment of such conditions for the emergence of GGE equations of states, or of hydrodynamcis, have been provided 12 . In this section we propose such conditions. We provide arguments to the effect that, under homogeneous time evolution, if densities and currents become, at long times, smooth enough in space-time, with a variation scale growing unboundedly with time, then the GGE equations of state and GHD emerge. In other words, we show that GGE equations of state hold in homogeneous, stationary states; and if the size of fluid cells, wherein uniform near-homogeneity and near-stationary hold, grow with time, then GGE equations of states are approached and GHD becomes increasingly accurate.
A free-particle model is characterized by the fact that ϕ(θ) = 0. For simplicity and clarity, in the following we specialize to the case of a single relativistic particle, but the derivation below can be generalized straightforwardly (to many particles, and to other dispersion relations). Let us therefore consider some initial state · · · , and let us evaluate in this state observables evolved in time: Of course, it cannot be expected in general that GGE equations of state emerge for any initial state, as cases where hydrodynamics fail certainly exist. Hence we need a condition which will guarantee that such "pathological" cases are avoided. A natural condition is the requirement that the long-time limit be smooth enough.
We first assume that everywhere in space-time (at positive times), average densities and currents O(x, t) stay uniformly finite. Let us also assume that, as time t becomes large, and uniformly within some region R of space-time that is unbounded in the positive time direction, average densities and currents display order-1 variations in space-time on lengths scales that diverge as t grows. We express this latter assumption more precisely by considering averages over Gaussian cells centered at x, t of extent T = T (t): Then the assumption is that there is a T = T (t) growing unboundedly with time, such that lim λ→0Ō (x, t; λ) = O(x, t) uniformly on (x, t) ∈ R.
(A3) For any finite (x, t), it is clear that the limit is as above; the assumption is that this holds uniformly in R, this being most nontrivial in the long-time subregion of R.
Then, under this assumption, we argue in below that the GGE equations of state emerge uniformly at long times in R 13 . In order to make this conclusion more precise, recall that the averages of conserved densities q[h] and currents j[h] associated to one-particle eigenvalue h(θ) are linear functionals of h as per (5): For a generic state and generic x, t, the densities ρ p (θ; x, t) and ρ c (θ; x, t) are functionally not related to each other. The emergence of the GGE equations of state is the statement of the emergence of the relation (7), or equivalently (6) with (8) (or (12)). In free relativistic particle models, this is particularly simple as the effective velocity is the group velocity v gr (θ) = tanh θ: the relation is ρ c (θ; x, t) − v gr (θ)ρ p (θ; x, t) = 0. We show that this relation emerges uniformly in the region R as t → ∞: A sketch of the proof is as follows. Let j[h](x, t) be the current associated to the GGE determined by quasiparticle density ρ p (θ; x, t). Then this implies that the difference j[h](x, t) − j[h](x, t)) goes to zero uniformly as above. This gives rise to the integral form of conservation equations, with uniform correction terms that are smaller than the total length of the path: We therefore find the emerging hydrodynamic conservation equations, in integral form, for hydrodynamic variables. Assuming differentiability, this implies the differential form (3).
We finally note that we may apply the above result to the case where the state is stationary and homogeneous. In this case, it is clear that the assumption is fulfilled, and we conclude that in such states, be them GGE states or not, averages of local densities and currents must be reproducible by a GGE.
In a free particle model, average densities and currents take are bilinears in terms of canonical annihilation and creation operators A(θ), A † (θ). Therefore, they take the following general form are linear functionals of h (here and below indices in A j , E j and p j represent the rapidity argument θ j , and E j is the energy and p j the momentum). Recall that · · · represents the initial state.
In specific models, it is a simple matter to evaluate the coefficients b[h](θ 1 , θ 2 ),b[h](θ 1 , θ 2 ), c[h](θ 1 , θ 2 ) and c[h](θ 1 , θ 2 ) explicitly. In some simple free-fermionic models, these coefficients may be simple enough to guarantee that, with Galilean invariance, the hydrodynamic equations hold exactly independently of the initial state and at all times [48]. However, here we leave these coefficients as general as possible, and impose only conditions that arise from general principles.
we will conclude that we must have b ≥ −1; the distribution A † 1 A 2 may also containt a delta-function term of the type f (θ 1 )δ(θ 1 − θ 2 ) with f (θ) decaying fast enough at infinity. We consider Gaussian-cell averages of densities, q[h](x, t; λ) (see (A2)) (the same conclusion is obtained using currents instead of densities). This should stay finite, in particular, with T = t, λ = 1 and x = 0, as t → ∞. We use and similarly for the y integral in (A2), as well as the mode expansion (A7). We see, from the fact that E 1 +E 2 is always positive, that all terms in (A7) involving A 1 A 2 and its hermitian conjugate will have exponentially decaying contributions as t → ∞. The remaining terms are where p 12 := p 1 − p 2 and E 12 := E 1 − E 2 . We recall that b[h](θ 1 , θ 2 ) is regular at θ 1 = θ 2 . We further assume that it behaves well enough at large rapidities, so that we do not worry about the large-rapidity region of the integrals.
The eventual delta-function term in A † 1 A 2 leads to a finite contribution to (A12) by our assumptions concerning behaviors at infinite rapidities. On the other hand, at large t, the algebraic contribution of A † 1 A 2 to (A12) can be analyzed by a stationary phase argument. Setting u := p 2 and w := p 2 2 + E 2 2 , the stationary phase occurs at θ 1 − θ 2 =: θ 12 = θ ⋆ := iu/(wt) + O(1/t 2 ). Keeping only up to the quadratic terms in θ 12 − θ ⋆ in the exponential and using A † (A13) Finiteness thus requires b ≥ −1.
Now consider δ δh(θ)j [h](x, t; λ) − tanh(θ) δ δh(θ)q [h](x, t; λ) (A14) for some T = T (t) in (A2) that grows unboundedly with t. Again, we see that all terms in (A7) involving A 1 A 2 and its hermitian conjugate will have exponentially decaying contributions in (A14) as t → ∞. Terms involving A † 1 A 2 , on the other hand, are of the form (A15) We may bound this integral by replacing the oscillatory factor e iE12−ip12x by 1. At large T this can then be analyzed by a stationary phase argument. The position of the stationary phase is exactly θ 1 = θ 2 , hence the main contribution occurs around θ 1 ≈ θ 2 . Thanks to (A8), we find Therefore, the delta-function part of A † 1 A 2 does not contribute to the integral (A15), and the algebraic contribution becomes, as t → ∞, (A17) This is clearly uniform on (x, t) ∈ R. Since b ≥ −1, as a consequence, we have found that lim t→∞ δ δh(θ)j [h](x, t; λ) − tanh(θ) δ δh(θ)q [h](x, t; λ) = 0 (A18) uniformly on R. By the assumption (A3), this is sufficient to show (A10). This is of course far from being a complete or rigorous proof. For instance, we have omitted the discussion of how the assumption (A3) is uniform with respect to the observables O themselves (allowing us to take h(θ)derivatives). We have also omitted the detailed dependencies on θ 1 , θ 2 in expressions of the form O(θ c 12 ), while these are important to make sure that the rapidity integrals are finite. In addition, of course, the stationary phase arguments, while treated with some care, would need to be developed in order to become rigorous. Nevertheless, we believe this provides the main arguments, and shows how GGE equations of state may indeed emerge.

Appendix B: Derivation of hydrodynamic equations within inhomogeneous fields
In order to describe the first part of the result, equation (41), consider the conservation law of the conserved density q i with respect to the time evolution generated by a conserved quantity Q k , GGE averages of the associated currents can be evaluated using (27) as j k,i = j[h k , h i ], which, thanks to (27), takes the explicit form Equation (46) (which implies (41)) is shown as follows. Locality of densities imply that there exists a field O j,i (x, y) supported at x = y (i.e. local at this position) such that Since q j (x) and q i (x) are local conserved densities, they are not affected by any nontrivial renormalization, and therefore O i (x, y) can be written as a finite sum of terms with increasing derivatives of the delta function, where O k,i;ℓ (x) are local fields. Integrating over y, by (B1) we find that On the other hand, integrating over x, we obtain and therefore comparing with (B1) we can make the following identification, using the fact that the only local fields whose derivative is zero are those proportional to the identity: It can be seen to vanish as follows. We write it as the following quantity, involving an averages · · · in any GGE: By symmetry, this constant is zero whenever q k and q i are both parity symmetric or parity anti-symmetric, or whenever their combined transformation under some internal symmetry is nontrivial. One can argue this constant should in fact be identically zero as follows. h(θ) and g(θ) that decay fast enough at infinite rapidities. Let us also consider the GGE · · · to be a thermal state in the limit of large temperatures. In this limit [36], the occupation number n(θ) has a large flat plateau, and decays to zero beyond this plateau. The regions where it starts decaying to zero are further and further away from θ = 0 as the large temperature limit is taken. Therefore, in (10) and (9), for h, g as above, we may consider n(θ) to be a constant, independent of the rapidity. Hence by integration by parts, we have (h ′ ) dr = (h dr ) ′ . Thus, using (27) and (9) (and its symmetry property), we have That is, in this limit j[h, g] + j[g, h] = 0. Further, in the infinite temperature limit the state is the trace state, which has the cyclic property AB = BA . As a consequence 14 , in this limit [q[g](x), q[h](0)] = 0. Therefore, since A[h, g] is independent of the state, we must have A[h, g] = 0. We thus conclude that this is the zero bilinear functional, and thus A i,k = 0 for all i and k. Note that one can further check that the result (B7) with A k,i = 0 agrees, in the case where q i and q k are either energy or momentum densities, with the firstderivative terms of the commutators of the stress-energy tensor calculated in [49].
We can then compute the time evolution within the inhomogeneous field as follows: 14 Taking the infinite-temperature limit in QFT is delicate, as large temperatures bring the system much beyond the quantum critical point. However, choosing h and g to decay at large rapidities amounts to a UV regularization of the fields q[h](x) and q[g](x) (which are therefore not local anymore). This UV regularization guarantees that the energy scale of the temperature, in the largetemperature limit, is beyond the UV scale of the observables, and thus the limit is indeed described by the microscopic formula, which is a trace state.
where W (x) = k V k (x)h k is the one-particle externalfield function (for every x, it is a function of θ). We have used integration by parts, assuming that boundary terms at infinity do not contribute. The terms omitted are "higher-derivative terms": they are composed of products of the first or higher derivative of the potentials V k (x) times local fields and their derivatives, with, in total, two or more space derivatives.
Integrating over a large space-time cell Ω, we obtain the integral form of a conservation equation, ∂Ω d x ∧ q( x) = S Ω , x = (x, t), q = (q, j), for the density q = q i and the modified current j = j i + j[W, h i ], with sources within the cell, S Ω = Ω dxdt j([h i , ∂ x W ](x, t). We may now make the hydrodynamic assumption that averages of local observables are evaluated in local GGEs, and reverting to the differential form of this conservation equation, this shows (46).
In a pure force field, i.e. with W (x) ′ = 0, the equation simplifies, as in this case j[W (x), h i ](x) = 0. For evolution within a pure force field, we are therefore left with (B10) which implies (41).
In order to show (47) and (48), Equation (B9) is written, using TBA and in particular using (B2) and the symmetry of the bilinear form (9), as (here for lightness of notation, we omit the explicit θ and x dependences, and recall that primes ( ′ ) indicate θ-derivatives). Using integration by parts for the last term in the square brackets, and using the fact that this holds for every function h i (assuming completeness of this space of functions), we obtain Let us use integral-operator notations, with measure dθ/(2π). Consider the diagonal operator N with kernel N (θ, α) = 2π n(θ)δ(θ − α)δ a,b , the vectors p ′ , E ′ and h 0 with elements p ′ (θ), E ′ (θ) and h 0 (θ) respectively, and the operator ϕ with kernel ϕ(θ, α). Then Using the first, second and last of these relations, as well as (29), we see that we can combine the term ∂ x (v eff ρ p ) with ∂ x (V k n (h ′ k ) dr ) into ∂ x (v eff [E+W ]ρ p ). On the other hand, using the first and the third, as well as (49), we see that k ∂ x V k (nh dr k ) ′ = 2π∂ θ (a eff ρ p ). Therefore, this indeed reproduces (48).