Generalized hydrodynamics with space-time inhomogeneous interactions

TL;DR

This paper introduces a new hydrodynamic framework extending Generalized Hydrodynamics to handle systems with space-time varying interactions, enabling analytical descriptions of complex out-of-equilibrium dynamics.

Contribution

It develops a novel analytical approach for inhomogeneous interactions in integrable systems, applicable to both quantum and classical models, advancing the understanding of non-equilibrium phenomena.

Findings

Successfully describes interaction quenches in trapped Bose gases

Benchmarks results in XXZ spin chain and Sinh-Gordon model

Enables analytical treatment of space-time dependent interaction modulations

Abstract

We provide a new hydrodynamic framework to describe out-of-equilibrium integrable systems with space-time inhomogeneous interactions. Our result builds up on the recently-introduced Generalized Hydrodynamics (GHD). The method allows to analytically describe the dynamics during generic space-time-dependent smooth modulations of the interactions. As a proof of concept, we study experimentally-motivated interaction quenches in the trapped interacting Bose gas, which cannot be treated with current analytical or numerical methods. We also benchmark our results in the XXZ spin chain and in the classical Sinh-Gordon model.

Equations

$$\underset{\text{TL}}{\lim }\frac{1}{L}⟨\{\lambda {\}}_{i=1}^N|{\hat{\mathcal{Q}}}_j|\{\lambda {\}}_{i=1}^N⟩=\int \text{d}\lambda q_j(\lambda )\rho (\lambda ),$$

$${\tau }^{\text{dr}}(\lambda )=\tau (\lambda )-\int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }\Theta (\lambda -\mu )\vartheta (\mu ){\tau }^{\text{dr}}(\mu ),$$

$$\hat{H}(\alpha )=\int \text{d}x\hat{\mathbf{h}}(x,\alpha (t,x)).$$

$${\partial }_t\rho +{\partial }_x(v^{\text{eff}}\rho )+{\partial }_{\lambda }\left(\frac{{\partial }_t\alpha f^{\text{dr}}+{\partial }_x\alpha {\Lambda }^{\text{dr}}}{({\partial }_{\lambda }p{)}^{\text{dr}}}\rho \right)=0$$

$$f(\lambda )=-{\partial }_{\alpha }p(\lambda )+\int \frac{\text{d}\mu }{2\pi }{\partial }_{\alpha }\Theta (\lambda -\mu )({\partial }_{\mu }p{)}^{\text{dr}}\vartheta (\mu ),$$

$$\Lambda (\lambda )=-{\partial }_{\alpha }ϵ(\lambda )+\int \frac{\text{d}\mu }{2\pi }{\partial }_{\alpha }\Theta (\lambda -\mu )({\partial }_{\mu }ϵ{)}^{\text{dr}}\vartheta (\mu ).$$

$${\hat{\mathcal{Q}}}_j=\int \text{d}x{\hat{\mathbf{q}}}_j(x,\alpha (t,x)).$$

$${\hat{\mathbf{q}}}_j^{\text{H}}(t+\text{d}t,x,\alpha (t+\text{d}t,x))=e^{i\text{d}t\hat{H}(t)}{\hat{\mathbf{q}}}_j(x,\alpha (t+\text{d}t,x))e^{-i\text{d}t\hat{H}(t)}.$$

$${\partial }_t{\hat{\mathbf{q}}}_j^{\text{H}}(t,x,\alpha (t,x))={\partial }_t\alpha (t,x){\partial }_{\alpha }{\hat{\mathbf{q}}}_j(x,\alpha (t,x)+i\int \text{d}y[\hat{\mathbf{h}}(y,\alpha (t,y)),{\hat{\mathbf{q}}}_j(x,\alpha (t,x))]).$$

$$\begin{gathered} i\int \text{d}y[\hat{\mathbf{h}}(y,\alpha (t,y)),{\hat{\mathbf{q}}}_j(x,\alpha (t,x))]=i\int \text{d}y[\hat{\mathbf{h}}(y,\alpha (t,x)),{\hat{\mathbf{q}}}_j(x,\alpha (t,x))]+ \\ i\int \text{d}y(x-y){\partial }_x\alpha (t,x)[{\partial }_{\alpha }\hat{\mathbf{h}}(y,\alpha (t,x)),{\hat{\mathbf{q}}}_j(x,\alpha (t,x))]+.... \end{gathered}$$

$${\partial }_t{\hat{\mathbf{q}}}_j(x,\alpha (t,x))+{\partial }_x{\hat{\mathbf{j}}}_j(x,\alpha (t,x))-{\partial }_t\alpha (t,x){\partial }_{\alpha }{\hat{\mathbf{q}}}_j(x,\alpha (t,x)-{\partial }_x\alpha \hat{\Phi }(t,x)=0.$$

$${\hat{\Phi }}_j(t,x)=i\int \text{d}y(x-y)[{\partial }_{\alpha }\hat{\mathbf{h}}(y,\alpha (t,x)),{\hat{\mathbf{q}}}_j(x,\alpha (t,x))].$$

$$⟨{\partial }_t{\hat{\mathbf{q}}}_j^{\text{H}}(t,x,\alpha (t,x))⟩\simeq {\partial }_t⟨{\hat{\mathbf{q}}}_j(x,\alpha (t,x)){⟩}_{\text{GGE}(t,x)},⟨{\partial }_x{\hat{\mathbf{j}}}_j^{\text{H}}(t,x,\alpha (t,x))⟩\simeq {\partial }_x⟨{\hat{\mathbf{j}}}_j(x,\alpha (t,x)){⟩}_{\text{GGE}(t,x)}.$$

$${\partial }_t⟨{\hat{\mathbf{q}}}_j⟩+{\partial }_x⟨{\hat{\mathbf{j}}}_j⟩-{\partial }_t\alpha ⟨{\partial }_{\alpha }{\hat{\mathbf{q}}}_j⟩-{\partial }_x\alpha ⟨{\hat{\Phi }}_j⟩=0.$$

$$⟨{\hat{\mathbf{q}}}_j(x,\alpha (t,x)){⟩}_{\text{GGE}(t,x)}=\int \text{d}\lambda q_j(\lambda ,\alpha (t,x))\rho (t,x,\lambda ).$$

$${\partial }_t⟨{\hat{\mathbf{q}}}_j⟩=\int \text{d}\lambda {\partial }_t\alpha {\partial }_{\alpha }{q}_j(\lambda )\rho (\lambda )+q_j(\lambda ){\partial }_t\rho (\lambda ).$$

$$⟨{\hat{\mathbf{j}}}_j(x,\alpha (t,x)){⟩}_{\text{GGE}(t,x)}=\int \text{d}\lambda q_j(\lambda ,\alpha (t,x))v^{\text{eff}}(t,x,\lambda )\rho (t,x,\lambda ).$$

$${\partial }_x⟨{\hat{\mathbf{j}}}_j(x,\alpha (t,x)){⟩}_{\text{GGE}(t,x)}=\int \text{d}\lambda {\partial }_x\alpha {\partial }_{\alpha }{q}_j(\lambda )v^{\text{eff}}(\lambda )\rho (\lambda )+q_j(\lambda ){\partial }_x\left(v^{\text{eff}}(\lambda )\rho (\lambda )\right).$$

$$⟨{\partial }_{\alpha }{\hat{\mathbf{q}}}_j(x,\alpha )⟩=\int \text{d}\lambda {\partial }_{\alpha }{q}_j(\lambda ,\alpha ){\rho }_j(\lambda )+\frac{1}{2\pi }{f}^{\text{dr}}(\lambda ){\partial }_{\lambda }{q}_j(\lambda ,\alpha )\vartheta (\lambda ),$$

$$⟨{\partial }_{\alpha }{\hat{\mathbf{q}}}_j⟩=\int \text{d}\lambda {\partial }_{\alpha }{q}_j(\lambda ){\rho }_j(\lambda )-\frac{1}{2\pi }{q}_j(\lambda ){\partial }_{\lambda }\left(f^{\text{dr}}(\lambda )\vartheta (\lambda )\right).$$

$$\int \text{d}\lambda q_j(\lambda )\left[{\partial }_t\rho (\lambda )+{\partial }_x\left(v^{\text{eff}}(\lambda )\rho (\lambda )\right)+\frac{{\partial }_t\alpha }{2\pi }{\partial }_{\lambda }\left(f^{\text{dr}}(\lambda )\vartheta (\lambda )\right)\right]+{\partial }_x\alpha \left[-⟨{\hat{\Phi }}_j⟩+\int \text{d}\lambda {\partial }_{\alpha }{q}_j(\lambda )v^{\text{eff}}(\lambda )\rho (\lambda )\right]=0.$$

$${\partial }_t\rho +{\partial }_x\left(v^{\text{eff}}\rho \right)+{\partial }_t\alpha \frac{1}{2\pi }{\partial }_{\lambda }\left(f^{\text{dr}}\vartheta \right)+{\partial }_x\alpha \chi =0.$$

$$\vartheta (\lambda )=2\pi \frac{\rho (\lambda )}{({\partial }_{\lambda }p{)}^{\text{dr}}}.$$

$${\partial }_t\left[\frac{({\partial }_{\lambda }p{)}^{\text{dr}}}{2\pi }\right]=\frac{{\partial }_t\alpha }{2\pi }{\partial }_{\lambda }\left(p(\lambda )-\int \text{d}\mu {\partial }_{\alpha }\Theta (\lambda -\mu )\rho (\mu )\right)-\int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }\Theta (\lambda -\mu ){\partial }_t\rho (\mu ).$$

$$\begin{gathered} \int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }\Theta (\lambda -\mu ){\partial }_x(v^{\text{eff}}(\mu )\rho (\mu ))={\partial }_x\left[\int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }\Theta (\lambda -\mu )v^{\text{eff}}(\mu )\rho (\mu )\right]-{\partial }_x\alpha {\partial }_{\lambda }\left[\int \frac{\text{d}\mu }{2\pi }{\partial }_{\alpha }\Theta (\lambda -\mu )v^{\text{eff}}(\mu )\rho (\mu )\right]= \\ =-\frac{1}{2\pi }{\partial }_x\left[({\partial }_{\lambda }ϵ{)}^{\text{dr}}-({\partial }_{\lambda }ϵ)\right]-{\partial }_x\alpha {\partial }_{\lambda }\left[\int \frac{\text{d}\mu }{2\pi }{\partial }_{\alpha }\Theta (\lambda -\mu )v^{\text{eff}}(\mu )\rho (\mu )\right] \end{gathered}$$

$$\begin{gathered} \int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }\Theta (\lambda -\mu ){\partial }_{\mu }\left(f^{\text{dr}}(\mu )\vartheta (\mu )\right)=-\int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }{\partial }_{\mu }\Theta (\lambda -\mu )f^{\text{dr}}(\mu )\vartheta (\mu )= \\ ={\partial }_{\lambda }\left[\int \frac{\text{d}\mu }{2\pi }{\partial }_{\lambda }\Theta (\lambda -\mu )f^{\text{dr}}(\mu )\vartheta (\mu )\right]=-{\partial }_{\lambda }\left[f^{\text{dr}}(\lambda )-f(\lambda )\right]. \end{gathered}$$

$$\begin{gathered} {\partial }_t\left[\frac{({\partial }_{\lambda }p{)}^{\text{dr}}}{2\pi }\right]=-\frac{{\partial }_t\alpha }{2\pi }{\partial }_{\lambda }{f}^{\text{dr}}(\lambda )-\frac{1}{2\pi }{\partial }_x({\partial }_{\lambda }ϵ{)}^{\text{dr}} \\ -\frac{1}{2\pi }{\partial }_x\alpha \left\{-{\partial }_{\lambda }\left[{\partial }_{\alpha }ϵ(\lambda )-\int \frac{\text{d}\mu }{2\pi }{\partial }_{\alpha }\Theta (\lambda -\mu )\vartheta (\mu )({\partial }_{\mu }ϵ{)}^{\text{dr}}\right]-\int \text{d}\mu {\partial }_{\lambda }\Theta (\lambda -\mu )\chi (\mu )\right\}. \end{gathered}$$

$$\partial_{t}\vartheta(\lambda)+v^{\text{eff}}(\lambda)\partial_{x}\vartheta(\lambda)+\frac{\partial_{t}\alpha f^{\text{dr}}(\lambda)}{(\partial_{\lambda}p)^{\text{dr}}}\partial_{\lambda}\vartheta(\lambda)+\\ +\frac{\partial_{x}\alpha}{(\partial_{\lambda}p)^{\text{dr}}}\Bigg{\{}\vartheta(\lambda)\partial_{\lambda}\left[\partial_{\alpha}\epsilon(\lambda)-\int\frac{\text{d}\mu}{2\pi}\partial_{\alpha}\Theta(\lambda-\mu)\vartheta(\mu)(\partial_{\mu}\epsilon)^{\text{dr}}\right]+2\pi\chi(\lambda)+\vartheta(\lambda)\int\text{d}\mu\,\partial_{\lambda}\Theta(\lambda-\mu)\chi(\mu)\Bigg{\}}=0$$

$$\chi (\lambda )=\frac{1}{2\pi }{\partial }_{\lambda }({\Lambda }^{\text{dr}}(\lambda )\vartheta (\lambda ))\text{ansatz},$$

$$\begin{gathered} {\partial }_t\vartheta (\lambda )+v^{\text{eff}}(\lambda ){\partial }_x\vartheta (\lambda )+\frac{{\partial }_t\alpha f^{\text{dr}}(\lambda )+{\partial }_x\alpha {\Lambda }^{\text{dr}}(\lambda )}{({\partial }_{\lambda }p{)}^{\text{dr}}}{\partial }_{\lambda }\vartheta (\lambda )+ \\ +\frac{{\partial }_x\alpha \vartheta (\lambda )}{({\partial }_{\lambda }p{)}^{\text{dr}}}{\partial }_{\lambda }\left[\Lambda (\lambda )+{\partial }_{\alpha }ϵ(\lambda )-\int \frac{\text{d}\mu }{2\pi }{\partial }_{\alpha }\Theta (\lambda -\mu )\vartheta (\mu )({\partial }_{\mu }ϵ{)}^{\text{dr}}\right]=0. \end{gathered}$$

Full text

Introduction. —

Exploring the out-of-equilibrium behavior of quantum many-body systems is nowadays among the most active research areas in physics, due to a successful synergy between theoretical and experimental advances Bloch_rev ; Trot12 ; Lan13 ; Mein13 .

How, and in what sense, does a coarse-grained thermodynamic description emerge through dynamical evolution in isolated out-of-equilibrium many-body systems? One-dimensional systems represent an ideal playground to address this question: there, remarkably powerful tools exist, both theoretical (such as conformal field theory difrancesco and integrability Korepin ; smirnov ) and computational (such as Matrix Product States methods Sch05 ). Integrability is ubiquitous in the low-dimensional world (and experimentally realized Giam ), with applications ranging from spin chains Korepin to continuum models, the latter having Lorentz smirnov or Galilean LiLi63 ; LiLi63_bis invariance, or neither Po01 . Integrable models are characterized by the presence of infinitely many conserved charges $\hat{\mathcal{Q}}_{j}$, which can be used to exactly determine their thermodynamics taka . In recent times, the importance of quasi-local charges has moreover been underlined IlMePrZa16 .

The last decade has witnessed exact results reaching out-of-equilibrium protocols as well: great attention has been devoted to the homogeneous sudden quantum quench calabrese-cardy (see also Ref. special_issue and reference therein). Because of the conserved quantities, the system exhibits local relaxation to a state that is not thermal PozMeWe14 ; MePoz15 ; DeNaWo14 ; WoDeNa14 ; Poz14 ; PiPoVe17 ; Po18 ; PoPiVe18 ; SoTaMu15 , but rather emerges from a Quench Action CaEs13 ; Ca16 or (where applicable) a Generalized Gibbs Ensemble (GGE) Rigol07 ; Rigol09 which accounts for all the relevant charges.

More recently, the focus has been on quenches from spatially inhomogeneous systems. A new theoretical toolbox, dubbed Generalized Hydrodynamics (GHD) transportbertini ; hydrodoyon1 allows us to address this problem. In Refs. transportbertini ; hydrodoyon1 GHD dealt with inhomogeneous states evolving under a homogeneous Hamiltonian. Several applications have been explored GHD3 ; GHD6 ; GHD7 ; GHD8 ; GHD10 ; F17 ; DS ; ID117 ; DDKY17 ; DSY17 ; ID217 ; CDV17 ; Bas_Deluca_defhop ; Bas_Deluca_defising ; PeGa17 ; CDDK17 ; mazza2018 ; BePC18 ; Bas2018 ; Kormos2018 ; DoyonSphon17 ; Doyon17 ; BasDeLu18 ; BiCoRoDeLuMa19 ; BePiKo19 ; DeLuCoDeNa17 , for instance, including diffusive corrections DeBD18 ; vasseurdiff ; GoVa18 ; AgGoVa19 or applying it to classical models BDWY17 ; Sp19 ; Do19 ; BuCaMo19 ; CaoBuMo18 . Also, combined with the quasiparticle picture for integrable systems AlCal17 ; AlCal17a , GHD allows to describe the entanglement spreading after inhomogeneous quenches Al18 ; BeFAPiCa18 ; Al19 ; AlBeFa18 . Very recently, it has been shown that GHD provides the correct theoretical framework to describe atom-chip experiments SCBoDoDu19 .

When comparing with actual experiments, inhomogeneities, for instance due to external trapping potentials, should ideally be kept into account. Strictly speaking, inhomogeneities break integrability, but smooth variations can still be captured by invoking local relaxation to a (locally homogeneous) integrable model.

Inhomogeneities in the dynamics have already been studied with some limitations for either spatial GHD3 ; BuCaMo19 or temporal changes BasDeLu18 , opening the possibility, for example, of studying the famous Quantum Newton Cradle experiment exp2 through GHD CDDK17 .

However, the current state of the art cannot capture changes in the interparticle interactions, leaving perturbative methods EckKo10 ; MoKe10 or bosonization techniques PoKo11 ; DzTy11 ; ChSc16 ; DoHaZa11 ; BeCiKoOr14 ; PolGR08 ; BiDeLuViRoMAFA16 as the only methods to tackle this experimentally relevant situation.

In this Letter, we present a complete GHD approach that allows to treat the dynamics under integrable Hamiltonians with space-time inhomogeneous interactions. Our results significantly extend the current GHD framework exhausting all the possible inhomogeneities which can be considered on a pure hydrodynamic level, disclosing the full power of GHD in describing experimentally relevant protocols. We discuss the potential applications of our result to interaction changes in the Lieb-Liniger model LiLi63 ; LiLi63_bis (Fig. 1-2), which is of primary experimental interest. So far, the primary analytical tool used in dealing with time-dependent interactions has been the Luttinger Liquid approach Ha81 ; Ha81_bis , recently generalized to include spatial inhomogeneities DuStViCa17 ; DuStCa17 ; BrDu18 ; BrDu17 ; EiBa17 ; RuBrDu19 ; MuRuCa19 . In contrast with GHD, this method is nevertheless confined to the low-energy excitations. We numerically benchmark the GHD predictions both in the quantum and classical realms, considering the XXZ spin chain and the classical sinh-Gordon field theory, showing once again the wide applicability of our results. Furthermore, we improve the numerical method proposed in Ref. GHD7 to solve GHD equations, promoting it from a first order to a second order algorithm in the time step, providing a great stability enhancement.

Thermodynamics of integrable models. —

The Thermodynamic Bethe Ansatz (TBA) technique is nowadays a textbook topic taka : here we present the basic concepts for the sake of a self-contained exposition. The Hilbert space of integrable models can be understood in terms of multiparticle states $|\{\lambda\}_{i=1}^{N}\rangle$, labeled by suitable parameters $\lambda$ called rapidities Korepin ; smirnov . Quasiparticles undergo pairwise elastic scatterings, which are described by an interaction-dependent scattering matrix $S(\lambda)$. These states are common eigenstates of the full set of (quasi-)local charges. In the thermodynamic limit (TDL), we switch to a coarse-grained description through a rapidity (root) density $\rho(\lambda)$ taka , which gives the density of rapidities within the interval $(\lambda,\lambda+\text{d}\lambda)$. The root densities are in a one-to-one correspondence with the possible thermodynamic states of the system, such as GGEs IlQuNaBr16 or thermal states and fix the (extensive part of) the expectation value of the local charges

[TABLE]

together with any other local (in real space) property of the system, according to the Quench Action approach CaEs13 ; Ca16 ; CuPa19 . The function $q_{j}(\lambda)$ in Eq. (1) is called the charge eigenvalue.

Nontrivial interactions induce collective effects. For example, the group velocity of the quasiparticles, which is defined as $v(\lambda)=\partial_{\lambda}\epsilon/\partial_{\lambda}p$, with $\epsilon(\lambda)$, $p(\lambda)$ the energy and momentum eigenvalues respectively, is “dressed” as $v^{\text{eff}}(\lambda)=(\partial_{\lambda}\epsilon)^{\text{dr}}/(\partial_{\lambda}p)^{\text{dr}}$, with $\epsilon^{\text{dr}}$ and $p^{\text{dr}}$ the dressed quasiparticle energy and momentum. These are obtained by using that an arbitrary dressed quantity $\tau^{\text{dr}}(\lambda)$ is defined through the integral equation

[TABLE]

with $\tau(\lambda)$ the “bare” quantity. Here $\Theta(\lambda)=-i\log S(\lambda)$, with $S(\lambda)$ the two-body scattering matrix encoding the interaction, while $\vartheta(\lambda)=2\pi\rho(\lambda)/(\partial_{\lambda}p)^{\text{dr}}$ is the so-called filling function. We summarized the TBA considering a single particle species, but the construction is easily generalized to several species of excitations and bound states.

Emergent hydrodynamics with space-time inhomogeneous interactions. —

TBA describes homogeneous stationary states. Instead, we now consider smooth space-time inhomogeneities, both in the initial state and in the Hamiltonian. We imagine a family of integrable models parametrized by a coupling $\alpha$, with Hamiltonians

[TABLE]

Crucially, in Eq. (3) $\alpha$ is a function of both space and time. We consider models in the continuum for simplicity, but the same construction can be repeated on the lattice.

Spatial inhomogeneities of the initial state on the same typical length-scale of the variation of $\alpha$ are allowed. We are then interested in describing the system at the Eulerian scales $(\Delta t,\Delta x)\sim((\partial_{t}\alpha)^{-1},(\partial_{x}\alpha)^{-1})$, considering at the same time the limit of infinitely smooth variations $\partial_{t}\alpha\sim\partial_{x}\alpha\to 0$.

Closely following the same argument presented in Refs. transportbertini ; hydrodoyon1 , in this limit we can invoke local relaxation to an inhomogeneous GGE, associated with a weakly inhomogeneous root density $\rho(t,x,\lambda)$. We report the details of the derivation of the GHD equations in the Supplemental Material (SM) suppl . Here, we rather present the result, discussing its physical interpretation and validity regime, together with possible applications.

Our main result is that $\rho(t,x,\lambda)$ satisfies the following hydrodynamic equations as

[TABLE]

where we dropped the space-time dependence to lighten the notation. In Eq. (4) $v^{\text{eff}}$ is the dressed velocity of the quasiparticles. Only first-order derivatives appear, implying that the equation is invariant under the rescaling $(t,x)\to(At,Ax)$, with $A\in\mathbb{R}^{+}$. For a space-time homogeneous dynamics ($\partial_{x}\alpha=\partial_{t}\alpha=0$), the standard GHD equations are obtained transportbertini ; hydrodoyon1 . The forces $f$ and $\Lambda$ are obtained by solving

[TABLE]

[TABLE]

Here $\vartheta=2\pi\rho/(\partial_{\lambda}p)^{\text{dr}}$ is the filling function. As usual in GHD, Eq. (4) has a clear semiclassical interpretation: $\rho(t,x,\lambda)$ locally describes the phase-space density of a collection of quasiparticles, moving with velocity $v^{\text{eff}}$ and subjected to force terms induced by the inhomogeneities, which can change the quasiparticles’ rapidity. The force terms account for both single particle as well as collective effects. The former are contained in the terms $\partial_{\alpha}p$ and $\partial_{\alpha}\epsilon$ in Eqs. (5-6). These are the energy-momentum changes of a single excitation of rapidity $\lambda$ induced by the inhomogeneities: the change in the dispersion relation causes the excitation to accelerate. Force terms due to inhomogeneities have been previously derived in Ref. GHD3 , for spatially inhomogeneous potentials linearly coupled to the charge densities, which nevertheless cannot induce any inhomogeneity in the scattering data of the model. In Ref. BasDeLu18 , slow magnetic flux changes in the XXZ spin chain have been studied. In both cases, only single-particle effects arise in the GHD equation and can now be regarded as a particular case of our more general findings.

The integrals in Eqs. (5)-(6) are entirely due to collective behaviors and have never been derived in previous studies. Because of modifications in the interparticle interactions, encoded in the scattering phase $\Theta$, the excitations experience force fields caused by the surrounding particles.

For spatial-homogeneous interactions, i.e. $\partial_{x}\alpha=0$, we are able to derive Eq. (4) for rather generic integrable models suppl . In the presence of spatial inhomogeneities, thus $\partial_{x}\alpha\neq 0$, Eq. (4) is derived in the presence of Lorentz invariance suppl and in Galilean invariant models through a non relativistic limit KoMuTr10 ; KoMuTr09 ; BaMuDeLu16 ; BaMuDeLu17 (see SM suppl ). Outside of the mentioned cases, we present Eq. (6) as a conjecture, although well supported by numerical evidence (see Fig. 3). As a further nontrivial check, thermal states are shown to be steady states of the GHD equation (4) with $\partial_{x}\alpha\neq 0$ suppl .

We stress that, in order to have a weakly varying (locally) integrable model, a smooth dependence of $\hat{\bf h}(x,\alpha)$ (3) on the coupling does not suffice: the whole set of (quasi-)local charges must be smooth as a function of $\alpha$. For example, our method cannot be applied to interaction changes in the XXZ spin chain with $|\Delta|<1$, which has a fractal dependence on the coupling taka .

Applications and numerical checks. —

We now show the wide applicability of our results. GHD equations are numerically solved according with the method described in the SM suppl , where we also present a short summary of the TBA of the models here investigated. In Fig. 2 we show a possible application to an experimentally relevant setup, namely a (slow) interaction quench in the interacting Bose gas LiLi63 ; LiLi63_bis . We mention that there are no alternative analytical and numerical methods to address this type of protocols. Closely related setups have already been experimentally addressed FaCaJoBo16 .

The Hamiltonian of the Lieb-Liniger model reads $\hat{H}=\int\text{d}x\,\{\partial_{x}\hat{\psi}^{\dagger}\partial_{x}\hat{\psi}+c(t)(\hat{\psi}^{\dagger})^{2}(\hat{\psi})^{2}+V(x)\hat{\psi}^{\dagger}\hat{\psi}\}$, with $[\hat{\psi}(x),\psi^{\dagger}(y)]=\delta(x-y)$. The gas is loaded in a harmonic trap in a low-temperature state, the interaction $c(t)>0$ is then slowly increased. This induces a non trivial evolution of the quasiparticle densities, which are reported in Fig. 2 $(a.1)$. As the interparticle repulsion is increased quasiparticles increase their rapidity $\lambda$ (reflected in the stretching of the initial blob along the vertical direction) and escape from the center of the trap. The local total density $n(t,x)=\int\text{d}\lambda\,\rho(\lambda)$ of the quasiparticles is shown in Fig. 2 $(a.2)$. Interestingly, the quench induces a breathing mode, which is long-lived in harmonic potentials FaCaJoBo16 . This is clear from Fig. 2 $(a.3)$, where we show the density $n$ in the center of the trap as a function of time. In Fig. 2 $(b.1)\,(b.3)$ we focus on the slow quench in an anharmonic trap. As it is clear from Fig. 2 $(b.1)$ the anharmonicity causes a spiral motion in the filling which develops a fractal structure as times passes CDDK17 ; CaoBuMo18 , which is smoothed due to the discretization used to solve the GHD equation (4). Now a much faster relaxation is observed as compared with the harmonic case, due to dephasing CDDK17 ; BaCoSo17 .

In Fig. 3 we focus on the XXZ spin chain with Hamiltonian $\hat{H}=\sum_{j=-L}^{L}\{\hat{S}^{x}_{j}\hat{S}^{x}_{j+1}+\hat{S}^{y}_{j}\hat{S}^{y}_{j+1}+\Delta_{j}(t)\hat{S}_{j}^{z}\hat{S}^{z}_{j+1}+B_{j}\hat{S}^{z}_{j}\}$, where $\hat{S}_{j}^{\alpha}$ are standard spin$-1/2$ operators. The system is initialized in a confining magnetic field and in a low temperature thermal state, with a uniform interaction $\Delta_{j}>1$. Then, $\Delta_{j}$ is slowly changed with time in the form of a traveling wave (see Fig. 3). In Fig. 3 we compare the GHD predictions for the local magnetization and the local energy density with tDMRG simulations itensor ; greg ; metts , finding excellent agreement. tDMRG methods suffer strong limitations in the accessible time scale, therefore we revert to the classical world to explore longer time scales. In Fig. 4 we benchmark the GHD in the classical sinh-Gordon model (see Ref. DeLucaMuss16 and BDWY17 for the TBA). The model describes a scalar field $\phi$ with Lagrangian $\mathcal{L}=\int\text{d}x\,\{\partial_{\mu}\phi\partial^{\mu}\phi/2-(m/g)^{2}[\cosh(g\phi)-1]\}$, with $m$ the mass and $g$ the interaction parameter. We consider two protocols. In the first case (panels $(a.1)-(a.2)$) the system is initially prepared in a thermal state with an inhomogeneous temperature profile. The system is then evolved with the sinh-Gordon Hamiltonian with inhomogeneous coupling $g\to g(x)$. Such an inhomogeneity does not affect the single-particle dispersion law (see SM suppl ), providing an ideal benchmark to test the collective effects in Eq.(5-6). Symbols are Monte Carlo data BDWY17 , whereas the lines are the GHD results. In $(a.1)$ we show the expectation value of the vertex operator $e^{g\phi}$ as a function of position $x/L$ for different times, in $(a.2)$ we plot $e^{g\phi}$ at $x=L/2$ as a function of time. The agreement with the GHD is spectacular. The dashed line is the GHD result neglecting the collective effects, i.e. the integrals in the right hand side in Eq. (5), which clearly have a crucial role. In the second protocol (panel $(b)$), we start from a homogeneous thermal ensemble, then a mass inhomogeneity is slowly activated with a linear ramp in time up to $t=L/4$ and then kept constant. Mass changes do not affect the scattering data suppl , implying that only single-particle effects in Eq.(5-6) play a role. Despite the form of the force terms, mass inhomogeneities do not belong to the class of inhomogeneities described in Ref. GHD3 and thus provide a nontrivial benchmark of our findings.

Conclusions and outlook.—

The success of hydrodynamic approaches is hard to overestimate. GHD merges the hydrodynamic framework with integrability, providing unprecedented levels of accuracy in describing out-of-equilibrium systems. In this Letter we extended the reach of this program, providing hydrodynamic equations which account for arbitrary (smooth) inhomogeneities in the couplings and state. Several interesting questions are left out for the future. Our analysis holds true when the model has a smooth dependence on the inhomogeneous coupling, but there could be special points (or regions) where this hypothesis breaks down. Understanding the behavior of protocols overcoming such special points is surely a compelling quest, which can unveil a rich phenomenology (see Ref. BasDeLu18 for a closely related problem). Including higher order corrections in the derivative expansion at the root of Eq. (4) is another important direction. Finally, it is important to devise numerical schemes based on molecular dynamics, such as the flea gas GHD6 , to simulate the GHD equations (4).

Acknowledgments. —

A.B. is indebted with A. De Luca for useful discussions and work on related subjects. The authors acknowledge the support from the European Research Council under ERC Advanced Grant No. 743032 DYNAMINT.

A Derivation of the GHD equations

Following the original references Ref. S_transportbertini ; S_hydrodoyon1 , we assume local relaxation to a (weakly) inhomogeneous GGE, yet to be determined. The (local) GGE is unambiguously fixed by the expectation value of all the (quasi-) local charges. Let us consider the family of integrable models described by the parameter-dependent Hamiltonian $\hat{H}(\alpha)$ out of which we constructed the inhomogeneous Hamiltonian. Similarly, we consider the parameter-dependent (quasi-)local charges $\hat{\mathcal{Q}}_{j}(\alpha)=\int\text{d}x\,\hat{\bf q}_{j}(x,\alpha)$ of the homogeneous system, then construct

[TABLE]

The operator (S1) pointwise resembles a local charge, but it is not conserved anymore due to the inhomogeneity, which breaks integrability on a large scale. Nevertheless, the knowledge of $\langle\hat{\bf q}_{j}(x,\alpha(t,x))\rangle$ for any $j$ fixes the local GGE. Let us write the Heisenberg equation of motion for the local density after an infinitesimal evolution $t\to t+\text{d}t$. Let $\hat{\bf q}^{\text{H}}_{j}(t,x,\alpha(t,x))$ be the local charge density in the Heisenberg picture, its time variation receives a contribution from the Hamiltonian evolution and one from the parametric change

[TABLE]

Above, we added a label ”$t$” to the Hamiltonian to stress its explicit time dependence. Expanding at $O(\text{d}t)$ we find

[TABLE]

Above, $\hat{\text{\bf h}}$ is the local Hamiltonian density. We further manipulate the integral expanding $\alpha(t,y)\simeq\alpha(t,x)+(y-x)\partial_{x}\alpha(t,x)+...$ . Higher derivatives can be neglected in the limit of smooth variations.

[TABLE]

In the first term of the r.h.s., the parameter $\alpha$ is constant and we can compute the expression within the homogeneous case, resulting in the divergence of the proper current operator of the homogeneous model. Therefore, we get

[TABLE]

Above, we drop further orders in the derivative expansion of $\alpha$ (negligible at first order in the infinitely smooth limit) and defined

[TABLE]

From the Heisenberg equation of motion, we want now to move to expectation values and invoke local relaxation to the inhomogeneous GGE. To this aim, we approximate the expectation values of space-time derivatives of charges and currents with the derivatives of the expectation values on the inhomogeneous GGE S_transportbertini ; S_hydrodoyon1 , i.e.

[TABLE]

Enforcing this approximation, we are finally lead to the (infinite set of) equations

[TABLE]

For simplicity, we drop the explicit dependence of the various operators and the expectation values are meant to be taken over the inhomogeneous GGE at the point of interest. Enforcing these equations on the complete set of charges, we aim for an equation for the root density: in this perspective, we need the expectation value of the various operators. The charge expectation value is the simplest: again, we focus on a single type of excitation, but everything we say is readily generalized to multiparticle species.

[TABLE]

Above, $q_{j}$ is the charge eigenvalue in which we made explicit the dependence on the inhomogeneous coupling. Taking the time derivative we have

[TABLE]

Above, we suppress the explicit space-time dependence for the seek of a lighter notation. The expectation value of the current is less trivial and it has been only recently computed S_transportbertini ; S_hydrodoyon1 , making possible the first formulation of GHD (which lacks the terms $\propto\partial\alpha$ in Eq. (S8))

[TABLE]

Above, the effective velocity has a space-time dependence due both to the dressing and to the parametric dependence on the coupling $\alpha$. We write its spatial derivative as it follows (below, the $(t,x)$ dependence is neglected in the notation since no ambiguities arise)

[TABLE]

Computing the remaining terms is an open problem, which we managed to partially solve. Indeed, $\langle\partial_{\alpha}\hat{\bf q}_{j}\rangle$ can be exactly computed through a generalization of the Hellmann-Feynman theorem: we postpone the derivation to the end of this Section and, for the time being, just quote the result. On arbitrary GGEs we have

[TABLE]

where $\vartheta$ is the filling and the function $f$ is defined in the main text Eq. (5) ($f$ has a parametric dependence on $\alpha$ which we drop for the seek of simplicity). It is useful to perform an integration by parts: assuming that the boundary terms in the integral vanish (it is usually the case, see however Ref. S_BasDeLu18 ) we have

[TABLE]

The knowledge of $\langle\partial_{\alpha}\hat{\bf q}_{j}\rangle$ is enough to deal with those protocols where the dynamics is time-dependent, but homogeneous, i.e. $\partial_{x}\alpha=0$. Notice that we do not need to require the homogeneity of the state. However, we want to provide an answer for arbitrary inhomogeneities, thus $\partial_{x}\alpha\neq 0$.

Computing $\langle\hat{\Phi}_{j}\rangle$ is much more complicated and we did not succeed in providing a first-principle derivation. However, invoking some reasonable assumptions which we discuss later on, the natural ansatz for the GHD equation presented in the main text emerges. Lately, the ansatz can be proven in presence of Lorentz invariance.

Let us plug (S10-S12-S14) into (S8) leaving $\langle\hat{\Phi}_{j}\rangle$ implicit

[TABLE]

It is convenient to stop for a moment and consider $\partial_{x}\alpha=0$. In this case, following S_transportbertini ; S_hydrodoyon1 , we invoke the completeness of the charges and replace the infinite set of integral equations (holding true for any charge $\hat{\bf q}_{j}$) with a differential equation for $\rho$, obtained posing to [math] the term in Eq. (S15) proportional to $q_{j}(\lambda)$. The presence of the unknown term $\langle\hat{\Phi}_{j}\rangle$ prevents us from straightforwardly apply the same reasoning to the case $\partial_{x}\alpha\neq 0$. However, we assume the existence of a GHD equation for the root density which, compared with the case $\partial_{x}\alpha=0$, adds a yet unknown contribution

[TABLE]

Above, $\chi(t,x,\lambda)$ is due to the second term in (S15). Invoking the locality of the GHD equation, $\chi(t,x,\lambda)$ must be completely determined by the model at $(t,x)$, i.e. by $\rho(t,x)$ and $\alpha(t,x)$. It cannot contain space or time derivatives neither of the root density or of the coupling, since these terms would be next-to-leading order in the weakly-inhomogeneity approximation. The problem is now reduced to the determination of $\chi$. To this aim, we convert Eq. (S16) into an equation for the filling function $\vartheta$, namely

[TABLE]

The experience gained from the previous literature (see e.g. Ref. S_transportbertini ; S_hydrodoyon1 ; S_GHD3 ; S_BasDeLu18 ) teaches us that simple equations for the filling should be expected: this rewriting leads us to a very natural ansatz for $\chi$. In order to reach the desired equation, we start computing the time derivative of the dressed momentum derivative $\partial_{t}(\partial_{\lambda}p)^{\text{dr}}$. From the definition of the dressing we have (again, we suppress the explicit $(t,x)$ dependence if no ambiguities arise)

[TABLE]

The first term in round brackets is readily identified with $-f(\lambda)$, as defined in the main text Eq. (5). In the term where $\partial_{t}\rho$ appears, we take advantage of the hydrodynamic equation (S16). Furthermore, we use the following identities

[TABLE]

and

[TABLE]

Above, we integrated by parts assuming zero contribution from the boundary terms, use the symmetry of the kernel and finally the definition of the dressing. Collecting the various terms we can write

[TABLE]

Using now this last result and the definition of the filling Eq. (S17), we finally reach the following hydrodynamic equation

[TABLE]

This is how further we can go without any additional hypothesis on $\chi$ or symmetries of the system. Notice that the contribution proportional to $\partial_{t}\alpha$, passing from Eq. (S16) to Eq. (S22), retains a very simple form, while the term $\partial_{x}\alpha$ looks strangely complicated. Inspired by the $\propto\partial_{t}\alpha$ term, we make the following ansatz

[TABLE]

with $\Lambda$ a still unknown function, Eq. (S22) is then greatly simplified

[TABLE]

At this point, it is very tempting to assume the identification

[TABLE]

i.e. Eq. (6) of the main text. This immediately enforces the hydrodynamic equation for the filling

[TABLE]

which is equivalent to Eq. (4) of the main text. Apart from the appealing formal structure, nontrivial checks can be performed. In the main text we provided numerical benchmarks of our result in a variety of contexts, finding excellent agreement. Furthermore, we can explicitly check that thermal states in the local density approximation are steady states of the hydrodynamic equation, as it should be. This check is performed in the next short subsection.

Lastly, we provide a derivation of our ansatz in Lorentz invariant models, from which we can assess galilean invariant models through proper non relativistic limits.

Check: thermal states are steady states of the GHD equation

As long as we are interested in GGEs described by thermal states, their filling is best parametrized in terms of the effective energy $\varepsilon$ S_taka as it follows

[TABLE]

The effective energy satisfies the following integral equation

[TABLE]

Or, equivalently

[TABLE]

We now consider an inhomogeneous system which is in a thermal state with inverse temperature $\beta$: within the local density approximation, the local GGE is fixed as per above where the energy eigenvalue has a parametric dependence on the position. Of course, such a state must be a steady state for the GHD equation. Thus, we plug Eq. (S27) in Eq. (S26) assuming $\partial_{t}\alpha=0$ (but keeping $\partial_{x}\alpha\neq 0$) and imposing $\partial_{t}\vartheta=0$. We then reach the following equation for the effective energy

[TABLE]

Deriving the defining equation of the effective energy Eq. (S29) in the rapidities we readily get $\partial_{\lambda}\varepsilon=(\partial_{\lambda}\epsilon)^{\text{dr}}$. Instead, deriving with respect to the position we find $\partial_{x}\varepsilon=-\partial_{x}\alpha\Lambda^{\text{dr}}$. Thus, Eq. (S30) is satisfied.

Derivation of the ansatz in the relativistic invariant case

In addition to integrability we now assume the system to be relativistic invariant (we set the speed of light equal to unity, for simplicity). We start from the hydrodynamic equation in terms of the filling (S22), but no hypothesis on the $\chi$ functions are made. We rewrite (S22) in a more compact way, collecting into an unknown function $w(\lambda)$ the $\propto\partial_{x}\alpha$ term

[TABLE]

We now enforce relativistic invariance on the dispersion law, having $\epsilon(\lambda)=m\cosh\lambda$, $p(\lambda)=m\sinh\lambda$, with $m$ the mass of the fundamental excitation. Therefore, it holds true

[TABLE]

We now construct the contravariant momentum $P^{\mu}(\lambda)=(\epsilon(\lambda),p(\lambda))$, furthermore we collect in an unique two component vector the force terms $F^{\mu}=(f^{\text{dr}}(\lambda)\partial_{\lambda}\vartheta(\lambda),w(\lambda))$. The hydrodynamic equation can be then rewritten as (sum over repeated indexes)

[TABLE]

Since $\vartheta$ is a scalar under Lorentz boosts, $(P^{\mu})^{\text{dr}}$ inherits the same transformation properties of $P^{\mu}$. Therefore, $(P^{\mu})^{\text{dr}}\partial_{\mu}\vartheta$ is a Lorentz scalar. In order to complete the hydrodynamic equation to a Lorentz scalar, we are forced to require $F^{\mu}$ to be contravariant.

We can now use a Lorentz boost to fix $F^{1}$, using the knowledge of $F^{0}$ and the transformation properties under Lorentz boosts. Let us consider a boost of velocity $v$, then $F^{\mu}\to(F^{\mu})^{\prime}$, in particular the first component

[TABLE]

with $\gamma=1/\sqrt{1+v^{2}}$. Using the definition of $f(\lambda)$, the identities (S32) and the tranformation properties of energy and momentum, from the above equation we can read $F^{1}$, which turns out to be

[TABLE]

with

[TABLE]

i.e. Eq. (S25) specialized to the Lorentz-invariant case.

Expectation value of the derivative of charges

During the derivation of the GHD equations, we postponed the proof of Eq. (S13), i.e. $\langle\partial_{\alpha}\hat{\bf q}_{j}(x,\alpha)\rangle$ computed on an arbitrary GGE. We now provide the proof through a suitable generalization of the Hellmann-Feynman theorem. Firstly, we should take a step back from the thermodynamic limit and consider the system at finite size $L$, periodic boundary conditions (PBC) are assumed. Let us consider a state $|\{\lambda_{i}\}_{i=1}^{N}\rangle$: due to the PBC, the rapidities must satisfy the Bethe-Gaudin equations S_taka

[TABLE]

Above, $I_{i}$ are suitable integers. Of course, we are lastly interested in the thermodynamic limit $N,L\to\infty$ ($N/L$ constant). In view of the representative state approach S_CaEs13 ; S_Ca16 , in the thermodynamic limit we can equivalently compute $\langle\partial_{\alpha}\hat{\bf q}_{j}(x,\alpha)\rangle$ on a single state rather than the whole GGE ensemble, provided that the root density associated with the representative state equals the GGE root density.

Rather than labeling the state with the rapidities, we use the Bethe integers. Moreover, we take advantage of the homogeneity of the GGE and compute the derivative of the whole charge, rather than its density

[TABLE]

The expectation value $\langle\{I_{i}\}_{i=1}^{N}|\partial_{\alpha}\hat{\mathcal{Q}}_{j}(\alpha)|\{I_{i}\}_{i=1}^{N}\rangle$ can be computed using the fact that $|\{I_{i}\}_{i=1}^{N}\rangle$ is an eigenstate of the charge

[TABLE]

and generalizing the Hellman-Feynman theorem

[TABLE]

Above, the derivative is taken keeping the Bethe integers fixed. Since the norm of the state is constant, we get the identity $\partial_{\alpha}\left(\langle\{I_{i}\}_{i=1}^{N}|\hat{\mathcal{Q}}_{j}(\alpha)|\{I_{i}\}_{i=1}^{N}\rangle\right)=\langle\{I_{i}\}_{i=1}^{N}|\partial_{\alpha}\hat{\mathcal{Q}}_{j}(\alpha)|\{I_{i}\}_{i=1}^{N}\rangle$. Taking the derivative of the expectation value of the charge, we get two effects: one due to the parametric change of the charge eigenvalues, the other due to a rearrangement of the rapidities caused by a modification of the scattering phase shift in Eq. (S37)

[TABLE]

When the thermodynamic limit is enforced, the first term above simply becomes

[TABLE]

Instead, the second term requires extra manipulations. Indeed, deriving the Bethe Gaudin equations (S37) we get

[TABLE]

In the thermodynamic limit, the above equation becomes

[TABLE]

with

[TABLE]

i.e. Eq. (5) presented in the main text. Replacing the last finding into Eq. (S41) we finally get

[TABLE]

i.e. Eq. (S13), as we desired.

B Numerical solution of the GHD equation

This section is dedicated to numerical methods for solving the GHD equation. It is convenient to look at the equation in terms of the filling (S26). Interestingly, it admits the following implicit solution

[TABLE]

where

[TABLE]

Above, the effective velocity and the forces must be computed at the integration time, i.e. using the root density and the coupling at that time. Furthermore, they must be computed along the trajectories $(x(\tau,t),\lambda(\tau,t))$. Checking that Eq. (S47) satisfies the GHD equation is immediate, however the solution is only implicit, since Eq. (S48) depends on $\vartheta$ through the dressing operations. Nevertheless, the implicit solution is very useful in constructing numerical algorithms. We introduce a time step $\text{d}t$ and we are interested in updating the filling from $t$ to $t\to t+\text{d}t$, thus we write

[TABLE]

Eq. (S49) is in principle exact for any $\text{d}t$: errors are introduced only when we approximate $x(t+\text{d}t,t)$ and $\lambda(t+\text{d}t,t)$.

A first order method

A first order method is readily obtained with the crude approximation

[TABLE]

This provides a $O(\text{d}t)$ algorithm: indeed, at any update of the fillings an error $O(\text{d}t^{2})$ is introduced and, in order to reach a time $t$, we need $t/\text{d}t$ steps. Therefore, at time $t$ we accumulate an error $\sim t\text{d}t$. This method has already been proposed in Ref. S_GHD7 (albeit in absence of force terms). In Ref. S_GHD7 it has been observed that estimating the integrals appearing in $x(t+\text{d}t,t)$ $\lambda(t+\text{d}t,t)$ with the endpoints of the leap (instead of the starting one as per above) makes the algorithm more stable, which however remains first order in time. Here, we further improve the algorithm providing a $O(\text{d}t^{2})$ method.

A second order method

A better approximation for $x(t+\text{d}t,t)$ and $\lambda(t+\text{d}t,t)$ can be obtained taking the middle points in the integrals (S48) rather than the extrema. Therefore

[TABLE]

with $x^{\prime}=x(t+\text{d}t/2,t)$ and $\lambda^{\prime}=\lambda(t+\text{d}t/2,t)$. The exact expressions for $x^{\prime}$ and $\lambda^{\prime}$ are unknown, but we can estimate them at first order

[TABLE]

This approximation contributes with a $O(\text{d}t^{3})$ correction to (S52-S53). Using middle points in the time leap requires computing the fillings at times $n\text{d}t$ and $(n+1/2)\text{d}t$, where $n$ is an integer. Therefore, this method is two time slower than the first order one if the same time step $\text{d}t$ is used, but the global error grows as $\sim t\text{d}t^{2}$. So, in general cases this algorithm outclasses the previous one, since much larger time steps can be considered.

There is a subtlety in this algorithm that needs to be mentioned, i.e. we need the fillings at times $t=0$ and $t=\text{d}t/2$ as a starting points. While the filling at time $t=0$ is simply the initial condition and provided by the TBA solution, the filling at $t=\text{d}t/2$ is not. In order to determine it, we choose a second time step $\text{d}t^{\prime}\ll\text{d}t/2$ and approximate $\vartheta(t+\text{d}t^{\prime}/2,x,\lambda)$ according to the first order algorithm, then the filling is evolved with the second order algorithm up to time $t+\text{d}t/2$ using a time step $\text{d}t^{\prime}$. At this point, both $\vartheta(t,x,\lambda)$ and $\vartheta(t+\text{d}t/2,x,\lambda)$ are known and we can proceed with the second order algorithm with time step $\text{d}t$.

C The TBA of the models of interest

In this short Section, for the sake of completeness, we briefly review the TBA description of the models we looked at. For a more detailed presentation of the TBA method, the reader can refer to Ref. S_taka . For any model, we also shortly mention the details of the numerical methods used in making the plots presented in the main text.

The interacting Bose igas

The interacting Bose gas describes bosons with contact interaction and it is known to be integrable since a long time S_LiLi63 ; S_LiLi63_bis . Within the second quantization formalism, its Hamiltonian reads

[TABLE]

The fields $\hat{\psi}^{\dagger}(x)$,$\hat{\psi}(x)$ are bosonic creation and annihilation operators $\left[\hat{\psi}(x),\hat{\psi}^{\dagger}(y)\right]=\delta(x-y)$. The interaction strength is assumed to be positive $c>0$ and we explicitly introduced the chemical potential $\mu$ which, once it is made inhomogeneous, can describe external traps.

Within the repulsive regime, the model does not have bound states, therefore its TBA is formulated in terms of a single species of particle with bare energy and momentum given by

[TABLE]

The rapidity lives on the whole real line $\lambda\in(-\infty,\infty)$. The expectation value of energy and density of particles, which are the observables on which we focus on, are (the thermodynamic limit is always enforced)

[TABLE]

Analyzing the model by mean of coordinate Bethe Ansatz, the following scattering matrix can be derived

[TABLE]

Thermal states can be described according to Eq. (S27) and Eq. (S28).

Details of the numerical simulations

In Fig. 2 we numerically simulated an interaction quench for a trapped interacting Bose gas. The GHD equations are solved with the second order method presented in Section B using a time step $\text{d}t=0.025$. The instantaneous TBA equations are solved by discretizing the integrals using Gaussian quadratures, thus converting the linear integral equations into finite-dimensional vector-matrix equations. The rapidity space has a cut off $|\lambda|\leq 3$ and its discretized on a lattice of $100$ points. The spatial coordinates are taken within the interval $x\in[-3,3]$ and are discretized on a lattice of $100$ points, with constant lattice space. In order to check the precision of the solution, we monitored the conservation of the total number of particles which is constant with $\lesssim 0.5\%$ fluctuations over the explored time scales.

The XXZ spin chain

The XXZ spin chain is governed by the Hamiltonian

[TABLE]

Above, $\hat{S}_{j}^{x,y,z}$ are usual spin$-\frac{1}{2}$ operators. Differently from the Lieb Liniger model, the XXZ spin chain always supports bound states and thus the TBA requires multiple root densities. The thermodynamics is greatly affected by the value of $\Delta$, in particular the cases $|\Delta|<1$ and $|\Delta|\geq 1$ require a different discussion. For $|\Delta|<1$ the TBA has a fractal dependence on the value of $\Delta$ S_taka . For this reason, inhomogeneous space-time dependent $\Delta-$profiles within this phase lay outside of the applicability of our method, which requires a smooth dependence of the model on the coupling.

Instead, the $|\Delta|\geq 1$ case is not pathological: more specifically, we focus on $\Delta\geq 1$ and in the positive magnetization sector $B<0$ (which implies $\langle S^{z}_{j}\rangle>0$). The TBA description requires infinitely many root densities, usually called strings, $\{\rho_{j}(\lambda)\}_{j=1}^{\infty}$. Accounting for several strings in the TBA is straightforward.

In the $\Delta\geq 1$ case the rapidities are confined to a Brillouin zone $\lambda\in[-\pi/2,\pi/2]$. To each string are associated an energy $\epsilon_{j}(\lambda)$ and a momentum $p_{j}(\lambda)$ ($j\in\{1,2,3,...\}$)

[TABLE]

where the angle $\theta$ parametrizes the coupling $\Delta=\cosh\theta$. Among the possible relevant observables, the expectation values of the Hamiltonian and local magnetization are of outmost simplicity

[TABLE]

The scattering phase is promoted to be a matrix with indexes running over all the possible strings

[TABLE]

The dressing operation now keeps in account the presence of multiple strings, therefore a function $\tau_{j}(\lambda)$ is now dressed according to

[TABLE]

The thermal states are now described by the set of equations

[TABLE]

[TABLE]

Notice that for $B<0$ the fillings are exponentially vanishing while increasing the string index $j$, thus the infinite set of strings can be truncated only to the first ones, the quality of the approximation being decided by the magnetic field $B$ and the inverse temperature $\beta$. We mention that the ground state, i.e. $\beta\to\infty$, is such that $\vartheta_{j>2}(\lambda)=0$, thus we can use only the first string to describe it.

Details of the numerical simulations

In Fig. 3 we provide a benchmark of the GHD equations against tDMRG S_itensor simulations. For what it concerns the GHD simulations, with the parameters we choose (i.e. low temperature) we found that retaining only the first two strings gives a satisfactory precision. The rapidity space is discretized into $50$ points and integral equations are solved by means of Gauss quadratures. The position lives on an interval $[-1,1]$ which is discretized into $100$ equally spaced lattice points. The time evolution is solved according to the second order algorithm with time step $\text{d}t=0.0125$.

For the tDMRG simulations we employed the standard purification method S_Sch05 to represent the initial density matrix. The time evolution was implemented by using the MPO representation for the evolution operators $e^{-iHt}$. To mitigate the error associated with the time discretization we employed the scheme presented in Ref. S_greg , which allows one to obtaine an accuracy ${\mathcal{O}}(dt^{5})$. The application of the MPO evolution operator is implemented by using the fitting algorithm described in Ref.S_metts . In our simulations we used $dt=0.1$. The maximum bond dimension employed was $\chi\approx 500$.

The classical sinh-Gordon model

The sinh-Gordon model is a relativistic field theory of a scalar field $\phi$, governed by the Lagrangian

[TABLE]

The model is integrable both at classical and quantum level. The reader could be more familiar with the thermodynamics of the quantum system, but a proper semiclassical limit of the latter readily gives access to the GGE S_DeLucaMuss16 and GHD S_BDWY17 of the classical version. We leave to the original references the details and present here the relevant results.

The classical shG model can be described in terms of a single species of particle, having energy and momentum eigenvalues given by

[TABLE]

Notice that, in contrast with the quantum case, no renormalization of the mass occurs and the single particle eigenvalues are independent from the interaction $g$. The scattering phase $\Theta$ is singular and defined as

[TABLE]

where the limit $\gamma\to 0^{+}$ must be enforced after the integrations have been carried out. For example, the dressing operation is actually defined as

[TABLE]

The filling of thermal states is written in terms of the effective energy $\varepsilon(\lambda)$ as

[TABLE]

which satisfies the following integral equation

[TABLE]

The expectation value of the energy is UV divergent on thermal states, similarly to what it happens in the famous black-body catastrophe. For this reason, we revert to other local operators with well-defined UV properties, namely the vertex operators $e^{kg\phi}$. Their expectation values on arbitrary GGEs are recursively fixed by the following set of integral equations S_BDWY17

[TABLE]

where

[TABLE]

Above $\mathcal{P}$ stands for the principal value regularization of the singular integral. Eq. (S73) allows for a recursive determination of $\langle e^{kg\phi}\rangle$ for $k=1,2,...$ using the fact that, for $k=0$, the vertex operator becomes the identity \langle e^{kg\phi}\rangle\Big{|}_{k=0}=\langle 1\rangle=1.

Details of the numerical simulations

In Fig. 4 we compare the GHD predictions against Monte Carlo simulations. The GHD is solved with the second order algorithm of Section B with time step $\text{d}t=0.025$. The singular nature of the integral equations requires a careful discretization whose details can be found in Ref. S_BDWY17 : for our purposes, we restricted the rapidities on a finite interval $[-10,10]$ which is discretized into $200$ equispaced lattice points. The spatial direction is restricted on the interval $[-1,1]$ which is discretized on a lattice of $200$ equispaced points.

The shG model is directly simulated through Metropolis-Hasting techniques presented in Ref. S_BDWY17 : the interval $[-L,L]$, together with the temporal direction, is discretized on a tilted squared lattice (lattice space $a=0.025$ and length $L=30$, i.e. $600$ points in the spatial direction). The initial configurations are sampled from a (inhomogeneous) thermal ensemble generated through a Metropolis-Hasting algorithm. Subsequently, each initial configuration is then deterministically evolved in time: observables are then averaged on the initial conditions. We took roughly $3.5\times 10^{5}$ realizations.

D The interacting Bose gas as non relativistic limit of the sinh-Gordon model

The derivation of the GHD equations can be transferred, without the need of any ansatz, from the relativistic world to the non relativistic one, through proper non relativistic limits. In this short section, for the sake of completeness, we mention how this operation can be performed on the interacting Bose gas, viewed as the non relativistic limit of the quantum sinh-Gordon model. We have already briefly presented the interacting Bose gas and the classical sinh-Gordon model in the previous section, here we must now revert to the quantum sinh-Gordon model explicitly restoring the speed of light $c_{\text{light}}$, which will be then send to infinity. The interacting Bose gas has been identified as the NR limit of the shG model in Ref. S_KoMuTr10 ; S_KoMuTr09 , then the same approach has been extended to a larger class of models in Ref. S_BaMuDeLu16 ; S_BaMuDeLu17 . Here, we leave to the original references a careful treatment of the limit, presenting only the most important steps. The Lagrangian of the model is

[TABLE]

Above, we restored the speed of light and redefine the interaction $g$ in such a way to make a direct contact with the interacting Bose Gas. The identification is achieved in the limit $c_{\text{light}}\to\infty$ that, in terms of the rescaled interaction in the shG model, corresponds also to the weakly interacting limit. The identification is readily seen at the level of scattering matrix. Indeed, the shG quantum model possesses an unique species of excitation with scattering matrix

[TABLE]

where the parameter $\alpha$ is

[TABLE]

Since the relativistic momentum is $p(\theta)=c_{\text{light}}M\sinh\theta$ (with $M$ the renormalized mass), posing $\theta\sim\lambda/(mc_{\text{light}})$ (where we use that in the weakly interacting limit the renormalized mass tends to the bare one and that the rapidity in the interacting Bose gas is simply the momentum) we readily find that the shG scattering matrix collapses to the interacting Bose gas’s one

[TABLE]

with the r.h.s. being the scattering matrix of the interacting Bose gas Eq. (S59). The mapping can be extended to the whole Thermodynamic Bethe Ansatz S_KoMuTr10 ; S_KoMuTr09 , where it has been understood that the filling of the shG model simply becomes that of the interacting Bose gas, once the NR limit has been taken

[TABLE]

This makes very easy to consider the non relativistic limit of dressed quantities. For example, let us consider a test function $\tau_{\text{shG}}(\theta)$ in the shG model and the relative dressing

[TABLE]

Replacing $\theta=\lambda/(mc_{\text{light}})$ and using Eq. (S79), we simply have, for example, for the derivative of the momentum

[TABLE]

where on the l.h.s we consider the momentum on the shG model dressed according to the shG TBA, while on the r.h.s we have the interacting Bose gas momentum and the dressing operation according to that model. It is now a matter of a simple exercise to extend the above identity to all the terms appearing in the GHD equation, obtaining the hydrodynamic of the interacting Bose gas as the non relativistic limit of that of the quantum sinh-Gordon model.