Quantum shock waves and population inversion in collisions of ultracold atomic clouds

Using a time-dependent density matrix renormalization group (TDMRG) approach we study the collision of one-dimensional atomic clouds confined in a harmonic trap and evolving with the Lieb-Liniger Hamiltonian. It is observed that the motion is essentially periodic with the clouds bouncing elastically, at least on the time scale of the first few oscillations that can be resolved with high accuracy. This is in agreement with the results of the “quantum Newton cradle” experiment of Kinoshita et al. [Nature (London) 440, 900 (2006)]. We compare the results for the density profile against a hydrodynamic description, or generalized nonlinear Schrödinger equation, with the pressure term taken from the Bethe ansatz solution of the Lieb-Liniger model. We find that hydrodynamics can describe the breathing mode of a harmonically trapped cloud for arbitrary long times while it breaks down almost immediately for the collision of two clouds due to the formation of shock waves (gradient catastrophe). In the case of the clouds' collision TDMRG alone allows one to extract the oscillation period which is found to be measurably different from the breathing mode period. Concomitantly with the shock waves formation we observe a local energy distribution typical of population inversion, i.e., an effective negative temperature. Our results are an important step towards understanding the hydrodynamics of quantum many-body systems out of equilibrium and the role of integrability in their dynamics.

Figure 1

Comparison between TDMRG data (black line) and the solutions of the hydrodynamic equation (4) [red (gray) line]. The top panels refer to the breathing mode quench while on the bottom ones the collision of clouds is illustrated. In each panel several different values of $g_B$ are shown, from free bosons $g_B=0$ (bottom) up to the Tonks-Girardeau limit $g_B=+\infty$ (top). Density profiles for different interaction strengths have been shifted vertically for clarity. Time is in units of the half period defined as $\tau =\pi /{\omega }_1$ for the collision quench and as $\tau =\sqrt{3}\pi /{\omega }_1$ for the breathing quench where the value of ${\omega }_1$ is given in Table 1. In the case of the breathing mode the perfect match between TDMRG and hydrodynamics (the density profiles are overlapping) is a strong indicator of the accuracy of the simulations. In the case of the collision quench, shock wave formation is signaled by steep density profiles at $t\sim 0.35\tau$ and the growth of oscillations present only in the hydrodynamic profiles. The amplitude of these oscillations increases with $g_B$ and they are very evident for $g_B=2.0Ja$ and in the Tonks-Girardeau limit at $t=0.5\tau$. This leads to the breakdown of the hydrodynamic approximation, Eq. (4), which is unable to capture the true quantum dynamics where the ripples are smoothed out by quantum fluctuations.

Figure 2

Comparison between TDMRG (black line) and hydrodynamic [red (gray) line] density profiles at $t=0$ (bottom), after one oscillation period ${\tau }^{*}\left(g_B\right)$ (middle), and after two periods (top). The density profiles at different times have been shifted in the vertical direction. How the renormalized oscillation period has been extracted from the TDMRG data is explained in Sec. 3b and Fig. 3. Note that the hydrodynamic simulations match the TDMRG results for several oscillation periods in the case of the breathing mode while in the collision of clouds the approximation breaks down before a single oscillation is completed due to shock wave formation (see Fig. 1). In both cases the full quantum dynamics exhibit (quasi-)periodicity for any interaction strength (the data shown here refer to $g_B=1.0Ja$).

Figure 3

Upper plot: Deviation $\Delta \left(t\right)$ as a function of time [Eq. (7)]. The black dashed lines are the data relative to the Bose-Hubbard discretization and the light gray solid ones are relative to the $XXZ$ spin chain discretization. The arrows indicate the instants where the density is closest to the initial one [minima of $\Delta \left(t\right)$]. In the lower plot the renormalized oscillation period ${\tau }^{*}\left(g_B\right)/\tau$ is shown as a function of the interaction strength $g_B$ extracted from the minima of $\Delta \left(t\right)$ (see the upper plot). The light gray triangles (black circles) refer to the $XXZ$ (Bose-Hubbard) discretization. The gray squares are relative to the breathing period for which hydrodynamics and TDMRG agree. The light gray (gray) arrow at the bottom left represents the frequency extracted from the $g_B=0$ data for the collision (breathing) quench. They deviate from the exact result ${\tau }^{*}\left(0\right)=\tau$ in the continuum limit since lattice effects distort the density profile in time. As can be seen in Fig. 1, higher densities are explored in the $g_B=0$ case; thus the continuum limit approximation is less accurate.

Figure 4

(a) Snapshots of the density profile $\rho (x,t)$ for the clouds' collision as in Fig. 1 ($g_B=0.2Ja$, $\gamma =2$). (b) Energy distribution function $f(x,\varepsilon ,t)$ (9) calculated from the Wigner function (8). The color plots show the values of $f(x,\varepsilon ,t)$ in the $(x,\varepsilon )$ plane at the corresponding times in panel (a). The vertical black, gray, and light gray lines correspond to the sections for fixed $x$ of $f(x,\varepsilon ,t)$ shown in Fig. 5. Note that shock-wave formation is characterized by a highly nonequilibrium distribution. (c) and (d) Same as panels (a) and (b) but for $g_B=2.0Ja$ $(\gamma =20)$.

Figure 5

Energy distribution function $f(x,\varepsilon ,t)$ corresponding to the first ($t=0\tau$) and third ($t\sim 0.3\tau$) snapshots in Figs. 4 and 4(d). The black, gray, and light gray lines correspond to $x/a=25$, 40, and 60 for the two upper quadrants ($g_B=0.2Ja$) [Fig. 4], and $x/a=15$, 30, and 50 for the two lower quadrants ($g_B=2.0Ja$) [Fig. 4], respectively. Note that neglecting oscillations for large $\varepsilon$—due to the finite number of particles—the initial distribution decreases monotonically, while a maximum develops for finite $\varepsilon$ after the shock-wave formation at $t\sim 0.3\tau$, i.e., a population inversion.

Figure 6

Inverse information compressibility as a function of time for various interaction strengths $g_B/\left(Ja\right)=0.2$, 1.0, and 2.0. Note the divergence in correspondence to a fully developed shock wave at $t=0.4$–0.5$\tau$. See Fig. 4, snapshots at $t=0.38\tau$ in panel (a) and at time $t=0.32\tau$ in panel (c).