Relaxation of Bosons in One Dimension and the Onset of Dimensional Crossover

We study ultra-cold bosons out of equilibrium in a one-dimensional (1D) setting and probe the breaking of integrability and the resulting relaxation at the onset of the crossover from one to three dimensions. In a quantum Newton’s cradle type experiment, we excite the atoms to oscillate and collide in an array of 1D tubes and observe the evolution for up to 4.8 seconds (400 oscillations) with minimal heating and loss. By investigating the dynamics of the longitudinal momentum distribution function and the transverse excitation, we observe and quantify a two-stage relaxation process. In the initial stage single-body dephasing reduces the 1D densities, thus rapidly drives the 1D gas out of the quantum degenerate regime. The momentum distribution function asymptotically approaches the distribution of rapidities, which are conserved in an integrable system. In the subsequent long time evolution, the 1D gas slowly relaxes towards thermal equilibrium through the collisions with transversely excited atoms. Moreover, we tune the dynamics in the dimensional crossover by initializing the evolution with different imprinted longitudinal momenta (energies). The dynamical evolution towards the relaxed state is quantitatively described by a semiclassical molecular dynamics simulation. 1 ar X iv :1 80 4. 01 96 9v 2 [ co nd -m at .q ua nt -g as ] 2 0 Ju l 2 02 0 SciPost Physics Submission


Introduction
The study of relaxation, thermalization and equilibration in an isolated many-body quantum system [1][2][3][4] has a long history starting with von Neumann [5]. An integrable system will not fully thermalize [6,7], but dephase towards a generalized Gibbs ensemble (GGE) [8][9][10], reflecting its many conserved quantities. Bosons in one dimension (1D) [11,12] are a model system to study these fundamental questions at the interface between microscopic quantum evolution and statistical physics. If the interactions are short range, the physics is described by the integrable Lieb-Liniger model [13][14][15].
However, in physical realization in the laboratory, "integrability" is never perfectly maintained. If the system is only approximately integrable, the dephased state is pre-thermal and is expected to reach thermal equilibrium at a much later time [16,17]. The integrability of identical bosons in 1D with short-range interactions can be broken by numerous effects [18][19][20][21][22][23][24][25][26][27]. Moreover, in a real experimental implementation, the strict 1D condition is maintained only within a limited time frame. For example, inevitable imperfections and noise lead to heating, which drives the system into the dimensional crossover regime, breaking integrability to some extent [28,29].
The pioneering work of Kinoshita et al. [6] introduced the Newton's cradle into the microscopic world as a typical near-integrable model. Since then it has attracted broad interests in the quantum gas community [24][25][26][27]. In this "quantum Newton's cradle", opposite longitudinal momenta are imprinted on the atomic ensembles, which then oscillate in the 1D-trap and collide. If the imprinted momenta are much smaller than required to excite transverse excitations, the oscillations persist for many atomic collisions. In Kinoshita et al. [6] the 1D condition is guaranteed by applying a very shallow longitudinal trap, so that high energy particles can leave the 1D traps at both ends. In return, the system exhibits significant loss [30][31][32].
In our work, we revisit the quantum Newton's cradle setting and study the dynamics of the longitudinal momentum distribution function (MDF) and the transverse excitation at the onset of 1D-3D crossover [18,33,34]. Opposite to Ref. [6], we retain nearly all the atoms during the dynamical evolution, thus significantly extend the observation time up to 400 oscillation periods. The non-equilibrium dynamics occur in two stages, dominated by different mechanisms. At a short time scale, the single-body dephasing reduces the 1D density, thus rapidly drives the 1D system out of the degenerate regime and asymptotically approaches the non-degenerate limit. During this time, the oscillation period averaged distribution of quasimomenta (also referred to as rapidities) is nearly conserved, while the MDF deforms significantly and finally in the non-degenerate limit tends to coincide with the quasimomentum distribution. At longer times, the system smoothly evolves into a dimensional crossover from 1D to 3D with a small fraction of the atoms transversely excited. The system finally relaxes to an equilibrium with a Gaussian momentum distribution suggesting a thermal final state.
Within the non-degenerate regime, accurate predictions for the system can be obtained very efficiently using a simple molecular dynamics (MD) calculation [35][36][37][38]. MD is a wellestablished numerical tool of simulation of many-body dynamics in various fields such as physics of fluids or chemical physics. In MD, the center-of-mass dynamics of particles (molecules) are classical, subject to the Newtonian mechanics. The internal degrees of freedom can be discretized (quantized). The main prerequisite to use the MD is therefore the nondegeneracy of the system, which is necessary to allow for the classical description of centerof-mass trajectories of individual molecules in a simulation. Alternatively, one could describe the system using the recently developed theory of generalized hydrodynamics (GHD), which has been proposed to describe dynamics of 1D quantum gases at or close to the integrable point [27,[39][40][41][42][43]. Unlike other methods, GHD takes into account interactions between particles and remains across all phases of the Lieb-Liniger model. Recently, we have extended the applicability of GHD to the dimensional crossover regime by including transversely excited states [44]. Although that approach displayed very promising results on describing realistic systems in short to intermediate time scales, employing the full machinery of GHD to describe a primarily non-degenerate gas seems disproportionate. Therefore, within the present paper's scope, we opt for the much simpler MD simulations, whose efficiency allows us to account for the whole spectrum of transverse excitations in an easy and efficient way and makes it more suitable for describing the intermediate to long-term dynamics. The simulation results show an excellent quantitative agreement with the experimental observations. The paper is organized as follows. In Sec. 2, we introduce the experimental setup, methods, 2d optical lattice Bragg pulses imaging beam (for measuring MDF) imaging beam (for band mapping) TOF Figure 1: Experimental setup. The 1D Bose gas is prepared by adiabatically loading a BEC into a 2D optical lattice (green arrows). A sequence of Bragg pulses (blue arrows) is applied to excite the atoms to oscillate and collide in the 1D traps. After some duration, the atoms are released from the 1D traps and detected after TOF (yellow arrows), providing us the longitudinal MDF (horizontal imaging) and the fractional population of atoms in transverse states (vertical imaging). and conditions. In Sec. 3, we present the experimental observations of MDF. In Sec. 4, we evaluate the dephasing process and observe that it is independent of the 1D density. In Sec. 5, we study the relaxation of MDF under the near-integrable condition, as the 1D systems evolve from the degenerate to non-degenerate regimes. In Sec. 6-7, we introduce the MD simulation and study the dynamical process towards thermalization at the onset of 1D-3D crossover. We compare the predictions of our simulation to experimental measurements in both longitudinal and transverse degrees of freedom. In Sec. 8, we further explore the dynamics initialized with a significant number of atoms having larger collision energies than the threshold of transverse excitation. In Sec. 9, we discuss the effect of virtual excitations as an alternative mechanism of breaking integrability in 1D. In Sec. 10, we draw the conclusions.

Preparation and characterization of 1D gases
We start our experiment by preparing a 87 Rb BEC in an all-optical trap. The BEC is then adiabatically loaded into a square array of 1D traps formed by a 2D optical lattice created with two retro-reflected perpendicular laser beams (see Fig. 1). The atoms are strongly confined in transverse directions ( y and z) in 1D traps with a trap frequency of ω ⊥ /2π = 31.0(3) kHz and weakly confined in the longitudinal direction (x) with ω /2π = 83.3 (8) Hz (see Appendix A.1 for more details). The prepared 1D gases are characterized by the following key parameters: Atom number. In the experiment, the atom number per tube varies across the 1D trap array. To balance the contributions from each tube, we calculate the weighted-average atom number.
The effect of the lattice loading procedure on the atom-number is considered by following the method proposed in Ref. [45]. The number of atoms per tube strongly depends on the total atom number, N t ot , of the BEC. By tuning N t ot between 1 × 10 4 and 1 × 10 5 , the atoms are distributed in 400 to 1200 tubes with a weighted-average atom number per tube, N , ranging from 40 to 130.
Interaction strength. The interaction strength is characterized by the dimensionless Lieb-Liniger parameter γ = c/n 1D (we use the standard notation for the Lieb-Liniger model: c = mg 1D /ħ h 2 and n 1D is the mean 1D density) [12][13][14]46]. To account for the inhomogeneous distribution of atoms over the tubes, we calculate the weighted average for γ. At the beginning of the evolution, γ 0 = 1.5 and 0.7 for the systems with N = 40 and 130, respectively. As the dephasing proceeds, the atoms spread along the tubes and the 1D density n 1D drops, thus increasing γ. When the system is fully dephased, γ ∼ 6 for N = 40, and γ ∼ 2.5 for N = 130.
Temperature and degeneracy. The temperature of the 1D gases is evaluated with the halfwidth at half-maximum (HWHM) of the Lorentzian-like longitudinal MDF before the Bragg pulses are employed. The inhomogeneous density profile is considered via the local-density approximation [47][48][49][50], and the significantly strong interaction is taken into account in a quantum Monte Carlo calculation [51]. We obtain the temperatures T = 34 nK ( T = 1.6) and T = 94 nK ( T = 4.4) for N = 40 and N = 130, respectively. The reduced temperature T = 2ħ h 2 k B T /(mg 2 1D ), together with γ, characterize the degeneracy of 1D gases. For T 1, the crossover from degenerate to non-degenerate regime occurs at γ T −1/2 [52]. As such, our Newton's cradle experiments start in the intermediate regime of degeneracy and approach the non-degenerate limit as the dephasing occurs.
Heating rate. Heating is evaluated by holding a BEC in the identical trapping potential without a Bragg-pulse excitation. Within the time scale of our experiment (4.8 s), we observe excitation of about 3% of the atoms to the first transversely excited state and hardly any to the second state (see Appendix A.2).
Tunneling between 1D tubes. Under our experimental conditions, the residual single-atom tunnel coupling (J < 0.01 s −1 ) is too small to influence the longitudinal dynamics. The Josephson frequency in our system is ω J = 2J(2J + 2µ/ħ h) ħ hk 2 L /m, which holds even for higher transversely excited states n r = 1 or 2 which show an increased J. Here µ is the chemical potential. As such, we can safely neglect the effect of tunneling.

Excitation of longitudinal motion
To start the non-equilibrium dynamics in the longitudinal direction, a sequence of two Bragg pulses using retro-reflected λ = 852 nm light is employed on the 1D gases (Fig. 1). Neutral atoms initially at rest are diffracted to a series of momentum states e 2niħ hk Br ag g x |ψ 0 〉 when exposed to the standing wave light pulses, where k Br ag g = 2π/852 nm is the recoil momentum and |ψ 0 〉 is the atom's (low-momentum) state before scattering [53][54][55][56][57][58]. Following the method described in Ref. [59] we can accurately control the resulting momentum distribution function (MDF) by tuning the pulse intensity V 0 and the time sequence of pulse lengths t 1 , t 3 , and separation t 2 (illustrated in Fig. 2). In this way we prepare a large variety of initial MDFs from where our study starts. The ones used in this paper are presented in Fig. 3    is the threshold momentum, and two head-on colliding atoms with momenta ±k th have the required collision energy (E th = 2ħ hω ⊥ ) for populating the transversely excited states [18]. The momenta ±2k Br a g g correspond to ∼ 40% of the excitation energy E th , while ±4k Br ag g and ±6k Br a g g are beyond the excitation threshold. In Sec. 3-7, we will first investigate the dynamics initialized with the distributions shown in Fig. 3 (a) and (b), where hardly any atom has the energy to be transversely excited at the beginning of dynamical evolution 1 . In Sec. 8, we will further show the results starting with higher-energy distributions shown in Fig. 3 (c)-(h). 1 The quantities that characterize the dynamics during a collision process are actually the quasimomenta rather than the bosonic momenta of two particles. Owing to the interaction, the quasimomentum distribution is in general broader than the MDF. Thus, we expect < 0.1% and ∼ 2% of atoms obtain the quasimomenta beyond the excitation threshold for N = 40 and N = 130, respectively. These values are estimated according to the method described in Sec. 6. Figure 4: Oscillations in the first period (time resolution of measurements ∆t = 0.2 ms). The distributions are shown in log scale with the maximum normalized to unity. The oscillation periods for ±2k Br ag g , ±4k Br ag g , ±6k Br ag g are measured to be 12 ms, 13.2 ms, 15.2 ms, respectively.

Post-pulse evolution and detection
After the Bragg pulses, we keep the lattice on for some duration t, during which the atoms oscillate in the tubes with a period of T = 12 ms and collide with each other. At the end of evolution, the atoms are released from the optical lattice and expand in 3D for a long time-of-flight (TOF). The image taken in the horizontal plane provides us a density profile given by the MDF in the longitudinal direction of 1D gases. Alternatively, we apply the imaging beam vertically and detect the population of atoms in transverse states (see Fig. 1). See Appendix A.3-A.4 for more technical details on the methods of detection and data analysis.

Basic observations: two-stage relaxation process
To study the non-equilibrium dynamics after the Bragg pulses, we explore the MDF f (t, k) for an evolution time up to 4.8 s (400 oscillation periods). The early stage of evolution is dominated by the oscillation of the kicked momentum components in the longitudinal trap, as shown in Fig. 4 for strong excitation pulses ( Fig. 3 (g)). The momentum peaks ±2k Br ag g , ±4k Br a g g , and ±6k Br a g g exhibit different oscillation periods, which stem from the anharmonicity of the longitudinal confinement in the tubes, caused by the Gaussian profile of the lattice beams. The measurements are consistent with the expectations based on our trap parameters within reasonable experimental imperfections, for example, the nonideal lattice beam quality, the imperfect overlap between lattice beams, etc. Fig. 5 shows the dynamics of MDF for the case of excitation with only ±2k Br ag g momentum components ( Fig. 3 (a) and (b)). From these measurements, one clearly observe two distinct stages: Stage I, single-body dephasing: With increasing time, the clear oscillations in the early stage become more and more blurry and finally completely dephase after about 60 oscillations (720 ms). Even though the MDFs behave markedly different, the dephasing seems to be happening over a similar time scale in both cases (atom numbers per tube: N = 40 and N = 130). We will quantitatively compare the dephasing processes in these two cases and discuss them in more detail in Sec. 4 (Stage Ia). Accompanied by the dephasing process, the Bragg peaks in the MDF broaden and become rounded over time at a near-integrable point. This deformation is also observed on the oscillation period averaged profile, and it is expected to be induced by the 1D density decrease caused by the dephasing effect. We will investigate this observation in Sec. 5 (Stage Ib). Stage II, relaxation towards thermal equilibrium: At long times (t > 720 ms, 60 periods), the MDF is completely dephased and does not show short time variations within one oscillation period. The MDF further evolves towards a Gaussian distribution. We conjecture that the observed long time relaxation is dominated by the scattering processes mainly between atoms in the transverse ground state and a small population in the transversely excited state. The twobody collisions involving transversely excited atoms allow redistribution of the longitudinal momenta, thus breaking the integrability and constituting the onset of the crossover between 1D and 3D. As we will show in Sec. 6, this conjecture is supported by the excellent agreement between the experimental observations and the results of MD simulations implemented on a semiclassical model using experimentally determined parameters as input.

Stage Ia: single-body dephasing
The simplest and most direct explanation of the observed blurring of the oscillations in Fig. 5 is the single-body dephasing. It does not require interactions and is caused by a diffusion of the relative oscillation phase of each particle. In our experimental setup, the single-body dephasing has two distinct contributions: (i) The anharmonicity of the longitudinal Gaussian confinement within each tube makes the atom with larger kinetic energy oscillate at a slightly larger period as can be observed in Fig. 4. This effect leads to the dephasing of atoms with different energies within each tube. (ii) The inhomogeneity among tubes, which is also a result of the Gaussian profile of the lattice beams, leads to different oscillation frequencies in different tubes and to dephasing of the oscillations among the tubes.
Compared to the inhomogeneity of the traps, we are concerned more about the anharmonicity of longitudinal confinement, which broadens the longitudinal spatial distribution of atoms in each tube. The resulting decrease of atomic density drives the 1D gases out of the degenerate regime. We will demonstrate this effect in Sec. 5.
The dephasing takes effect until the MDF stops varying during one oscillation period. The process can be characterized by the deviation of the MDFs from their oscillation period averaged profile, where F (t, k) is the average profile of MDF over one oscillation period. For the convenience of comparison, the prefactor C D = 2 × 10 7 /N 2 rescales the calculation according to the atom number and lifts the results close to 1. In Fig. 6 we plot D(t) as a function of time for N = 40 and N = 130. D(t) decreases as the variation of MDF in one period decays over time. The dephasing we observe is independent of the atomic density. In both cases, the MDFs are fully dephased at around 720 ms (marked with the vertical dashed line), even though the dephased distributions are distinct due to the different relaxation rates (see Sec. 6). Beyond 720 ms D(t) reaches a plateau, which is dominated by the imaging noise floor. The density-independent dephasing rate can be explained by a collisionless classical calculation with our experimental conditions [25].

Stage Ib: dynamics of the momentum distribution function at a near-integrable point
To study the long time evolution, we average out the effect of single-body dephasing by calculating the oscillation period averaged MDF F (t, k). As the relaxation occurs, F (t, k) evolves from a double-peak profile to a peak-rounded profile and finally approaches a Gaussian distribution. The evolution is characterized by the deviation of F (t, k) from its closest thermal distribution. More specifically, we calculate the summed square of residuals between F (t, k) and its best fit to a Gaussian distributionF (t, k) where C R = 10 4 /N 2 . This approach has emerged as the most robust way to characterize the distance from a relaxed Gaussian (thermal) equilibrium state under our experimental conditions with inevitable noise and fluctuations 2 . In Stage I of evolution, the dynamics are well described by the evolution of the estimated MDF (dotted curves). As the system approaches the non-degenerate limit, the estimated MDF tends to coincide with the QMD (solid curves). In Stage II of evolution, the 1D gases relax towards a Gaussian MDF, which can be interpreted as a thermal equilibrium state. The molecular dynamics simulations agree nicely with the experimental data in Stage II. The vertical dashed line indicates the end of single-body dephasing and demarcates the two stages of relaxation. The shaded area indicates the noise floor.
We clarify that the rapid decay of R(t) in Stage I of evolution stems from the evolution of MDF within the integrable model. In the degenerate limit, the MDF of a 1D Bose gas is described by a Lorentzian-like profile and dominated by the phase fluctuations (see Appendix C or Ref. [49,60]). The MDF is in general much narrower than the distribution of quasimomenta [27,[61][62][63][64][65][66][67], which are the conserved quantities characterizing the integrable many-body system. As a result of the single-body dephasing derived from the anharmonicity of confinement, the 1D densities decrease and asymptotically approach one-fourth to one-third of the initial values. Thus, at the end of Stage I, γ T 1/2 is greater than 7.6 (N = 40) and 5.2 (N = 130) 3 , suggesting that the dephased gases are in the non-degenerate regimes. The change of the interparticle interaction during this process deforms the MDF and brings it close to the quasimomentum distribution. As we observe in our experiments, the narrow momentum peaks present at the beginning of evolution are washed out and become rounded, resulting in the rapid decay of R(t) in the first 720 ms. When the single-body dephasing comes to its end, the decrease of density stops, and the evolution towards the non-degenerate limit is greatly slowed down. The process of the MDF approaching the quasimomentum distribution is illustrated in Fig. 8. Since the description of this figure involves the MD simulation that will be introduced in Sec. 6, we will explain the details later in the relevant paragraphs. Note that the integrability is (nearly) preserved during Stage I, especially for the N = 40 case. The heating process affects the dynamics very little at short to intermediate time scales. Due to the parity conservation, an inelastic collision occurs when at least two atoms appear in the first transversely excited state (or at least one atom in the second excited state) in the same tube. When the atom number per tube is relatively low as the case in our experiments, the effect of heating on the longitudinal motion is postponed. We will demonstrate this statement in more depth in Sec. 6.

Stage II: relaxation towards a Gaussian (thermal) momentum distribution and molecular dynamics simulation
At longer times, the gases relax towards thermal equilibrium through inelastic collisions that involve the transversely excited modes. This process is triggered by the minimal heating effect and the tiny population of atoms with momenta |k| > k th , marking the onset of dimensional crossover. As shown in Fig. 7, R(t) in Stage II decreases until it reaches a plateau, which is dominated by the imaging noise floor, indicating that the MDF of the 1D gas is indistinguishable from a Gaussian (thermal) distribution.
To further study the relaxation in the dimensional crossover, we implement a molecular dynamics (MD) simulation in a semiclassical framework. The dynamical evolution is described by quasiparticles characterized by their spatial coordinates, quasimomenta, and transverse modes' occupations. In the longitudinal direction, the quasiparticles are initialized according to the thermodynamic Bethe ansatz [68] but move as classical point-like objects. The transverse degrees of freedom are treated as discrete quantum levels, which is the new key (quantum) ingredient of our model. This model is valid in our Newton's cradle experiments because the quantum correlations are strongly suppressed in the first ten oscillation periods (5% of the evolution time), as the chemical potentials of the 1D gases turn from positive to negative resulting from the decreases of 1D densities. Note that this time is too short to appreciably change the oscillation period averaged quasimomentum distribution due to the interplay between the effects of atomic scattering and the longitudinal trapping, which is known to break down the integrability and to induce relaxation [27] (see also [23]). As the dephasing in each tube happens, the 1D gases enter the non-degenerate regime, where γ T −1/2 is fulfilled (see Sec. 2.1 and 5). In the non-degenerate regime, the filling of states is much smaller than unity, whereby the effects of quantum statistics become negligible. When close to the non-degenerate limit, quasiparticles become individual atoms, and quasimomenta can be interpreted as usual momenta.
Scattering and transitions between these transverse states are calculated. The collisions between quasiparticles follow the parity selection rules, and the lowest excitation energy is E th . The transition matrix elements determining the transition probabilities are obtained from quantum-mechanical calculations similar to those of Ref. [46]. The transverse state of a quasiparticle is specified by a number n of transverse excitation quanta, with n = 0 corresponding to the ground state of the transverse motion. We do not resolve the degenerate sublevels but invoke the statistical weight (i.e., degeneracy) w n = n + 1 of the corresponding state of an isotropic 2D harmonic oscillator.
A harmonic longitudinal potential is assumed in order to make calculations simple and fast. We performed two tests to accept this assumption. (i) We tested an anharmonic potential U 0 tanh 2 (x/∆x) that admits analytic integration of the equations of motion. The parameter U 0 was taken equal to the lattice depth and the typical length scale ∆x was chosen to provide the harmonic potential 1 2 mω 2 x 2 for |x| ∆x. Using these parameters, the effect of the anharmonicity of the potential was found to be small. (ii) We initialized the simulation of dynamics in a harmonic trap with a fully dephased distribution, where the two Bragg components are indistinguishable and overlap during the whole oscillation period. Compared to the normal case starting with two distinct Bragg peaks, we observed a very small effect on the calculated quasimomentum distributions within our model. In other words, the anharmonicity is important during the initial relaxation stage only and we restricted ourselves to the harmonic model.
Each numerical realization corresponds to a single tube. The number of quasiparticles per tube N as an input parameter is set to the weighted-average value measured in experiments. The initial distribution of quasiparticles is sampled according to a thermal distribution calculated by solving the Bethe-ansatz equations at the measured temperature in an experimentally defined harmonic trap. Afterwards, each quasiparticle obtains a boost of quasimomentum −2k Br a g g or +2k Br a g g with equal probability (1 − η)/2. As such, we kick the quasiparticles with Bragg momenta, leaving a fraction of η of quasiparticles at the trap center. To match with the experimentally measured initial MDFs, η is set to be 9% and 20% for N = 40 and N = 130, respectively. Subsequently, we propagate the quasiparticles over time according to the functions described in Appendix B.
The change of transverse states of quasiparticles due to heating in the optical lattice is included in the model. Its probability per unit time per quasiparticle is denoted by Γ . We assume N Γτ 1, whereτ = 2π/(ω N 2 ) is the typical time between two subsequent atomic collisions. Within our simulation, we set Γ = 0.0375 s −1 for both N = 40 and N = 130, which is determined by the heating rates measured in experiments (as shown in Appendix A.2) . Whenever a collision occurs (i.e., when coordinates of two neighboring atoms coincide), we check both the possibility of the change of transverse states for the colliding pair of quasiparticles and the possibility of the change of transverse states for all the quasiparticles. We generate a pseudorandom number ζ uniformly distributed between 0 and 1. If ζ < exp(−N Γ τ), where τ is the time elapsed since the previous collision, then no state change occurs. Otherwise, we pseudorandomly select one of the N quasiparticles and change its transverse excitation num- For comparison purposes, we estimate the MDFs (green dotted curves) from the quasimomentum distributions (QMD) calculated in MD simulations (red solid curves) and observe good agreement with the experimentally measured profiles (black dashed curves). As the 1D gases approach the non-degenerate limit, the MDF and QMD tend to coincide.
ber n j to |n j + 1| or |n j − 1|, each channel having the probability of 50%. Since quasiparticles are predominantly in the ground state, the most probable process n j = 0 → n j = 1 leads to the energy supply to the system (heating).
The quasimomentum distribution derived from the MD simulation is distinct from the MDF (of the real physical bosons) when the system is far from the non-degenerate limit. For comparison with the experimentally measured profiles, we need to estimate the corresponding bosonic MDFs for the MD results. The relation between both distributions is not straightforward. There is no general analytic approach to calculate the MDF in the Lieb-Liniger model, and only numerical methods were for example applied in [27,[61][62][63][64][65][66][67]. Within the scope of this paper, we use an estimation of the MDF as outlined in Appendix C instead of an exact numerical calculation. Fig. 8 presents a comparison between the calculated quasimomentum distributions from the MD calculation (red), the estimated MDF (green) and the experimentally measured MDF (black dashed curves), which are all shown in oscillation period averaged profiles. The estimated MDFs are in good agreement with the experimental measurements throughout the entire dynamical evolution. Meanwhile, we observe the increasing similarity between the MDF and the quasimomentum distribution during Stage I of evolution.
The relaxation processes of the MD results are characterized by following Eq. (2), in the same way as in processing the experimentally measured profiles (see Fig. 7). R(t) calculated from both MDFs (R MDF , dotted curves) and quasimomentum distributions (R QMD , solid curves) are compared to the experimental results (R EXP , open circles). In Stage I of evolution, R MDF displays similar features with R EXP . The relaxation of MDF due to the transition from the degenerate to non-degenerate regimes is captured by the estimation of MDF. After the dephasing ends, R MDF asymptotically approaches R QMD , and at longer times the two results tend to coincide.
Furthermore, it is demonstrated by the plot of R QMD that the system relaxes towards a Gaussian (thermal) distribution at an accelerating rate. As has been mentioned at the end of Sec. 5, the 1D condition of the cradle system is preserved until the transversely excited states are populated with a threshold of atom number. For a 1D system with N = 40, it at least needs 5% of atoms in the first transversely excited state to break the integrability with inelastic collisions. These atoms may come from the minimal residual heating. From the plot of R QMD for N = 40, we observe that it takes about 1 s to start the relaxation towards a Gaussian (thermal) momentum distribution. The onset of relaxation leads to an energy transfer from the longitudinal to transverse degrees of freedom (an increase of the population in the transversely excited states), which in return intensifies the relaxation itself.
In contrast to N = 40, the relaxation towards a Gaussian (thermal) momentum distribution occurs earlier and faster for N = 130. It is because that the higher atom number offers a larger chance for the atoms to be excited transversely, so that it equivalently lower the threshold. Additionally, owing to the higher collision rate, the excited atoms spend less time in the upper state before returning back to the ground state through collisions. During this process, the excitation energy is deposited to the ground state. While, in the lower atom number case, the excited atoms stay longer in the upper state, having a small but non-zero probability of being de-excited through external disturbances. Moreover, in the initial distribution of quasiparticles for the MD simulation, we expect about 2% of quasiparticles obtaining quasimomenta larger than k th for N = 130. This fraction of atoms is not observed in Fig. 3(b) because of the narrower profile of MDF compared to the quasimomentum distribution. In contrast, this fraction is expected to be almost zero for N = 40. The very rare high-energy quasiparticles speed up the relaxation to some extent. For the above reasons, we observe faster relaxation for N = 130.  To further illuminate the relaxation process driven by the inelastic collisions, we study the dynamics of atoms in the transverse states by band mapping in the deep-lattice limit [69]. The fractional populations of atoms in the first and second transversely excited states (notated by η 1 and η 2 , respectively) are extracted throughout the evolution from the experimental band mapping images (see Appendix A.4). We observe that approximately 11% and 1% of the atoms are excited to the first and the second transversely excited states in 4.8 s for N = 130. These fractions are larger than the ones introduced by heating, meaning that energy is transferred from longitudinal to transverse directions in the cradle experiment.
In Fig. 9, we compare the experimental data with the MD simulations and observe very good agreement also in the transverse degrees of freedom. The minimal discrepancy from simulations is expected to be caused mainly by the additional excitation during the lattice unloading procedure for band-mapping. This conjecture is demonstrated by measuring η 1 in the first half oscillation period after the Bragg pulses (shown in the inset of Fig. 9). The additional excitation is clearly observed when the atoms are at the oscillation phase with large momenta at the end of the unloading procedure. For this reason, we accept the local minimum of η 1 in a half oscillation period as the most credible measurement of the in-trap transverse excitation; for example, the result at t = 3 ms represents the situation in the first half period. This effect becomes inevitable as the oscillations of the particle are diffused and becomes weak at longer times as most of the atoms are scattered to the low-energy regime. In this section, we go deeper into the crossover regime by initializing the atoms in the 1D tubes with considerably higher longitudinal momenta (energies) (Fig. 3(c)-(h)). Atoms that are kicked to large momenta ±4k Br ag g and ±6k Br ag g obtain enough energy to be transversely excited through subsequent collisions. Compared to the dynamics discussed in Sec. 3-7, these high energy collisions rapidly drive the system out of 1D integrability and towards thermalization. In Fig. 10, we present the evolution of R(t) for the dynamics initialized by the distributions shown in Fig. 3. For comparison, the branches (a) and (b) are replicas of the results shown in Fig. 7. The dynamics starting with higher imprinted energies exhibit faster relaxation.

Starting with stronger excitations of the longitudinal motion
The data is again compared to the MD simulations implemented for describing the dynamics. The initial quasimomentum distribution for the simulation is obtained by assuming the same respective particle number in each momentum peak with the experimentally measured initial MDF. The discrepancy between the MDF and the quasimomentum distribution reduces faster in the case with higher energies. Apart from the density reduction due to the single-body dephasing, the increase of temperature also drives the system out of the degenerate regime, especially in the high-energy cases. The MD simulations show good agreement with experimental observations in Stage II (after the quasimomentum distribution coincides with the MDF). The time scale for reaching an MDF that is indistinguishable from a Gaussian (the time when R(t) drops below the noise floor) is accurately predicted by the simulation. The validity of the semi-classical model is confirmed throughout the regimes from 1D to the dimensional crossover.

Integrability breaking through virtual excitations
Up to now, we considered only collisions with enough energy to excite the transverse degrees of freedom. However, even when the two-body collision energy is below the threshold of populating the transverse states, the latter can be virtually excited. A collision with a third atom can return the system on the energy shell and simultaneously redistribute momenta of the threeatoms involved in such an effective three-body collision [18,21]. However, this integrabilitybreaking mechanism does not contribute much to the relaxation of MDF in our experiment. Indeed, the rate of the velocity-changing collisions per atom due to this mechanism is given by [22] Here, l ⊥ = ħ h/mω ⊥ , a s is the 3D s-wave scattering length, and g 2 (0) is the density-density correlation function at zero distance. g 2 (0) is equal to 2 in a non-degenerate gas, to 1 in a weakly interacting quasicondensate and rapidly tends to zero if γ → ∞ [70]. Since the 1D density substantially decreases after the rapid dephasing stage, the typical time 1/Γ ∞ of the MDF relaxation due to virtual transverse excitation extends well beyond the time scale of our experiment.

Conclusion
We have investigated the relaxation processes of bosons at the onset of the dimensional crossover from 1D to 3D. We demonstrate that the system relaxes in two stages under different mechanisms. At short times, the single-body dephasing rapidly drives the 1D gases into the nondegenerate regime, during which the momentum distribution function deforms and asymptotically approaches the quasimomentum distribution of Lieb-Liniger model. At longer times, a tiny fraction of atoms in the transversely excited states triggers the transition from 1D to 1D-3D crossover. Subsequently, the system relaxes towards an equilibrium with a Gaussian momentum distribution through inelastic two-body collisions at an accelerating rate. A molecular dynamics simulation was implemented for efficiently modeling the non-equilibrium dynamics. Meanwhile, we proposed a simple method of estimating the momentum distribution functions in the whole regimes of quantum degeneracy. The numerical results quantitatively fit the experimental observations from short to long time scales in all three dimensions. Moreover, the long-term dynamics of a Newton's cradle with minimal heating and loss as can be obtained in a red-detuned lattice offers a model system to test theoretical methods for describing the complex dynamics in many-body systems at the point of breaking integrability. Future prospects include detailed studies of integral dynamics and its breakdown in the framework of the recently developed generalized hydrodynamics (GHD) [39,40]. In a first step, we have recently extended the applicability of GHD to the dimensional crossover regime [44] and tested it with the data at short to intermediate time scales. But still, many open questions remain, such as for example the many-body dephasing induced by non-trivial interactions [27], or the effect of atom losses [71]. We hope our investigations presented here will pave the way towards a more comprehensive understanding of the non-equilibrium quantum physics in the dimensional crossover regime and the influence of the effectively compactified dimensions.

A.1 Preparation of 1D gases
The 87 Rb BEC is produced in the Zeeman sublevel F = 1, m F = −1 by evaporative cooling in a crossed optical dipole trap. In the final stage of the evaporative cooling, the atomic cloud is levitated by switching on a magnetic field gradient in the vertical direction and decompressed by reducing the trap laser power. The total atom number, N t ot , is tuned between 1 × 10 4 and 1 × 10 5 by holding the BEC for different times in a shallow trap, where the BEC is overcooled due to the low trap depth. Afterwards, the BEC cloud is adiabatically transferred from the optical dipole trap to a 2D square optical lattice located in the horizontal plane. To avoid interference, the two lattice beams derived from a fiber laser are detuned 220 MHz from each other and have orthogonal polarization. The beam waist (w ol = 145 µm) of the optical trapping beam is much larger than the BEC.
During the lattice loading procedure, the lattice depth is exponentially ramped to the max-

A.2 Heating process and atomic loss
The heating in an optical lattice is mainly caused by two reasons (i) the spontaneous scattering of lattice laser photons; (ii) the trap fluctuations (including the laser intensity fluctuations and the pointing stabilities of lattice laser beams) at specific frequencies. The former mechanism heats an atomic system by transferring the photon recoil momenta to atoms. The latter excites atoms to higher transverse states, and the energy may be deposited into the longitudinal kinetic energy through the subsequent inelastic collisions. The heating effect is in general stronger for 1D gases with higher atomic densities because of the larger collision rates.
In our experiments, the heating process is studied by observing the evolution of 1D gases held in the identical lattice setup without the Bragg-pulse excitation. The MDFs for both N = 40 and N = 130 exhibit the expansion of the momentum peaks (see Fig. 11). By summing up the contribution on each pixel of the MDF measurement, we obtain the increase of kinetic energy of 0.06ħ hω ⊥ /s and 0.09ħ hω ⊥ /s for N = 40 and N = 130, respectively. On the other hand, we also estimate the heating rate by evaluating the energy growth in the transverse dimensions and observe a rate 0.006ħ hω ⊥ /s for N = 130 (see Fig. 12). Most of the transversely excited atoms appear in the first state, while the signal in the second excited state is so weak that it is submerged in the imaging noise.
Although the heating rate in the ground state is usually higher in a red-detuned lattice [28][29][30][31], we suppress the heating in our system to a minimal value by the large detuning of the lattice laser and the carefully controlled environment. The heating rates achieved in our experiments are at least twice as low as observed in Ref. [6] in a blue-detuned lattice.
The atom loss observed in our experiments is between ∼ 4%/s (N = 40) and ∼ 7%/s (N = 130), which are about one order of magnitude lower than observed in Ref. [6]. Such low loss rates enable us to study the long-term dynamics of bosons out of equilibrium. The loss rates are nearly constant throughout the evolution. We do not observe any significant three-body loss as seen in Ref. [6,32].  Figure 13: An example of the horizontal imaging. A region of interest (ROI) is arranged to contain all atoms, and two background regions BG1 and BG2 are selected next to ROI. The pixel size is 6.45µm×6.45µm in the object plane. To obtain the longitudinal MDF and estimate the noise floor, we integrate the images in these regions over the transverse direction, respectively. The MDF (blue) and noise level (green and purple) are shown on the right.

A.3 Detection of longitudinal momentum distribution function
The 1D gases are detected by standard absorption imaging after being released from the lattice trap and expanding in 3D. To keep the signal-to-noise ratio of images at a comparable level, the expansion time for N = 40 and 130 are set to 10 ms and 30 ms, respectively. In the horizontal plane, the image is taken along the bisector of two lattice beams. The lattice is turned off in 500 µs, during which the interparticle interaction vanishes. It is fast compared to the longitudinal dynamics but slow enough to prevent the atomic cloud from spreading too much in transverse directions. By integrating the image over the transverse direction in the region of interest (ROI), we obtain the longitudinal distribution profile. This profile approaches MDF  f (k) after a long TOF. To assess the impact of the imaging noise, which mainly stems from the photon and atom shot noise on the CCD camera, we choose two background regions (BG1 and BG2) beside the ROI with the same size (see Fig. 13) and consider them as the noise floor in the data analysis.

A.4 Detection of populations in transverse states
In the vertical direction, band mapping is applied to obtain information of the population in the respective energy-band. The lattice depth is exponentially turned off in 2 ms. During the turnoff, the crystal momentum is mapped to the free particle momentum, and afterwards, the Brillouin zones are imaged.
Two examples are shown in Fig. 14 to explain the method of evaluating the fractional population of atoms in each state. White boxes separate the 2D band-mapped distributions into Brillouin zones corresponding to the ground state, the first and second excited states. The biggest issue to overcome for achieving an accurate fraction in each state is the overlap between adjacent Brillouin zones due to the broadening of the quasi-free-particle momentum distribution. Firstly, we integrate the band-mapped distribution in the central region (the unshaded area) in z direction over the width of the first Brillouin zone. Secondly, we fit the integrated distribution K( y) to a summation of 180 Gaussian curves arranged with equal spacing and identical r.m.s. width between y = −3ħ hk ol R and y = +3ħ hk ol where A( y c ) is the amplitude of the Gaussian curve centered on y c 4 . By integrating A( y c ) in the corresponding regions, we get the fractional population of atoms in each state. The same calculation is repeated in the z direction, and η 1 , η 2 are calculated according to the results from both y and z dimensions. Since the atoms in the second excited state located in the four corners are not included in the calculations in both dimensions, η 2 is multiplied by 1.5 under the assumption of a uniform distribution.

B Model of the atomic collision in molecular dynamics calculation
In this section, we describe how atomic collisions are modeled. Since the system is 1D and quasiparticles are indistinguishable, we can always consider an ordered array of quasiparticles, In what follows, it is convenient to introduce the scaled coordinatesx j = x j /l and quasimomentaq j = q j l , where l = ħ h/(mω ). For a given configuration of N quasiparticles in the phase space we calculate the time of the first collision, i.e. the first (smallest) time when the coordinate of any two neighboring quasiparticles coincide. The oscillatory motion of the jth atom is described bȳ Then we calculate the collision time τ j for the jth and ( j + 1)th quasiparticles: and find the smallest one, We propagate the quasiparticles until the time t + τ according to Eq. (5) and then decide, according to the probabilities (see below) and using a pseudorandom number generator, what happens to the transverse states of the involved quasiparticles. The probabilities of the change of the transverse state are based on the standard quantum mechanical expressions, which can be easily derived for a pair of colliding quasiparticles with the initial state of their relative motion in the ( y, z)-plane as the transverse ground state [46]. However, for the sake of simplicity, we neglect any dependence of the transverse transition probabilities on the transverse quantum states of colliding quasiparticles. As such, the scheme of transverse transitions is simplified.
The collisions are assumed to be instantaneous, in other words, the Wigner delay time is neglected. This is justified by the observation that, even though the interplay between Wigner delay time and the longitudinal trapping potential will lead to appreciable thermalization, in our non-degenerate system this time scale is very long [23] and exceeds the duration of the experiment.
To be definite, consider a collision of the quasiparticles 1 and 2. Their coordinates at the collision time are x 1 = x 2 and the respective quasimomenta are ħ hq 1 and ħ hq 2 . The quasimomentum of the relative motion is canonically conjugate to the interatomic distance x 2 − x 1 and defined as The total quasimomentum of the pair is denoted by Concerning the transverse quantum numbers, we begin with the option Because of the parity conservation, the transverse energy of a pair of quasiparticles in the course of a collision can change by a multiple of 2ħ hω ⊥ . If the kinetic energy of the relative motion is less than 2ħ hω ⊥ , then the increase of the transverse energy is impossible. In the opposite case, the increase of the transverse energy by 2ħ hω ⊥ is possible. The probability of such an event is where α 1D = c/2. A pseudorandom number ζ uniformly distributed between 0 and 1 is generated. If ζ < P ↑ then we raise the transverse excitation energy by 2 quanta. To preserve the ordering of quasiparticles in the course of the subsequent evolution, we assign the new (primed) quasimomenta to them as follows: With the help of a new pseudorandom number we assign the new transverse quantum numbers with the following probabilities: The detailed balance condition should be satisfied: the number of transitions up and down per unit time must be the same on average. Therefore, if n 1 = n 2 > 0 then we allow for the transition to the state characterized by where The probability of this process is P n 1 ,n 1 →n 1 −1,n 1 −1 = n 1 n 1 + 1 where The prefactor in front of P ↓ in Eq. (8) ensures the detailed balance. The condition of the downward transverse transition corresponds to the pseudorandom number ζ falling between P ↑ and P ↑ + P n 1 ,n 1 →n 1 −1,n 1 −1 .
This is always the case when two quasiparticles in the ground transverse states collide with the energy insufficient for excitation by two transverse quanta. Consider now another possibility n 1 = n 2 .
Here an important simplification of the model comes into play. If the transverse states of colliding quasiparticles are different, we neglect, except of a special case described below, the change of the set of the transverse excitation numbers, allowing for the exchange process only, when the transverse excitation numbers associated with the two quasimomenta ħ hq 1 and ħ hq 2 are interchanged: The probability of the exchange process is given by A special case is given by In this case, in addition to the exchange process, the decrease of the larger of the transverse excitation numbers by 2 can happen, as it is required by the detailed balance: n 1 = n 2 = min(n 1 , n 2 ) and the quasimomenta after collision are given by Eq. (7). The respective probability is given by P |n 1 −n 2 |=2→n 1 =n 2 = 1 2 min(n 1 , n 2 ) + 1 min(n 1 , n 2 ) + 3 P ↓ , where P ↓ is again given by Eq. (9).

C Estimation of bosonic momentum distribution function
For a degenerate 1D Bose gas, the MDF, defined as the Fourier transform of the correlation function, is distinct from the quasimomentum distribution within the Lieb-Liniger model. In the regime where the temperature is low (below the chemical potential) and the Lieb-Liniger parameter γ is not excessively large, the correlation function for bosonic Luttinger liquid at x ħ h/mc s is written as where C 0 ∼ 1, c s is the speed of sound, k T = k B T /(ħ hc s ) and K = πħ hn 1D /(mc s ) is the Luttinger liquid parameter. The MDF W (k) = +∞ −∞ d x 2π g 1 (x)e ikx , and it is expressed via Euler's betafunction [60] W (k) = C 0 2 1 2K where B(x, y) = Γ (x)Γ ( y)/Γ (x + y). W (k) consists of a Lorentzian profile on top of a pedestal. The central Lorentzian is restricted to k < k T , and k T is the momentum where the Bose-Einstein distribution starts to deviate from the Rayleigh-Jeans classical limit and to decrease exponentially. For larger momenta k k T , W (k) decreases ∝ const/k, slower than Lorentzian. For much larger momenta k k C , W (k) is determined by Tan's contact and decreases ∝ C k −4 , where C is the Tan's contact and the momentum k C is approximately equal to the maximum quasimomentum at zero temperature. There are known approaches to precisely calculate the value of Tan's contact, see, for example, Ref. [51]. Considering the experimental uncertainties, we use an asymptotic approximation for large momenta k ξ −1 h , where ξ h is the healing length, for a weakly interacting quasicondensate. For stronger interactions, the qualitative picture is similar. The modified MDF is For k ξ −1 h , W (k) ∝ k −4 . In general, the W (k) expressed by Eq. (13) and (14) is expected to be much narrower than the corresponding quasimomentum distribution.
In the non-degenerate limit, the MDF coincides with the quasimomentum distribution. Furthermore, when the temperature T is relatively high, the MDF approaches the Maxwell-Boltzman distribution, Let us consider that we try to derive the MDF corresponding to a distribution of quasiparticles ρ target (x, q), namely the target distribution. The basic idea to estimate the MDF is done by fitting the target distribution with a sum of multiple thermal distributions ρ(x, q) = i ρ i (x, q), which are calculated by solving the Bethe-ansatz equations in a harmonic trap defined with experimentally measured parameters. In a general case of quantum Newton's cradle experiments, the ρ(x, q) consists of three components of quasiparticles, described by thermal distributions ρ 1 , ρ 2 and ρ 3 , respectively. The central component ρ 1 centers at the origin of the phase space with zero mean quasimomentum (〈q 1 〉 = 0). The symmetric Bragg components ρ 2 and ρ 3 are derived from the Bragg pulses, and they are shifted by the mean quasimomentum boosts 〈q 2 〉 and 〈q 3 〉 (〈q 2 〉 = −〈q 3 〉 = 2k Br ag g ). All of ρ i are normalized to the respective quasiparticle number N i , subject to the restriction i N i = N .
To imitate the effect of single-body dephasing, we drop the densities of ρ 2 and ρ 3 by a scale factor α n , which is calculated by assuming a linear expansion of the Bragg components in the direction of motion. α n decreases from 1 at t = 0 and stops decreasing when two Bragg components merge until fully dephased. The expanding rate is determined by evaluating D(t) for the calculated distributions and comparing the results with experimental observations as shown in Fig 6. To make the fitting procedure simple and fast, we integrate the distributions ρ target (x, q) and ρ(x, q) over the direction of motion and seek for the minimum discrepancy between the distributions in the radial coordinate. We accept this simplification in a quasi-harmonic potential when the interaction is not excessively large. The distribution in the radial coordinate hardly changes within one period of oscillation. The best-fit distribution returns us the chemical potential µ i and temperature T i for each component. Following Eq. (13)(14)(15), the MDFs in the degenerate and non-degenerate limits are derived.
As we discussed in the main text, the cradle system evolves from the degenerate regime to the non-degenerate regime during the dynamical evolution. To interpolate the crossover between two limits, we propose an empirical formula: convolving the MDF for the degenerate limit W i (k) with its counterpart for the non-degenerate limit M i (k) Here we modify Eq. (15) via M i (k) = M i (k/β) so that we rescale the width of the profile. β is tuned from 0 to 1 in the crossover regime and it follows β ∝ (µ i /T i ) 2 when µ 2 < 0. For the central component ρ 1 , µ 1 −T 1 is obtained throughout the entire evolution of time, resulting in an MDF very close to its corresponding quasimomentum distribution. While for the Bragg component ρ 2 and ρ 3 , µ 2 and µ 3 evolve from positive to negative and in longer times becomes much smaller than −T 2 and −T 3 . Thus, we obtain peaked MDFs in the early stage of evolution and gradually rounded MDFs that asymptotically approach the quasimomentum distributions as the system evolves towards the non-degenerate limit.
Since the mean quasimomentum equals to the mean momentum, f i (k) is shifted to be centered at 〈q i 〉. The MDF for the entire cloud f (k) = i f i (k).