Large-scale description of interacting one-dimensional Bose gases: generalized hydrodynamics supersedes conventional hydrodynamics

TL;DR

This paper demonstrates that generalized hydrodynamics (GHD) provides a more accurate and comprehensive description of one-dimensional Bose gases than conventional hydrodynamics (CHD), especially in the presence of sharp profiles and at nonzero temperatures.

Contribution

The authors establish that GHD supersedes CHD for 1D Bose gases, justifies GHD from hydrodynamic principles, and numerically confirm its accuracy across the full interaction range.

Findings

GHD reduces to CHD at zero temperature without shocks.

Sharp profiles in CHD dissolve into GHD's hierarchy, avoiding shocks.

Numerical simulations confirm GHD's accuracy at zero temperature.

Abstract

The theory of generalized hydrodynamics (GHD) was recently developed as a new tool for the study of inhomogeneous time evolution in many-body interacting systems with infinitely many conserved charges. In this letter, we show that it supersedes the widely used conventional hydrodynamics (CHD) of one-dimensional Bose gases. We illustrate this by studying "nonlinear sound waves" emanating from initial density accumulations in the Lieb-Liniger model. We show that, at zero temperature and in the absence of shocks, GHD reduces to CHD, thus for the first time justifying its use from purely hydrodynamic principles. We show that sharp profiles, which appear in finite times in CHD, immediately dissolve into a higher hierarchy of reductions of GHD, with no sustained shock. CHD thereon fails to capture the correct hydrodynamics. We establish the correct hydrodynamic equations, which are finite-dimensional reductions of GHD characterized by multiple, disjoint Fermi seas. We further verify that at nonzero temperature, CHD fails at all nonzero times. Finally, we numerically confirm the emergence of hydrodynamics at zero temperature by comparing its predictions with a full quantum simulation performed using the NRG-TSA-ABACUS algorithm. The analysis is performed in the full interaction range, and is not restricted to either weak- or strong-repulsion regimes.

Equations

$$H=\int \mathrm{d}x\left(\frac{1}{2m}{\partial }_x{\psi }^{†}{\partial }_x\psi +\frac{c}{2}{\psi }^{†}{\psi }^{†}\psi \psi \right),c>0$$

$${\partial }_t{\mathtt{q}}_i+{\partial }_x{\mathtt{j}}_i=0.$$

$$\phi (\theta )=-\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}\theta }\log S(\theta )=2c/({\theta }^2+c^2).$$

$${\mathtt{q}}_i=m\int \frac{\mathrm{d}\theta }{2\pi }n(\theta )h_i^{\mathrm{dr}}(\theta ),{\mathtt{j}}_i=m\int \frac{\mathrm{d}\theta \theta }{2\pi }n(\theta )h_i^{\mathrm{dr}}(\theta )$$

$$f^{\mathrm{dr}}(\theta )=f(\theta )+\int \frac{\mathrm{d}\alpha }{2\pi }\phi (\theta -\alpha )n(\alpha )f^{\mathrm{dr}}(\alpha ).$$

$${\partial }_{t}n(\theta )+v^{\mathrm{eff}}(\theta ){\partial }_{x}n(\theta )=0,$$

$$v^{\mathrm{eff}}(\theta )=\frac{(E^{\prime }{)}^{\mathrm{dr}}(\theta )}{(p^{\prime }{)}^{\mathrm{dr}}(\theta )}=\frac{{\mathrm{id}}^{\mathrm{dr}}(\theta )}{1^{\mathrm{dr}}(\theta )}$$

$$H_V=\int \mathrm{d}x\left(\frac{1}{2m}{\partial }_x{\psi }^{†}{\partial }_x\psi +\frac{c}{2}{\psi }^{†}{\psi }^{†}\psi \psi +V(x){\psi }^{†}\psi \right).$$

$${\partial }_t{\theta }_j^{\pm }+v_{\{\theta \}}^{\mathrm{eff}}({\theta }_j^{\pm }){\partial }_x{\theta }_j^{\pm }=0.$$

$${\partial }_t\rho +{\partial }_x(v\rho )=0,{\partial }_{t}v+v{\partial }_{x}v=-\frac{1}{m\rho }{\partial }_x\mathcal{P}$$

$$v_{\mathrm{F}}(\theta )=v_{\{-\theta ,\theta \}}^{\mathrm{eff}}(\theta )$$

$$\partial_{t}\theta^{\pm}+\Big{(}(\theta^{+}+\theta^{-})/2\pm v_{\rm F}\big{(}(\theta^{+}-\theta^{-})/2\big{)}\Big{)}\partial_{x}\theta^{\pm}=0.$$

$$\mathcal{P}(\rho )=m^2{\int }_{-{\Lambda }_{\rho }}^{{\Lambda }_{\rho }}\frac{\mathrm{d}\theta \theta }{2\pi }U(\theta )=-m^2{\int }_{{\Lambda }_{\rho }}^{{\Lambda }_{\rho }}\frac{\mathrm{d}\theta }{2\pi }\mathcal{E}(\theta )$$

$$U(\theta )=\theta +{\int }_{-{\Lambda }_{\rho }}^{{\Lambda }_{\rho }}\frac{\mathrm{d}\alpha }{2\pi }\phi (\theta -\alpha )U(\alpha ),$$

$$\mathcal{E}(\theta )=\frac{{\theta }^2}{2}-\frac{\mu ({\Lambda }_{\rho })}{m}+{\int }_{-{\Lambda }_{\rho }}^{{\Lambda }_{\rho }}\frac{\mathrm{d}\alpha }{2\pi }\phi (\theta -\alpha )\mathcal{E}(\alpha ),$$

$$\rho =m{\int }_{-{\Lambda }_{\rho }}^{{\Lambda }_{\rho }}\frac{\mathrm{d}\theta }{2\pi }V(\theta )$$

$$V(\theta )=1+{\int }_{-{\Lambda }_{\rho }}^{{\Lambda }_{\rho }}\frac{\mathrm{d}\alpha }{2\pi }\phi (\theta -\alpha )V(\alpha ).$$

$${\partial }_i\rho (\Lambda )={\partial }_i\Lambda \frac{\partial \mu (\Lambda )}{\partial \Lambda }\frac{\partial \rho (\Lambda )}{\partial \mu (\Lambda )}={\partial }_i\Lambda m{v}_{\mathrm{F}}(\Lambda )\kappa (\Lambda ),$$

$${\partial }_i\rho (\Lambda )=\frac{{\partial }_i\Lambda m{V}^2(\Lambda )}{\pi }.$$

$${\partial }_i\mathcal{P}(\Lambda )={\partial }_i\Lambda m{v}_{\mathrm{F}}(\Lambda )\rho (\Lambda ).$$

$${\partial }_t\Lambda +\eta {\partial }_x\Lambda +v_{\mathrm{F}}(\Lambda ){\partial }_x\eta$$

$${\partial }_t\eta +\eta {\partial }_x\eta +v_{\mathrm{F}}(\Lambda ){\partial }_x\Lambda$$

$${\partial }_t\delta {\theta }^{\pm }(x,t)\pm v_{\mathrm{F}}({\theta }_{\mathrm{F}}){\partial }_x\delta {\theta }^{\pm }(x,t)=0.$$

$$\begin{array}[]{rcl}\partial_{t}n+\partial_{x}(nu)&=&0\\ (\partial_{t}+u\partial_{x})u&=&-\frac{1}{\rho}\partial_{x}P\\ (\partial_{t}+u\partial_{x})\tau+\frac{2}{3}(\partial_{x}u)\tau&=&0\end{array}$$

$${\partial }_t{\mathtt{q}}_i+{\partial }_x{\mathtt{j}}_i=0,i=0,1,2,$$

$$\frac{\epsilon (\theta )}{T}$$

$$n(\theta )$$

$$\frac{\partial {\mathtt{q}}_j}{\partial {\beta }_i}$$

$$\frac{\partial {\mathtt{j}}_j}{\partial {\beta }_i}$$

$$\sum _{j=0}^2\left[C_{ij}{\partial }_t{\beta }_j(x,t)+B_{ij}{\partial }_x{\beta }_j(x,t)\right]=0,$$

$$C_{ij}$$

$$B_{ij}$$

Full text

Introduction. Modern experiments with ultracold atoms confined in “cigar-shaped” traps coldatoms ; qnc or in atom chips chips provide real-world implementations of one-dimensional (1d) many-body systems giamarchi , and represent an important challenge for theoretical physics. Even though it is widely accepted that 1d clouds of bosonic atoms are described, at the microscopic scale, by the paradigmatic Lieb-Liniger (LL) model olshanii ; review_giamarchi , solving this model in experimentally relevant out-of-equilibrium inhomogeneous situations for more than a few dozens of atoms is a task that is out of the reach of modern theoretical methods, including state-of-the-art numerical ones.

Yet, it is a classic result of XXth century mathematical physics that, in its homogeneous, translation-invariant version at equilibrium, the LL model is exactly solvable by means of the Bethe ansatz LiebLiniger , and its equation of state can be calculated exactly yangyang . It then seems reasonable to use this equation of state as the basic input into a coarse-grained, hydrodynamic, approach, that is expected to be applicable as soon as typical lengths of variations of the local density are large enough as compared to inter-particle and scattering distances (the Euler scale) – in much the same way that classical hydrodynamics describes water waves. Such a ”conventional hydrodynamic” (CHD) approach – defined in Eqs. (10) below –, has been used extensively in the cold atoms literature over the past decade chdpapers ; damski ; pv14 , and has sometimes been viewed as a consequence of the Gross-Pitaevskii equation review_giamarchi ; pv14 in the regime of small interaction strength.

However, a key physical feature of the LL model is overlooked in CHD: the fact that it admits infinitely many conservation laws. Indeed, CHD focuses only on a few quantities, like the particle density, the momentum density or the energy density. But the LL model possesses infinitely more conserved quantities. Those are not just a mathematical curiosity: they can have dramatic physical consequences, as illustrated by the quantum Newton cradle experiment qnc : the crucial observation of undamped oscillations in this experiment is connected with the lack of conventional thermalization thermalization .

The full connection between generalized thermalization and many-body dynamics was only recently uncovered ghd ; bertini1 . The fundamental precepts of hydrodynamics – the emergence of local entropy maximization – were used in systems with an infinite number of conservation laws in order to form the theory of generalized hydrodynamics (GHD). It is a type of Euler-scale hydrodynamics, but with an infinite-dimensional space of fluid states accounting for the large state manifold accessible by generalized thermalization. In practice, GHD consists in an infinite set of coupled continuity equations (one for each conserved charge), that can, at least in principle, be worked out with numerical solvers for non-linear partial differential equations.

In this letter, we show that GHD supersedes CHD. For this purpose, we focus on far-from-equilibrium waves emanating from a density accumulation in the LL model. The density waves are a good illustration for our purposes, but the main results and mechanisms are general. They have been studied in several ways in the past decade in the free Fermi gas bumps_free_fermions , in the effective theory of the non-linear Luttinger liquid bumps_NLLL , and in the Calogero-Sutherland model bumps_Calogero , quantum Hall edges bumps_FQH , and the LL model damski using CHD. In particular, all these references – see also pv14 – pointed out that the applicability of CHD was limited by the appearance of shocks. In this Letter, we show that GHD is the proper hydrodynamic framework to go beyond the latter.

We demonstrate that, only at zero temperature and for finite evolution times does CHD coincide with GHD. CHD being a finite-component hydrodynamics, it inevitably leads to “gradient catastrophes” and shock propagations thereon. In contrast, we show that at zero entropy, GHD decomposes into a hierarchy of instantaneously invariant finite-dimensional subspaces, whose exact hydrodynamic equations we establish. These are described by a multitude of Fermi seas, the stability of which is a consequence of integrability. We show that shocks dissolve as the system leaves the CHD subspace into a higher-dimensional reduction of GHD. No shock propagates in this process, as instead smoothness is re-established. We note that an important practical consequence of the zero-entropy reduction is that the infinite system of coupled non-linear equations of GHD collapses to a finite number of equations that are computationally easy to solve, taking typically a few minutes on a laptop. We also numerically verify that at nonzero temperature, the GHD evolution, which necessitates the full infinite-dimensional space, is different from CHD at all times. In the density wave problem, a stark difference is that no sharp profile develops in GHD, while CHD based on the finite-temperature LL equations of state has gradient catastrophes. Finally, at zero temperature and using a local density approximation for the initial fluid state, we compare the hydrodynamic prediction for the space-time density profile with a simulation of the full quantum model obtained from the NRG-TSA-ABACUS algorithm abacus ; konik ; caux_konik , and find perfect agreement.

GHD. The Hamiltonian of the repulsive LL model is

[TABLE]

for the complex bosonic field $\psi(x)$, where $m$ is the mass (we set $\hbar=1$ throughout the manuscript). An inhomogeneous initial state $\langle\cdots\rangle$, to be specified below, is set to evolve unitarily with $H$.

Since the LL model is integrable, it admits infinitely many conservation laws $\partial_{t}q_{i}+\partial_{x}j_{i}=0$. This includes the gas density $q_{0}=\psi^{\dagger}\psi$, the momentum density $q_{1}=-{\rm i}\psi^{\dagger}\partial_{x}\psi+h.c.$, and the energy density $q_{2}$ (the integrand in (1)). According to the principles of hydrodynamics, if averages of conserved densities $\langle e^{{\rm i}Ht}q_{i}(x)e^{-{\rm i}Ht}\rangle$ and currents $\langle e^{{\rm i}Ht}j_{i}(x)e^{-{\rm i}Ht}\rangle$ have smooth enough space-time profiles, they can be described by space-time dependent local states that have maximized entropy with respect to the conserved charges afforded by the dynamics. Eulerian hydrodynamics, which neglects viscosity effects and is valid at large scales, is formed of the ensuing macroscopic conservation laws. In integrable systems, entropy maximization leads to generalized Gibbs ensembles (GGEs) thermalization ; GGEreviews with (formal) density matrix $\rho_{\rm GGE}=e^{-\sum_{i}\beta_{i}Q_{i}}$, $Q_{i}=\int q_{i}(x){\rm d}x$. Therefore, $\langle e^{{\rm i}Ht}{\cal O}(x)e^{-{\rm i}Ht}\rangle\approx{\rm tr}[\rho_{\rm GGE}(x,t){\cal O}]$, where the only space-time dependence is in $\rho_{\rm GGE}(x,t)$. The macroscopic conservation laws of generalized hydrodyamics (GHD) are the infinite number of equations for the density averages ${\tt q}_{i}(x,t)={\rm tr}[\rho_{\rm GGE}(x,t)q_{i}]$ and the current averages ${\tt j}_{i}(x,t)={\rm tr}[\rho_{\rm GGE}(x,t)j_{i}]$:

[TABLE]

The set of ${\tt q}_{i}$ fixes the GGE state, and thus can be seen as a set of fluid variables for GHD. In the manifold of GGE states, the currents ${\tt j}_{i}$ have a fixed functional form in terms of the densities ${\tt q}_{i}$: these are the equations of state, which fully determine the GHD model at hand.

An efficient treatment of hydrodynamics requires an appropriate choice of fluid variables. Instead of the ${\tt q}_{i}$, the most powerful fluid variables are obtained in terms of the emerging quasi-particles of the integrable model. In the repulsive LL model, there is a single quasi-particle species. The interaction in the LL model is fully described by the two-particle scattering matrix $S(\theta-\theta^{\prime})$, a function of velocity differences. The object of importance is the differential scattering phase LiebLiniger ,

[TABLE]

The quasi-particle can be seen as a spinless real fermion, which is free in the Tonks-Girardeau (TG) limit $c=\infty$ (hard-core repulsion).

States $|\theta_{1},\ldots,\theta_{N}\rangle$ are described by the velocities $\theta_{k}$ of the quasi-particles. Each conserved charge $Q_{i}$ is characterized by its one-particle eigenvalue $h_{i}(\theta)\propto\theta^{i}$, with $Q_{i}|\theta,\ldots,\theta_{N}\rangle=\sum_{k}h_{i}(\theta_{k})|\theta_{1},\ldots,\theta_{N}\rangle$. For instance, the particle number has eigenvalue $h_{0}(\theta)=1$, the momentum $h_{1}(\theta)=p(\theta)=m\theta$, and the energy $h_{2}(\theta)=E(\theta)=m\theta^{2}/2$. In the thermodynamic limit, the eigenstates are expressed in terms of $\rho_{\rm p}(x,\theta){\rm d}x{\rm d}\theta$, the number of quasi-particles in the phase-space region $[x,x+{\rm d}x]\times[m\theta,m(\theta+{\rm d}\theta)]$, leading to average densities ${\tt q}_{i}=\int{\rm d}\theta\,\rho_{\rm p}(\theta)h_{i}(\theta)$. The most convenient fluid variable is the occupation function $n(\theta)=\rho_{\rm p}(\theta)/\rho_{\rm s}(\theta)$, where $\rho_{\rm s}$ is the state density, $2\pi\rho_{\rm s}(\theta)=m+\int{\rm d}\alpha\,\varphi(\theta-\alpha)\rho_{\rm p}(\alpha)$. The density and current averages take the form ghd ,

[TABLE]

where the dressing operation is defined by

[TABLE]

These establish a relation between the ${\tt j}_{i}$’s and the ${\tt q}_{i}$’s, and thus the equations of state. For the LL model they were derived in ghd by extending the theory of (generalized) TBA caux_konik ; GTBA .

It was realized in ghd ; bertini1 that demanding the continuity equations (2) together with the averages (4) implies the continuity equation at the level of quasi-particles:

[TABLE]

where the effective velocity $v^{\rm eff}(\theta)$ is the velocity of elementary excitations vexcitations

[TABLE]

with ${\rm id}(\theta)=\theta$. These are the GHD equations in terms of quasi-particle fluid variables in the LL model. Since $v^{\rm eff}(\theta)$ depends on the fluid state through the function $n$, these are nonlinear equations for an infinity of functions of space-time (one for each velocity $\theta$).

Some intuition into Eqs. (6), (7) can be gained by looking at the TG limit $c=\infty$. In this case, $n(\theta)$ is the fermion occupation number at each momentum $m\theta$, at position $x$. This is the Wigner function of the state wigner (the partial Fourier transform of the fermion-fermion correlator). The effective velocity is equal to the particle velocity, $v^{\rm eff}(\theta)=\theta$, and (6) simply reproduces the exact evolution equation for the Wigner function, a direct consequence of the Schrödinger equation, as exploited in bumps_free_fermions ; Wigner_oscillations (see also the Supplementary Material (SM)). The quasi-particle occupation $n(\theta)$ may thus be viewed as the generalization of the Wigner function to non-free-fermion systems, with time-evolution governed by GHD (6)-(7).

Zero-entropy GHD. Natural initial conditions are ground states within inhomogeneous potentials $V(x)$,

[TABLE]

With a slowly varying potential, local averages are well described by a local-density approximation (LDA) review_giamarchi . LDA provides a GHD initial condition, a fluid of local zero-temperature states, which at every point $x$ is the ground state of $H+V(x)Q_{0}$. In this section, we observe that GHD equations give rise to finite-dimensional hydrodynamics when one restricts to the subspace of fluid states with zero entropy such as those. An analogous observation was made previously for free fermions bumps_free_fermions and for the Calogero-Sutherland model bumps_Calogero .

Recall that the occupation function at zero temperature is $n_{T=0}(\theta)=\chi(\theta\in[-\theta_{\rm F},\theta_{\rm F}])$ (where $\chi$ is the indicator function) where $\theta_{\rm F}$ is the Fermi pseudo-velocity, which depends on the chemical potential. Let us consider the space of zero-entropy occupation functions which have exactly $2k$ jumps, characterized by $2k$ velocities $\cdots<\theta_{j-1}^{+}<\theta_{j}^{-}<\theta_{j}^{+}<\theta_{j+1}^{-}<\cdots$ bounding separate Fermi seas: $n(\theta)=\chi(\theta\in\cup_{j=1}^{k}[\theta_{j}^{-},\theta_{j}^{+}]))$. We show that under GHD evolution, any smooth fluid whose state lies in such a space at all positions $x$, stays so for short enough times. Time evolution leads to displacements of Fermi points. Thus at zero entropy, GHD is reduced to hydrodynamics with a finite number of fluid variables.

Indeed we have $\partial_{x}n(\theta)=-\sum_{\epsilon=\pm}\sum_{j=1}^{k}\epsilon\partial_{x}\theta_{j}^{\epsilon}\,\delta(\theta-\theta_{j}^{\epsilon})$, and thus the time derivative $\partial_{t}n(\theta)$ is supported on the finite set of velocities $\theta_{j}^{\pm}$. A solution to (6) is therefore provided by setting $\theta^{\pm}_{j}=\theta^{\pm}_{j}(x,t)$ with

[TABLE]

We expect the solution to (6) in the space of smooth fluid space-time functions to be unique, based on such rigorous results in related classical gases dobrods . Thus it is given by solving (9) as long as no shock develops. Here, more explicitly, the effective velocity is $v^{\rm eff}_{\{\theta\}}(\alpha)={\rm id}^{\rm dr}_{\{\theta\}}(\alpha)/1^{\rm dr}_{\{\theta\}}(\alpha)$ with the dressing operation $f^{\rm dr}_{\{\theta\}}(\alpha)=f(\alpha)+\sum_{j=1}^{k}\int_{\theta_{j}^{-}}^{\theta_{j}^{+}}{\rm d}\gamma\,\varphi(\alpha-\gamma)f^{\rm dr}_{\{\theta\}}(\gamma)$. The resulting equations (9) will be referred to as $2k$-hydrodynamics ($2k$HD).

Eqs (9) are Euler-type hydrodynamic equations for a fluid with finitely-many components. Finite-component fluids are expected to develop shocks. Therefore, Eqs (9) are expected to hold only for finite times. However, contrary to true conventional finite-component fluids, where viscosity effects, present beyond the Euler scale, dominate and produce entropy at shocks, the presence of infinitely-many conservation laws forbids sustained entropy production in GHD. Any shock instantaneously dissolves into the higher-dimensional solution space of GHD. More precisely, $2k$HD solutions become multivalued at the time of the appearance of the shock, but here this multivaluedness is physically meaningful, representing a higher number of Fermi seas. Thus shocks in $2k$HD resolve by increasing the number of Fermi seas, passing to $(2k+2)$HD. We exemplify this in Fig. 1a, where after $t_{shock}=37$ we begin to simulate, effectively, 4HD equations. This has previously been observed in free fermion models bumps_free_fermions , thanks to an analysis based on the Wigner function; here it is generalized to the fully interacting LL model. We have also confirmed this shock dissolution mechanism of GHD in a nontrivial classical gas with the same hydrodynamic equations as those of the LL model, see the SM.

GHD and conventional hydrodynamics (CHD). Starting with a smooth fluid of local zero-temperature states, GHD reduces to 2HD, where every local fluid cell is the Galilean boost of a zero-temperature state. As a consequence, 2HD is in fact equivalent to the conventional hydrodynamics (CHD) of Galilean fluids,

[TABLE]

where $\rho={\tt q}_{0}$ is the fluid density and $\mathcal{P}$ is the pressure (for a proof of the equivalence, see the SM). The first equation is conservation of mass, the second, of momentum. The pressure $\mathcal{P}=\mathcal{P}(\rho)$ gives the equations of state of the fluid, and here equals the momentum current ${\tt j}_{1}$ in the zero-temperature state with density $\rho$. The explicit equations of state are obtained from (4) and (5) (see SM).

CHD has been used as an important tool in analyzing the dynamics of 1d Bose gases chdpapers ; damski ; pv14 . It has sometimes been presented as a consequence of the Gross-Pitaevskii equation review_giamarchi —itself valid only in the limit of small interaction strength $c$—, and it was never quite clear what exactly the range of validity of CHD was in the full interaction range of the LL model. Our analysis clearly shows that CHD is valid —in the sense that it coincides with GHD— only at zero temperature, and before the first shock. We conclude that, in any other situation, CHD is not applicable and leads to quantitatively wrong results. To illustrate this, a comparison of CHD at finite temperature and GHD is shown in Fig.1b; the initial state is the same in both cases (obtained from LDA at finite temperature), but one sees that the density profiles differ significantly at finite time; moreover, CHD has solutions up to a finite shock time, while GHD has no shocks and has solutions at arbitrarily long times.

Comparison with microscopic simulation of the LL model. We consider evolution from the ground state of (8) with a background chemical potential $\mu_{\infty}$ perturbed by a Gaussian, $V(x)=-\mu_{\infty}-Ue^{-ax^{2}}$. The initial density profile accumulates around $x=0$, and is asymptotically nonzero. Two procedures are compared: (1) the ground state is exactly evaluated using the NRG-TSA-ABACUS algorithm konik ; caux_konik (see the SM for details), and then evolved unitarily; and (2) the ground state is approximated using LDA, and this initial fluid state is evolved using 2HD (see the SM for a review of standard conditions for the hydrodynamic regime, which are fulfilled by the choice of parameters below). Fig. 2 provides the result for a choice of parameters corresponding to the local dimensionless coupling $\gamma(x)=mc/\rho(x)$ of the order of 1, thus the system is in an intermediate regime with nontrivial interactions being important. We observe that GHD is in excellent agreement with NRG-TSA-ABACUS numerics at almost all times except near the right and left boundaries at $t=72$, and provides a substantially better approximation than linear sound waves (see the SM). It can also be seen that, since $v^{\rm eff}_{\{\theta\}}(\theta^{+})$ is always greater than the background Fermi velocity corresponding to $\mu_{\infty}$, the propagation speed of 2HD is larger than that of the sound wave. It is remarkable that the complex (zero-temperature) dynamics of the LL gas is exactly described by 2HD, a simple set of differential equations.

At large times, discrepancies, though very tiny in the above graph, are expected to emerge. One reason is that, as explained above, shocks attempt to form, and as variations become larger, conditions for the hydrodynamic regime break down. Higher-derivative effects, such as viscosity, become more important. Recent observations in the related hard rod gas ds suggest however that such higher-derivative effects play only a small role. Another cause for late-time discrepancy is that LDA is not exact. As higher-order charges are more sensitive to the large-scale variations of the potential, LDA gives an extremely good approximation to the particle density, but describes poorly densities of higher-order charges $Q_{i}$. As time passes, the effects of the latter under the full GHD evolution eventually breaks 2HD. An analysis of the Wigner function $n(\theta)$ at the free-fermion point $c=\infty$, where GHD is exact (no viscosity is neglected), gives further insight (see the SM).

Conclusion. We showed that widely used conventional hydrodynamics of the Lieb-Liniger model correctly describes interacting Bose gases, but that this holds only at zero temperatures and for finite times. We provided exact, simple hydrodynamic equations valid beyond gradient catastrophes, where no shocks are sustained. These are zero-entropy reductions of GHD, which are finite-component fluid equations easily solvable on a laptop. This suggests that fluids of integrable models avoid entropy production thanks to the large space of fluid states. We provided compelling evidence for the emergence of GHD in the LL gas in the limit of slow variations of the density profile. This provides a crucial dynamical extension of LDA that is valid beyond previously existing frameworks. As a future direction, it would be very interesting to apply our method to more experimentally relevant situations such as “Quantum Newton’s cradle”-type protocol qnc , known to be beyond the reach of CHD pv14 .

Acknowledgements: We are extremely grateful to J.-S. Caux for access to the ABACUS software package without which the Lieb-Liniger simulations reported here would not have been possible. We thank B. Bertini, M. Fagotti, and H. Spohn for stimulating discussions, and the anonymous referees for useful comments. TY is grateful for the support from the Takenaka Scholarship Foundation and the hospitality at the Tokyo Institute of Technology. BD and JD thank the IESC Cargèse for hospitality. RMK’s research effort here was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract No. DE- AC02-98CH10886.

I Conventional hydrodynamics and sound wave

I.1 2HD = CHD

In the set of zero-entropy states, GHD boils down to 2HD, the hydrodynamics of two conserved quantities described by Eq. (9) in the main text, with $k=2$. Here we show that 2HD is equivalent to the conventional hydrodynamics (CHD) of Galilean fluids defined by Eqs. (10) in the main text.

First, Eq. (9) simplifies thanks to Galilean invariance. Let

[TABLE]

be the equilibrium Fermi velocity, with Fermi pseudo-velocity $\theta$. Let $\Lambda=(\theta^{+}-\theta^{-})/2$ and $\eta=(\theta^{+}+\theta^{-})/2$. Shifting the integration variable by $\eta$ in the dressing operation, we obtain f^{\rm dr}_{\{\theta^{-},\theta^{+}\}}=\big{(}f\circ\iota_{\eta}\big{)}^{\rm dr}_{\{-\Lambda,\Lambda\}}\circ\iota_{-\eta} where $\iota_{\eta}(\theta)=\theta+\eta$ is the shift function. Using ${\rm id}\circ\iota_{\eta}={\rm id}+\eta$ and $1\circ\iota_{\eta}=1$, we find $v^{\rm eff}_{\theta^{-},\theta^{+}}(\alpha)=v_{\rm F}(\alpha-\eta)+\eta$, and with $v_{\rm F}(-\alpha)=-v_{\rm F}(\alpha)$ it follows that $v^{\rm eff}_{\{\theta^{-},\theta^{+}\}}(\theta^{\pm})=\pm v_{\rm F}(\Lambda)+\eta$. This leads to

[TABLE]

That is, the hydrodynamics is completely determined by the functional form of the equilibrium Fermi velocity. Notice that the propagation velocity $\pm v_{\rm F}(\Lambda)+\eta$ of $\theta^{\pm}$ equals $\pm v_{\rm F}(\Lambda)+v$, where $v={\tt j}_{0}/{\tt q}_{0}={\tt q}_{1}/{\tt q}_{0}=\eta$ is the fluid velocity. This is a simple consequence of the Galilean velocity-addition formula. We note that this 2HD equation (12) was in fact already derived in bump_ll_A , although the authors did not provide the interpretation that $\theta^{\pm}$ are“dynamical” Fermi pseudo-velocities.

Next, we observe that the equation of state of CHD – the relation between the pressure and the density in Eq. (10) in the main text – is obtained as follows. The pressure is given in two equivalent ways

[TABLE]

where the function $U(\theta)$ and $\mathcal{E}(\theta)$ solve

[TABLE]

[TABLE]

where $\mathcal{E}(\theta)$ is defined so that it satisfies $\mathcal{E}(\Lambda_{\rho})=0$ from which $\mu(\Lambda_{\rho})$ is also fixed. The Fermi pseudo-velocity $\Lambda_{\rho}$ is determined by

[TABLE]

where the function $V(\theta)$ solves

[TABLE]

Observe that $U(\theta),V(\theta)/2\pi,$ and $\mathcal{E}(\theta)$ are nothing but the dressed momentum, the density of state, and the pseudo-energy, respectively, at $T=0$ with the Fermi pseudo-velocity $\Lambda_{\rho}$.

Using these expression, we can explicitly demonstrate that the CHD equations are recast into Eqs. (12). The derivative of the density $\rho(x,t)=\rho(\Lambda(x,t))$ with respect to $i$ ($i=x,t$) can be written in the following way

[TABLE]

where $\kappa(\Lambda):=\partial\rho(\Lambda)/\partial\mu(\Lambda)$ is the compressibility. Here we also used $\partial\mu(\Lambda)/\partial\Lambda=mv_{\rm F}(\Lambda)=mU(\Lambda)/V(\Lambda)$ where the first equality can be shown by differentiating $\mathcal{E}(\Lambda)=0$ with respect to $\Lambda$. Using an identity $\kappa(\Lambda)=V^{2}(\Lambda)/(\pi v_{\rm F}(\Lambda))$ KBI_A , we obtain

[TABLE]

The density itself also bears a similar expression $\rho(\Lambda)=mv_{{\rm F}}(\Lambda)\rho_{\rm s}^{2}(\Lambda)/\pi$ in terms of $v_{\rm F}(\Lambda)$ and $V(\Lambda)$ KBI_A . An analogous manipulation on derivatives of the pressure $\mathcal{P}(\Lambda)=\mathcal{P}(\rho(\Lambda))$ is done using a thermodynamic relation $\partial\mathcal{P}/\partial\mu=\rho$, which readily follows from (13), yielding

[TABLE]

Necessary ingredients are now ready; combining them and $v={\tt j}_{0}/{\tt q}_{0}=\eta$, the continuity equation and the Euler equation in (10) become, respectively,

[TABLE]

which are clearly equivalent to (12).

I.2 Sound wave

When the fluid of the LL gas is in the linear-response regime, we can describe the fluid as a sound wave: expanding $\theta^{\pm}(x,t)=\pm\theta_{\rm F}+\delta\theta^{\pm}(x,t)$, where $\theta_{\rm F}$ is the Fermi-pseudo velocity of the unperturbed LL gas, the 2HD equation (12) becomes

[TABLE]

This is the sound wave propagating with the Fermi velocity $v_{\rm F}(\theta_{\rm F})$, and describes the dynamics of the linear Luttinger liquid. Note that, in the LL gas, the Fermi velocity is same as the sound velocity $v_{\rm s}=\sqrt{\partial\mathcal{P}/\partial(m\rho)}$, which is readily confirmed.

I.3 Conventional hydrodynamics at finite temperature

Here we explain how to generalize the CHD Eqs. (10) in the main text to the finite temperature case, needed to produce the data in Fig. (1.b) in the main text, see also Fig. 3 below.

To do CHD at finite temperature, one needs to keep track of three conserved quantities: the number of particles, the momentum, and the energy. These correspond to the three densities ${\tt q}_{0}$, ${\tt q}_{1}$ and ${\tt q_{2}}$ respectively. In principle, one can write three continuity equations for these densities, which can be recast in the more conventional form of Euler hydrodynamics

[TABLE]

where $n$ is the particle density, $u$ is the local velocity of the fluid, $\tau$ is the kinetic energy, and $P$ is the pressure.

We find it more convenient, however, to work directly with the continuity equations

[TABLE]

and with the Lagrange multipliers $\beta_{0}$, $\beta_{1}$ and $\beta_{2}$, such that the reduced density matrix at a point $(x,t)$ is approximated by the Generalized Gibbs Ensemble with three charges, $e^{-\beta_{0}Q_{0}-\beta_{1}Q_{1}-\beta_{2}Q_{2}}$ where $Q_{0}=N$ is the number of particles, $Q_{1}=P$ is the momentum, $Q_{2}=H$ is the Hamiltonian, and $\beta_{2}$ is the inverse temperature while $-\beta_{0}/\beta_{2}$ is the chemical potential. The occupation function $n(\theta)$ in this Generalized Gibbs Ensemble is obtained by solving the (Generalized) Thermodynamic Bethe Ansatz equations,

[TABLE]

where $q_{0}(\theta)=1$, $q_{1}(\theta)=m\theta$, $q_{2}(\theta)=m\theta^{2}/2$. Then, using the results of Ref. doyon_spohn_A , one gets

[TABLE]

The CHD equations at finite temperature can then be written as

[TABLE]

where

[TABLE]

These are the equations we solve to produce the data for CHD at finite temperature in Fig. (1.b) in the main text and in Fig. 3.

II The free fermion case

Hydrodynamics of 1d free fermions has been studied thoroughly in the past decade bumps_free_fermions_A ; Wigner_oscillations_A , with the conclusion that the leading role is played by the Wigner function. For the convenience of the reader, we now discuss the aspects of the free fermion case that shed light on GHD, in relation with the discussion in the main text; for full details about the free fermion case, the reader is referred to the original references bumps_free_fermions_A ; Wigner_oscillations_A .

II.1 GHD for free fermions, and the role of the Wigner function

In the Tonks-Girardeau limit $c\rightarrow\infty$, the bosons are impenetrable. Impenetrable bosons can be mapped to non-interacting fermions through the Jordan-Wigner mapping,

[TABLE]

and a similar relation for the annihilation operator $\Psi_{\rm F}$, such that they satisfy the canonical anticommutation relations $\{\Psi^{\dagger}(x),\Psi(x^{\prime})\}=\delta(x-x^{\prime})$ and $\{\Psi(x),\Psi(x^{\prime})\}=\{\Psi^{\dagger}(x),\Psi^{\dagger}(x^{\prime})\}=0$. In terms of these operators, the Hamiltonian of the Lieb-Liniger model with $c\rightarrow+\infty$, or Tonks-Girardeau gas, is simply

[TABLE]

What is crucial here is that the hamiltonian is quadratic, so the particles created/annihilated by $\Psi^{\dagger}_{\rm F}$/$\Psi_{\rm F}$ propagate freely, and interact with other particles only through the Pauli principle—as they have fermionic statistics, unlike the original bosonic particles which entered the definition of the model—. There is no interaction energy, and this is of course a dramatic simplification. Thanks to translation invariance $H_{\rm TG}$ is diagonal in Fourier space,

[TABLE]

where $\Psi_{\rm F}^{\dagger}(k)\,:=\,\int e^{ikx}\,\Psi_{\rm F}^{\dagger}(x)\,dx$.

Because the Hamiltonian is quadratic and diagonal in $k$-space, we see that for any function $h(.)$, one can construct a conserved charge

[TABLE]

This provides a (continuous) family of commuting conserved charges,

[TABLE]

The associated charge density is

[TABLE]

where $h(y):=\int\frac{dk}{2\pi}e^{iky}h(k)$ is the Fourier transform of $h(k)$. For an analytic function $h(k)$, the Fourier transform decays exponentially, so the charge density $q[h](x)$ is local in the sense that it has exponentially decaying tails.

If we allow ourselves to take $h(k)=\delta_{k_{0}}(k)\equiv\delta(k-k_{0})$, then we obtain a charge density with tails that do not decay, but apart from that, the conserved charge $Q[\delta_{k_{0}}]=\frac{1}{2\pi}\Psi_{\rm F}(k_{0})\Psi^{\dagger}_{\rm F}(k_{0})$ is a perfectly decent conserved quantity. Of course, we see that it is nothing but the mode occupation at momentum $k_{0}$. The associated charge density reads [from now we drop the subscript ’[math]’ in ’$k_{0}$’]

[TABLE]

The charge density $q[\delta_{k}]$ obeys the following evolution equation, with $v(k)=k/m$,

[TABLE]

which follows from a straightforward calculation using the Hamiltonian (28). Thus, we can identify the current operator associated to the density $q[\delta_{k}](x)$ as

[TABLE]

such that the continuity equation $\partial_{t}q[\delta_{k}]+\partial_{x}j[\delta_{k}]=0$ holds. Notice also that, by linearity, this gives a simple closed formula for the current associated to any charge density $q[h]$,

[TABLE]

We stress that Eq. (32)-(33) is an exact formula that holds for the operators themselves; we have not taken any expectation values yet. It is also important to emphasize that this particularly simple form of the evolution equation is possible only because the problem is galilean invariant. If we were dealing with a lattice gas, the velocity would not simply be $v(k)=k/m$, and there would be additional terms generated by the bracket $[H,.]$ in Eq. (32).

In order to connect the simple exact result (33) for currents in the Tonks-Girardeau gas to the more general discussion of hydrodynamics in the main text, we need to turn to expectation values. Namely, we take expectation values $\left<.\right>$ in some initial state $\left|\psi_{0}\right>$ (or some initial density matrix), such that the exact relation (32)-(33) becomes

[TABLE]

In the language of the Thermodynamic Bethe Ansatz used in the main text, the quantity ${\rm q}[\delta_{k}]$ is nothing but the density of particles $\rho_{\rm p}$. In the main text, we chose to work with the occupation number $n(k)=\rho_{\rm p}(k)/\rho_{\rm s}(k)$, where $\rho_{\rm s}(k)$ is the state density. But in the Tonks-Girardeau limit, the latter is just a numerical constant, $\rho_{\rm s}(k)=1/(2\pi)$, so we end up with

[TABLE]

The function $n(x,k)$ is an extremely well-known object in quantum mechanics: it is the Wigner function. Roughly speaking, the Wigner function represents the probability of finding a particle at position $(x,k)$ in phase space. It is an extremely useful tool because it allows to get a simple (semi-classical) picture of a quantum particle, or, as we are doing here, of a large collection of non-interacting particles.

In conclusion, in the Tonks-Girardeau limit, the GHD equations are nothing but the well-known evolution equation for the Wigner function,

[TABLE]

which can be solved trivially in terms of the Wigner function of the initial state at $t=0$,

[TABLE]

We are looking at this well-known equation in a slightly non-standard way though: for each $k$ we view $n(x,k)$, or $\rho_{\rm p}(x,k)$, as the expectation value of a single charge density $q[\delta_{k}]$ at position $x$. In that sense, Eq. (36) is an infinite set of continuity equations, one for each mode $k$.

II.2 Initial Wigner function: deviation from LDA

This section is adapted from Bettelheim and Wiegmann Wigner_oscillations_A .

In the main text, we focus on an initial state that is the ground state of the Hamiltonian with an external potential

[TABLE]

It is the Wigner function of this ground state that enters the formalism of GHD. In the main text, it is obtained from LDA: in the limit where the variation of $V(x)$ is small at the interparticle distance, the system is locally translation invariant, so the Wigner function should be well approximated by the one of a Fermi sea with position-dependent Fermi momentum $k_{\rm F}(x)=\sqrt{-2mV(x)}$,

[TABLE]

where $\Theta(.)$ is the Heaviside function. Thus, the initial occupation number $n(x,k)$ obtained from LDA is indeed of the form discussed in the main text, $n_{\rm LDA}(x,k)=1$ if $k^{-}(x)<k<k^{+}(x)$, and $n_{\rm LDA}(x,k)=0$ otherwise, for $k^{+}=-k^{-}=k_{\rm F}$. The evolution from such an initial state is then governed by 2HD until a shock appears, which dissolves into 4HD, which itself dissolves into 6HD, and so on, as discussed in the main text.

An important question that needs to be addressed is the accuracy of the approximation $n(x,k)=n_{\rm LDA}(x,k)$ in the initial state. An estimate of the ground-state propagator in real-space is Wigner_oscillations_A

[TABLE]

Expanding the argument of the exponential, one gets

[TABLE]

From this expression, one can calculate the Wigner function, and get a better result than the crude approximation $n_{\rm LDA}(x,k)$. For notational simplicity, we focus on the behavior of $n(x,k)$ for $k$ close to $k_{\rm F}(x)$, which is determined only by the first term in the propagator (so we drop the c.c. term),

[TABLE]

where ${\rm Ai}(.)$ is the Airy function, and $C(x):=(k_{\rm F}^{\prime\prime}(x)/8)^{-1/3}$ is a new characteristic length scale. The convolution with the Airy function induces oscillations of $n(x,k)$ as a function of $k$, close to the Fermi point $k_{\rm F}(x)$, see Fig. 4.

This discussion shows that replacing the true initial Wigner function $n(x,k)$ by $n_{\rm LDA}(x,k)$ is a valid approximation only if the typical period of the oscillations is small compared to $k_{\rm F}(x)$, such that, when calculating average values of observables, integration over $k$ in the interval $[-k_{\rm F},k_{\rm F}]$ kills the oscillatory terms. We thus arrive at the criterion $1/C(x)\ll k_{\rm F}(x)$ for the validity of LDA, or in terms of the local density $\rho(x)$,

[TABLE]

The deviation of the true Wigner function $n(x,k)$ from the LDA result $n_{\rm LDA}(x,k)$ affects the expectation values of charge densities ${\tt q}[h]$. For the charges $q_{j}$ associated to $h(k)=k^{j}$, the direct calculation from Eq. (39), including the c.c. term, shows that

[TABLE]

For the first few $j$’s, this gives

[TABLE]

Interestingly, since $1/C(x)^{3}$ is a total derivative, when we integrate over $x$, we see that the expectation values of the total charges $\left<Q_{j}\right>=\int dx{\tt q}_{j}(x)$ remain unaffected by the oscillations of the Wigner function for $j\leq 3$. At the level of the expectation values of total charges, the deviation from the LDA result is seen only for $j\geq 4$.

III The hydrodynamic regime

The hydrodynamic regime. In classical systems, one assumes that local entropy maximization has occurred because of the fast motion of classical particles and small mean free path (separation of scales). Common wisdom then suggests that hydrodynamics emerges at large times. However, in quantum systems, the eigenstate thermalization hypothesis suggests that near-homogeneous, near-steady states are well approximated, from the viewpoint of local observables, by (generalized) Gibbs ensembles, and thus their dynamics by GHD. In effect, “molecular chaos” is already present within homogeneous and stationary states. One would therefore expect that time evolution from near-homogeneous, near-steady initial states is instantaneously well represented by GHD. It is possible to suggest, from physical intuition, quantitative criteria. Consider a state with density $\rho(x)$. The total variation, over an interval $d$, of the inter-particle length $l_{\rm int}(x)=1/\rho(x)$ is $\Delta(x,d)={\rm max}(|l_{\rm int}(x+y)-l_{\rm int}(x)|:y\in[0,d])$. On one hand, it is natural to require the relative total variation over $d$ to be small whenever $l_{\rm int}(x)$ is a significant proportion of $d$. That is, $\Delta(x,d)/l_{\rm int}(x)$ times $l_{\rm int}(x)/d$ is $\ll 1$; taking $d\to 0$, this is boils down to $|\partial_{x}l_{\rm int}(x)|\ll 1$, or equivalently $l_{\rm int}(x)\ll\rho(x)/|\partial_{x}\rho(x)|$. On the other hand, the scattering length $\varphi(\theta)/p^{\prime}(\theta)$ should also be small. It is velocity dependent, and an upper bound is $l_{\rm scat}=\max(\varphi(\theta)/p^{\prime}(\theta):\theta\in{\mathbb{R}})=2/(mc)$. Scattering will occur in a relatively homogeneous region if $\Delta(x,l_{\rm scat})/l_{\rm int}\ll 1$, which boils down to $l_{\rm scat}\ll\rho(x)/|\partial_{x}\rho(x)|$. Thus, sufficient conditions for the validity of GHD should be

[TABLE]

IV The NRG-TSA-ABACUS algorithm

In order to describe the behavior of the Bose gas pre- and post-quench, we must deal with two problems: one, first finding the ground state of the gas in a one-body Gaussian potential, and then two, determining its post-quench evolution in time. As we will see, it is only the first that is non-trivial.

To find the pre-quench ground state, we employ the truncated spectrum approach, an approach first introduced by V. Yurov and A. Zamolodchikov to study simple perturbations of the conformal minimal models TCSA_A . This approach in general enables the study of models of the form,

[TABLE]

where $H_{integrable}$ will for us be the Lieb-Liniger model and $V_{pert}$ the one-body potential,

[TABLE]

where $\hat{\rho}(x)$ is the density operator. The method employs the eigenstates of $H_{integrable}$ as a computational basis. The method presumes that we have complete control of this basis, knowing both the energies of the eigenstates as well the matrix elements of any relevant operator with respect to this basis. And in fact we do as the spectrum of the Lieb-Liniger as well as matrix elements of the density operator are computable. Because we are using the eigenstates of an interacting theory as a starting point for our numerics, we, at the very start, incorporate strong correlations into the problem. This, in principle, alleviates the numerical burden of solving the full Hamiltonian, at least in comparison of starting with a non-interacting basis. For a comprehensive review of truncated spectrum methods, see Ref. tsa_review_A .

The essential idea behind the truncated spectrum approach is to order the eigenstates of $H_{integrable}$ by importance as determined by some metric. In its simplest form, this metric is their unperturbed energy. This (infinite) list of states is then truncated, keeping the most (N, say) important states (as measured by the metric). This list of states, $\{|E_{i}\rangle\}^{N}_{i=1}$, where state, $|E_{i}\rangle$, has energy, $E_{i}$, is then used to form the Hamiltonian matrix,

[TABLE]

This finite dimensional matrix can then be easily diagonalized and both the spectrum and matrix elements of the full theory determined.

To compute the energies, $E_{i}$, as well as the matrix elements, $\langle E_{i}|V_{\rm pert}|E_{j}\rangle$, for our perturbed Lieb-Liniger model, we employ ABACUS abacus_A . ABACUS is a software package that enables the efficient computation of these quantities. These are non-trivial because the determination of energy of an eigenstate, $|E_{i}\rangle$, involving $N$ particles necessitates the solution of N-coupled non-linear equations. These are known as the Bethe ansatz equations:

[TABLE]

Here $p_{i}$ is the momentum associated with the $i-$th particle. To determine the matrix element $\langle E_{i},\{p_{i}k\}^{N}_{k=1}|V_{\rm pert}|E_{j},\{p_{j}k\}^{N}_{k=1}\rangle$, ABACUS uses a determinantal expression for this quantity first developed by N. Slavov slavnov_A .

The use of a simple truncation of the spectrum works exceedingly well for perturbations of simple conformal minimal models. However for more complicated models, there are strong truncation effects. To lessen these effects, a variety of strategies are available, both analytical tcsa_anal_rg_A and numerical konik_A . The approach we employ here is to adapt the numerical renormalization group nrg_A first developed by Kenneth Wilson for the study of quantum impurity problems. The basic idea is to perform a set of iterative exact diagonalizations where the most important states are taken into account in the first set of diagonalizations and successively less important states are taken into account in later diagonalizations. In this way the truncation level of the eigenspace of $H_{Lieb-Liniger}$ can be made much greater. Rather than employ a basis of $\sim 10^{4}$ states, one can operate with bases of sizes on the order of $\sim 10^{5}-10^{6}$.

The metric that determines which states in the computational basis are important need not be energy. In fact for the perturbed Lieb-Liniger model, the energy is not the most robust of metrics. We have instead developed an adaptive metric that uses the magnitude of matrix elements of $V_{pert}$ of an unperturbed state with the ground state of the full model (as computed using a small exact diagonalization) caux_konik_A ; scnt_A . This metric is much more efficient at determining which unperturbed states will be important for the full problem and enables the equivalent study of bases with $\sim 10^{6}-10^{7}$ states.

Using this approach we then arrive at the ground state of the gas in the Gaussian potential in the form

[TABLE]

We are confident that this representation is an accurate one for several reasons. For our purposes here we did not need to consider a very strong Gaussian potential, nominally because we wished to mimic the zero temperature hydrodynamic limit. In past work, i.e. Refs. caux_konik_A ; kam_A , we considered stronger one-body potentials with larger number of particles and were able there to establish the robustness of the technique. Furthermore we are solely interested in the computing the ground state of the gas in the one-body potential. If we were in a position to need excited states, we would have to be more cautious.

Having the ground state in hand, Eqn. 46, we can then turn to the second problem, determining it’s post-quench evolution in time. This in fact is trivial because the post-quench Hamiltonian is the Lieb-Liniger model itself. Thus the time evolved state, $|\psi(t)\rangle$, is simply

[TABLE]

that is, because our computational basis is equal to the post-quench eigenbasis, the coefficients of expansion simply pick up phase factors that depend on the energies, $\{E_{i}\}$. It is worthwhile to stress that because we can compute $E_{i}$ to arbitrary accuracy, we can then evolve the state to arbitrary times without incurring phase errors.

While it is more difficult, we can also consider quenches where we quench from one one-body potential to another non-trivial potential. In such a case, we must compute using our NRG-TSA-ABACUS algorithm a large number of excited states post-quench so that we can expand the pre-quench ground state into the post-quench eigenbasis. This is however a tractable problem as demonstrated in Ref. kam_A where quenches from parabolic to cosine potentials were studied. However it does place limits on the strength of the quench as well as the number of particles that can be treated.

V Additional graph

In the main text, we showed a graph for curves with $\gamma$ being of order 1, which is an intermediate regime. It is then natural to ask how curves would look like in two extreme regimes: the TG regime ($\gamma\gg 1$) and the Gross-Pitaevskii regime ($\gamma\ll 1$). It turns out, however, that the qualitative behavior basically does not differ whatever the coupling constant is. Namely, with any coupling constant, an initial density accumulation splits into two bumps and each bump moves towards right/left with a same speed. Nonetheless, it is instructive to see how a density wave evolves when a trap potential is not smooth enough.

For this purpose, we use a LL gas with $c=1$ as in Fig.2, and a potential $V(x)=-\mu_{\infty}-Ue^{-ax^{2}}$ with $U=0.05$, $a=1/64$, and $\mu_{\infty}$ is determined so that $N=24$. In Fig.5 we observe that the initial density prepared by LDA is already not agreeing with NRG-TSA-ABACUS numerics, that is, the initial state is not in the hydrodynamic regime defined by (41). This discrepancy is then amplified as time goes, and we also see that the height of the dots by NRG-TSA-ABACUS numerics is lowering. We expect that these effects are generated by the neglected higher conserved charges in preparing the LDA.

VI Shock dissolution of 2HD into 4HD: a comparison with the classical flea gas

The mechanism of shock-dissolution elaborated in the main text can be further supported by a comparison with a molecular dynamics of GHD (classical flea gas) flea_A . The flea gas algorithm was recently developed by making use of a surprising equivalence between dynamics of the hard rod-like model and GHD. A protocol with the same parameters as in Fig. 1 in the main text is compared against a simulation using the flea gas in Fig. 6, showing a perfect agreement. This confirms, in a nontrivial gas (although classical instead of quantum), that the shock dissolution mechanism proposed in the main text is correct.

VII Sketch of the algorithm for zero-entropy GHD

In this section, we briefly sketch the algorithm that we use to solve zero-entropy GHD. It is the one we used, for instance, to produce the data plotted in Fig. 1.a in the main text.

The idea is to work with a finite set of curves in the $(x,\theta)$ plane, each of these curves being a set of points where the occupation function $n(x,\theta)$ jumps from [math] to $1$. We call $\Gamma_{j}$, $j=1,\dots,p$ these different curves. For example, in Fig. 1.a, there are $p=2$ curves, and they are both plotted in red.

Then we view the zero-GHD entropy equation (Eq. (9) in the main text) as an evolution equation for this set of curves. Namely, if $(x_{j},\theta_{j})$ is a point on the curve $\Gamma_{j}$, then after a small time $\delta t$, its new position is $(x_{j}+\delta x_{j},\theta_{j})$ where

[TABLE]

The velocity $v^{\rm eff}_{\{\theta\}}(\theta_{j})$ is calculated from the set $\{\theta\}$ of all points that lie at the intersection between the vertical line passing through $(x_{j},\theta_{j})$ and all the curves $\Gamma_{j^{\prime}}$. It is obtained by solving the integral equation given in the main text, in the paragraph following Eq. (9).

This leads to the following algorithm. Numerically, each curve $\Gamma_{j}$ is encoded as the interpolation of some discrete set of points $(x_{j,a},\theta_{j,a})$ for $a=1,2,\dots$. The most basic version of the algorithm uses linear interpolation, but one can also use more refined interpolation schemes such as splines. To go from the configuration at time $t$ to the configuration at time $t+\delta t$, we do a loop over all the points labelled by $(j,a)$. For each point:

• •

we start by finding all the intersections between the vertical line at $x=x_{j,a}$ and all the curves $\Gamma_{j^{\prime}}$ for $j^{\prime}=1,\dots,p$. This gives us a set $\{\theta\}$.

• •

we order the elements of the set $\{\theta\}$ and label them as $\theta_{1}^{-}<\theta_{1}^{+}<\theta_{2}^{-}<\dots<\theta_{k}^{-}<\theta_{k}^{+}$.

• •

we calculate ${\rm id}^{\rm dr}_{\{\theta\}}(\theta_{j,a})$ and $1^{\rm dr}_{\{\theta\}}(\theta_{j,a})$. To do this, one needs to solve the linear integral problem $f^{\rm dr}_{\{\theta\}}(\alpha)=f(\alpha)+\sum_{q=1}^{k}\int_{\theta_{q}^{-}}^{\theta_{q}^{+}}d\gamma\varphi(\alpha-\gamma)f^{\rm dr}_{\{\theta\}}(\gamma)$. This is done by discretizing the integral, which leads to a matrix formulation of the problem, of the form $f^{\rm dr}=f+M\cdot f^{\rm dr}$ where $f$ and $f^{\rm dr}$ are finite-dimensional vectors.

• •

we evaluate the ratio $v^{\rm eff}_{\{\theta\}}(\theta_{j,a})={\rm id}^{\rm dr}_{\{\theta\}}(\theta_{j,a})/1^{\rm dr}_{\{\theta\}}(\theta_{j,a})$, and use it to calculate the new position of the point $(x_{j,a}+\delta x_{j,a},\theta_{j,a})$ after $\delta t$, according to Eq. (48).

We find that the algorithm performs better if one reparametrizes the curves $\Gamma_{j}$ from time to time. Namely, it is possible that the distance between successive points $(x_{j,a},\theta_{j,a})$ and $(x_{j,a+1},\theta_{j,a+1})$ increases during the evolution, so that the initial discrete set of points does not provide a good description of the curve $\Gamma_{j}$ any longer. To avoid that, it is useful to chose a new discrete set points on the curve $\Gamma_{j}$ every once in a while, and then carry on with the evolution.