Prethermalization Revealed by the Relaxation Dynamics of Full Distribution Functions

We detail the experimental observation of the non-equilibrium many-body phenomenon prethermalization. We study the dynamics of a rapidly and coherently split one-dimensional Bose gas. An analysis based on the use of full quantum mechanical probability distributions of matter wave interference contrast reveals that the system evolves towards a quasi-steady state. This state, which can be characterized by an effective temperature, is not the final thermal equilibrium state. We compare the evolution of the system to an integrable Tomonaga-Luttinger liquid model and show that the system dephases to a prethermalized state rather than undergoing thermalization towards a final thermal equilibrium state.


Introduction
A general understanding of the dynamics of non-equilibrium quantum many-body systems is an important unsolved problem impacting many areas of physics [1,2]. In particular, the question of why and how isolated quantum systems relax toward equilibrium states has only been studied for a very limited number of special systems [3][4][5]. Moreover, in many cases, such as when dealing with integrable systems exhibiting many constants of motion, efficient relaxation can be completely absent [6] or strongly inhibited.
An intriguing phenomenon that has recently been experimentally demonstrated in this context is prethermalization [7,8]. The concept of prethermalization was first introduced in order to explain the remarkable success in applying thermodynamic models to heavy-ion collision experiments on time scales much shorter than the expected thermalization time of the scattered particles [7,[9][10][11]. In particular, the successful application of hydrodynamics [12] relies on the assumption of being at least approximately close to thermal equilibrium at every point in space and the presence of an equation of state [13,14].
Taking a low-energy quark-meson model as a model system, Berges et al [7] suggested an explanation in 2004 for this unexpected, apparent early thermalization. In their theoretical analysis, Berges et al were able to show that after a very rapid initial evolution, their model system relaxes to a non-thermal quasi-steady state. This intermediate state is very robust with respect to changing initial conditions and already shows many of the bulk characteristics of the final thermal state. For example, the ratio of kinetic energy to pressure (i.e. the equation of state) or a properly defined kinetic temperature is already practically indistinguishable from their values in thermal equilibrium. In contrast, mode quantities such as the occupation numbers of different momentum modes are still far away from a true thermal equilibrium even though during the initial rapid evolution they also reach quasi-steady values. Their work showed that although the model system had not yet fully thermalized, it could in many aspects be described using a thermal model. Accordingly, the above described phenomenon was termed prethermalization.
In the present understanding, prethermalization is characterized by the rapid establishment of a quasi-stationary state on time scales much shorter than the expected thermalization time. This state is long lived and already exhibits some thermal-like properties, but can still be very different from the true thermal equilibrium state of the system. A full relaxation to thermal equilibrium, if present at all, is then expected to happen on a further much-longer time scale. While thermalization can be interpreted as a full loss of the memory of the system about its initial state, prethermalization describes only a partial loss of information about the initial state.
In the following, we discuss the direct observation of the prethermalized state presented in [8] in more detail. We give, in particular, the details of the experimental system and analysis methods used to reveal the prethermalization and present additional evaluations that support the findings of [8]. An even more complete description can be found in the PhD Thesis of Gring [25].

Experimental system
The experimental study is performed using trapped 1D Bose gases. Such systems offer two unique advantages for non-equilibrium experiments. First, on the experimental side, realizing them with ultracold atoms facilitates a precise preparation and probing of the system. Secondly, on the theoretical side, 1D Bose gases offer a model system that contains complex manybody physics, but can still be captured with reasonable theoretical effort, particularly due to the existence of effective models that allow one to describe the essential physics in a relatively simple way [26]. Furthermore, the homogeneous 1D Bose gas with repulsive contact interactions is an example of a fully integrable quantum system [27,28]. The approximate realization of such a system in experiments thus allows the study of thermalization in the vicinity of multiple conserved quantities and hence the study of the interplay between integrability, many-body dynamics and thermalization. Figure 1 summarizes the main idea of our experimental study of prethermalization using 1D Bose gases in the quasi-condensate regime [29][30][31]. In this regime, density fluctuations are strongly suppressed and the gas is characterized by strong phase fluctuations. The properties of these phase fluctuations are determined by the temperature and the density of the system. In the  [8]. (a), (b) A phase fluctuating 1D quasi-condensate is coherently split creating two 1D gases with almost identical phase profiles φ L (z) and φ R (z) (represented by the black solid lines). The gases are then allowed to evolve in the doublewell potential for some time t e which leads to strong fluctuations in the local phase difference φ(z) and a decrease of the phase correlation length λ φ . The question at which this experiment aims is whether and how this dynamical state reaches the thermal equilibrium state corresponding to two independently created quasi-condensates. In this thermal equilibrium state the phase difference between the 1D gases fluctuates strongly along its length and the correlation length λ T is determined by the temperature and density of each cloud. (c) The phase difference φ(z) between the two 1D gases is probed through time-offlight matter wave interference of the two gases where the local relative phase is directly transformed to a local phase shift of the interference pattern. The contrast C(L) of the axially integrated interference pattern can then be used as a direct measure of the strength of the relative phase fluctuations. Figure adapted from Gring et al [8]. experiment such a single trapped quasi-condensate is rapidly and coherently split, producing two uncoupled 1D gases with identical phase profiles. The aim of our study is to probe how these initially almost perfect correlations of the relative phase become obscured over time and if the thermal equilibrium state corresponding to two completely independent gases is finally reached [32,33]. To this end, the two gases are allowed to evolve in the double-well potential for a varying evolution time t e before the relative phase correlations are probed via time-offlight matter wave interference (figure 1(c)). As differences in the relative phase lead to a locally displaced interference pattern, the contrast of the longitudinally integrated interference pattern is a direct probe for the dynamics of the system [34][35][36][37]. Example interference patterns after various evolution times, demonstrating the loss of the initial phase coherence, are shown in figure 1. Comparison of the initial non-equilibrium state prepared by the coherent splitting process and the thermal equilibrium state of two independent quasi-condensates. The system can be described by a combination of relative degrees of freedom (a) and common degrees of freedom (b). After the coherent splitting the common degrees of freedom contain all the thermal energy E ∼ k B T of the system, whereas the relative degrees of freedom contain only the quantum noise of the splitting process. Here, k B is Boltzmann's constant and T the temperature of the system. In thermal equilibrium all degrees of freedom contain the same amount of thermal energy and are thus characterized by the same temperature.

Initial non-equilibrium state and thermal equilibrium
A prerequisite for non-equilibrium experiments is the ability to precisely prepare and characterize both the initial non-equilibrium state and the expected thermal equilibrium state of the system. One of the key advantages of coherently split 1D Bose gases is that both these states can be prepared and described with high precision. For a general system of two spatially separated 1D Bose gases, the excitations can be described by anti-symmetric and symmetric longitudinal modes which relate to the relative and common degrees of freedom of the two halves of the system. They are given by for the phase, and for the density. Here φ L,R (y) describes the longitudinal phase profiles and n L,R (y) the density of the left and the right gas, respectively. After the coherent splitting, the individual phase profiles of the two halves are almost identical and hence the relative phase profile is almost flat (figure 2, left column). In terms of excitations, this means that all the thermal excitations are initially contained in the common degrees of freedom. As we will detail in section 6, the relative degrees of freedom, on the other hand, are initially populated only by quantum noise created in the splitting process.
In thermal equilibrium, we expect to find the system in a state where the energies in the relative and in the common degrees of freedom are equal (figure 2, right column). This thermal equilibrium state corresponds to the situation of two independently created quasi-condensates where the phase profiles of the clouds are uncorrelated down to the thermal phase correlation lengths, resulting in a relative phase profile that fluctuates strongly along the length of the system [32,33]. In experiment, the thermal equilibrium situation can be purposely created, by splitting a thermal gas into two, followed by the creation of two independent quasi-condensates through further cooling in the double well.
The transition between the initial non-equilibrium state and the thermal equilibrium state can therefore be directly probed by studying the relative phase fluctuations.
Considering the above, we note that any thermalization mechanism must necessarily redistribute the mode populations of the relative and common degrees of freedom to lead to the thermal equilibrium situation.

Full distribution functions of matter wave interference contrast
Previous experiments have also investigated the dynamics of split Bose gases [38,39] or equivalent systems [40]. In these works, though, due to the lack of theoretical or experimental tools for a full characterization of the transient states, the driving mechanism of the dynamics could not or only indirectly be revealed. Consequently, the state to which the system decayed remained elusive.
This directly illustrates the key challenge common for most experiments trying to observe non-equilibrium dynamics: the scarcity of experimental tools for characterizing the complex many-body states during the evolution. In this work, this problem is approached by generalizing tools developed for studying equilibrium systems of spatially separated 1D Bose gases [32,33,35,41] by applying them to the non-equilibrium case, as recently proposed in [23,24]. For this purpose, we measure the time evolution of full quantum mechanical probability distribution functions (FDFs).
When dealing with quantum many-body systems, valuable information about the underlying quantum state of the system can be gained through noise correlation measurements. Performing such measurements has already deepened our understanding of quantum mechanics as it led to the discovery of the Hanbury Brown-Twiss effect [42] which triggered the development of modern quantum optics [43]. Furthermore, the study of current fluctuations led to important observations in quantum-Hall systems [44,45]. Recently in atomic physics the analysis of noise correlations revealed the coherence properties of atom lasers [46] and enabled observations of the Hanbury Brown-Twiss effect for massive fermions and bosons [47]. It was further suggested [48] and experimentally demonstrated [49][50][51][52][53] that noise correlations in time of flight can be used to probe strongly correlated equilibrium states of quantum many-body systems.
In 1D systems, fluctuations play a much more pronounced role than in their threedimensional (3D) counterparts [54,55]. The interference pattern of two expanding 1D quasicondensates therefore inherently contains strong noise and fluctuations that can be used to study the many-body state of the system. When performing interference experiments with spatially separated 1D Bose gases, one is furthermore in the advantageous position that the interference pattern itself is directly related to the correlation functions of the system. For example, the mean squared value of the integrated matter wave interference contrast |C| 2 is a measure of the integrated two-point phase correlation function. Similarly, by exploiting the strong shot-toshot fluctuations of the contrast, higher moments of the contrast can be obtained which contain information about higher-order correlation functions [35,56]. An experimentally accessible quantity capturing all these higher moments of the contrast and therefore also the higher-order correlation functions is given by the full distribution function of the contrast P(C). This function describes all the fluctuations of C in the measurement, as P(C) dC measures the probability of observing a contrast between C and C + dC. Higher moments of C k of the contrast can be obtained by integration C k P(C) dC.
It is important to point out that the method of using FDFs to characterize a system requires the detection of single realizations of the quantum system in question. If only ensemble averages can be measured in the experiment, the statistics of those values will always be Gaussian due to the central limit theorem and the characteristic higher moments of the observable will not be accessible. The method of analyzing quantum states through the FDFs has already been successfully used to characterize isolated [32] and tunnel-coupled [33] systems of two separated 1D Bose gases in equilibrium. Within the context of [23,24], this method was extended theoretically to non-equilibrium systems, whereas the experimental demonstration of this method is the topic of the work presented here and in [8].

Experimental details
The preparation of the atomic sample is performed following our standard procedure [57] using magnetically trapped 87 Rb in the 5S 1/2 F = 2, m F = 2 state on an atom chip [58]. The atom chip is a current-carrying gold structure micro-fabricated on a silicon substrate. A basic outline of the most important structures on the chip is given in figure 3(a).

Realizing a quench: the coherent splitting process
We first create a single-well trap using the static magnetic field of a 100 µm wide wire in combination with homogeneous bias fields. Longitudinal confinement is provided by additional wires. The trap is located at a distance of about 100 µm away from the atom chip surface (see figure 3(a)). This strongly anisotropic trap has measured trap frequencies of ω ⊥ = 2π × (2.1 ± 0.1) kHz (radial confinement) and ω = 2π × (11 ± 1) Hz (axial confinement).
To create the double-well potential, we use an adiabatic dressed-state potential [60,61]. By applying RF current to two 30 µm wide wires that are adjacent to the central 100 µm wire, the cloud is split along the radial direction of the trap and perpendicular to the direction of gravity. This ensures that any potential sag due to gravity is common for both wells and that the 1/r dependence of the RF fields emanating from the wires is the same for both wells. For the splitting process, the amplitude of the RF current is linearly increased from 0 to typically 22.5 mA within 17 ms and the frequency is 30 kHz detuned to the red of the F = 2, m F = 2 → F = 2, m F = 1 transition at the minimum of the initial static magnetic trap. This creates a double well with a separation of (2.75 ± 0.05) µm, a simulated barrier height of (2.9 ± 0.1) kHz and measured trap frequencies for each well of ω ⊥ = 2π × (1.4 ± 0.1) kHz and ω = 2π × (7 ± 1) Hz. Great care was taken to guarantee stability and repeatability of the  [8]. (b) Simulation of the longitudinal density profiles of trapped quasi-condensate gases following [59]. The black solid line shows the profile of the cloud before splitting for an atom number of 6000 atoms. The colored dotted, dotted dashed and dashed lines show the profile of the sum density of two quasi-condensates trapped in the double-well trap with a longitudinal trapping frequency of 6.5, 7 and 7.5 Hz. For the experimentally realized double-well trapping frequency of 7 Hz, the profile matches almost exactly to that of the single-well trap at 11 Hz. splitting process, especially to gain a symmetric splitting with almost zero mean atom number difference between the two halves. For the fine tuning of this quantity, the relative amplitude of the RF current applied to the two 30 µm wide wires can be slightly adjusted.
Following [33,62], any residual tunnel coupling between the two halves of the system was calculated to be J 2π × 0.1 Hz. In this regime, a residual tunnel coupling would have no significant effect on the results presented in this work [41]. The absence of tunnel coupling was also confirmed by experiments with independently created condensates in the same double well [33]. In the measurements presented in section 5, we further used a larger final RF amplitude which resulted in a larger splitting distance and larger barrier height corresponding to an even lower bound for the residual coupling and found no difference to the previous results.
Although the splitting RF current was ramped up in a time of typically 17 ms, the actual splitting process is much faster than the time scale of the longitudinal dynamics and happens close to the end of the ramp when the coupling between the two gases vanishes. From the evolution of the mean relative phase between the two gases for deliberately imbalanced double wells [63], we determined this decoupling to occur (15 ± 0.5) ms after the start of the RF ramp. All evolution times given in the following refer to this point in time. To capture also the evolution that happens after the splitting but still during the ramp up, it is possible to probe the atoms during the ramp up.
While the actual splitting is much faster than the ramp up of the RF current, the trap is nevertheless continuously deformed throughout the whole ramp up time. The choice of 17 ms total duration is a compromise between realizing a fast splitting and avoiding the excitation of collective oscillations of the cloud. The trapping frequencies of the final double well were chosen to be such that the longitudinal profile of the clouds before and after splitting are almost Splitting of the gases is performed along the x-axis, the y-axis corresponds to the axial direction of the 1D Bose gases, gravity points downwards along the z-axis. The atomic density distribution in time of flight can be imaged along each axis of the coordinate system using absorption imaging. Imaging along the x-axis is used to transversally probe single gases before the splitting or the incoherent sum of the two independent gases in the double well, the latter case being depicted in the figure. The corresponding absorption picture shows strong density fluctuations which are used for thermometry (for details see text). Imaging along the y-axis and the z-axis is used to probe the interference pattern of two 1D Bose gases. Imaging along the y-axis results in a projection of the interference pattern along its length, while imaging along the z-axis reveals the full undulating structure of the interference pattern. exactly matched, as can be seen from the simulations in figure 3. This minimizes the creation of longitudinal breathing oscillations of the clouds during splitting. As we will discuss in section 7 these breathing oscillations and the general instability of the splitting process severely limited previous studies of coherently split 1D Bose gases [38]. We find that longitudinal center-ofmass oscillations of the clouds after splitting cannot be completely avoided without the use of more complex splitting protocols. However, these small oscillations do not significantly alter the physics as they do not change any significant parameters of the cloud such as density or total size.

Probing the system
The system is studied using absorption imaging, which integrates the 3D structure of the expanding atomic density distribution along the one direction of space. Figure 4 shows schematic views of the integrated density distribution as seen from each of the three detection directions realized in the experiment.
Imaging along the x-direction is used to probe the system from its transversal direction. This can be applied both to extract information about a single quasi-condensate before the splitting process as well as about the incoherent sum of two quasi-condensates expanding from the double well. As the phase of a quasi-condensate fluctuates in space, density fluctuations will appear along its length when it is released from the trap and allowed to expand in time of flight [64][65][66]. These spatial density fluctuations are a direct consequence of matter wave interference between the different positions along the single quasi-condensate which have different phases. In equilibrium, the in situ phase profiles are determined by the temperature of the system, so that the measurement of the spectrum of density fluctuations in time of flight can be used for thermometry of the system before the splitting process, as experimentally demonstrated in [66]. All initial temperatures mentioned in the following, as well as the temperatures of pairs of independently created gases in the double well have been determined using this method.
Information about the relative phase between the two 1D clouds can be obtained from absorption images along the longitudinal axis of the two gases in the y-direction, or transversally along the vertical z-axis.
In an intuitive picture, the interference pattern can be thought of as an array of many thin interference patterns stacked up together along the longitudinal direction of the cloud, where the phase of each thin interference pattern depends on the local in situ phase difference between the clouds. The axially integrated contrast C(L) of the interference is therefore a direct measure for the relative phase fluctuations within the integration length L [34][35][36][37]. In the following, we will discuss the technical implementation of the detection systems for the longitudinal and vertical directions in detail and show how we extract the integrated contrast C(L) for different integration lengths L.

Optical slicing method.
In this scheme, the imaging beam propagates along the longitudinal axis of the 1D gases and directly performs an integration of the interference pattern along the length of the gases. In order to vary the integration length L, we perform the following slicing scheme using optical pumping (figure 5(a)). Shortly (1 ms) after the two halves of the split condensate are released from the trap, a shadow of width L is imaged on the atoms using a 50-125 µs pulse of optical pumping light. The frequency and the polarization of the light are chosen such that the atoms, initially being in the 5S 1/2 F = 2, m F = 2 state, are excited to a state in the 5P 3/2 F = 1 manifold ( figure 5(b)). From there, the atoms preferentially decay to a state in the 5S 1/2 F = 1 manifold which is dark to the imaging beam operating on the 5S 1/2 F = 2 → 5P 3/2 F = 3 transition. Atoms that are in the protected shadow area of length L are not optically pumped and can be imaged after whichever chosen time of flight using an absorption imaging system operating along the axial direction of the system.
The branching ratios of the 5S 1/2 F = 2 → 5P 3/2 F = 1 pumping transitions are such that, on average, only very few photons are needed for pumping. Already after one photon has been absorbed and remitted an atom has a probability above 80% of being in a dark state. This ensures that the atoms remain virtually undisturbed by the optical pumping process.
To create the sharp-edged shadow for the optical slice, we use a multi-lens high-resolution imaging objective in reverse. The imaging system has a maximum resolution of 2.5 µm, and a demagnification of approximately 1:4.6. The shadow is produced by a glass target with a series of opaque bars (produced at the ZMNS TU-Wien). In order to be able to change the integration length L in each experiment, the bar target is mounted on a computer-controlled motorized translation stage. The effect of the optical pumping light on the detected atomic density is presented in figures 5(c) and (d). Figure 5(c) shows examples of sliced atomic density distributions after optical pumping and a subsequent time of flight. The images were obtained using the same imaging system (x-direction) that is used in reverse for the optical pumping. Example interference images are shown in figure 5(d) for three integration lengths and a hold time of 7 ms in the double-well potential. Note that due to the evolution in the double-well potential, the contrast of the interference pattern which is integrated over the full cloud is already significantly reduced in most of the images. Furthermore, since for shorter L more and more atoms are cut away, the detected atom number becomes lower, giving a limit to the shortest integration length which can be investigated.

Direct imaging method.
In the vertical direction, we can directly image the undulating structure of the interference pattern and freely select the integration length via postprocessing. To this end an imaging system is employed which observes the atoms from below ( figure 6).
The implementation of such an imaging system is challenging since the chip surface blocks the line of sight, hindering standard absorption imaging schemes. One possibility to circumvent this problem is to use a light sheet fluorescence detector as in [33,68]. Here, we instead rely on absorption imaging, as it allows for a higher imaging resolution.
To this end, we reflect the imaging beam from the gold surface of the atom chip close to normal incidence. In this procedure the high quality of the gold layer is essential to achieve good images. For a collimated imaging beam the reflection means that the light passes the atom cloud twice [69]. This results in a second, virtual image behind the atom chip, which is out of focus and overlaps with the primary image. Since the imaging of the interference pattern requires the imaging beam to propagate almost parallel to the fringes, it is not possible to offset the two images by using a large angle of incidence when propagating the imaging beam onto the atom chip surface. In order to eliminate the second image we therefore focus the imaging beam close to the atoms, avoiding the first absorption process as shown in figure 6(a). This limits the field of view to the area of a 'spotlight' (figure 6(b)) and ensures that the light only interacts with the atoms once. Also, by aligning the slowly diverging beam onto the central 100 µm wire of the chip, interference effects caused by diffraction from the microscopic chip structures can be minimized. We note that a similar method was used in [70]. In figure 6(b) we show an experimentally obtained image, where the imaging system was set up for a time of flight of 16 ms. The image clearly reveals the full structure of the spatially varying interference pattern, which can then be used for further analysis.

Extraction of the matter wave interference contrast
The main observable used in this work is the contrast C(L) of the matter wave interference pattern integrated over a length L ( figure 1(c)). For the optical slicing, the integration length L is selected in the imaging setup by choosing a specific mask, and the integration is performed optically. In the direct imaging method the recorded image of the interference pattern is integrated over the length L during post-processing. In both cases this results in a line profile (see figure 1). This line profile is fitted with a sinusoidal function f L (x) having a Gaussian envelope where ω is the cloud size, is the fringe spacing and θ L is the global phase of the interference pattern integrated over L, i.e. the phase of the line profile. As interactions do not play a significant role in the radial time-of-flight expansion for our parameters, the fringe spacing can be directly related to the double-well-potential separation d = ht/m , where t is the time of flight, m is the mass of the atoms and h is Planck's constant [71]. Thus, the stability of the double-well splitting can be monitored with interferometric precision. The contrast C(L) is extracted from this fit and used for further analysis that is presented in the remainder of this paper. Repeated realizations of the experimental cycle then allow us to build the time-dependent FDFs of the interference contrast C(L).

Comparison of the two imaging methods
The use of two technically different imaging systems allows us to preclude systematic effects of the imaging system on the observed contrast evolution.
Example FDFs obtained for different integration lengths L from both imaging systems are shown in figure 7. The total area under the curves is normalized to 1 in order to represent a probability density. The upper panel shows FDFs of the squared contrast C 2 . While the overall form of the FDFs shows qualitative agreement, a small offset between the results of the two imaging systems remains. This is because the absolute value of C 2 is affected by the finite imaging resolution and any mis-alignment of the imaging system. These effects can be eliminated by analyzing contrast distributions that are normalized to the mean contrast C 2 . After this procedure the offset vanishes, as shown in the lower panel of figure 7.
This comparison confirms that the form of the C 2 / C 2 distributions is not influenced by any technical effects of the imaging systems and can therefore be used to study the dynamics of the system. Note that the normalization to the mean contrast also removes further systematic errors affecting both imaging systems, such as, for example, imperfections in the switch-off of the trapping fields.
In the following, the only and important practical difference between the two methods is thus that the total measurement time with the direct imaging is significantly shorter than the time needed to obtain the data with the optical slicing method. This is a consequence of the fact that with the direct imaging method only one experimental run is necessary to extract information on all length scales down to the imaging resolution. Systematic studies of large parameter spaces or long evolution times, as presented in section 5, are thus always performed with the direct imaging system. Comparison of distributions of the squared contrast obtained with the optical slicing method and the direct imaging system for the same doublewell configuration. The evolution time is 10 ms and each distribution contains 300 realizations, which minimizes statistical effects. In the upper row we plot exemplary distributions of the squared contrast, in the lower row the squared contrast is normalized to its mean. While the former shows slight differences for the results of both imaging systems, the latter shows very good agreement (see text for details). Figure 8(a) shows the measured time dependence of the FDFs of interference contrast C(L) for the double-well system which was initialized through a coherent splitting process, as described in the previous section. The data were obtained using the optical slicing detection scheme. The central line density is ρ = (32 ± 7) atoms µm −1 (error given by standard deviation) in each well at an initial temperature of the unsplit cloud of T = (120 ± 30) nK. For the initial gas before the splitting, this results in a chemical potential of µ ∼ (2/3)hω r and an energy k B T ∼hω r due to temperature. These parameters are similar to the parameters of other experiments with 1D gases in microtraps, e.g. [72,73], and give an effectively 1D Bose gas in the quasi-condensate regime.   After the splitting, a slight net mean imbalance n = (N L − N R )/(N L + N R ) of n = (1.5 ± 1)% is present. Here, N L and N R are the atom numbers of the left and right well, respectively. Owing to technical noise, the width of this imbalance distribution is around two to three times larger than the quantum shot noise limit of √ N , where N = N L + N R is the total atom number of the system. The imbalance has mainly two effects. The net difference in chemical potential resulting from the mean imbalance leads to a trivial common phase evolution of the two wells [37,74,75]. In addition to this, the width of the imbalance distribution leads to an additional broadening of the global phase distribution. It is a strength of our approach that this technical noise does not affect the dynamics of the contrast distributions, as its effect on the dynamics is negligible compared to the quantum noise associated with the splitting process for the probed integration lengths L.

Dynamics of the full quantum mechanical probability distribution functions of matter wave interference contrast
For the shortest evolution time t e = 2 ms, the C 2 / C 2 distributions are peaked on all lengths L, being almost all identical in form and having a low probability of obtaining a low contrast C(L). As time evolves, the distributions behave very differently on different length scales. For long integration lengths L, the FDFs evolve to become exponential in form, i.e. there is a high probability of observing a low contrast C(L) and the initial correlations of the system appear to be lost. For short integration lengths, however, the FDFs remain peaked and there is a very low probability of observing a low contrast C(L). This directly demonstrates the persisting memory in the system of the initial state. Furthermore, the strong length dependence is a direct signature of the multimode nature of 1D Bose gases [40,76]. Finally, the shape of the distributions evolves very quickly for the first 12 ms and very little between t e = 12-27 ms, showing that the system reaches a quasi-steady state. This imposes the question of whether this state corresponds to the thermal equilibrium state of the system.

Emergence of an effective temperature
We will now address the question of whether the observed quasi-steady state already corresponds to the thermal equilibrium state of the system.
If the system has reached thermal equilibrium, the measured distributions plotted in figure 8(a) should be consistent with equilibrium theory [35,41,56]. The obtained nonequilibrium distributions of C 2 / C 2 were therefore fitted with equilibrium distributions taking temperature as a free fit parameter. To this end, a χ 2 analysis [77] was performed. Figure 8(a) shows the experimentally obtained distributions together with the respective best fitting equilibrium curves. The agreement is very good for longer evolution times, but less good for early evolution times. Figure 8(b) shows the corresponding reduced χ 2 values of the fits. For short evolution times the reduced χ 2 differ significantly from 1, showing that the experimentally obtained distributions are inconsistent with equilibrium theory. For evolution times longer than 12 ms, the observed value of the reduced χ 2 approaches 1, which shows that the experimental data agree well with equilibrium theory. This indicates that the thermal-like nature of the distributions is established dynamically during the evolution.
It is thus possible to associate the observed steady state with a temperature T eff extracted from equilibrium theory. For the last three evolution times t e = 12, 17, 27 ms the fits give T eff = (14 ± 4), (17 ± 5), (14 ± 4) nK, respectively. This is, however, a factor of eight lower than the initial temperature of the unsplit system. Furthermore, the corresponding effective thermal phase coherence length of the relative phase between the two gases is λ φ =h 2 ρ/mk B T eff = 13 +5 −3 , 10 +5 −2 , 13 +5 −3 µm. This is much longer than the expected λ T = (1.5 ± 0.6) µm deduced from the initial temperature of the unsplit cloud. For a detailed study of the emergence and properties of this length scale see [76].
As a further illustration of the striking difference between the observed FDFs and the ones expected for the corresponding equilibrium system, we performed interference experiments with two independently created 1D gases, employing the same double-well potential as for the non-equilibrium situation. The measured equilibrium distributions are plotted in figure 9. Even for T = (60 ± 20) nK the measured C 2 / C 2 distributions are exponential in form for all length scales probed, which is in stark contrast to the observed non-equilibrium steady-state distributions but in perfect agreement with what is expected from equilibrium theory. A second, colder dataset at T = (27 ± 7) nK visualizes the crossover from an exponential to a Gumbel-like shape of the equilibrium distributions [32,56].
It is, however, a justified question as to whether the very low observed T eff of the quasisteady state could also be explained by a cooling effect involved in the splitting process or, stated differently, how does the initial temperature of the unsplit cloud relate to the final temperature of the gas after splitting and after a possible equilibration? In principle, the rapid splitting procedure is a rather violent process which is very likely to transfer some energy to the system. This would suggest that the final equilibrium temperature should be even higher than the initial temperature of the unsplit system, in contrast to the steady state that we observe.
On the other hand, the splitting process leads to a decompression of the gas, which, for a 3D ideal gas would lead to a temperature reduction. If we assume an adiabatic splitting a comparison of the trap frequencies leads to a lower bound on the decompressed temperature of (0.6 ± 0.1) × T [78]. This temperature reduction has been confirmed by preparing a thermal 3D cloud of atoms in the initial single-well trap and turning on the double-well potential using the same 17 ms ramp as in the non-equilibrium measurements. For this procedure, we observe a decrease in temperature of (0.59 ± 0.15) × T in agreement with decompression.
To our knowledge, no quantitative descriptions of this process exist for degenerate 1D gases. The exact evolution and distribution of the thermal energy during the splitting therefore remains a topic of ongoing research. The only quantities one can compare to are the initial temperature T = (120 ± 30) nK of the unsplit system and the lower bound of T ∼ (70 ± 20) nK calculated under the assumption of an adiabatic decompression. Figure 8 therefore also shows the obtained values for the reduced χ 2 if the experimentally obtained non-equilibrium distributions are fitted with equilibrium theory for T = 120 nK. The plot shows that these high temperatures are clearly rejected by the χ 2 test. In addition, in figure 8, the corresponding equilibrium distributions for T = 30 nK (green line) are plotted for comparison. This shows that even for T = 30 nK, which is much lower than the lower temperature bound obtained from decompression, the equilibrium distributions are still significantly different from what is observed in the experiment for the non-equilibrium system.
These observations directly demonstrate that the quasi-steady state we observe has thermallike properties, but is not the true thermal equilibrium of the system. We thus associate this quasi-steady state with prethermalization, as introduced in [7] and suggested for split 1D Bose gases in [24]. The system decays rapidly to a state whose contrast distributions are thermal like in form, exhibiting a temperature almost an order of magnitude lower than the temperature of the initial unsplit system. As we will see in section 6, the relaxation of the system to this prethermalized state is very well described by a dephasing of the momentum modes of the system [23,34,76]. We will first briefly discuss the persistence of this prethermalized state as the system continues to evolve.

Long-term evolution of the system
As shown in section 4, the non-equilibrium system decays rapidly to a prethermalized state. However, the question still remains as to any subsequent evolution of the system, particularly whether the system will eventually reach the true thermal equilibrium corresponding to two independent gases.
In the measurements presented in the following, the initial temperature of the unsplit system was (78 ± 10) nK and the peak atomic line density of the split system was (38 ± 9) atoms µm −1 in each of the two wells, i.e. approximately similar parameters to those presented in section 4 above. Furthermore, the final amplitude of the current for the dressing RF creating the double well was increased to 23.8 mA leading to a higher potential barrier, larger separation of the double well and therefore placing an even lower upper bound on any hypothetical residual tunnel coupling. To ensure that the actual splitting process was also the same, the same linear ramp speed as in the measurements in section 4 was used, but the total ramp duration was extended to 18 ms. Figures 10 and 11 summarize the main results of the measurements. Figure 10(a) shows the evolution up to a time of almost 200 ms of the mean squared contrast C 2 for the full length of the system. The fast initial decay to the prethermalized state, followed by a very slow further evolution, is clearly visible. As an inset, example interference pictures obtained using  [8,76] together with example pictures obtained using the direct imaging system. (b) FDFs of C 2 / C 2 extending to longer evolution times. Again, the experimentally obtained data (histograms) can be very well described by equilibrium theory where the distributions obtained at the best-fitting temperature are plotted using red solid lines. The blue dashed lines show theoretical equilibrium distributions for the initial temperature T = 78 ± 10 nK of the unsplit cloud. Figure adapted from Gring et al [8].  [8]. the direct imaging system are shown. For very short evolution times the fringes are almost straight, directly visualizing the coherence of the splitting process. For larger evolution times the fringes become more and more wiggly, leading to the observed strong decrease in the integrated interference contrast. Again, for evolution times after the initial decay, the distributions of C 2 / C 2 can be described by equilibrium theory, as shown for example in figure 10(b) for various evolution times t e = 12-182 ms. Figure 11(a) displays the minimal reduced χ 2 values for the inspected evolution times which shows that for evolution times after the initial rapid decay, the reduced χ 2 is similar to 1. The extracted effective temperatures are again very low compared to the initial temperature of the system, while higher temperatures clearly fail to describe the measured distributions, as demonstrated by the χ 2 analysis for larger temperatures in figure 11(a).
Observing the long slow decay of C 2 seen in figure 11(a) one may be tempted to assume that this is the system slowly thermalizing. In order to analyze this slow subsequent evolution of the quasi-steady state, the effective temperatures for all times after the initial decay were extracted, as shown in figure 11(b). As can be seen, T eff rises slowly over time at a rate of (0.14 ± 0.04) nK ms −1 . This is, however, comparable with the measured heating rate of the atom trap of (0.11 ± 0.06) nK ms −1 which was characterized independently using equilibrium quasicondensates. This indicates that either no thermalization is present, or, if it is present, that it is a very subtle process.

Luttinger liquid description of the dephasing of the system
In the previous sections it was shown that the system created through the coherent splitting of a single quasi-condensate rapidly relaxes to a quasi-steady state that has thermal characteristics but which does not correspond to the true thermal equilibrium state of the system. This observation can be very well understood by the theoretical model presented in [23,24] where the appearance of a thermal-like state is explained through the dephasing of the multimode system. Before presenting a quantitative comparison of this model to our experiment, we summarize the most important points of this theory in order to give an intuitive picture for an understanding of the prethermalized state. For a more detailed discussion the reader is referred to the original work [23,24]. For an alternative description of the rapidly split quasi-condensate in terms of an effective-tunneling model, see [79].

Summary of the theoretical model
The evolution of the interference pattern is directly determined by the evolution of the relative phase between the two uncoupled halves of the split gas. We describe this evolution using a Tomonaga-Luttinger liquid approach. In this low-energy approximation, the relative degrees of freedomφ(y) andn(y) perfectly decouple from the common degrees of freedom of the systemφ com (y) andn com (y), where the experimental observables are now described by operators. The resulting Hamiltonian for the relative degrees of freedom is of the form H =¯h c 2 L/2 −L/2 dy[ K π (∇φ(y)) 2 + π Kn 2 (y)], where K = πξ h ρ is the Luttinger parameter, c is the speed of sound, ξ h is the healing length and ρ is the atomic density.
The evolution ofφ(y) can be described in Fourier space by a set of decoupled harmonic oscillators of collective modes with momentum k. A collective mode with momentum k modulates the relative density and relative phaseφ(y) along the condensate in a sinusoidal fashion on a length scale of ∼1/k. This is visualized in figure 12.
In the experiment, the splitting is performed fast in comparison to the time scale for the spread of perturbation which is set by the inverse chemical potential t split < ξ h /c = 1/µ. Thus, there is no time for the atoms to develop correlations along the longitudinal trap axis such that the relative density fluctuations n 2 t e =0 = ρ/2 are completely random and not affected by atomic interactions. In particular, there is no correlation between modes with different momenta. Also, for each atom the decision of going to either half of the split system is random and uncorrelated with other atoms, leading to a binomial distribution of atom number fluctuations in each small segment of the 1D system. The respective width of the relative phase distribution follows from the Heisenberg uncertainty relation. Initially, at t e = 0, the relative phase fluctuations are governed completely by the atomic shot noise and are practically negligible at length scales > ξ h . In other words, φ 2 t e =0 ∼ 0 and its fluctuations are strongly suppressed. In contrast, initial density fluctuations exhibit for |k| ξ −1 h a large excess of noise compared to the zero-temperature state of two split quasi-condensates. In the latter stationary case, the relative density fluctuations are suppressed by atomic repulsion. As a result, almost all the energy is initially distributed in the density fluctuations, with each collective mode k containing equal energy but different populations that scale as 1/|k| as the energy per quanta in each mode is ∝ c|k|. In particular, this leads to the weaker effects of modes with higher momenta on the evolution of the system.
During the dynamics of the system the energy of each harmonic oscillator mode k oscillates between the fluctuations in density and fluctuations in phase, driven by interactions between the atoms. This results in a harmonic time dependence of the fluctuation amplitude with a period ∝ 1/c|k|. For short evolution times, all phase fluctuation amplitudes grow in magnitude, which leads to a scrambling of the relative phaseφ(y) along the axial direction. This results in an initial rapid decrease in the interference contrast C(L). For longer times, the oscillations in different k-modes dephase and the system reaches a quasi-steady state which is thermal in its appearance. In particular, correlation functions take a form that is algebraically equivalent to their respective counterparts in equilibrium. The effective temperature of the quasi-steady state corresponds to the energy that was equally introduced to all k-modes by the splitting process. In analogy to the equipartition theorem [80] one can deduce that the quasi-steady state can be described by an effective temperature T eff given by [24] where µ is the chemical potential, g is the 1D interaction strength and k B is Boltzmann's constant. As shown in [24] the detailed form of the FDFs of C(L) and C 2 (L) can then be obtained by sampling the amplitude of the fluctuations of the phase with Gaussian statistics. The choice of a particular integration length L during the detection process represents a filter for the effects of different k-modes ( figure 12). For an integration length L, modes with wavelengths 2π/k < L produce strong fluctuations ofφ(y) within L, while modes with wavelengths 2π/k > L uniformly affectφ(y) in L and just change the overall phase of the interference pattern. This means that the contrast C of the integrated interference pattern is governed by the dynamics of modes with wavelength 2π/k < L.
The higher the population of a particular k-mode, the greater the mean amplitude of its phase fluctuation. For long integration lengths L, there are many significantly populated modes with k > 2π/L which affect the integrated interference contrast C. With increasing evolution time t e , the dynamics of these large amplitude modes leads to a strong decrease in the probability of observing a high contrast C and the FDF of the contrast takes the form of an exponential decay. This is the contrast decay regime. For short integration lengths L, the contrast C is influenced only by sparsely populated modes with high momenta so the FDF of the contrast preserves a non-zero peak for the whole evolution time. This is the phase diffusion regime. These two regimes have been experimentally studied in detail in [76]. The observation of these two regimes is a direct and intuitive visualization of the multimode nature of 1D Bose gases [40,76].
Note that the above discussion is only valid in the case of a perfect symmetric splitting. As discussed in detail in [24], a density imbalance between the left and right wells leads to a coupling between the relative and common degrees of freedom. In this case, the initial temperature of the unsplit cloud, which defines the initial state of the common degrees of freedom, influences the dynamics of the relative degrees of freedom and hence the dynamics of the interference pattern. However, the treatment in terms of collective excitations of the system remains valid along with the distinction between phase diffusion and contrast decay dynamics depending on the integration length L. The decay of the contrast then happens in two steps. Firstly, there is an initial rapid decay analogous to the symmetric splitting case. Secondly, the quasi-steady state which is reached, very slowly evolves in time due to the coupling between relative and common modes [24]. This is, however, not expected to lead to true thermalization as the system remains integrable and the modes will not fully equilibrate.

Quantitative comparison between experiment and theory
A quantitative comparison between experiment and the Tomonaga-Luttinger liquid model is presented in figure 13. We find good agreement of the theoretical description with the experimental data without any free fitting parameter.
Considering that the theory neglects the effect of the longitudinal trapping potential and contains no fit parameter, only the experimentally measured input parameters, the agreement of experiment and theory is very good. For the first two evolution times the description is, however, less accurate for two reasons. Firstly, the extraction of the contrast is very sensitive for the early evolution times, leading to a broadening of the experimental distributions by fitting uncertainties. Secondly, the disagreements also indicate that the initial preparation of the system through coherently splitting a single cloud is not fast enough to be instantaneous and perfectly uniform as assumed in theory.
By preparing 1D Bose gases with a wide range of different parameters, we can further systematically check the predictions of our theoretical model. Figure 14 shows the predicted scaling (equation (4)) of the effective temperature with density, as well as the independence of T eff from the initial temperature before splitting.
These two observations provide further strong support for our theoretical model and the interpretation of the observations presented in this work as prethermalization. In particular, the observation that the properties of the quasi-steady state are independent of the initial temperature of the unsplit system is another indicator that the observed state is clearly different from the thermal equilibrium state. The observed linear scaling of T eff with density ρ further shows that the prethermalized state is solely defined by the quantum noise associated with the splitting process. The apparent small systematic offset of the experimentally obtained T eff and the theoretical description visible in figure 14(a) can be attributed to imperfections in the experimental splitting process.
Returning to the discussion of the theory in section 6.1, one sees that the relative degrees of freedom φ and n do indeed decouple from the common degrees of freedom φ com and n com , at least on the time scales on which we observe the system. This is particularly interesting, since there is no a priori reason to assume that the different degrees of freedom should not couple and equilibrate in our system. Figure 13. Comparison of distributions of the normalized squared contrast C 2 / C 2 obtained in section 4 to the theoretical description of the rapidly split quasi-condensate presented in the previous section 6.1. The experimental data are plotted using histograms and the theoretical simulations solid red lines. The light red shaded areas denote the errors resulting from the uncertainty of the experimentally measured theory input parameters. These input parameters were the experimentally measured values of density in a single well ρ = (32 ± 4) atoms µm −1 , trap frequency ω ⊥ = (1.4 ± 0.1) kHz, imbalance n given in per cent as (0.1 ± 0.7)%, temperature of the unsplit system of T = (78 ± 10) nK and the uncertainty of ±0.5 ms for the point in time of t e = 0.

Dephasing, not thermalization
In an earlier experiment by Hofferberth et al [38], the dynamics of a coherently split 1D Bose gas was analyzed through the decay of the coherence factor (t) = Re e iφ s (z,t) dz . The coherence factor was first introduced in [34] and is closely related to the mean contrast. Note that while the observable C(L) is more easily extracted from experiments, it is more Figure 14. (a) Dependence of T eff on density ρ and (b) independence of T eff from the initial temperature T of the system before splitting, corrected for the scaling of T eff with density. The (black) solid line corresponds to the theoretical prediction k B T eff = gρ/2. The black (green) data points correspond to the dataset presented in section 4 (5), respectively. Figure taken from Gring et al [8].
complicated to calculate theoretically. A characteristic sub-exponential decay of the coherence factor (t) ∝ exp[−(t/t 0 ) α ] (with α ∼ 2/3) was observed, as predicted by Burkov et al [81] for the thermalization of such a system, where t 0 is the characteristic time scale for the thermalization. This observation was interpreted as evidence for thermalization. The basic physical system investigated in [38] was very similar to that in the new experimental setup presented in this paper. However, at the time of [38], neither the experimental methods nor the theoretical understanding were sufficiently developed to make a full analysis and characterization of the many-body state to which the system decayed. As demonstrated in this work, this characterization was now possible through the new, more complete theoretical description of the non-equilibrium processes [23,24] and the ability to experimentally measure the FDFs which enabled us to characterize the state to which our system decays. We would like to revisit here this old work [38] in the light of our new results.
In the work of Hofferberth et al [38], the coherent splitting was done in the vertical direction, along the direction of gravity. This allowed the direct observation of the phase fluctuations along the length of the system (as with the direct imaging system in this work, but without having to deal with the problem of the chip surface being in the line-of-sight of the imaging beam). However, in general, this simplified imaging scenario comes at the price that a vertical splitting process is not very well controlled and non-symmetric due to the interplay between the pull of gravity and the 1/r dependence of the amplitude of the RF radiation emitted from the atom chip wires. As an illustration, a 1 µm height difference corresponds to a change in potential energy of more than 2 kHz, which is more than the chemical potentials in the 1D gases. This gravitational sag can be partially compensated by adjusting the polarization of the RF fields, by changing the distance of the trap to the chip surface and by taking advantage of the 1/r dependence of the trapping and RF fields used to generate the double-well potential [82]. However, this constrains the double-well splittings that can be implemented. Furthermore, the compensation is never perfect and variations of the potential depth during the splitting process remain that are on the order of the chemical potential of the 1D gases.
As a consequence, significant shot-to-shot variations in the atom number difference between the two wells were observed in the previous work [38]. These overall atom number fluctuations led to an additional phase evolution φ(t) = t h µ, where µ is the difference in chemical potential between the two wells in each experimental run [37,74,75].
Consequently, in [38], these common fluctuations were rejected by evaluating the coherence factor through (t) = Re e i(φ(y,t)− φ(t)) dy = e iφ(y,t) dy| . Here φ(t) was determined for each individual interference pattern via φ(t) = arg ( dy e iφ(y,t) ), which for an infinitely long system makes the expression for the coherence factor equivalent to the original formulation of Burkov et al [81]. This way of analyzing the coherence factor, however, also rejected the contribution to the dynamics coming from the k = 0 mode (i.e. the phase evolution which, in a semi-classical picture, can be interpreted as being induced by the quantum shot noise of the difference between the total number of particles in each well). For a finite size system, this k = 0 contribution to the phase evolution is suppressed only by 1/ √ N where N is the total number of particles in both wells [34].
In addition, in [38], a strong longitudinal breathing mode was observed after the splitting which limited measurements to the first 11 ms of the evolution. As discussed in section 3 we split horizontally, perpendicular to the direction of gravity in the new experimental setup presented here. In this way, the symmetry during the splitting process is easily conserved, which makes the initial non-equilibrium state much better defined and the final evolution more repeatable. The atom chip in the new setup is designed in such a way that the longitudinal confinement can be varied independently of the radial confinement. This makes it possible to choose combinations that allow the splitting of a single trap into a double well with a minimized excitation of longitudinal breathing modes. This made the probing of long evolution times possible.
Furthermore, in the new experimental setup, in comparison to [38], we split much further away from the atom chip surface, facilitating a double-well separation that is much greater for a given barrier height than in the experimental setup detailed in [38]. Simulations suggest that this leads to a ten times faster break down of the tunnel coupling between the wells during the splitting process because, due to the larger separation, the overlap of the wavefunctions is much smaller once the classical connection between the wells breaks down. Hence, the requirement of fast symmetric splitting to reach the initial state which allows the observation of prethermalization (discussed in [23,24]) is better met in this work.
Any comparison to current theoretical models of the fast initial evolution requires a clear starting point of the evolution as the exact value of the characteristic exponent α is very sensitive to the position of t e = 0. By analyzing the new data and from novel insight from recent experiments on tunnel-coupled systems [33], it became clear that it is very difficult to determine the exact point in time when the two clouds separate during the splitting, the dynamics of the different momentum modes during the final moments of the splitting, and when an independent evolution starts. In the present evaluation, we used similar criteria as in [38] to set the starting time of the evolution (simulated plasma frequency 10 ms and an observable phase drift). We find that the starting point of the evolution is very difficult to define to better than ±0.5 ms, given our current experimental precision and theoretical understanding. This is not a problem for observing and analyzing the FDFs and the long-term behavior, but it prevents us from reliably evaluating α from our data. However, for completeness, taking the first 9 ms of the evolution presented in section 5 and analyzing the coherence factor in the same manner as it was performed in [38], we extract a value for α that is compatible with the previous result. Yet, studying the system for different initial temperatures, we observe an independence of the characteristic relaxation time t 0 from temperature, in stark contrast to the prediction of Burkov et al [81]. For an ideal, sudden splitting, our model predicts an exponential decay with α = 1 [24,34]. In the experiment, the time for the splitting is finite. Thus, to settle the very interesting question of how the observed dynamics starts and to study its details will require extensive additional theoretical work on the dynamics of the splitting process itself, as well as advanced experimental techniques [83].
In view of the present experiment and the theoretical analyses, the interpretation of Hofferberth et al [38] should therefore be revised as showing the fast integrable dephasing of the 1D systems, and not thermalization, as was suggested in [38] by the comparison to Burkov et al [81]. This revised view is in line with the observation of the absence of thermalization in 1D gases in optical lattices [6]. We point out that for the present experiment, even independent of our theoretical model, the observed clear difference between the FDFs of normalized squared contrast obtained for an equilibrium system formed by cooling in the double well (figure 9) and the FDFs observed after dynamic splitting (figure 8) provides direct experimental evidence that no thermalization is observed in our experiment.

Conclusion
We have detailed the experimental observation of the non-equilibrium many-body phenomenon of prethermalization. This was accomplished by observing the evolution of a rapidly and coherently split 1D Bose gas. Through the use of FDFs of the contrast of matter wave interference patterns, the quasi-steady state to which the system decays was found not to be a thermal equilibrium state. Using the theory developed for equilibrium systems, we were able to show that the FDFs of the interference contrast in the quasi-steady state are thermal-like in form, but display an effective temperature that is independent of the initial temperature of the gas before the splitting process. The early-time evolution of our system is well described by a Tomonaga-Luttinger liquid model that describes the dephasing of the system to the prethermalized state. This evolution can be visualized by means of FDFs and then contrasted to the true thermalization expected at long times. Moreover, these results clearly illustrate the power of this new method of probing the dynamics of quantum many-body systems through FDFs. Beyond the present work, FDFs are a very general and widely applicable method that can be extended to other quantum many-body systems [84].
The nature of the true thermalization that manifests itself as a slow evolution toward a steady state after the initial rapid decay to the prethermalized state still remains an open question. A possibility would be the integrability violation through virtual [85][86][87] or even real [85] excitations of the radial degrees of freedom.
In fact, the question of thermalization in our system is closely connected to the question of how two quantum mechanically correlated but spatially separated objects can lose their memory of that initial correlation, i.e. how do classical properties emerge from a closed quantum system? To answer these questions of thermalization and the emergence of classicality is an ongoing experimental [88] and theoretical endeavor.