Double light-cone dynamics establish thermal states in integrable 1D Bose gases

We theoretically investigate the non-equilibrium dynamics in a quenched pair of 1D Bose gases with density imbalance. We describe the system using its low-energy effective theory, the Luttinger liquid model. In this framework the system shows strictly integrable relaxation dynamics via dephasing of its approximate many-body eigenstates. In the balanced case, this leads to the well-known light-cone-like establishment of a prethermalized state, which can be described by a generalized Gibbs ensemble. In the imbalanced case the integrable dephasing leads to a state that, counter-intuitively, closely resembles a thermal equilibrium state. The approach to this state is characterized by two separate light-cone dynamics with distinct characteristic velocities. This behavior is rooted in the fact that in the imbalanced case observables are not aligned with the conserved quantities of the integrable system. We discuss a concrete experimental realization to study this effect using matterwave interferometry and many-body revivals on an atom chip.


Introduction
Non-equilibrium dynamics of isolated quantum systems play a central role in many fields of physics [1]. An important question in this context is whether the unitary evolution can lead to the emergence of thermal properties. For example, the eigenstate thermalization hypothesis conjectures that dephasing can lead to thermalization in systems with a chaotic classical limit [2][3][4][5]. On the other hand, integrable or many-body localized systems are expected not to thermalize at all [6], but instead relax to generalized thermodynamical ensembles [7][8][9][10]. An important role in both cases is played by the signal propagation during the nonequilibrium dynamics. It has been shown for many systems [11][12][13] that this propagation follows a light-conelike linear evolution in time with a characteristic velocity. This behavior has important consequences for the growth of entanglement in such systems [5,14].
In this manuscript, we study the dynamics in a pair of bosonic 1D quantum gases with number imbalance using the Luttinger liquid formalism. The dynamics of such quantum wires has recently been studied in great detail using atom chips [10,13,[15][16][17] or in optical lattices [18][19][20][21][22]. In particular, 1D Bose gases have been established as a prime experimental realization of a nearly integrable system with strongly suppressed thermalization [6,15,21]. So far, the consequences of this integrable behavior have mainly been studied for individual gases or sets of nearly identical gases. Here, we show that for imbalanced pairs of gases, i.e. gases that differ in their mean density, integrable dephasing can establish a state that closely resembles thermal equilibrium. The dynamics towards this state are found to be exceptionally rich, including a metastable thermallike state described by a generalized Gibbs ensemble and two distinct light-cone dynamics, each exhibiting their individual characteristic velocity. Beyond the fundamental interest in thermalization dynamics, our study is of high relevance for ongoing experiments, where density imbalances are often unavoidable and can thus fundamentally affect the interpretation of the results.

Model
In the following we consider an initial 1D Bose gas in thermal equilibrium that is coherently split into two gases (figure 1). This scenario is motivated by experiments with 1D Bose gases on atom chips, which have recently been established as an important model system to study non-equilibrium dynamics [10,13,15,23,24]. In these experiments gases are confined in two radial directions using a strong potential characterized by the trapping frequency w^, such that  m w <k T , B and they behave effectively one-dimensional. Here, μ is their chemical potential and T their temperature.
The effective low-energy description for each of the gases is given by the Luttinger liquid Hamiltonian  ò r f =  + ⎡ Figure 1. Relaxation dynamics in a pair of quantum wires. A phase-fluctuating 1D Bose gas in thermal equilibrium is coherently split into a pair of gases using a double well potential. This leads to almost identical phase profiles f 1 (z) and f 2 (z) (indicated by the solid black lines) for the two resulting gases and to a binomial distribution of atoms between them. We study here the relaxation of this highly correlated state. If the two gases after the split have the same densities ρ 1 =ρ 2 , there is no coupling between the anti-symmetric or relative degrees of freedom (d.o.f.), which contain only small fluctuations after the quench, and the symmetric or common degrees of freedom of the total system, which contain large fluctuations after the quench. In the dynamics, these fluctuations of the relative and common d.o.f. reach thermal-like steady states with different temperatures. If, however, the two gases after the split have different densities, relative and common degrees of freedom are coupled and their fluctuations equilibrate over a second, much slower timescale.
where -H and + H are again LL Hamiltonians We note the change in pre-factors, i.e. the factor 1/2 that appears in addition to the mean density (ρ 1 +ρ 2 )/2 in the first term and the appearance of g instead of g/2 in the second term of H ± , as compared to the original Luttinger Hamiltonians describing the individual gases(equation (1)). These changes arise due to our specific definition of the symmetric and anti-symmetric degrees of freedom in equation (2). Coupling between the new degrees of freedom is mediated by , which leads to a coupling between symmetric and anti-symmetric degrees of freedom. In the latter case, the fluctuations in the two condensates evolve differently causing a relative dephasing of the two gases over time.
We note in passing that this scenario holds promise as a platform to study spin-charge physics within the Luttinger liquid framework [27,28], with  H describing the spin and charge degrees of freedom, respectively. In experiments the relevant parameters can be made fully tunable by imbalancing both the densities and also the coupling constants (see appendices A and B).
Previously, the corresponding dynamics were obtained by diagonalizing -H [26,[29][30][31]. However, this is not sufficient to capture the full dynamics if Using this Hamiltonian to solve the equation of motion for the phase operator we find Here, the speed of sound in the individual gases is given by denote the initial values for the phase and density fluctuations, respectively. We can transfer this result into the familiar symmetric/anti-symmetric basis using q 2, . Our aim is to investigate the dynamics after a single 1D gas in thermal equilibrium (with temperature T in ) is split into two parts. The Hamiltonian of this initial gas is of the form given by equation (7) with 1D density ρ=ρ 1 +ρ 2 . The thermal expectation values for the second moments in classical field approximation are therefore where we have used f f = - † k k and = - † n n k k . All first moments vanish. Note that experimentally the splitting often results in a change in radial trapping frequency, such that the gases before and after splitting will differ in their 1D interaction strength wĝ . In the following, we will denote the interaction strength before splitting with g in and after splitting with g f to include this possibility in our calculations.
In the limit of fast splitting no information can propagate along the axis of the gases and we can assume that they both have the same phase profile after the splitting, which is identical to the profile of the initial gas. Moreover, we want to assume that in this case the probability for each atom to go into well 1 is ρ 1 /(ρ 1 +ρ 2 ) independently from where the other atoms go (binomial splitting) [26]. With these assumptions, we find for the fluctuations of the individual gases right after splitting The first term in equation (10) represents the shot noise from the binomial splitting process, which is anticorrelated, as expressed by the factor 2δ i, j −1. The second term in equation (10), as well as equation (11), stem from the thermal fluctuations of the initial condensate, and describe correlated fluctuations. Again, all first moments vanish. Assuming Gaussian fluctuations, the second moments are sufficient to fully describe the system. This assumption is justified for long enough length scales containing a large number of particles. Note that the same assumption has to be made for the validity of the Luttinger liquid model and, also, typically only such length scales are accessible in experiments.

Results
We can now investigate the dynamics by combining equations (10)-(12) with equation (8). To clearly reveal the relevant dynamics we use the approximation of small imbalance (D  1) between the two gases, and obtain with the average velocity c +  = (c 1 +c 2 )/2 and and the velocity difference c − =(c 1 −c 2 )/2. For a derivation without the small imbalance approximation see appendix C.
In the following discussion of the dephasing dynamics we first focus on the case when the length scales under consideration are much smaller then the system size (large/infinite system limit). Finite system sizes and the occurrence of revivals are discussed in section 4.
If there is no imbalance (i.e. Δ=0 and c − =0, c +  = c 1 =c 2 =c) the terms proportional to -( ) c kt sin 2 vanish and we obtain the well known dephasing dynamics where thermal correlations are established with a light-cone [31]. In this process thermal correlations instantaneously emerge locally within a certain horizon, while they remain non-thermal outside of the horizon. This horizon spreads through the system with a characteristic velocity that is given by c=c 1 =c 2 [13,30,31]. One can see this from equation (14) by realizing that at a certain time t all modes down to a lower bound given by 2ct·k lower =2π have dephased (note that the factor 2 comes from the square of the sine). The bound k lower therefore corresponds to the length-scale 2ct. The result of the dephasing dynamics is a prethermalized state with a temperature This temperature can be identified directly from equation (13) in the dephased limit, i.e. by averaging over kt and comparing to the result for a pair of gases (each with 1D density (ρ 1 +ρ 2 )/2) in thermal equilibrium (classical field approximation). It corresponds to the energy -( ) k T B eff that is added to the relative degrees of freedom during the splitting quench [15,29].
The corresponding expression to equation (13) for the common phase variance q á ñ 2 is derived in appendix C. One observes that the symmetric degrees of freedom exhibit the temperature in which is coming from the initial thermal fluctuations. However, the initial temperature is decreased through an interaction quench (second term in the brackets). In the splitting not only the density but also the density fluctuations are halved [32][33][34]. As they enter quadratically in the Hamiltonian, this leads to a decrease of a factor 2 in energy, which can be further modified by the aforementioned change in the interaction constant g. While the individual correlation functions of relative and common degrees of freedom are thus thermal, the state of the total system is non-thermal and has to be described by a generalized Gibbs ensemble with the two temperatures  ( ) T eff , respective [10,15]. With imbalance ( D > | | 0) the dephasing dynamics changes significantly. The aforementioned light-cone dynamics to the prethermalized state still proceeds with the average velocity c + . In addition, examining the additional terms in equation (13) we identify a second dephasing timescale characterized by the slower velocity c − . After complete dephasing (with the fast as well as the slow velocity), we end up with a second thermal-like state. Following the same procedures as before we can identify the temperature of this state to be identical for both relative and common degrees of freedom and given by 4 Again, this result can be interpreted intuitively in terms of the corresponding energies The first term corresponds to half the energy that is initially contained in the common degrees of freedom (equation (16)), the second term to half the energy introduced to the relative degrees of freedom during the quench (equation (15)).
Equation (17) hence describes an equipartition of energy that is dynamically established by the coupling term H c .
Note that equation (17) remains true, even without the assumption of small imbalance 5 , which was used to obtain equation (13). For typical parameters in atom chip microtraps the change in confinement leads to To visualize the corresponding dephasing dynamics leading to this equipartition in detail, we calculate the two-point phase correlation function [13]  This function measures the correlation between the relative phases θ(z) at two arbitrary points z and z′ along the length of the system and can directly be measured in the experiments [13]. An alternative but related probe, the interference contrast, is discussed in appendix D.
As discussed in the previous section, the initial fluctuations in our model are Gaussian and of course remain so during the evolution with the quadratic Hamiltonian. Therefore, the phase correlation function can be rewritten in the form CF z z , e z z 1 2 2 . In the limit of an infinitely large system this leads to In figure 2 we plot equation (19) for increasing evolution times, revealing a double light-cone. First, the system relaxes to the prethermalized state with exponentially decaying (thermal) correlations. For longer evolution times, the system relaxes further to the second thermal-like steady state. As in the previous light-cone-like relaxation to the prethermalized state, the system reaches this new final state for a given time only up to a certain horizon, but then follows a different shape beyond that point. The position of this horizon moves with a second characteristic velocity that is given by the (typically small) velocity difference c − of the individual gases. While both symmetric and anti-symmetric degrees of freedom reach a thermal-like state with temperature  ( ) T f , the complete system still differs from the thermal equilibrium of two condensates with equal density (ρ 1 +ρ 2 )/2 in the aspect that cross-correlations q q á ñ +k k , , between symmetric and anti-symmetric degrees of freedom do not vanish (see appendix C). Note that for the thermal equilibrium of two gases with unequal densities, the cross-correlations between common and relative degrees of freedom also do not vanish. However, they are still of different magnitude than in the completely dephased case.
Another way to discuss the question of whether the system dephases to thermal equilibrium is to have a look at the quantities for the individual gases. These could e.g. be studied in experiments using density fluctuations in time of flight [35] or by probing the density fluctuations in situ [36]. The phase variance of the individual gases is given by 4 Note that this temperature is obtained by comparison with the thermal equilibrium of two gases of equal density (ρ 1 +ρ 2 )/2. When comparing to two gases with unequal density ρ 1 and ρ 2 , the effective temperature is given by r r r r +´ ( ) . 5 This can be best seen from the dephased state for the phase variances of the individual condensates(equation (20)) and the crossterms dephasing to 0(equation (25)).  This expression is different from the results for the symmetric/anti-symmetric basis, highlighting how the observed dynamics and, in particular, also their timescales, are indeed intimately connected to the choice of observable.
In detail, the timescale for the dynamics within a single gas is, as expected, given by their speed of sound c i . However, the cross-correlations of the form f f á ñ k k 1, 2, dephase to zero with the slow velocity c − (see appendix C). After complete dephasing we therefore end up with two independent gases, which independently appear to be in thermal equilibrium with their respective temperatures ( ) T i f . However, for all the dynamics described only dephasing and no true thermalization has taken place. The intuitive reason for this complex behavior is that the two imbalanced gases are non-identical and dephase with respect to each other. Therefore, their individual and the difference between the ( ) T i f tends to zero and they approach the final temperature of the symmetric and anti-symmetric degrees of freedom given by equation (17). Also, already in the approach of this limit the small difference in final temperatures can be challenging to measure in an experiment. In both cases the system would thus appear completely thermalized independent of the choice of basis, with symmetric, anti-symmetric and individual degrees of freedom all exhibiting the same temperature. However, with imbalance going to zero, this approach of the final temperature would become infinitely slow.

Influence on many-body revivals
While the double light-cone dynamics are clearly visible in the correlation functions calculated for an infinite system, observing them directly in an experiment with a finite size system, in particular when a typically harmonic longitudinal confinement is present, is challenging. In particular, due to the nonlinear excitation spectrum in harmonic traps [30] the effect is severely scrambled by highly irregular many-body revivals. Examples of this behavior are shown in the appendix E.
However, regular and well controlled many-body revivals have recently been observed for the first time in homogeneous trapping potentials [24], demonstrating their power to probe the dephasing and higher-order interactions of phonon modes. In the following we thus illustrate the influence of our effect on such many-body revivals.
To this end, we first repeat our calculations for a system with periodic boundary conditions, but with a typical experimental finite size of 100 μm x400 h . In figure 3(a) we show the corresponding results for the phase correlation function. Due to the finite number of momentum modes they show clear rephasing behavior, as experimentally observed in [24]. However, due to the imbalance the two velocities in the system can be observed through the presence of two different types of revivals-slow revivals resulting from c − and fast revivals coming from c + . The value of the phase correlation function reached in the slow revivals depends on how well slow and fast revivals coincide.
A scenario more relevant for an experimental realization is the one of fixed boundary conditions ( f ¶ ¶ = z 0), which corresponds to the case of a hard walled box. This boundary conditions guarantee that the particle current at the box walls vanishes. The corresponding results are shown in figures 3(c), (d). Again, a clear distinction between slow and fast revivals can be observed, which can directly be connected to the two characteristic velocities.
Our observations have important practical consequences for the experimental study of integrability breaking in 1D Bose gases [37][38][39][40][41][42]. As the effect described here and true thermalization through integrability breaking would essentially lead to the same experimental signatures (i.e. thermal correlations corresponding to a temperature given by equation (17)) they would be very challenging to disentangle from measurements of correlation functions alone. In particular, any experimental effort clearly has to take both effects into account simultaneously. The many-body revivals presented in figure 3 provide additional tools for such studies.

Discussion
We have observed how the dephasing of an imbalanced pair of 1D Bose gases can result in states which are, for all practical purposes, indistinguishable from thermal equilibrium. This is due to a coupling of the relative and common degrees of freedom that is mediated by the relative dephasing of the individual gases. It is important to note that this observation of an apparent thermalization relies on the thermal-like initial conditions that were imposed on the system by the coherent splitting process. The system always retains a strong memory of the initial conditions and thus has not truly reached global thermal equilibrium. For example, if the system was initialized with other non-thermal initial conditions like the ones demonstrated in [10], it would equilibrate, but never appear thermal in its correlation functions [43].
Interestingly, the observed dynamics are closely related to the measurement process. In experiments, fluctuations of the anti-symmetric degrees of freedom are probed. These degrees of freedom exhibit a rapid relaxation with a single timescale if there is no imbalance, and a relaxation with two distinct timescales if there is imbalance. The same timescales govern the relaxation of the symmetric degrees of freedom. In contrast to that, if the properties of a single gas were accessible in experiment, their individual correlations would already look completely relaxed after the first, rapid timescale. This highlights how even in integrable systems, observables need to be properly aligned with (i.e. chosen such that they are sensitive to) the integrals of motion to reveal the integrable nature of the complex many-body dynamics. We note that the experiment in [24] recently revealed related behavior, where many-body revivals could be observed in certain correlation functions but not in others. This points to a general connection between the choice of measurement basis and the observed relaxation dynamics and will thus be an interesting topic for future research. For small imbalances we can assume r r r where  1,2 are time-independent constants. This expression describes a dephasing of the cross terms to 0 with the velocity difference c − .

Appendix D. Interference contrast
In experiments the expectation value in the definition of the phase correlation function is realized through an average over many experimental runs. We have previously demonstrated [13] that the number of runs that is required for a statistically meaningful determination of the correlation functions can dramatically increase for evolution times t?10 ms. However, for the parameters used in this work, the fully thermal state is not expected beforet 100 ms or more. For longer evolution times it could thus be beneficial to probe the level of coherence between the two gases using the mean squared interference contrast á ñ C 2 of the matterwave interference pattern in time-of-flight [13,15,49]. It is well known from one of our previous experiments [49] that this procedure involves an unknown factor that describes the finite resolution and other spurious experimental effects. We thus suggest to extract the relative phase q from every longitudinal position z of the interference pattern and calculate the contrast via the identity  Here,  is the two-point phase correlation function discussed in the main text (equation (19)) and L is a length scale over which the interference pattern is integrated. Because of the integration this procedure can be more robust against statistical fluctuations than the phase correlation function alone. With the identity given above, it is straight forward to generalize our predictions for the dynamics to the contrast. An example is shown in figure D1.

Appendix E. Harmonic traps
The Luttinger Hamiltonian as written in equation (1) stays valid for inhomogeneous density profiles ρ i (z). For the calculation in the harmonic trap, we assume a Thomas-Fermi profile for the density distribution before splitting, and rescaled density distributions for the evolution of the fluctuations after the quench. For the latter, the initial density profile depending on the total atom number N 1 +N 2 and on the trap frequencies ω z (longitudinal) and w^, in (radial) is simply multiplied by the factor + ( ) N N N i 1 2 . With this density distributions the Hamiltonian can be diagonalized with the help of Legendre-Polynomials [30,51]. Note that due to the density dependence of the shot noise fluctuations introduced in the splitting process, the initial density fluctuations expanded in Legendre-Polynomials are not diagonal anymore [52].
The results of the calculation are shown in figure E1. The incommensurate excitation energies of the trapped system lead to very complex dephasing and rephasing dynamics. As already discussed in the main text, such complex dynamics make it very challenging to experimentally disentangle different competing integrable (such as the one presented here) and non-integrable (such as thermalization) mechanisms. Figure D1. Contrast as another experimental probe. The slow and fast revivals that have been identified for the phase correlation function are also observable using the mean squared contrast á ñ C 2 . Here, we have used the same boundary conditions ( f ¶ ¶ = z 0) and parameters as in figure 3(c), i.e. a system size of 100 μm, initial temperature of T in =50 nK, imbalance Δ=0.  , the perpendicular trap frequency is w p =2 1.4 ,f kHz. Note that we do not need to specify a longitudinal trap frequency for the evolution after the quench as we simply assume rescaled density profiles. The inset shows the corresponding evolution of the interference contrast. The coordinates z and ¢ z have been chosen symmetrically around the center of the trap in both (a) and (b).