Benchmarks of Generalized Hydrodynamics for 1D Bose Gases

Generalized hydrodynamics (GHD) is a recent theoretical approach that is becoming a go-to tool for characterizing out-of-equilibrium phenomena in integrable and near-integrable quantum many-body systems. Here, we benchmark its performance against an array of alternative theoretical methods, for an interacting one-dimensional Bose gas described by the Lieb-Liniger model. In particular, we study the evolution of both a localized density bump and dip, along with a quantum Newton's cradle setup, for various interaction strengths and initial equilibrium temperatures. We find that GHD generally performs very well at sufficiently high temperatures or strong interactions. For low temperatures and weak interactions, we highlight situations where GHD, while not capturing interference phenomena on short lengthscales, can describe a coarse-grained behaviour based on convolution averaging that mimics finite imaging resolution in ultracold atom experiments. In a quantum Newton's cradle setup based on a double-well to single-well trap quench, we find that GHD with diffusive corrections demonstrates excellent agreement with the predictions of a classical field approach.

Introduction.-Thestudy of dynamics of integrable and near-integrable quantum many-body systems has been a thriving area of research for more than a decade since the landmark experiments on relaxation in the quantum Newton's cradle setup [1] and in coherently split one-dimensional (1D) Bose gases [2].During this time, an in-depth understanding of the mechanisms of thermalization and emergent out-of-equilibrium phenomena within these systems has been developed [3][4][5][6][7][8].A recent breakthrough in this area has been the discovery of the theory of generalized hydrodynamics (GHD) [9,10] (for recent reviews, see [11][12][13]).This new theory is capable of simulating large-scale dynamics of integrable and nearintegrable systems across a significantly broader range of particle numbers and interaction strengths than those accessible using previous approaches [14][15][16].Because of its broad applicability, GHD is currently regarded as well on its way to becoming "a standard tool in the description of strongly interacting 1D quantum dynamics close to integrable points" [16].
In the years since its discovery, GHD has been rapidly developed to include diffusive terms [17][18][19][20][21][22], particle loss [23], calculations of quantum and Euler-scale correlations [24][25][26][27][28][29], as well as the incorporation of numerous beyond-Euler scale effects [30][31][32][33] (see also [34][35][36][37] in a special issue).Recently, GHD applied to a 1D Bose gas has been experimentally verified in a variant of the quantum Newton's cradle setup in the weakly interacting regime [15], and in a harmonic trap quench in the strongly interacting regime [16].In both cases, GHD provided an accurate coarse-grained model of the dynamics, exceeding conventional (classical) hydrodynamics.In addition to comparisons with experiments, GHD was benchmarked against other established theoretical approaches-most prominently for the 1D Bose gas and XXZ spin chain [9,10,14,15,24,27,33,[38][39][40][41][42].As the purpose of these initial benchmarks was to validate GHD, the typical dynamical scenarios considered were in regimes where GHD was expected to be a valid theory.In all such cases GHD demonstrated very good agreement with the alternative approaches.On the other hand, in scenarios involving, for example, short wavelength density oscillations due to interference phenomena (which are not captured by GHD), it was conjectured that GHD would nevertheless adequately describe spatial coarse-grained averages of the more accurate theories [14,15,32].More generally, it is of significant interest to scrutinize the performance of GHD by extending its benchmarks to a more challenging set of dynamical scenarios.This is important for understanding exactly how GHD breaks down when it is pushed towards and beyond the limits of its applicability.
In this Letter, we systematically benchmark the performance of GHD for the 1D Bose gas in several paradigmatic out-of-equilibrium scenarios.In particular, we focus on the regime of dispersive quantum shock waves emanating from a localized density bump of the type explored recently in Ref. [43].We use an array of theoretical approaches, including finite temperature c-field methods, the truncated Wigner approximation, and the numerically exact infinite matrix product state (iMPS) method, spanning the entire range of interaction strengths, from the nearly ideal Bose gas to the strongly interacting Tonks-Girardeau (TG) regime.We also analyse the dynamics of a localized density dip which sheds grey solitons, hence benchmarking GHD in scenarios not previously considered.In doing so we address the question of how well GHD predictions agree with coarse-grained averaging of the results of the more accurate theoretical approaches.Additionally, we explore the dynamics of a thermal quasicondensate in a quantum Newton's cradle setup [15,44,45] using Navier-Stokes type diffusive GHD [17,20,46], and address the question of characteristic thermalization rates [19,44].
Expansion from a localized density bump.-Webegin our analysis by considering dispersive quantum shock waves of the type studied recently in Ref. [43].

arXiv:2208.06614v3 [cond-mat.quant-gas] 13 Apr 2023
More specifically, we first focus on the weakly interacting regime of the 1D Bose gas of N particles, and consider the dynamics of the oscillatory shock wave train generated through a trap quench from an initially localized perturbation on top of a flat background to free propagation in a uniform box of length L with periodic boundary conditions [47,48].The weakly interacting regime is characterized by the Lieb-Liniger [49,50] dimensionless interaction parameter γ bg = mg/ 2 ρ bg 1, defined with respect to the background particle number density, ρ bg , where g > 0 is the strength of repulsive contact interaction and m is the mass of the particles.
In our first example, we consider the case of a large total number of particles, N = 2000, and γ bg = 0.01, so that the gas is in the Thomas-Fermi regime where the interaction energy per particle dominates the kinetic energy.We assume that the gas is initialized in the zero-temperature (T = 0) ground state of a dimple trap that results in the density profile of Eq. ( 1) given in Appendix A. At time τ = 0, the dimple trap is suddenly switched off, and we follow the evolution of the system in a uniform 1D trap.In Figs. 1 (a) and (b), we show snapshots of the density profiles at different times, and compare the GHD results with those obtained using the mean field Gross-Pitaevskii equation (GPE) and the truncated Wigner approximation (TWA) which incorporates the effect of quantum fluctuations ignored in the GPE [51].The snapshot at τ = 0.00014, which corresponds to the onset of a shock formation due to a large density gradient, shows excellent agreement between GHD and the more accurate microscopic approaches.Such an agreement at early times is remarkable given that GHD, which is derived here at Euler scale [52], becomes formally exact only in the limit of infinitely large length and time scales [11,14,27].
Past this time, the GPE and TWA show the formation of an oscillatory shock wave train, which has been identified in Ref. [43] as a result of self-interference of the expanding density bump with its own background.The interference contrast in this regime is generally large, even though the quantum fluctuations present in the TWA approach cause a visible reduction in contrast compared with the mean-field GPE result.The GHD prediction, on the other hand, completely fails to capture the oscillations, as these occur on a microscopic lengthscale.The characteristic period of oscillations here (which we note are chirped) is given approximately by the healing length l h = / √ mgρ bg (l h /L = 0.0057) which is smaller than the width σ (σ/L = 0.02) of the initial bump and hence represents the shortest lengthscale of the problem in the bulk of the shock wave train.Thus, even though the local density approximation (required for GHD to be applicable to an inhomogeneous system in the first place) is valid for the initial Thomas-Fermi density profile, the failure of GHD at later times is expected since it is not supposed to capture phenomena on microscopic lengthscales, which emerge here dynamically.
FIG. 1. Dimensionless density profiles ρ = ρL of quantum shock waves in the 1D Bose gas, as a function of the dimensionless coordinate ξ ≡ x/L at different times τ ≡ t/mL 2 .In (a) we show the initial (τ = 0) and time-evolved (τ = 0.00014) profiles of a weakly interacting gas at zero temperature, for γ bg = 0.01 and N = 2000 (with N bg 1761 the number of particle in the background).Due to the symmetry about the origin, we only show the densities for ξ > 0. In (b), the timeevolved profile is shown at τ = 0.0007.Panel (c) demonstrates the results of finite resolution averaging of both GPE and TWA data from (b) and compares them with the same GHD result.Panel (d) shows the same system as in (b), but at finite temperatures, simulated using the stochastic projected GPE (SPGPE) [43]; the dimensionless temperature T here is defined according to T = T /T d , where T d = 2 ρ 2 bg /2mkB [50].Panel (e) compares GHD predictions with exact diagonalization (ED) results in the TG regime (γ bg → ∞) for N = 1000 (N bg 884), at τ = 0.00004.In all examples, the initial profiles are characterized by the amplitude height β = 1 and dimensionless width of the bump σ = 0.02; see Appendix A for details.
Despite this failure, GHD clearly captures the average density of the oscillations for the fully formed shock wave train, similar to that shown in [13].This is consistent with the analysis of Bettelheim [53], who showed that the Whitham approach, which allows one to write equations for averaged quantities in the oscillatory shock wave train, is equivalent to GHD in the semiclassical limit (c = mg/ 2 → 0) of the Lieb-Liniger model [54,55].This is also consistent with the expectation that GHD in an interfering region would correspond to a coarse-grained average density [14,15].To quantitatively assess this expectation, we perform a type of convolution averag-FIG.2. Quantum shock waves at zero temperature for N = 50 particles (N bg 44.03), over the entire range of interaction strengths.In all examples, the initial density profiles (not shown) closely match Eq. ( 1) in Appendix A, with β = 1 and σ = 0.02.In all panels, we show the GHD (dashed lines) and iMPS (full lines) results for the evolved density profiles at two time instances.In (a) there is no phase coherence beyond the mean interparticle separation (1/ρ bg L 0.0227), whereas in (e) the shortest lengthscale that determines the characteristic period of oscillations is given by the width of the initial Gaussian bump σ (σ/L = 0.02), which is much smaller than the healing length l h (l h /L = 0.227).
ing that mimics the finite resolution of in-situ imaging systems used in quantum gas experiments (see Appendix B).As the imaging resolution is usually unable to resolve wavelengths on the order of the healing length (typically in the submicron range), one expects that such averaging will smear out the interference fringes seen in the GPE and TWA data-just as GHD implicitly does.In Fig. 1(c) we show the results of convolution-averaged density profiles performed on the GPE and TWA data of Fig. 1(b) and compare them with the same GHD curve.The level of agreement between all three curves is now remarkable-a result which was not a priori obvious for both GPE and TWA under this model of coarse-graining.This highlights the quantitative success of GHD in describing the dynamics on large scale despite interference or short-wavelength phenomena being present.
In our second set of examples, shown in Fig. 1(d), we consider the same shock wave scenario, except now for a phase fluctuating quasicondensate at finite temperatures.Here, the effect of thermal fluctuations is expected to lead to a smearing of the interference contrast due to a reduced thermal phase coherence length in the system, l T = 2 ρ bg /mk B T [56][57][58].A well-established theoretical approach to model this is a c-field stochastic projected GPE (SPGPE) approach [59, 60] (see also [44,[61][62][63]), and we indeed observe such smearing in Fig. 1(d) [64], in addition to seeing the expected very good agreement of GHD with these c-field results.
Our third example is shown in Fig. 1(e) and lies in the TG regime of infinitely strong interactions, γ bg → ∞.It further illustrates the same observation-that the performance of GHD improves with the loss of phase coherence in the system, wherein interference phenomena are suppressed.Here, we compare the predictions of GHD for the shock wave scenario at T = 0 with the results of exact diagonalization.In the TG regime, the system does not posses phase coherence beyond the mean interparticle separation 1/ρ bg , hence the absence of interference fringes in the evolution of a density bump whose initial width is larger than 1/ρ bg [43].Accordingly, we see very good agreement of GHD with exact diagonalization, ig-noring the small-amplitude density ripples that can be seen in the exact result.Such density ripples (which we note have different origin to Friedel oscillations) have been predicted to occur in the ideal Fermi gas by Bettelheim and Glazman [65] (see also [66]).By the Fermi-Bose mapping [67, 68], these same ripples should emerge in the TG gas, which we confirm here through exact diagonalization.However, their description lies beyond the scope of GHD as a large-scale theory [69].
The final set of examples for the evolution of a density bump is shown in Fig. 2. Here, we consider a range of interaction strengths, starting from very strong and going back [from (a) to (e)] to weak interactions, all at zero temperature and N = 50.We compare the GHD results with iMPS simulations, which are numerically exact at all interaction strengths [43].At this relatively low particle number, the strongly interacting regime displays Friedel oscillations which appear in the iMPS result and are, as expected, absent from the prediction of GHD.However, there is generally good agreement between GHD and iMPS at large scale.As the interaction strength is reduced, and hence the phase coherence of the gas increases, the Friedel oscillations disappear and interference fringes return, which now have period ∼ σ (with σ < l h ) since the gas is no longer in the Thomas-Fermi regime.The worst performance of GHD is observed for γ bg = 0.01, which lies in the nearly ideal (noninteracting) Bose gas regime for N = 50.In this regime, the local density approximation, intrinsic to GHD [14][15][16]50], is no longer valid even for the initial density profile, and we see that Euler-scale GHD breaks down both spatially and temporally, explaining the failure of GHD to agree with iMPS results even in the coarse-grained sense.
In addition to considering the dynamics of a localized density bump, we have also analyzed evolution of an initial density dip in a uniform background.This scenario is known to shed a train of grey solitons in the mean-field GPE treatment [47,48,70], and the results of comparison of GHD simulations with those of GPE and TWA are presented in Appendix C. The overall conclusions regarding the performance of GHD in this scenario are similar to those for a density bump, including good agreement of GHD with coarse-grained averages of GPA and TWA results in the soliton train region.
Quantum Newton's cradle in a thermal quasicondensate.-Ourfinal scenario for benchmarking GHD is in a variant of the quantum Newton's cradle setup for a weakly interacting 1D Bose gas in the quasicondensate regime.Namely, we analyze the release from a symmetric double-well trap to a single-well harmonic trap of frequency ω, similar to the type utilized in Ref. [15].Here, we use the SPGPE to simulate collisional dynamics and eventual thermalization, as in Ref. [44], and for the sake of one-to-one comparison, we also simulate the same system using the Navier-Stokes type of diffusive GHD [17,46], solved using a secondorder backwards-implicit algorithm [18,38,71].
Comparison of the results using the two methods are shown in Fig. 3, where we illustrate the evolution of the density distribution [(a) -for diffusive GHD, and (c) -for SPGPE] over the initial few oscillations, as well as after sufficiently long time, when the system has already thermalized.In Ref. [72] we give further details of how the final relaxed states were assessed within GHD and the SPGPE, whereas here, in Figs. 3 (b) and (d), we simply show the respective relaxed density profiles, along with their corresponding thermal equilibrium profiles from Yang-Yang thermodynamics [50,[72][73][74], as well as density profiles at earlier times illustrating their contrast to the relaxed state.The overall conclusion here is that GHD demonstrates excellent agreement with SPGPE in both short-and long-term dynamics, as well as in the characteristic thermalization rate [75].
We have also simulated the quantum Newton's cradle experiment in the original Bragg pulse scenario [1], except in a weakly interacting quasicondensate regime.In this scenario, we observe different thermalization rates in GHD and SPGPE simulations, and we discuss these results and the reasons behind the discrepancy in the Supplementary Material [72].
Summary.-We have benchmarked GHD in a variety of out-of-equilibrium scenarios in a 1D Bose gas against alternative theoretical approaches which are not limited to long-wavelength excitations.In particular, we have focused on systems supporting dispersive quantum shock waves and soliton trains, demonstrating that GHD generally agrees with the predictions of these approaches at sufficiently high temperatures and strong interactions.Here, the good agreement stems from a reduced phase coherence length of the gas, which in turn leads to a suppression of interference phenomena and therefore an absence of high-contrast short-wavelength interference fringes in the density.At low temperatures and weak interactions, where interference phenomena are more pronounced, the predictions of GHD only agree with a coarse-grained convolution averaging approximation.The effect of such averaging is similar to having finite imaging resolution in quantum gas experiments, and explains why GHD may perform well when compared to experiments, whilst departing from the predictions of theoretical approaches that are valid at short wavelengths.We have also benchmarked Navier-Stokes GHD within a quantum Newton's cradle setup for a doublewell to single-well trap quench of a weakly interacting quasicondensate, observing excellent agreement with the SPGPE in both transient dynamics and final relaxed state, as well as in the characteristic relaxation timescale.K. V. K. acknowledges stimulating discussions with I. Bouchoule, M. J. Davis, and D. M. Gangardt.This work was supported through Australian Research Council (ARC) Discovery Project Grants No. DP190101515.
Appendix A: Parametrization of the density bump.-Theinitial density profile in Fig. 1 (a), in dimensionless units, is set to where the dimensionless coordinate, time, and density are introduced, respectively, according to ξ ≡ x/L, τ ≡ t/mL 2 , and ρ(ξ, τ ) ≡ ρ(x, t)L, with ρ bg = ρ bg L = N bg being the dimensionless background density equivalent to the total number of particles in the background, 2σ )] from the normalization.In addition, the width and amplitude of the bump above the background are characterized by the dimensionless parameters σ ≡ σ/L and β > 0, respectively.
The associated trapping potential that is required for preparation of such a density profile as an initial ground or thermal equilibrium state of the 1D Bose gas in different regimes is discussed in Ref. [43].Within the meanfield approximation, described by the Gross-Pitaevskii equation, the density profile of Eq. ( 1) corresponds to the mean field amplitude being initialized as a simple Gaussian bump superimposed on a constant background, Appendix B: Finite resolution averaging -Finite resolution averaging procedure implemented in Fig. 1(c) emulates the finite spatial resolution of experimental absorption imaging systems.Following Ref. [76], we denote the impulse response function of the imaging system by A(x), which we here assume to be a normalized Gaussian.The impulse response for a pixel of width ∆ centered at x p is then, The measured atom number in the given pixel is then given by where N provides the correct normalization for the total particle number in the limit of zero pixel width.
In our particular example of such averaging, the density profile ρ(x) (at any given time step, with the time argument t being omitted here for notational simplicity) is convoluted with a Gaussian resolution function of width w = 1 µm and then averaged over a finite pixel size ∆ = 4.5 µm, as in Ref. [76].These absolute values translate to dimensionless values of w/L = 0.01 and ∆/L = 0.045, assuming L ∼ 100 µm, with results being generally insensitive to the exact values of these parameters around these typical values.For comparison, the healing length in this example is equal to l h /L = 0.0057.Considering 87 Rb atoms, which have a scattering length of a 5.3 nm, in a system of size L = 100 µm, this corresponds to an absolute healing length of l h = 0.57 µm.These choices of dimensionless parameters, and γ bg = 0.01, can be realized at a background density of ρ bg 1.8 × 10 7 m −1 , with an interaction parameter g 2 ω ⊥ a 1.4 × 10 −38 J•m [77], where ω ⊥ /2π 1.9 kHz is the frequency of the transverse harmonic trapping potential.
Appendix C: Dynamics of a localized density dip.-In this Appendix, we present the results of evo-FIG.4. Evolution of a density dip in a 1D Bose gas.Panel (a) shows the initial (τ = 0) and time-evolved (τ = 0.0005) density profiles from GPE, TWA and GHD simulations, for γ bg = 0.01 and N = 1688 (N bg 1761); panel (b) shows a timeevolved density profile at a later time (τ = 0.002), where we can see a fully formed train of three grey solitons in the meanfield GPE (full yellow) curve.Panels (c) and (d) compare the same GHD results (notice the different scale of vertical axis) at τ = 0.0005 and τ = 0.002 with the outcomes of finite resolution averaging of both GPE and TWA curves.In panel (e), we show a time-evolved snapshot of the density profile in the TG regime (γ bg → ∞) for N = 844 (N bg 880.5), and compare the GHD result with that of exact diagonalization (ED).Panel (f) is in the nearly ideal Bose gas regime, with γ bg = 0.01, N = 42 (N bg 44).In all examples, the initial density profile is given by Eq. ( 1) with β = −0.5 and σ = 0.02.lution of a localized density depression, after quenching (at time τ = 0) the initial trap potential with a localized barrier to uniform.We assume that the initial density profile is given by the same Eq.( 1), except with β being negative and satisfying −1 < β < 0.
In Figs. 4 (a) and (b), we consider the weakly interacting regime (with γ bg = 0.01) and show the results of the GPE, TWA, and GHD simulations, for a gas with N = 1688 atoms and the same N bg 1761 as in Fig. 1 (a).In this scenario, the steep gradient of the shock front forms as the background fluid flows inward and tries to fill the density depression.As a result, one first observes the emergence of large-amplitude structures, forming multiple density troughs, which then evolve into a train of grey solitons propagating away from the origin [47,48,55,70,78,79].The differences between the TWA and pure mean-field GPE results, seen in Figs. 4 (b), are consistent with previous observations [80][81][82] that quantum fluctuations lower the mean soliton speed and fill in the soliton core.The GHD result, on the other hand, fails to capture the solitonic structures, whose characteristic width (on the order of the microscopic healing length) lies beyond the intended range of applicability of GHD.
However, GHD still manages to adequately capture the coarse-grained description of the density across the soliton train, which is rather remarkable.This is seen in Fig. 4 (c) and (d), where we demonstrate the outcomes of finite resolution averaging applied to GPE and TWA results of panels (a) and (b), respectively.Similarly to Fig. 1(c), here we used the same normalized Gaussian resolution function of width 1 µm and adopted 87 Rb atoms as an example species for the relevant parameter values (see Appendix B).For panel (c) we used the same pixel size (∆ = 4.5 µm) as before, whereas for panel (d), due to the presence of fully formed grey solitons whose width is on the order of (2 − 4)l h , we used a twice larger pixel size (∆ = 9.0 µm).A larger pixel size here results in ∆/l h 16 1, which is required in order to comply with the large-scale framework of GHD.
The last two examples, shown in Figs. 4 (c) and (d), correspond, respectively, to the strongly interacting TG and nearly ideal Bose gas regimes.The overall behaviour and conclusions about the performance of GHD in these examples are the same as in the equivalent scenario of the density bump discussed earlier in Figs. 1 (e) and 2 (a).
Appendix D: Parametrization of the doublewell trap-The initial (pre-quench) double-well trap potential is set to V ( x) 2.16×10 −3 x 4 −5.27×10−1 x 2 in dimensionless form, where x = x/l ho , V = V / ω, and l ho = /mω, where ω is the post-quench single-well harmonic trap frequency.The initial dimensional temperature of the cloud, in harmonic oscillator units, is set to T = T /( ω/k B ) 205.In this configuration, the initial density profile for a total of N = 3340 atoms is double peaked, with the dimensionless interaction strength at either of the peaks given by γ max 0.0138.We briefly review the various formulations of generalized hydrodynamics (GHD) used in the main text and its relation to the Lieb-Liniger model of the 1D Bose gas.In particular we cover Euler-scale GHD, its zero-entropy subspace method, and its relation to the Wigner function formulation of the free Fermi gas.We also give a brief introduction to Navier-Stokes GHD and numerical tools used to simulate the dynamics of thermalization in Newton's cradle setups.

A. The one-dimensional Bose gas
The Lieb-Liniger model of a 1D Bose gas with repulsive contact interaction is a paradigmatic integrable quantum model, described by the following second-quantized Hamiltonian [S1, S2] Here, Ψ † (x) and Ψ(x) are the boson creation and annihilation field operators, obeying the canonical bosonic commutation relations [ Ψ(x), Ψ † (x ′ )] = δ(x − x ′ ), m is the bosonic mass, g is the one-dimensional interaction strength which is taken to be positive for repulsive interactions.For a uniform system of 1D density ρ, the exact eigenstates of the Lieb-Liniger model and its ground state properties can be found through the Bethe ansatz [S1], whereas the finite temperature equilibrium properties can be treated using Yang-Yang thermodynamic solutions [S3].The dimensionless interaction strength that characterises the system is introduced via γ = mg/ℏ 2 ρ.
In the presence of an external trapping potential, V (x), the Hamiltonian acquires an additional term, ´dx V (x) Ψ † Ψ.Through this, the model becomes generically non-integrable outside of the ideal (noninteracting) Bose gas limit or the Tonks-Girardeau limit of infinitely strong interactions.However, the exact solutions for a uniform gas can be still utilized for finding thermodynamic properties of inhomogeneous gases in the local density approximation (LDA) [S4].In the LDA, the Yang-Yang equations are solved using a local chemical potential µ(x) = µ 0 − V (x), where µ 0 is the global chemical potential of the system; as a result, the local dimensionless interaction strength γ(x) = mg/ℏ 2 ρ(x) acquires position-dependence through the inhomogeneity of the density profile ρ(x).

B. Euler-scale generalized hydrodynamics
GHD for the 1D Bose gas is expressed using the language of Yang and Yang's thermodynamic Bethe ansatz [S3, S5-S7].Here we present a brief formulation of first order (or Euler-scale) GHD in terms of the density of quasiparticles, f p (λ; x, t) [S3].Through f p (λ; x, t) we express the core 'Bethe-Boltzmann' equation of Euler-scale GHD in an inhomogeneous external potential [S8, S9], where λ is the rapidity or quasi-velocity variable of the quasiparticle excitations.
Through the occupation number function, n(λ; x, t), related to the density of quasiparticle and quasi-hole excitations, n = f p /(f p + f h ) [S3], one may express this hydrodynamic equation of motion in the form [S8, S9] Though they are equivalent, Eq. ( S3) is more convenient for simulating Euler-scale GHD, as it can be solved through the method of characteristics [S10, S11].The effective velocity, v eff [n], is a functional of the occupation number function, and is written in terms of the quasiparticle energy E and momentum p [S6, S7], and is the group velocity, v gr = ∂ λ E/∂ λ p, of the excitations renormalized by interactions through the dressing operation [S7, S11], where we suppress dependence on space and time variables for clarity.The scattering kernel, θ(λ), in Eq. ( S5) is the first model-dependent quantity utilized in GHD.For the Lieb-Liniger model it is given by θ(λ) = 2ℏg/(g 2 + (ℏλ) 2 ) [S1, S3, S12].The dressing equation may be seen as a generalization of the type of equation present throughout the thermodynamic Bethe ansatz; in particular, one may express the total state density, ].The other model-dependent quantities are the singleparticle eigenvalue functions, h i (λ) (i = 0, 1, 2, . . .), given in the Lieb-Liniger model by polynomials, h i (λ) ∝ λ i , with h 0 = 1, h 1 (λ) = p(λ) = mλ being the quasiparticle momentum, and h 2 (λ) = E(λ) = mλ 2 /2 -the quasiparticle energy, as examples [S1, S3].Correspondingly, the effective velocity simplifies to v eff (λ) = id dr (λ)/1 dr (λ), where id dr (λ) is the identity function id(λ) = λ, dressed according to Eq. ( S5), likewise 1 dr (λ) is the dressed unit function 1(λ) = 1 [S13].The evaluation of average charge densities takes a simple integral form [S3, S6, S7], where, for example, ⟨q 0 ⟩(x, t) = ρ(x, t) is the particle number density, ⟨q 1 ⟩(x, t) is the momentum or masscurrent density (where, in conventional hydrodynamics, ⟨q 1 ⟩(x, t)/mρ(x, t) = v(x, t) would correspond to the velocity field), and ⟨q 2 ⟩(x, t) = e(x, t) is the energy density.Evaluation of the initial equilibrium state in an arbitrary external potential requires knowledge of the local chemical potential, µ(x) = µ 0 − V (x).Under the LDA, the gas is treated as locally uniform across mesoscopic 'fluid cells', a process which coarse-grains over microscopic lengthscales [S14].Calculation of the equilibrium quasiparticle density on each fluid cell is achieved through an iterative numerical procedure utilizing the local chemical potential [S3, S15], with the GHD evolution of this density according to Eq. ( S3) occurring on macroscopic scales between these fluid cells.
In our GHD simulations of systems evolving from finite-temperature thermal equilibrium states, we used the iFluid software package [S15], which is an efficient, easily expandable open-source numerical framework based in Matlab.Simulations of systems evolving from zero-temperature ground states, on the other hand, were carried out using zero-entropy subspace methods found in Ref. [S13] and detailed below.
C. Zero-entropy subspace methods and the hard-core limit At zero temperature the occupation number simplifies to an indicator function, and 0 outside this region, where λ F is the Fermi rapidity and is generally dependent on the local effective chemical potential [S13].Thus, dynamics of the occupation number is restricted to its edges, or its 'Fermi contour' in systems of finite size [S16, S17].Numerically solving the GHD equation of motion of these systems may be accomplished through the zero-entropy algorithm presented in the supplementary material of Ref. [S13].This algorithm works with the finite set of curves defined by the edges of the zero-temperature occupation number, with each point on these curves given by a position and rapidity coordinate (x i , λ i ).For a single time-step, δt, the position x i is shifted by an amount v eff {λ} (λ i )δt, where calculating the effective velocity given in Eq. (S4) through the dressing operation in Eq. ( S5) is simplified due to the form of the zero-temperature occupation number described above.
In the Tonks-Girardeau limit of hard-core bosons, γ → ∞, at zero temperature there is an equivalence between the occupation number and the Wigner function in the semiclassical formulation for the free Fermi gas [S18].
Here, the Fermi rapidity is equivalent to the Fermi pseudovelocity λ F (x) = 2µ(x)/m [S13].In this limit, the GHD equation of motion is equivalent to the evolution of the Wigner function, and is given in the form of a Boltzmann equation [S13, S18], Quantum corrections to the semiclassical Wigner function description were analysed by Bettelheim et al., demonstrating the presence of long-lived 'quantum ripples' in the dynamics [S18-S20].Analysis presented in the supplementary material of [S13] demonstrated that these corrections to the local density approximation are negligible at large scales in the regime where In this limit, the Wigner function oscillates rapidly as a function of λ such that, when calculating observables using Eq.(S6), the integration over rapidity averages over these fast oscillations, and eliminates the oscillatory terms [S13].

D. Navier-Stokes GHD
Higher order corrections to the Euler-scale hydrodynamics presented above have recently been formulated, extending GHD to the Navier-Stokes scale and incorporating diffusive effects [S11, S21-S25].At this level, the hydrodynamic equation of motion is modified through the incorporation of a diffusion operator, D, arising through two-body scattering processes among quasiparticles [S26].Subsequently, it was shown in Ref. [S23], and later fully justified in Ref. [S24], that the diffusive hydrodynamic equation for the quasiparticle density, f p , in the presence of an external trapping potential is given by where (D∂ x f p )(λ; x, t) = ´dλ ′ D(λ, λ ′ ; x, t)∂ x f p (λ ′ ; x, t), and D(λ, λ ′ ; x, t) is the diffusion kernel, calculable from elements already present in the Euler scale theory [S21, S26].Through a scaling analysis performed at second order, it was shown in Ref. [S23] that corrections to Eq. (S9) arising from the Euler-scale force term do not contribute to the dynamics.It was thus demonstrated that, when an external trapping potential breaks the integrability of a gas, diffusive dynamics inevitably lead the system to thermalize at late times [S11, S23].Diffusive dynamics are soluble through a second-order backwards implicit algorithm, first demonstrated in Ref. [S27], and available within the iFluid package [S15].

II. QUANTUM NEWTON'S CRADLE
We provide further analysis for the thermalization of the double-well to single-well quantum Newton's cradle discussed in the main text.Additionally, we provide comparison between GHD and SPGPE dynamics for a quantum Newton's cradle under a Bragg pulse protocol, thermalization of which was recently studied in Ref. [S28].

A. Additional results for double-well to harmonic trap quench
Here, we further investigate dynamical thermalization of a double-well to single-well trap quench, demonstrated in Fig. 3 of the main text.To quantify the thermalization rate observed in Navier-Stokes GHD, we have access to the thermodynamic entropy density [S3] s which may be integrated for the total entropy per particle, S/k B N = N −1 ´dxs/k B .The total entropy per particle is known to plateau upon the system reaching a thermalized state [S23], and is demonstrated for the Navier-Stokes GHD simulations in Fig. S1(a), showing thermalization at time t ≃ 100/(ω/2π).For direct comparison between GHD and SPGPE results (as entropy is not accessible in SPGPE simulations), we also estimate the rate of thermalization through monitoring the peak density averaged over the region x/l ho ∈ [−2, 2], where l ho = ℏ/mω is the harmonic oscillator length.This quantity is plotted Figs.S1(a) and (c) for Navier-Stokes GHD and SPGPE simulation, respectively.Under non-equilibrium Newton's cradle dynamics, the peak density undergoes oscillations at twice the longitudinal trap frequency, eventually relaxing to a final thermal state at time t ≃ 100/(ω/2π) for both simulation methods.This thermalization time agrees with that extracted from the plateauing of GHD entropy per particle.As noted in Ref. [S23], observed thermalization times are generically observable-dependent, however the results presented here demonstrate that, for this system, the thermalization times estimated through relaxation of peak and entropy agree with each other.The temperature, T = T /(ℏω/k B ), and global chemical potential, µ 0 , of the relaxed state for the Navier-Stokes GHD result are fixed by the initial total energy and number of particles [S23].This determines the final temperature of the system to be T ≃ 213, with the corresponding Yang-Yang thermodynamic density profile shown as the cyan dotted line in Fig. 3(b) of the main text.Temperature estimation of the relaxed state for the SPGPE result is achieved via Yang-Yang thermometry of the relaxed density profile [S29, S30], and is shown as the cyan dotted line in Fig. 3(d) of the main text at a temperature of T ≃ 216.The small difference in the temperatures extracted for the relaxed states in GHD and SPGPE simulations comes from a small difference between the respective final density profiles, so that the best-fit density profiles from Yang-Yang thermodynamics return slightly different values of the respective temperatures.
Additionally, we calculate the proximity of evolving density distributions in both GHD and SPGPE simulations to the respective Yang-Yang thermal density distributions through the Bhattacharyya statistical distance [S31] where B(P, P ′ ) is the Bhattacharyya coefficient of two normalized probability density functions P (x i , t) and P ′ (x i ) of the same discrete variable x i , Here, P (x i , t) = ρ(x i , t)/ i ρ(x i , t) is the normalized evolving density profile of either GHD, shown in Fig. 3 (a), or SPGPE, shown in Fig. 3 (c) of the main text, with x i being the position on our computational lattice.P ′ (x i ) = ρ YY (x i )/ i ρ YY (x i ), on the other hand, is taken to be the normalized density profile of the Yang-Yang thermal state, ρ YY (x i ), fitted to the final relaxed density profile of GHD or SPGPE evolution.At any given time t, the Bhattacharyya distance serves as a measure of similarity between the instantaneous distribution P (x i , t) and P ′ (x i ).For P → P ′ , the Bhattacharyya coefficient becomes B(P, P ′ ) → 1 (due to the normalization condition) and hence D B → 0, implying complete overlap of the two distributions.As we see from Figs. S1 (b) and (d), as the evolving density profile approaches the relaxed state, the Bhattacharyya distance, D B , tends to a vanishingly small constant, indicating near perfect overlap between the relaxed and Yang-Yang thermal states.
We further point out here that as the SPGPE is a classical field method [S33, S34], it has an inherent dependency on the high-energy cutoff ϵ cut separating out the classical field region of relatively highly occupied modes from the sparsely occupied modes treated as an effective reservoir.More specifically, the results of SPGPE simulations can be cutoff dependent subject to the observable in question.In the simulations performed here, the highenergy cutoff was chosen according to the prescription of Ref. [S35], namely, such that the mode occupancy of the highest energy mode in the harmonic oscillator basis (of Hermite-Gauss polynomials) was on the order of ∼ 0.3.Such a cutoff mode occupancy, which is somewhat lower than a more conventional cutoff choice at mode occupancy of the order of ∼ 1 (see, e.g., [S36-S38] and references therein), must be taken in order to faithfully reproduce the tails of the initial thermal equilibrium density distribution.In doing so, we have also observed that the simulation of the quantum Newton's cradle system under variation of the high-energy cutoff shows a dependence of the thermalization time on this cutoff mode occupancy.In particular, as the cutoff mode occupancy is reduced (corresponding to an increase in the high-energy cutoff ϵ cut ), not only does the initial density profile match better to that of Yang-Yang thermodynamics, but the thermalization times are additionally shortened.As increasing the high-energy cutoff corresponds to including a larger number of thermally occupied modes, which speed up collisional relaxation, such shortening of thermalization times is expected.Ultimately, as we have seen from the results presented in Fig. S1, the characteristic thermalization times extracted from GHD and SPGPE agreed well with each other under the cutoff mode occupancy of 0.3.Accordingly, all SPGPE simulations reported here, including the ones corresponding to quantum Newton's cradle with a Bragg pulse (see next section) were carried out with ϵ cut corresponding to this optimal cutoff mode occupancy of 0.3.

B. Bragg pulse quantum Newton's cradle
Collisional dynamics and thermalization of a quantum Newton's cradle under a Bragg pulse protocol in the finite-temperature quasicondensate regime has recently been studied in Ref. [S28].Here, we provide an additional point of comparison between GHD and SPGPE dynamics through simulation of this system, studying the rate of convergence to thermalization within the two simulation methods.Dynamics are instigated from a thermal state in a harmonic trap of frequency ω, characterized by dimensionless interaction parameter γ 0 = 0.01 in the trap centre, with temperature T ≃ 152, and total atom number N ≃ 1960.A Bragg pulse is then applied to initiate Newton's cradle dynamics; this splits the initial atomic wavepacket at rest into two counter-propagating halves corresponding to ±2ℏk 0 diffraction orders of Bragg scattering [S28, S39], where k 0 is the Bragg momentum in wave-number units.
Implementing the Bragg pulse within the SPGPE method consists of replacing each stochastic realisation of the initial thermal equilibrium state, ψ(x, t = 0 − ) by a coherent superposition, ψ(x, t = 0 + ) = 1 √ 2 (e i2k0x + e −i2k0x )ψ(x, t = 0 − ), which subsequently evolve in time according to the projected GPE [S28, S39].Such a superposition of wavefunctions with positive and negative momentum boosts of ±2k 0 is known to be an excellent approximation to experimental implementation of the Bragg pulse [S28, S39, S40] via an external periodic lattice (Bragg) potential formed by two counterpropagating laser beams [S41].The microscopic effect of the realistic Bragg pulse in the SPGPE simulation can be seen as the fast oscillating interference fringes in the spatial density profile at early times, during the periods when the two counter-propagating halves of the density profile overlap spatially.
Simulating a Bragg pulse protocol in GHD, by contrast, is done through the quasiparticle density distribution, f p (λ; x, t), which does not directly depend on the bare atomic momenta, instead relying on the rapidity of quasiparticles, and is generally not equivalent to the momentum distribution, except in free models [S16, S42, S43].Implementing a Bragg pulse in GHD requires one to work with the quasiparticle density distribution of the initial thermal state, f t=0 − p (λ; x, t), adding positive and negative momenta to the quasiparticles with equal probability.Correspondingly, the post-pulse quasiparticle distribution may be modeled as Thus, in contrast to the quantum Newton's cradle presented in the main text, which is instigated in exactly the same way-through a sudden quench of the actual external trap potential V (x) from double to single well-in both the SPGPE and GHD simulation, the implementation of the Bragg pulse quantum Newton's cradle is achieved in SPGPE and GHD via two different approximations of the post-Bragg pulse state.Additionally, as the GHD is a large-scale theory, its post-Bragg pulse state fails to capture the fast-oscillating interference fringes observed in the SPGPE at early times.
In GHD implementation, one typically chooses the dimensionless quasimomentum of the Bragg pulse, λ Bragg = (m/ℏ)λ Bragg l ho (given here in harmonic oscillator units), to be equal to the Bragg momentum, q 0 ≡ k 0 l ho , however this is an approximation that assumes a large momentum difference between the clouds [S43].Here, we instead choose the Bragg pulse quasimomentum such that the total energy imparted in the GHD simulation is equal to the energy difference between the initial thermal state and the post-Bragg pulse state calculated in the SPGPE.Using this method, we observe a large difference between the physical Bragg momentum and the corresponding GHD quasimomentum under a low momentum Bragg pulse (∼219% difference at a Bragg momentum of q 0 = 1), with this difference decreasing for larger momentum Bragg pulses (∼ 16.8% difference at a Bragg momentum of q 0 = 5).
We first investigate a low momentum Bragg pulse of q 0 = 1 and illustrate the dynamics of the density profile for GHD and SPGPE simulations in Fig. S2(a) and (d), respectively.We observe a qualitative similarity in the oscillations of collisional dynamics between the two simulation methods, however, there is a clear disagreement in the rate of thermalization.The observed thermalization times of the Navier-Stokes GHD simulation is observed to be t f ≃ 80/(ω/2π), after approximately 160 oscillation periods, relaxing to a state of temperature T ≃ 214.In comparison, the SPGPE simulation of this q 0 = 1 Bragg pulse Newton's cradle system thermalizes at an earlier time of t f ≃ 30/(ω/2π), after approximately 60 oscillation periods, to a relaxed state at T ≃ 261.The observed difference in thermalization rates between Navier-Stokes GHD and the SPGPE methods of simulating this system likely stems from the approximate method utilized for simulating the Bragg pulse within GHD, described above.
Simulation of the same system, but with a larger Bragg momentum of q 0 = 5, is illustrated in Fig. S3(a)-(c) for GHD, and in Fig. S3(d)-(f) for SPGPE, respectively.Here, we observe a greater disparity in the thermalization rates with SPGPE simulations thermalizing at around t f ≃ 40/(ω/2π) (approximately 80 oscillation periods), and GHD at t f ≃ 350/(ω/2π) (approximately 700 oscillation periods).This disparity is confirmed when comparing the average peak density between GHD and SPGPE, shown in Fig. S3(b) and (e), respectively.
The density profile of the final relaxed state for the GHD simulation is shown in Fig. S3 (c), with its matching thermal density profile of temperature T ≃ 336.Likewise, the relaxed density profile for the SPGPE simula-tion is shown in Fig. S3 (f), with the density profile of its matching Yang-Yang thermodynamic state of temperature T ≃ 348.We thus observe that, for the Bragg pulse quantum Newton's cradle, the discrepancy in thermalization rates is present regardless of the momentum of the simulated Bragg pulse.However, short and intermediate time dynamics, along with the final relaxed states of both simulation methods, are seen to be qualitatively similar.
To summarise, we emphasize that the quantum Newton's cradle in double-well to single-well quench is implemented in exactly the same way in both the GHD and SPGPE approaches-via a sudden quench of the external trapping potential at time t = 0. Accordingly, as we discussed in the main text, the GHD and SPGPE thermalization times agree with each other in this setup.In contrast to this, the original Bragg pulse quantum Newton's cradle scenario of Ref. [S41], which employs a splitting of the initial quasicondensate wavefunction into two counter-propagating halves, is implemented in GHD in an approximate and qualitatively different way to the SPGPE [S28, S43].Because of this, we observe significantly different thermalization rates (especially for large Bragg momenta) using GHD and SPGPE simulations, even though the overall dynamics of collisional oscillations are still qualitatively very similar.

FIG. 3 .
FIG. 3. Evolution and thermalization of the density distribution ρ(x, t) in a quantum Newton's cradle setup initialized from a double-well to single-well trap quench, simulated using (a)-(b) Navier-Stokes GHD, and (c)-(d) SPGPE.The initial cloud of N = 3340 atoms at temperature T = 205 (in harmonic oscillator units) is prepared in a thermal equilibrium state of a symmetric double-well trap potential (see Appendix D for details).Panel (b) demonstrates the relaxed density profile of the Navier-Stokes GHD evolution at t = 100/(ω/2π) (black solid line), alongside a best-fit thermal equilibrium profile from Yang-Yang thermodynamics at T 213 (cyan dotted line), and an additional GHD density profile at earlier time t = 6.79/(ω/2π) (red dashed line).Panel (d) is the same as (b), but for the SPGPE, with the relaxed density profile at t 100/(ω/2π), Yang-Yang thermodyanmic density profile of T 216, and an additional density profile at t = 6.81/(ω/2π).
FIG. S1.Convergence to thermalization for the double-well to single-well trap quench demonstrated in Fig. 3 of the main text.Panel (a) shows the evolution of the peak density for the Navier-Stokes GHD simulation (averaged over the region x/l ho ∈ [−2, 2], where l ho = ℏ/mω is the harmonic oscillator length), alongside the respective total entropy per particle, both of which plateau upon reaching the final thermal state.Panel (b) demonstrates the Bhattacharyya statistical distance of the GHD evolving density profile to the corresponding Yang-Yang thermodynamic density profile (fitted to the final relaxed GHD profile) of temperature T = T /(ℏω/kB) ≃ 213 shown in Fig. 3 (b) of the main text.The final observed Bhattacharyya distance here is DB = 2.23 × 10 −4 , which is rather small and hence indicates near-perfect overlap of the two distributions.Similarly, panel (c) shows the evolution of the peak density for the SPGPE simulation (again averaged over the region x/l ho ∈ [−2, 2]) as it approaches thermalization and plateaus; (d) demonstrates the Bhattacharyya distance of the evolving SPGPE density profile to the Yang-Yang thermodynamic density profile of temperature T ≃ 216 shown in Fig. 3 (d) of the main text, with a final observed distance of DB = 4.28 × 10 −3 .

FIG. S2 .
FIG. S2.Dynamics of thermalization of a harmonically trapped quasiondensate for GHD and SPGPE evolution in a quantum Newton's cradle setup under a Bragg pulse of momentum q0 = k0l ho = 1.The corresponding value of Bragg quasimomentum for this system is λBragg = 3.19.Initial thermal states are characterised by a dimensionless interaction parameter at the trap centre of γ0 = 0.01, and dimensionless temperature T ≃ 152, with a total atom number of N ≃ 1960 [S32].Panel (a) shows Navier-Stokes GHD evolution of the density profile, ρ(x, t)l ho ; (b) demonstrates the evolution of the respective peak density, averaged over the region x/l ho ∈ [−2, 2] (blue, fast-oscillating curve), alongside the total entropy per particle of the GHD simulation (red, smooth curve), which both plateau upon reaching the final thermal state; (c) shows the density profile of the relaxed GHD state after dynamics at time t = 80/(ω/2π) (black solid line), as well as the Yang-Yang thermodynamic density profile (cyan dotted line) of temperature T ≃ 214 that fits best to that of the relaxed GHD state.Also shown in (c) is a GHD density profile at an earlier time t ≃ 6.56/(ω/2π) (red dashed line), chosen to illustrate its deviation from the density profile of the final relaxed state.Density profile evolution under SPGPE simulation is shown in (d); panel (e) demonstrates the evolution of the respective peak density (averaged over the region x/l ho ∈ [−2, 2]) as it approaches thermalization and plateaus; in (f) we show the relaxed density profile after dynamics at time t = 30/(ω/2π) (black solid line), as well a best fit Yang-yang thermodynamic density profile of temperature T ≃ 261 (cyan dotted line), and the SPGPE density profile at an earlier time t ≃ 6.54/(ω/2π) (red dashed line), which is already much closer to the final relaxed density profile compared to that of GHD simulation around the same time.

FIG. S3 .
FIG. S3.Same as in Fig. S2, but for q0 = 5.The corresponding Bragg quasimomentum for this system is λBragg = 5.84, whereas the total atom number, and the initial values of γ0 and T are the same.The final relaxed density profile from GHD in (c) is at time t = 350/(ω/2π) (black solid line), and is compared with a corresponding Yang-Yang thermodynamic density profile at temperature T = 336 (cyan dotted line), and a respective density profile taken at time t ≃ 21.5/(ω/2π) (red dashed line).The relaxed SPGPE state in (f) is given at a time t = 40/(ω/2π) (black solid line), with a respective Yang-Yang thermodynamic density profile of temperature T = 348 (cyan dotted line), and an additional density profile demonstrated at a time of t ≃ 22.1/(ω/2π) (red dashed line).
weakly interacting one-dimensional Bose gas, Phys.Rev.A 86, 033626 (2012).[59]Y.Castin, R. Dum, E. Mandonnet, A. Minguzzi, and I. Carusotto, Coherence properties of a continuous atom laser, Journal of Modern Optics 47, 2671 (2000).[60]P.Blakie, A. Bradley, M. Davis, R. Ballagh, and C. Gardiner, Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques, Advances in Physics 57, 363 (2008).[61]I.Bouchoule, S. S. Szigeti, M. J. Davis, and K. V.For the examples considered in Fig. 1(d), the thermal phase coherence lengths are lT /L 0.1 for T = 0.01 and lT /L 0.01 for T = 0.1.These values, in turn, are comparable or smaller than the width of the initial density bump σ = 0.02, implying a reduced phase coherence over the extent of the bump and hence loss of interference contrast upon expansion of the bump into the background [43].[65]E.Bettelheim and L. Glazman, Quantum Ripples Over a Semiclassical Shock, Phys.Rev.Lett.109,260602(2012).[66]In the example of Fig.1(e), which is for N = 1000 particles, we have been able to discriminate between the