UvA-DARE

Describing and understanding the motion of quantum gases out of equilibrium is one of the most important modern challenges for theorists. In the groundbreaking Quantum Newton Cradle experiment [ Kinoshita, Wenger and Weiss, Nature 440, 900 (2006) ] , quasi-one-dimensional cold atom gases were observed with unprecedented accuracy, providing impetus for many developments on the effects of low dimensionality in out-of-equilibrium physics. But it is only recently that the theory of generalized hydrodynamics has provided the adequate tools for a numerically efﬁcient description. Using it, we give a complete numerical study of the time evolution of an ultracold atomic gas in this setup, in an interacting parameter regime close to that of the original experiment. We evaluate the full evolving phase-space distribution of particles. We simulate oscillations due to the harmonic trap, the collision of clouds without thermalization, and observe a small elongation of the actual oscillation period and cloud deformations due to many-body dephasing. We also analyze the effects of weak anharmonicity. In the experiment, measurements are made after release from the one-dimensional trap. We evaluate the gas density curves after such a release, characterizing the actual time necessary for reaching the asymptotic state where the integrable quasi-particle momentum distribution


Introduction
In 2006, the pioneering experiment of the "Quantum Newton Cradle" [1] (QNC) provided a groundbreaking demonstration of the fundamental importance of the large number of conservation laws in the description of out-of-equilibrium one-dimensional (1d) quantum systems. In this experiment, an ultracold gas of Rubidium atoms is confined to one dimension by a strong transverse optical trap, and weakly confined to a finite region by a longitudinal quasi-harmonic potential. A sequence of Bragg pulses imparts a linear combination of oppositely-directed momenta to the initial cloud that lies at the center of the trap. After a short dephasing period, two independent clouds emerge, which oscillate within the trap and meet twice every period. Surprisingly, upon meeting and interacting, the clouds do not thermalize to a single zeromomentum cloud, as would happen in an ordinary gas. Instead, the two clouds re-emerge and continue their oscillations. This is likened to the Newton cradle, the famous desktop toy in which, upon collision, momentum is transferred between the end beads.
To a good approximation [2], the dynamics of a 1d gas of N identical bosonic atoms with mass m at positions x i is described by a hamiltonian with "delta" repulsion called the Lieb-Liniger model [3,4], where g > 0 is the repulsion strength and V (x) the longitudinal trapping potential. In the absence of a potential V (x), this model is integrable -it has an extensive number of conserved quantities-, a property that was conjectured to be at the root of the lack of thermalization in the QNC experiment. This has given rise to a wealth of theoretical developments on the generalization of thermalization in integrable models, following Refs. [5,6]. It was understood that, even after very long times, integrable systems fail to converge to a thermal Gibbs state. Instead, they reach a macrostate that is entirely characterized by the distribution of quasiparticles [7], in a way that parallels and generalizes the early work of Yang and Yang on the thermodynamics of the 1d Bose gas [8].
Realistic theoretical modeling of the QNC experiment requires to deal with N ∼ 10 2 -10 3 particles at finite repulsion strength -i.e. away from the hard-core limit g → +∞ that can be treated with exact determinantal methods [9]-and with an inhomogeneous potential V (x). This has remained completely out of reach of modern theoretical tools, and has represented one of the most prominent challenges in quantum many-body theory in the past decade. On the analytical side, the difficulty is twofold. First, it is known that the external potential breaks integrability. How, then, is the physics of integrability coming into play? Second, the setup is highly inhomogeneous, which is a major issue for most analytical techniques available in 1d [10]. On the numerical side, the situation is not better: modern tools for out-of-equilibrium quantum many-body physics like tDMRG [11] are limited to small numbers of particles and small times in QNC-like setups [12], while numerical methods based on integrability [13][14][15][16][17] break down because of the strong inhomogeneity.
A new set of theoretical tools and ideas have come to the fore in 2016 that, as we argue here, provide such a realistic modeling. These pertain to the theory of "generalized hydrodynamics" (GHD) [18,19], a hydrodynamic approach to 1d integrable models that is able to account for a wide variety of inhomogeneous situations [20], including states obtained from domain-wall initial boundary conditions [18,19,[21][22][23], the effects of external potentials [24] and the propagation of waves [25,26]. The last two of these works show excellent agreement with full numerical simulations of the Lieb-Liniger model for N ∼ 40 particles [25] based on NRG+ABACUS [13][14][15][16][17]27], and with tDMRG in the XXZ chain [26]. The goal of this Letter is to apply the newly developed GHD framework to provide a quantitatively accurate modeling of the QNC experiment that remains easily tractable numerically. large enough such that it contains a thermodynamically large number of bosons, and small enough such that the gas is homogeneous throughout the cell. As in any other hydrodynamic theory [28], at this scale GHD assumes that local maximization of entropy has occurred. In standard approaches [12,[29][30][31][32], this would imply that each fluid cell is locally in a Galilean boost of a thermal Gibbs state. Instead, GHD keeps track of the infinite set of conserved charges by representing the local macrostate by its distribution ρ p (θ ) of quasi-particles with velocity θ . In GHD, this distribution is position-and time-dependent, and is denoted ρ p (θ , x, t). In terms of the quasi-particles, the density of bosons n(x, t) = e iH t N j=1 δ(x − x j ) e −iH t is recovered, in the thermodynamic limit by integrating locally over all the quasi-particles: Similarly, all other densities of conserved charges in the Lieb-Liniger model, such as momentum, energy, or others, may be expressed as integrals dθ h(θ )ρ p (θ , x, t) for suitably chosen functions h(θ ), see e.g. Ref. [18].
The conservation of quasi-particle densities exchanged between neighbouring fluid cells fully fixes the dynamics [18,19]. Taking into account the trapping potential V (x), the GHD equation reads [24] where the effective velocity v eff itself depends on the local distribution of quasi-particles Here id(θ ) = θ and 1(θ ) = 1, and the "dressing" is defined for any function f (θ ) as the solution of the linear integral equation f dr (θ ) = f (θ )+ dθ Thus, GHD is a large-scale approach to the (inhomogeneous) Lieb-Liniger model (1) that requires to (i) specify an initial condition ρ p (θ , x, 0) at t = 0, and (ii) solve the partial differential equation (3).
A number of techniques are available to numerically solve Eq. (3). Numerical methods for solving partial differential equations [33] can be adapted to GHD -a discussion about this can be found in Ref. [26]-. Other methods have been discovered in the past months, that can be easier to implement, and in certain simple cases, possibly more efficient (an example is an exact solution expressed as a system of integral equations [34], which can be solved recursively on a computer very quickly; but until now this applies only to the case without an external potential V (x)).
An efficient technique [25] is that of the zero-entropy subspace. It uses the simplification of GHD for initial states which have zero Yang-Yang entropy, such as zero-temperature states. In this case, the solution is expressed via space-time dependent Fermi points, whose equations are simple enough to be directly amenable to numerical solution. Zero-entropy GHD works in the presence of external potentials, and it also sheds light on phenomena such as shock formation and dissolution [25]. Zero-temperature initial states are often a good first approximation for low-temperature experimental setups. It is possible to analyze such a zero-temperature version of the QNC by taking the initial state as the ground state in a double-well potential that splits the gas into two. By releasing the two clouds in a single harmonic well, the main phenomenon -the lack of thermalization upon cloud collisions -of the QNC experiment is observed (see Appendix I).
Finally, the technique which we chose to use here is a classical molecular simulator [35]. This is essentially a Monte Carlo technique: a gas of classical particles [36] is initially sampled Figure 1: Evolution of the density of quasi-particles ρ p (θ , x, t) -here plotted in the (x, θ )-planein the QNC setup, with parameters given in the text. The solution of the GHD equations are obtained from the flea gas [35]. The results are displayed for the harmonic trap (top row) and the weakly anharmonic one (middle row), on one period of the (quasi-)harmonic trap. (Bottom row) Corresponding density of particles n(x, t), for the harmonic trap (blue) and the anharmonic one (red).
according to the distribution of integrable quasi-particles determined, via appropriate integral equations, by the initial quantum state, and then let to evolve in a deterministic fashion by a specific dynamics that encodes the interaction of the quantum gas. The hydrodynamic description of this classical gas is, at the Euler scale, exactly the same GHD equation (3) as the one found in the quantum gas [35]. This quantum-classical equivalence provides an extremely efficient method for simulating solutions to the GHD equation, that is able to account simultaneously for arbitrary initial conditions and for external potentials.

Modeling the Bragg pulse sequence
We start from the thermal Gibbs state in a (quasi-)harmonic potential V (x). The exact experimental distribution pre-pulse is not known, and is not expected to be exactly thermal as cooling methods deplete the large-momentum region. However a thermal distribution is expected to be a good approximation, accounting well enough for the remaining energy that brings the system away from absolute zero temperature. The pre-pulse density ρ Gibbs p (θ , x) is obtained by searching for the finite-temperature hydrostatic solution of (3), which can be shown [24] to be equivalent to a local density approximation (LDA) [37], obtained using the Thermodynamic Bethe Ansatz [8]. This technique is known to give an accurate description of atomic gases at equilibrium in quasi-harmonic traps [38,39]. We model the post-pulse distribution ρ p (θ , x, 0) by imparting, in a random fashion, positive and negative momenta to the quasi-particles. More precisely, corresponding to a pulse where the momentum of a quasi-particle is kicked by ±mθ Bragg with equal probability 1 2 . In fact, results of [40] show that the momentum distribution function (MDF) of particles (both the real bosonic particles, and their fermionic Jordan-Wigner transform) are affected in this way by a Bragg pulse. To a good approximation the same holds In the anharmonic case, the distribution ρ p (θ , x) is strongly stirred up after a few trap periods, and it goes to stationary state that looks rotationally invariant in the (x, θ )-plane. for the quasi-particle MDF because of the large momentum difference between clouds: the inter-cloud interaction is therefore effectively hard-core, hence screened for the fermions (the intra-cloud uniform momentum kicks by ±θ Bragg are Galilean transformations). We have also considered more refined models for the Bragg pulse (see Ref. [40] for a detailed discussion), but we find that it does not affect drastically the evolution at later times, and this simple version is already in good agreement with what is seen in the experiment.

Results for the (harmonic) QNC
We work with the following parameters, which are close to the ones given in Ref. [1]. We take N = 350 bosonic atoms with mass m = 142.9 × 10 −27 kg in a harmonic trap V (x) = mω 2 x 2 /2 with period τ = 2π ω = 13ms, and with repulsion strength c = mg ħ h 2 = 12.3µm −1 . The pre-pulse state is a thermal Gibbs state at temperature T = 10µm −2 × ħ h 2 /m 57nK. We then apply the Bragg pulse sequence (random kick of the particle velocities by ±θ Bragg = ±ħ hk Bragg /m) with 2k Bragg = 25µm −1 . In the initial state, the density at the center of the trap is n(0) = 11.5µm −1 , which gives a dimensionless interaction strength γ = c/n(0) 1.07, close to the experimental value [1]. Notice that we are far from both the hard-core (or Tonks-Girardeau) limit γ 1 and the weakly interacting (or Gross-Pitaevskii) limit γ 1.
After the Bragg pulse, the dynamics of the Bose gas is given by Eq. (3), which we solve with the molecular dynamics simulation, see Fig. 1. We observe that the two blobs, originally separated in momentum space and symmetric with respect to θ → −θ , evolve by performing a deformed rotation-like movement around the origin of phase space. At time t τ/4, the two blobs have zero spatial overlap, corresponding to two well separated clouds in real space. At time t τ/2, they have again overlapping spatial support. Their evolution is not drastically affected by this overlapping, and it is clear, in the phase space picture, how the two atomic clouds can pass through each other. The actual gas density -bottom row of Fig. 1is obtained by integrating ρ p (θ , x, t) over θ according to Eq. (2).
We observe that the blobs are slightly slowed down when they overlap: it takes them slightly longer than a period τ to come back to their original position along the vertical axis. The inter-cloud interaction at spatial overlap is weak because of large momentum separation, but nonzero, with perceptible effect. Further, after several periods, the blobs elongate transversally, roughly towards the center of rotation in phase space, see Fig. 2. This slow "many-body dephasing" effect, controlled by how far the effective velocity is from the bare velocity 1 , is also due to inter-cloud interactions. Notice how, in Fig. 1, the blobs' shape stay relatively unchanged until just before half a period, t = 0.4τ, while after the clouds have passed through each other, t = 0.6τ, slight modifications have occurred. Particles of one cloud, going through the other cloud, are scattered on the scale of the scattering length. We observe that the particle density spreads mostly towards lower energies. The spreading must indeed be stronger at lower energies than at higher energies, because total energy is conserved, while the change in energy per distance (per momentum) is greater at higher energies than it is at lower energies (as is clear from the form of the potential, and of the kinetic energy as function of momentum). Intra-cloud scattering is also present but its effect is weaker as fewer scattering events occur.

Effects of weak anharmonicity
In the QNC experiment, the trapping potential is not exactly harmonic. To study the effects of anharmonicity, we now replace the harmonic trap by V (x) = 1 π 2 ω 2 2 (1−cos( πx )), both before and after the Bragg pulse. This form is chosen in order to mimic actual experimental setups, where potentials are often close to trigonometric functions. In particular, the anharmonicity has the property that the potential is smaller than harmonic further away from the centre (i.e. it becomes flatter). We chose = 100µm, and all other parameters are identical to the ones of the harmonic case. In Fig. 1 (second row), we see that the two blobs get deformed much more quickly than in the harmonic case; the distribution ρ p (θ , x, t) gets more and more stirred up after few periods (Fig. 2, second row). A similar effect was recently observed for a single particle in an anharmonic trap in Ref. [41]. In particular, we observe in Figs. 1 and 2 (second row) that the blobs elongate longitudinally, roughly in the direction of their motion in phase space. This is an effect of the anharmonicity < ∞: because the potential becomes flatter far from the origin, particles with higher velocities take longer to come back to their original position. This "single-body dephasing" effect is well captured by the spreading of blobs of noninteracting particles, and can be quantified by evaluating the nonzero difference ∆t between the periods of independent particles of different velocities. With minimal (maximal) velocity of 10µm/ms (40µm/ms) within one initial blob, one finds (∆t) max ≈ 1.2ms. The anharmonic mixing time is approximately τ · τ/(∆t) max ≈ 11τ, in agreement with Fig. 2. Many-body dephasing is also present in the anharmonic case. Noticeably, without interactions the original blobs would simply disintegrate into long spiraling filaments. We see instead at 10τ new structures forming. These appear as many-body dephasing takes effect, causing filaments to merge and high-energy (longer-period) tails to scatter to lower energies, opening the way for the quicker single-body effects to reform new blobs.

Discussion: thermalization?
We now turn to the fundamental question that was raised in 2006 in Ref. [1]: does the gas thermalize after a sufficiently large number of oscillations? We first consider this question within pure GHD. Within this context, the answer is negative. The reason is that, in addition to its particle number and its global energy GHD keeps track of many initial features. Indeed, we found that even in the presence of a trapping potential V (x), the GHD equation (3) possesses infinitely many conserved quantities Q(η), with a continuous parameter η ∈ [0, 1] -see Appendix III, that give rise to conservation of the Yang-Yang entropy [8] and generalizations thereof. This is incompatible with the system converging to a Gibbs state, even at infinite time. Much like phase-space preservation in classical mechanics may lead to fractal trajectories, in GHD these constraints give rise to fine structures developing at ever decreasing scales; see for instance the times of order t 6τ − 10τ in the anharmonic case of Fig. 2, Ref. [41] for a similar study in the hard-core (γ → ∞) limit, and the very recent work on the dynamics of a classical hard rod gas in a trap [42]. However, there is a sense in which GHD converges to a smooth, stationary state: this is through coarse-graining.

Microscopics and coarse-graining
GHD is an idealized hydrodynamic description with no UV cutoff. In contrast, the 1d Bose gas (1) is a microscopic model described by GHD only at larger scales, and is therefore only an imperfect realization of GHD. Now GHD predicts the appearance of fine structure in phase space, i.e. strong variations of ρ p (θ , x). The way the microscopic model treats the appearance of these UV degrees of freedom in GHD amounts to coarse graining: the density ρ p (θ , x) predicted from GHD gets essentially replaced by its average over fluid cells [x, x +d x] of size larger than the inter-particle distance, so eliminating the UV degrees of freedom from the problem. As a result of this coarse-graining, the entropy of the gas increases and the quantities Q(η) are not conserved. Interestingly, this effect was studied in Ref. [42] for the classical hard-rod gas [28], and was marked as the consequence of an interim chaotic regime.
In an attempt to understand better coarse-graining effects, we have discovered that under certain hydrodynamic conditions, the GHD equations are invariant under the coarse-graining procedure where quasi-particle phase-space densities are replaced by their local averagessee Appendix IV. In particular, GHD time evolution commutes with this coarse-graining procedure, as long as the coarse-gaining cells are mesoscopic and certain dephasing conditions hold. This suggests an amount of universality to the GHD solutions; for instance, as fine structures emerge, one may coarse-grain, still correctly describing the solution over time.
We have not been able to study coarse-graining effects directly in the 1d Bose gas (1), as it would require a full quantum simulation which is currently beyond reach. However, the loss of fine structures at large times is made clear as we vary the UV cutoff in our molecular simulation (Appendix IV). The UV cutoff can be taken as the number of classical particles used in the simulation, and this provides an indication for the relation between loss of fine structures and actual particle numbers in the 1d Bose gas. We also observed that at a coarse-grain level, a stationary state emerges -see Appendix V.
Although not a Gibbs state, as argued in [24] such a stationary state must be a "Generalized HydroStatics" solution: it satisfies ∂ x [v eff ρ p ] = (∂ x V /m)∂ θ ρ p (this was explicitly checked in the hard-rod case in [42]). We furthermore have partial numerical evidence that such a stationary state is a universal property of coarse-grained GHD, and not sensitive to particular microscopic realizations of GHD.

Profiles after Time-Of-Flight (TOF) in 1d
In Ref. [1], it is not the in situ density that is being measured, but the density profile after a trap release in 1d: in order to have a cloud that is large enough compared to the resolution of the camera, the longitudinal potential is suddenly released, and the atomic cloud expands for a time t TOF , before a picture is taken. When t TOF is sufficiently large, the real-space density profile n TOF (x) is directly related to the momentum distribution function (MDF) of integrable quasi-particles ρ p (θ , x ) just before the expansion, as (see Refs. [43][44][45][46][47][48] and Appendix VI). However, importantly, it is not only the asymptotic distribution for t TOF → ∞ that is accessible thanks to GHD. Instead, the full expansion of the cloud can be simulated by suddenly setting V = 0 in Eq. (3) at the expansion time. Since, during the whole expansion, the typical length scale of density variations grows proportionally to the typical inter-particle distance, the hydrodynamic approximation remains valid throughout. Hence, GHD is particularly well-suited for predicting the outcome of such measurements, even at finite Time-Of-Flights. To our knowledge, this is not accessible by other techniques. We compare the profiles for t TOF = 10ms (red), t TOF = 25ms (magenta), t TOF = 50ms (cyan) and t TOF = ∞ (dashed black, corresponding to Eq. (7)).
We have evaluated the expansion curves obtained after 10ms, 25ms and 50ms, as released from an evolution within the anharmonic potential for times from 0 to 2τ, and at time 10τ, see Fig. 3. The time necessary to reach the asymptotic state (7) depends strongly on the initial distribution. It is very short for distributions with high momenta such as just after the Bragg pulses (less than 5ms for an expansion at t = 0τ), but it is much longer when particles are slower and mostly lie at the edges of the potential (requiring more than 60 ms at t = 0.25τ). This corresponds to an expanded cloud that is 10 (for initially fast particles) to 30 (for initially slow particles) times larger than the original one.

3d expansion and momentum distribution
Longitudinal expansions are in strong contrast with three-dimensional trap releases. In the latter case, the profile after a long TOF is given by the MDF of the real physical bosons, as opposed to that of quasi-particles. It is also possible to calculate the bosonic MDF, but this is more difficult. An approximation scheme combining GHD data with the ABACUS algorithm is proposed in Appendix VII, where we evaluate the bosonic MDF at 10τ in the anharmonic case, showing that it is significantly different from the quasi-particle MDF.

Integrability breaking: GHD and a collision term
In the quantum Newton's cradle experiment, there are numerous sources of integrability breaking. These include heating and losses, the presence of a trap, interactions between neighbouring tubes as well as hopping events between transverse levels of a given tube [49][50][51].
We can estimate the timescales for these processes in a variety of ways. For integrability breaking due to a trap, an extensive analysis has been provided in [42] for the hard rod gas.
Here, for simplicity, we provide an estimate for the timescale using the notion of energy loss due to a "quantum jump" in our MD simulations. We elaborate on this in Appendix VIII A.
To estimate the timescales associated with other processes, we exploit the fact that GHD can readily take into account at least some integrability breaking processes through the addition of a collision term to the GHD equations. With such a collision term, the GHD equations, in terms of the space-time dependent filling fraction n(θ ) (the GHD normal coordinates) [18,19], read In lowest order perturbation theory, f collision can be computed in terms of the exact matrix elements of the Lieb-Liniger model, for instance using the density matrix elements computed in [52]. As a simple demonstration, we sketch in Appendix VIII B f collision due to intertube interactions. We show that for the parameters governing the original quantum Newton's cradle experiment [1], f collision leads to changes in the energy at roughly the same rate as our estimate for the timescale of integrability breaking due to the trap.

Conclusion
The recently discovered Generalized HydroDynamics (GHD) is an ideal tool in order to provide a full account of the quantum Newton cradle experiment within the Lieb-Liniger model at the Euler scale, fully in the interacting regime. We have observed non-thermalization at cloud collisions, many-body elongation of the oscillation period, and many-body and single-body dephasing. We have also shown that GHD is the ideal tool to study trap expansions.

A Double-well initial state
As a warm-up exercise, we have analyzed a setup that is a rather crude approximation of the QNC experiment [1], see Fig. 4. The initial state is the zero-temperature ground state of (1) in a double-well potential, which splits the gas into two well separated clouds. It is then evolved within a single harmonic well. The main phenomenon -the lack of thermalization upon cloud collisions -of the QNC experiment is observed. We compare both the zero-entropy GHD and molecular dynamics, with excellent agreement. In this setup, two important aspects of the QNC experiment are overlooked: the initial state is not at zero temperature, and the sequence of Bragg pulses produces an initial state with two sets of particles that are separated in momentum space rather than in real space. In the main text, we develop a more realistic approximation. The initial state is not a zero-entropy state, and we rely on the molecular dynamics simulation to solve the GHD equation. We take N = 250 particles, and (in this paragraph only) we work in units where ħ h = m = 1 and g = 2. We take the initial state as the zero-temperature ground state of (1) in a doublewell potential V (x) = V init (x) = 20 ((x/100) 4 − (x/100) 2 ), which splits the gas into two well separated clouds, each containing 125 particles. We construct the corresponding initial density of quasi-particles ρ p (θ , x, 0) by searching for the zero-temperature hydrostatic solution of (3), equivalent to a local density approximation (LDA) [37]. Then, at time t > 0, the double-well is switched off and replaced by a harmonic trap V (x) = mω 2 x 2 /2 = x 2 /800. The GHD equation (3) is integrated using the second method -zero-entropy GHD-and third method -molecular dynamics-above. The results are presented in Fig. 4. The two methods are compared with perfect match: this is compelling evidence for the reliability and robustness of both methods for solving the GHD equation in a trap. Paralleling the QNC experiment, we see that the two clouds initially propagate towards each other due to the harmonic potential, collide, and emerge to continue their quasi-harmonic motion. The interaction deforms the clouds substantially, but thermalization does not occur.

B Many-body dephasing and spreading in phase space
A free particle moving in a potential V (x) preserves, at all times, its total energy θ 2 /2 + V (x), where θ is its velocity (and here and below we take the particles' mass to be unity). Because of interactions, the particles of the Lieb-Liniger Bose gas of course do not preserve this energy. A related question is whether the quasi-particles of the GHD description of the gas do conserve it. We derive here explicitly the fact that the number N (E) of particles within the energy region θ 2 /2 + V (x) < E is not conserved (for any finite E), unless the effective velocity v eff (θ ) equals the bare velocity θ . Non-conservation is thus an effect of the interaction, and can be interpreted as a many-body dephasing effect. Conservation happens when the effective interaction is very weak: either in the Tonks-Girardeau limit, or in the free boson limit. In particular, the accuracy of the conservation of N (E) -the distance between effective velocity and bare velocity -is a nontrivial (and nonlinear) function of the particle density. This helps explain the difference between the strength of the many-body dephasing effect in the harmonic and anharmonic cases, as explained in the main text. Using Performing integration by part, the first term on the right-hand side gives and the second term Clearly, ∂ x θ ± (x) = ∓V (x)/θ ± . Thus we find Thus the change of N (E) is bounded by which is controlled by the relative difference between v eff (θ , x) and θ at the region's boundaries θ ± (x).
A similar calculation gives the change of the energy within this region, x) . (14) Further, by using the defining integral equation for the effective velocity, one has Therefore, (16) Again we see the same velocity difference playing an important role.
It is worth mentioning that for E → ∞, we have θ ± (x) → ±∞, and recall that v eff (θ , x) → θ as θ → ±∞. In this limit it is clear that both N (∞) and E(∞) are invariant, as they should as the system preserves the total number of particles and the total energy. It is also clear that in the free case, where v eff (θ , x) = θ , both quantities are preserved for all Es.

C A continuous family of conserved quantities for the GHD equation in a trap
Given a quasi-particle distribution function ρ p (θ ), one defines the occupation number n(θ ) = 2πρ p (θ )/1 dr (θ ), where the dressing is defined as in the main text. The occupation number n(θ ) is always between 0 and 1. As per the theory of GHD, this satisfies ∂ t n + v eff ∂ x n − (∂ x V /m)∂ θ n = 0. We find that, under GHD evolution in a trap -see Eq.
(3) in the main text- is a conserved quantity for any function f , as long as the quasi-particle density ρ p (θ , x, t) does not have discontinuities in θ or x. To see this, notice that Upon integrating over d x dθ , this gives zero using Stokes theorem, assuming zero quasiparticle density at infinity. Thus, ∂ t Q[ f ] = 0. In fact, this can be understood as the generalization of the fact that the total Yang-Yang entropy is conserved in perfect fluids, hence in GHD. Namely, the total Yang-Yang entropy is so it can be recast into the form of Q[ f ] with f (n) = 1 n (−(1 − n) log(1 − n) − n log n). The precise structure of the infinite number of dynamical symmetries generated by these conserved quantities is still not known, and deserves further studies.
Finally, we note that a particularly convenient choice of basis for these conserved quantities corresponds to the choice f (.) = δ(.−η) for η ∈ [0, 1], thus leading to the conserved quantities mentioned in the main text. Q(η)dη is the number of quasi-particles whose local occupation number lies between η and η + dη.  Figure 5: Diagram illustrating the fact that coarse-graining commutes with GHD time evolution. Starting from the state at t = 8τ in the anharmonic case (same data as in Fig. 2 in the main text), we first coarse-grain the system (i.e. we chose larger bins, and re-sample the phase-space distribution ρ p (θ , x) in our molecular simulator accordingly) and let it evolve during a time ∆t = 2τ. We compare the resulting distribution with the one obtained from evolving the system first, and then coarsegraining it. The two resulting distribution are almost identical.

D.1 Invariance of GHD under coarse-graining
We consider the equation where a is the acceleration. The proof below is specialized to the simple case of the Lieb-Liniger model, but it is straightforward to extend it to the general context of GHD. Consider coarse graining GHD, with coarse cells of area × in phase space. That is, denote and letρ dγd y ρ p (γ, y). (22) We make the following assumptions: (a) the acceleration is essentially constant on the scale , and (b) the velocity θ and (c) the differential scattering phase ϕ(θ ) are essentially constant on the scale . We also assume one of the following: either (d) the rapidity integral of quasi-particle densities and currents, on scale , are essentially constant on scale in the position variable; or (e) cells are uncorrelated, From the evolution equation, dγd y a( y)ρ p (γ, y) .

Assuming (a) we have
The above can also be written in integral form, so that the derivation holds in the space of weak solutions as well. Now definē v eff (θ , x) = C(θ ,x) dγd y v eff (γ, y)ρ p (γ, y) Then clearly Thus We now derive the integral equation forv eff (θ , x) that shows that it is determined bȳ Also from (c) we find in general This gives, for the second terms in (29), and similarly for the third term. Now, on the one hand, assuming (d), we have and, similarly, Putting these together, we find On the other hand, assuming (e) we directly obtain (33) and (34), and the result (35) follows again.
That is, we conclude thatv eff (θ , x) is the effective velocity associated to the coarse-grained densityv eff (θ , x) = v eff [ρ p ](θ , x). This is the GGE equation of state leading to GHD, and thus the coarse-grained equation (28) is GHD again.

D.2 Numerical analysis
Numerical simulations have been performed with the molecular dynamics simulator proposed in [35] (using a standard desktop computer, 3.8GHz, quad core). We have performed simulations with the exact parameters described in the main text (giving approx. 350 particles), Figure 6: Comparison of the density ρ p (θ , x) obtained from three different microscopic realizations of the same GHD equation (anharmonic case): the lower one ("less refined") is our molecular simulator with 350 classical particles, the middle one ("more refined") is with 32 × 350 classical particles, and the top one ("more and more refined") is with 256 × 350 classical particles. The latter is the one shown in the main text in Figs. 1 and 2. The corresponding density profiles are plotted below (green: 350 class. part., red: 32 × 350 class. part., blue: 256 × 350 class. part.). As the discretization is refined, the density ρ p (θ , x) converges to the true GHD solution. A given discretization with a finite UV cutoff cannot resolve the fine structures that appear in GHD at scales smaller than the cutoff. However, GHD remains valid at larger scales, and this then amounts to coarse-graining. as well as after rescaling all lengths by factors of 2 n for n = 1, . . . , 8. Since GHD is manifestly invariant under scaling of lengths, these represent different choices of microscopy, with different numerical precision for the solution to the GHD equations. The equivalent of a sampling of 2000 has been used (that is, approx. 2000/2 n samples). In the harmonic case, this was observed to give a noise level (as calculated by the relative L 1 distance between two equivalent sampling) of the order of 5% throughout the evolution, on spectral densities binned on a 70 × 70 lattice covering the range of Figs. 1 and 2 (main text). In this case, we have found it sufficient to take n = 4 (approx. 5600 particles): we observed a Yang-Yang entropy production of approx. 6% over 10τ, and no significant change of the GHD-conserved function Q(η) (defined in Section III).
The anharmonic case is much more delicate, and we have performed a more detailed numerical analysis. We show here results n = 0 (approx. 350 classical particles), n = 5 (approx. 11000 particles), and n = 8 (approx. 90000 particles; this with about half the sampling, using 4 samples only). The latter provides the results presented in the main text, see Fig. 2, reproduced (in the anharmonic case) in Fig. 6 for convenience, where we show the results at t = 0, 2τ, 4τ, 6τ, 8τ, 10τ. It is apparent that agreement between the choices n = 5 and n = 8 is relatively good, although small-scale structures are more clearly discerned in the latter case. We have observed a Yang-Yang entropy production of approx. 20% over the evolution time of 10τ for n = 8. More striking, however, is the analysis of the conserved quantity Q(η). For the less refined discretization (350 classical particles), we see that it is conserved up to 4τ − 7τ, then fine structures develop and are progressively erased by coarse (a) 350 classical particles (b) 32×350 classical particles (c) 256×350 classical particles Figure 7: Time-evolution of the quantities Q(η), that are exactly conserved in pure GHD, but not in microscopic realizations of GHD. We see that, as the discretization is refined, and more and more fine structures of GHD are probed by the microscopic model, the quantities are conserved on longer times. Differences of thickness between the various cases are due to noise level differences.
graining, and after that Q(η) is conserved again. For more refined discretizations, Q(η) is conserved for a longer time, see Fig. 7. The distribution at 10τ is relatively stable under change of the microscopy, as long as the number of particles is such that the corresponding coarsegraining is fine enough, in phase-space, for variations of the potential and scattering length to be small from cell to cell, yet large enough so that each cell contains a large number of particles (this happens to good approximation for n ≥ 6 on a binning of 70 × 70, for instance). This lay support to the idea that coarse-grained GHD leads to the large-scale evolution, independently from the microscopy.

E Numerically obtained stationary state in the anharmonic case, and evidence that it is not thermal
We have investigated the stationary state obtained at large time. The setup is the same as in the main-text, however in order to speed up the many-body dephasing, we take a slightly stronger anharmonicity. We take the confining potential as V (x) = (1+4(x/ ) 2 ) mω 2 x 2 2 with = 120µm and ω = 0.314ms −1 . We call τ = 2π ω . We observe that, at times than t > 12τ, the distribution ρ p (θ , x) looks stationary, see Fig. 8. We have compared this stationary distribution to the thermal distribution that has the same particle number and the same total energy. The two distributions are obviously different, as can be seen in the last plot of Fig. 8. Figure 8: Comparison of the "stationary state" obtained numerically at t = 20τ, and of the thermal distribution with the same number of particles and the same total energy.

F Trap release in 1d, and measurement of the momentum distribution of quasi-particles
Assume that we have a 1d Bose gas described, at some given time t (which we fix to zero in this appendix), by a distribution of quasi-particles ρ p (x, θ ), as in GHD. Assume that this density has a support that is contained in [−∆x, ∆x] × [−∆θ , ∆θ ] so that no particle is outside the box [−∆x, ∆x], and no particle has a velocity larger than ∆θ . Then we release the longitudinal confinement, and let the gas expand in 1d. In this appendix, we are going to derive the following result: for a sufficiently long time of flight T , the spatial density of bosons n(X , T ) = Ψ † (X )Ψ(X ) is given by the momentum distribution function of the quasi-particles before the release, n(θ ) = d xρ p (x, θ ): This result has long been known for the Tonks-Girardeau gas, where it has sometimes been dubbed "dynamical fermionization" [44,45]. For the interacting case, it seems to have been pointed out only recently [46,47]. Here, for the convenience of the reader, we provide a fully detailed proof of this result. The derivation of formula (36) consists of two steps. The first step is to note that there must exist a time T 1 that is large enough such that the local density n(x) is sufficiently low everywhere in the system, so that γ(x) = c/n(x) 1 for all x. At that time T 1 , the quasi-particle distribution is some function ρ p,1 (x, θ ), with support in [−∆x 1 , ∆x 1 ] × [−∆θ , ∆θ ] for some ∆x 1 . Since, by construction in the Bethe ansatz method, quasiparticle spectral densities are exactly conserved under quantum evolution, we have namely n(θ ) was conserved during the evolution from t = 0 to t = T 1 . Importantly, on the right-hand side, although the quantity ρ p (x, θ ) is meaningfully defined only if the state is weakly varying in space (so that we can approximate it by a collection of homogeneous fluid cells), its spatial integral makes sense beyond this regime. Indeed, it simply encodes, as a function of θ , the values of all extensive conserved quantities in the inhomogeneous initial state. Of course, in the application considered in the present work, ρ p (x, θ ) is obtained after time evolution within an inhomogeneous external potential using the hydrodynamic approximation, and thus the values of all extensive conserved quantities it encodes are likewise subject to the hydrodynamic accuracy. The second step goes as follows. At times T > T 1 , because γ is uniformly very large, the dynamics of the gas is captured by the Tonks-Girardeau hamiltonian, where Thus, we are back to the case of the Tonks-Girardeau gas, and we simply apply the results of [44,45]. For completeness, here we give a fully detailed calculation that leads to the wanted result.
The density of bosons at point X and time T is In the last line, we used the fact that is nothing but the Wigner function, so it is exactly the number of fermions at position x with momentum mθ , therefore it has to be equal to 2π m ρ p,1 (x, θ ).
[The factor 2π m simply comes from a difference in normalization convention between ρ p,1 for the Tonks-Girardeau gas and the Wigner function: for instance, the total number of particles is d , this gives (with the change of variables K = k+k 2 , q = k − k ): The starting point for the calculation of the bosonic momentum distribution function in the anharmonic trap is the full spatial density profile of the gas at time t = 10τ obtained from GHD and plotted in the left panel of Fig. 9. This density is subsequently divided into three separate regions, as illustrated again in the left panel of Fig. 9. For each region, a root density ρ i (θ ) is extracted. These are plotted on the right panel of Fig. 9. They satisfy the sum rule

G Computing the bosonic MDF
These distributions clearly show that the gas in each of the three regions is significantly excited away from the ground state. The next step is to compute the bosonic MDF separately on each of these three individual constituent representative states. To do this, we rescale the ρ i (θ )'s viã so that we are working at unit density, i.e. dθρ i (θ ) = 1.
We then use ABACUS [13][14][15][16][17] to compute the MDF on each representative state. This involves the following steps: starting from each individualρ i (θ ), a best-fitting discretized Bethe state |i〉 N is constructed at a chosen particle number N ABAC US (setting this equal to L ABAC US to stay at unit filling) by choosing a set of quantum numbers generated from the state's counting function, namely: adding a rapidity whenever L λ −∞ dλ ρ i (λ ) crosses a half-odd integer, and setting the quantum numbers to those giving the closest-matching set of rapidities; N is chosen even, and as large as practically possible. ABACUS is then run for the one-body correlation N 〈i|ψ † (x, t)ψ(0, 0)|i〉 N . The quality of the result is quantified by the saturation of the integrated intensity sum rule. On such highly-excited states, a large number of intermediate states must be summed over (for N ABAC US = 32, these were 74307322 (i = 1), 101334549 (i = 2) and 87195380 (i = 3), yielding saturations of 0.964412, 0.933457 and 0.979004).
Having these three MDFs (ñ i (k), i = 1, 2, 3), we now rescale back to MDFs (n i (k), i = 1, 2, 3) corresponding to the original three regions characterized by N i , L i . The relevant relation here is Here c i is a constant that can be determined by insisting that The results are displayed in Fig. 10. We then average over the three n i (k) to obtain what would be the MDF measured in the actual experiment. If we compare the r.h.s. of Fig.9 and Fig. 10, we see the bosonic MDF and ρ(θ ) are considerably different. However both have a double humped feature characteristic of post-Bragg pulse states.

H Timescales for Integrability Breaking
We consider in this section estimates for timescales due to integrability breaking. We do this in two parts. In the first part we consider the timescale associated with integrability breaking due to the trap. And in the second part, we consider the derivation of a GHD collision term coming from intertube interactions.

H.1 Estimate of Integrability Breaking due to Trap
We have argued in the main text that integrability breaking due to the trapping potential is small. We provide here an argument for this 2 .
One way to parameterize the integrability breaking is to consider how the potential energy of a particle changes when two particles collide in the presence of the trap. We can estimate this energy by using the molecular dynamics representation of the GHD equations. In the molecular dynamics, when two particles collide, they experience a quantum jump, ∆x, where the particles are displaced according to their relative momentum, p 1 − p 2 , and the strength of interactions: The quantum displacement gives us a scale for integrability breaking because upon displacement in the molecular dynamics simulation, the potential energy of a particle experiencing a quantum displacement changes by an amount (so violating energy conservation): Rather than trying to estimate this directly (it is difficult to provide even a back of the envelope estimate of this quantity as it requires accounting for both inter-and intra-cloud collisions), we read off the change of energy directly from our numerical simulations. We find it to be approximately 0.1% for each oscillation of the clouds. This rate is smaller than what we estimate in the next section for the intertube interactions present in the QNC experiment.

H.2 Collision Term due to Density-Density Couplings Between Tubes
In the main text we have discussed the possibility of adding a collision term to the GHD equations. Here we elaborate on how to compute the collision term and from it provide an estimate of the time scale for integrability breaking.
For this exercise, we are going to consider a system composed of two Lieb-Liniger models that are coupled by a density-density interaction. As we will be deriving the collision term in lowest order perturbation theory, the effects of having more tubes coupling to one another, either because the tubes are in an array of a given coordination number (as is typical) or because of long range dipolar forces (as in Ref. [49]) are additive. The case of two coupled Lieb-Liniger models is then sufficiently general.
The Hamiltonian we will then consider where ρ i,k , the Fourier component of the density operator in the i-th tube, is defined as where ψ i,q is the q-th Fourier component of the i-th chain field operator.
To compute the collision term, i.e. the rate of change of the quantum numbers I r in a state, we imagine that we have an initial state |i〉, characterized by a set of occupied quantum numbers These integers I r,s are the the Bethe integers for which one solves the Bethe ansatz equations describing the uncoupled chains. Now let n s (I r,s ) be the occupation of quantum number I r,s on chain s. By the Fermi golden rule, the rate of change of n s (I r,s ) from its initial value is where we are summing over all final states, f , and To develop this expression forṅ s (I r,s ) further, we write the states |i〉, | f 〉 explicitly as a product state of states belonging to the two chains: We, for simplicity, will take the initial states on the two chains to be equal, i.e. |i 1 〉 = |i 2 〉. We will also only consider the first set of final states that can lead to thermalization. Such states involve 2-particle-hole excitations on one chain, and 1-particle hole excitation on the other, i.e. | f s=1,2 〉 are given by | f 1 〉 = |i 1 ,ĥ 1 ,ĥ 1 , p 1 , p 1 〉, | f 2 〉 = |i 2 ,ĥ 2 , p 2 〉, where here a state |i s ,ĥ s ,ĥ s · · · , p s , p s , · · · 〉 (54) is defined as the state |i s 〉 on chain s with quantum numbers (holes) h s , h s , · · · removed and quantum numbers (particles) p s , p s , · · · added. States involving only 1-particle-hole excitation on each chain can lead to equilibriation between the chains, but will not thermalize non-Gibbsian distributions. As such, these final states will not be considered.
Here the first term in the above corresponds to the case where the 2-particle-hole excitation takes place on chain 1 and the 1-particle-hole excitation on chain 2 while the second term exchanges the chains where these two processes occur. We are not going to evaluate this expression in detail as the 2-particle hole matrix elements F p 1 ,p 1 ,h 1 ,h 1 are highly non-trivial. We can however provide an estimate of the time scale. We know that matrix elements F p 1 ,p 1 ,h 1 ,h 1 scale as 1/c (and so in the c = ∞ limit this process would be suppressed and thermalization would not occur). In general, density matrix elements involving n-particles and n-holes scale as 1/c n−1 . The contribution then from the above sum goes as N 3 /c 2 (N is the number of particles) with the resulṫ where here m is the mass of rubidium atom. We can write the density-density coefficient A coupling the tubes together as γ int er tube ρ 0 where ρ 0 is the background density in the tubes.
We know in the context of Ref.
[1] that γ int er tube γ int r atube . If, for example, γ int r atube ∼ 10 −2 , we see that the energy change per oscillation of the gas is on the order of 1%, similar to our estimates of the energy change due to integrability breaking arising from the trap. However this rough estimate shows that if γ int r atube ∼ γ int er tube as in Ref. [49], the fractional change in energy per oscillation cycle will be O(1).