Separation of Timescales in a Quantum Newton’s Cradle

Separation of Timescales in a Quantum Newton’s Cradle

R. van den Berg Institute for Theoretical Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands
   B. Wouters Institute for Theoretical Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands
   S. Eliëns Institute for Theoretical Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands
   J. De Nardis Institute for Theoretical Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands
   R.M. Konik CMPMS Dept. Bldg 734 Brookhaven National Laboratory, Upton NY 11973, USA    J.-S. Caux Institute for Theoretical Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands
July 25, 2019

For strongly repulsive bosons in one dimension, we provide detailed modeling of the Bragg pulse used in quantum Newton’s cradle-like settings or in Bragg spectroscopy experiments. By employing the Fermi-Bose mapping for a finite harmonically trapped gas and the Quench Action approach for a thermodynamic system on a ring, we reconstruct the exact post-pulse many-body time evolution of Lieb-Liniger gases in the Tonks-Girardeau limit, together with their changing local density profile and momentum distribution. Our results display a clear separation of timescales between rapid and trap-insensitive relaxation immediately after the pulse, followed by slow in-trap periodic behaviour.


The study of many-body quantum physics has in recent years been transformed by the progress achieved in experiments on ultracold atoms Bloch et al. (2008). The context of one-dimensional (1D) bosonic gases in particular provides a fertile ground for investigating physics beyond traditional paradigms Cazalilla et al. (2011), with concepts such as the Luttinger liquid and exact solvability Giamarchi (2004) playing a primary role.

One of the main experimental probes of cold gases is Bragg spectroscopy Martin et al. (1988); Stenger et al. (1999); Ozeri et al. (2005), which consists in applying a pulsed monochromatic laser grating onto the gas, thereby creating excitations at (multiples of) the recoil momentum , where is the wavevector of the laser. The precise time-dependent form of the pulse can be chosen to optimize (de)population of specific quantum states. In Wang et al. (2005); Wu et al. (2005), a two-pulse sequence was used to optimize the population of the first satellites as compared to the zero-momentum ground state of a Bose-Einstein condensate. The theoretical description of such a sequence relied on a two-state model in which many-body dynamics are not included. In one dimension however, many-body effects are inescapable. One of the fundamental models in this context is the Lieb-Liniger gas Lieb and Liniger (1963) of bosons in a 1D continuum interacting with a -function potential, giving a proper description of atoms in tight transverse confinement Olshanii (1998). This model is relevant to the description of experiments, most prominently the famous quantum Newton’s cradle experiment Kinoshita et al. (2006), in which a Bragg pulse is used to initiate the oscillations. Bragg spectroscopy has also recently been used to investigate correlated 1D Bose gases of rubidium Fabbri et al. (2015) and cesium Meinert et al. (2015). In these, the heating of the gas resulting from the Bragg pulse was measured and matched using linear response to theoretical calculations of the dynamical structure factor of the Lieb-Liniger gas at finite temperature Panfil and Caux (2014).

Our main objective here is to model the effects of Bragg pulses theoretically for correlated 1D Bose gases, from first principles and without approximation (and thus beyond linear response), for experimentally relevant setups. We focus on the Tonks-Girardeau limit Tonks (1936); Girardeau (1960); Haller et al. (2009) of strongly repulsive Lieb-Liniger bosons both on a periodic interval and in confining traps Minguzzi and Gangardt (2005); Pezer and Buljan (2007); Gangardt and Pustilnik (2008); Muth et al. (2010); Schenke et al. (2011); Astrakharchik and Pitaevskii (2013); Collura et al. (2013); Quinn and Haque (2014); Cartarius et al. (). Instantaneous Bragg pulses of varying amplitude and wavevector are studied via their effect on physical observables: the time-dependent local density of the gas, and the experimentally more easily accessible momentum distribution function (MDF).

We start by modeling the Bragg pulse as a standing wave forming a one-body potential that couples to the density , where the Bose fields obey the canonical equal-time commutation relations, . For a general Bragg pulse the gas is perturbed for a finite duration . We will however consider the regime where the motion of the particles during the pulse can be neglected (the Raman-Nath limit), in which case the Bragg pulse is also referred to as a Kapitza-Dirac pulse  Kapitza and Dirac (1933); Freimund et al. (2001). Taking the limit such that is kept finite, the Bragg pulse operator is given by


where we have used the convention . Applying such an instantaneous Bragg pulse to a ground state yields a post-pulse state , which can be interpreted as the initial state of a quantum quench Calabrese and Cardy (2006); Rigol et al. (2007); Polkovnikov et al. (2011). Typical experimental pulses Kinoshita et al. (2006); Sapiro et al. (2009); Fabbri et al. (2015); Meinert et al. (2015) correspond to Bragg momentum and , where is the mean density.

The post-pulse time evolution is driven by the Lieb-Liniger (LL) model of interacting bosons


In what follows we will focus on the hard-core Tonks-Girardeau (TG) limit Tonks (1936); Girardeau (1960) and we will consider the model both on a ring (periodic boundary conditions) and in a parabolic trapping potential.

In the hard-core limit, the bosonic many-body wavefunction can be related through the Fermi-Bose (FB) mapping Girardeau (1960) to the many-body wavefunction of free fermions , where and the fermionic wavefunction is the usual Slater determinant of the free single-particle (SP) wavefunctions, . Following Pezer and Buljan (2007); Cartarius et al. (); Collura et al. (2013), the bosonic one-body density matrix, defined as , is given in terms of a single determinant involving the time-dependent fermionic SP states. This allows for an efficient computation of the MDF .

Starting with the ring geometry, our ground state consists of SP plane waves, on which the Bragg pulse imprints a cosine phase due to the one-body potential,


with ground-state rapidities forming a Fermi sea with Fermi momentum . Note that the Bragg momentum is quantized due to the periodic boundary conditions of the ring, with . Expanding Eq. (3) in plane waves, the time-dependent SP wavefunctions after the Bragg pulse yield


where is the modified Bessel function of the first kind.

Contrary to the finite-size Fermi-Bose mapping, the Generalized Gibbs Ensemble (GGE) Rigol et al. (2007, 2008) and the Quench Action (QA) approach Caux and Essler (2013); De Nardis et al. (2014) enable the study of the Bragg pulsed system (on a ring) in the thermodynamic limit ( with fixed). The GGE can be constructed using the infinite number of conserved charges provided by the integrability of the LL model, with , and eigenvalues associated to a Bethe state . The expectation values of the charges on the initial post-pulse state can be computed using the overlaps , which can be derived from the matrix elements for the Bragg pulse between two Bethe states and Sup ()


where we defined . Taking the thermodynamic limit, the energy density of the system after the Bragg pulse yields . The GGE logic Rigol et al. (2007, 2008) then requires the expectation values of all charges to be reproduced by the equilibrated post-pulse system, described by a density of rapidites , i.e.


for all . This leads to the following stationary-state distribution Sup (); Bevilacqua et al. (2011),


where is the Heaviside step function. The saddle point distribution is a sum of copies of the ground-state density of rapidities, , shifted by multiples of and weighted by the modified Bessel functions.

This form of the stationary state is consistent with the QA approach Sup (), which furthermore states that the time evolution of local observables after a quantum quench is given by a sum over particle-hole excitations in the vicinity of  Caux and Essler (2013); De Nardis et al. (2014, ). One easily obtains the time evolution of the density of the gas,


with the Bessel function of the first kind. The time evolution of the density is compared to the FB result for in Fig. 1 and shows excellent agreement, with relative differences of order due to finite-size effects in the FB calculations. Throughout the paper, all data is produced setting . As a consequence of the Raman-Nath limit, the initial post-pulse density (at ) is unaltered from the flat ground state profile. A sharp density profile then develops, mimicking the one-body cosine potential that was instantaneously turned on and off, followed by relaxation back to a flat profile at time scales , with the sound velocity.

Figure 1: Time evolution of the density after a Bragg pulse with and , computed by (a) the FB mapping and (b) the QA approach. The relative differences between the two results due to finite-size effects are less than .

The QA approach also provides access to the time evolution of the MDF De Nardis and Caux (2014); Sup (). The result is plotted in Fig. 2 along with the FB result for . Except for minor disagreements in the sharp peaks due to finite-size effects, the large-system-size dynamics after the Bragg pulse is again well captured by a Fermi-Bose mapping for particles. At , using the commutation relations between the Bragg operator and the Bose fields, one can show that the momentum distribution is simply a sum of copies of the ground-state MDF Sup (), with a small- divergence , centered around multiples of . The sharply peaked MDF then relaxes to a characteristic ghost-like shape Kinoshita et al. (2006), with the mixing of particles with different momenta causing a substantially broadened stationary MDF.

In Fig. 3 the equilibrated MDF is shown for different values of and . Similar to the initial MDF, the late-time distribution behaves like a superposition of independent peaks shifted to multiples of . The width of each satellite shows no dependence on the value of , and is only influenced by the choice of . Since in the limit of the resulting MDF is just that of the ground state and would stay constant as time progresses, the broadening can be ascribed to interactions between particles belonging to different copies of the ground state density of rapidities in .

Figure 2: Time evolution of the MDF after a Bragg pulse with and , computed with the QA approach (left half) and the FB mapping (right half). Because the FB mapping treats a finite system () the momenta are quantized, causing less pronounced peaks for short times. All other results are in excellent agreement with the QA computations.
Figure 3: The relaxed MDF function for particles as a function of the Bragg momentum . In the upper panel , and in the lower panel is fixed to .

Next, we will use the FB mapping to investigate how these observations translate to the more experimentally relevant geometry of a harmonic trapping potential, with the Hamiltonian and the trapping frequency. The ground state SP harmonic oscillator wavefunctions are given by


for , with denoting the Hermite polynomials. Similar to Eq. (3), acting on these states with the Bragg operator leads to an additional phase of the one-body cosine potential in the SP wavefunctions. Using the propagator for the quantum harmonic oscillator (the Mehler kernel) Mehler (1866), we compute the time evolution of the SP wavefunctions Sup ():


As a consequence of the Tonks-Girardeau limit, the SP wavefunctions are periodic in time with period , which is reflected in observables such as the density and the MDF. This periodicity is expected to be broken by finite- interactions and anharmonicities in the trapping potential. The time evolution of the density and the MDF during one period is shown in Fig. 4, where the contributions from particles belonging to different momentum satellites are clearly distinguishable. During the initial stages of relaxation (and around multiples of ) the density shows strong oscillations and the initially sharply peaked MDF relaxes rapidly to a more broadened shape. This prerelaxation is well separated from the collective periodic motion due to the trap, suggesting that it is governed by the same physics as relaxation on a ring.

Figure 4: The time evolution of the density (top) and MDF (bottom) in the trap, computed with the FB mapping for , , and .

In Fig. 5 the density at early stages in the oscillation cycle is compared with the density dynamics on a ring, where the latter was supplemented by a local density approximation (LDA) to account for the classical expansion of the gas in the trap  Campbell et al. (2015); Sup (). The initial density profile is accurately reproduced by the LDA, except for small differences near the edges originating from the gradient in the local density not accounted for by the LDA Kheruntsyan et al. (2005); van Amerongen et al. (2008); Armijo et al. (2011); Jacqmin et al. (2011). Note however that these differences do not stay confined to the edges and propagate towards the center as time progresses.

The short-time MDF in the trap and ring geometry is shown in Fig. 6 up to . The initial distributions are nearly identical, after which the MDFs dephase in a similar fashion to a (pre)relaxed ghost-like shape. The strong similarities can be attributed to the short-range correlations characterizing the post-quench steady state, for which the one-body density matrix decays exponentially as a function of the distance between the particles. Large-distance effects due to the trap geometry lead to discrepancies only at low momenta . We conclude that the short-time dynamics in a trap closely resembles the dephasing on a ring and is thus governed by the physics of hard-core interactions. The time scale associated to this (pre)relaxation is much shorter than the collective oscillations in the trap. Considering conditions similar to the Newton’s cradle experiment, we estimate the short time scale to be of the order of s. Since this estimate is of the same order of magnitude as the pulse durations used in Kinoshita et al. (2006), an interesting next step would be to extend the study of the effects of interactions to longer pulses in the Bragg regime.

Figure 5: Time evolution of the density in a trap, computed with: the FB mapping for particles (left), the QA approach on a ring with an LDA accounting for the trap (middle). The difference between the two results is shown on the right. The Bragg pulse parameters are set to and with a trapping frequency .
Figure 6: Time evolution of the MDF for the trap geometry (solid lines) and the ring geometry (dashed lines), obtained with the FB mapping for particles. The trapping frequency is set to , and the Bragg parameters are given by and . The time step is set to .

An interesting open question is how these observations extend to finite- interactions. Away from the TG limit, the initial post-pulse MDF remains a weighted sum of copies of the ground-state MDF Sup (), leading to a decreasing spread of momenta in the satellites centered around multiples of , as one goes from the hardcore limit to the BEC limit (). Although finite- dynamics is currently unattainable, we expect the time scale of the rapid relaxation to grow for smaller , proportional to the inverse of the sound velocity, .


In summary, we have developed a theoretical description of the Bragg pulse for one-dimensional Bose gases and shown that the time evolution of physical observables for a Bragg pulsed Tonks-Girardeau gas in a trap is characterized by two well-separated time scales. The shortest time scale is dominated by the trap-insensitive hard-core interactions and causes a substantial broadening of the momentum distribution function well before the collective motion due to the presence of the trap sets in. Our work opens up the possibility to study the influence of interactions on more general pulse protocols and their detailed effects on experimentally relevant observables.

We thank M. Brockmann, N.J. van Druten, V. Gritsev, F. Meinert, H.-C. Nägerl, J. Schmiedmayer, F.E. Schreck, D. Weiss and J. van Wezel for useful discussions. This work was supported by the Netherlands Organisation for Scientific Research (NWO) and the Foundation for Fundamental Research on Matter (FOM), and forms part of the activities of the Delta-Institute for Theoretical Physics (D-ITP). This research was done in part under the auspices of the CMPMS Dept. at Brookhaven National Laboratory, which in turn is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, under Contract No. DE-AC02-98CH10886. We are grateful for support from the Centre de Recherches Mathématiques of the U. de Montréal, where this work was completed.


Supplemental Material

.1 Matrix elements and initial state on the ring

The initial state after the Bragg pulse is easily obtained from the matrix elements of the operator


In the Tonks-Girardeau (TG) limit on a ring geometry we use the eigenstates


with wavefunctions given by


By commuting through the creation operators, the matrix elements can be expressed as


The expectation value of the bosonic operators conspires with the signs in the Tonks-Girardeau wavefunctions leading to a determinant of -functions. Treating the coordinates as dummy variables under the integral sign, it is easy to rewrite the integral in factorized form as a determinant of integrals of the form


where we define and where we used that and lie on the momentum lattice with . This results in the matrix elements


The initial state


is easily expressed in the Tonks-Girardeau eigenbasis using the matrix elements .

.2 The stationary state on a ring from a GGE and the Quench Action approach

In order to implement the GGE logic Rigol et al. (2007, 2008), one starts with computing the conserved charges on the initial state. Let us focus on the case , for which the overlaps coming from Eq. (16) reduce to a simple product of modified Bessel functions. While odd charges are trivially zero, for the even charges we find at finite system size


where we defined and where the coefficients come from the sum over the order of the Bessel functions and are known recursively Bevilacqua et al. (2011). The sum over particles can be performed, after which the thermodynamic limit can be taken,


where is the average particle density. For example, the energy density pumped into the system by an instantaneous Bragg pulse is given by


One can show that the saddle-point density


reproduces these values of the conserved charges, i.e. for all , by performing the integral and recasting the infinite sum into the coefficients . Since the local conserved charges (if well defined) uniquely determine the saddle point, we have thus found the saddle-point density after a Bragg pulse for . For smaller Bragg momenta the computation becomes considerably more difficult due to the determinant structure of the overlaps, but one can show that the saddle-point density given in Eq. (21) is still correct.

The Quench Action (QA) approach Caux and Essler (2013); De Nardis et al. (2014) reproduces this saddle-point density for . As a consequence of working in the Tonks-Girardeau regime, there are many microstates with exactly the same overlap. We can rephrase the overlaps as


where is the number of rapidities with and . In the thermodynamic limit these numbers are given by


The normalized overlap coefficients have a well-defined thermodynamic limit,


where is the Heaviside step function and we used that . Furthermore, in the thermodynamic limit , where is the average particle density. The noncontinuous integrand will serve as the driving term of the GTBA equations. Note that in the second line we implicitely assume that when for any . The reason is that for Bethe states that do not obey this condition, the overlap is exactly zero (rapidities will never end up in those regions) and therefore . These states are therefore infinitely suppressed in the Quench Action saddle-point equations. Another way of seeing this is that originally the functional integral in the quench action approach is a sum over states with non-zero overlaps and these states are not in that sum.

Even when you restrict the support of the density function to these intervals, in this ensemble of states there are still many microstates that have zero overlap with the Bragg-pulsed ground state. The reason is that when a rapidity has moved to an interval , it is not in the other intervals and therefore leaves a hole there. This alters the usual form of the Yang-Yang entropy significantly. Given the fillings , the finite-size entropy is


which leads to a modified Yang-Yang entropy for the Bragg pulse from the ground state,


where we used the Tonks-Girardeau Bethe equation . The variation of the quench action should be restricted to densities for which when for any . Also, a Lagrange multiplier is added to fix the particle density to . The resulting GTBA equation is not an integral equation because of the Tonks-Girardeau limit,


This is solved by the normalized saddle-point density of Eq. (21), with . For smaller Bragg momenta the derivation of the saddle-point distribution using the QA approach remains an open problem, since the determinant structure of the overlaps prevents obtaining a straightforward thermodynamic limit of the overlap coefficients that is expressible in terms of a root density .

Moreover, one could question whether the time evolution of simple observables obtained from the QA approach using the saddle-point density of Eq. (21) is valid also for small Bragg momenta . However, the analysis of the FB mapping does not show any qualitative differences for Bragg momenta smaller or bigger than and the agreement with the time evolution of the QA approach is excellent for all . It therefore seems safe to assume that the time evolution from the QA approach is valid for all Bragg momenta.

.3 Time evolution of the momentum distribution function on a ring in the thermodynamic limit

In Ref. De Nardis and Caux (2014) the thermodynamic limit of the matrix elements of the the one-body density matrix between states with a countable number of particle-hole excitations on a thermodynamic state was computed. The result is obtained by decomposing the elements into a Fredholm determinant and a finite-size determinant accounting for the excitations. The result is


where is the density of rapidities of the thermodynamic state. The Fredholm determinants are denoted by Det, where . The kernels are given by


where is the density. The function is defined as


or equivalently via the integral equations


The QA approach yields the following expression for the time evolution of the one-body density matrix in the thermodynamic limit


where the effective densities of holes and particles are given by


with the density of holes of the saddle point given by . The differential overlap for a single particle-hole is obtained by taking a finite-size realization of the saddle point state and the modified state obtained by performing a single particle-hole on . The ratio of the two overlaps in the thermodynamic limit gives the differential overlap


The same argument gives the energy of the single particle-hole excitations


Since the overlaps only couple states such that the difference of their rapidities are multiples of , the sum over the excitations on the saddle point state reduces to


After this substitution the sum in Eq. (.3) becomes the definition of a Fredholm determinant, and using the definition of the saddle-point distribution in Eq. (21), the expression for the time evolution of the one-body density matrix can be rewritten as the difference of two Fredholm determinants of two infinite block matrices where each block and for any is an operator acting on ,


with the operators given by


Here and is the identity operator. The coefficients and the function are given by


In order to obtain the time evolution of the momentum distribution one needs to restrict the sum in Eq. (.3) to excitations with zero total momentum, namely


where we denoted the Fourier transform as . Using the identity


we obtain


.4 Time-evolved single-particle states in the trap

The propagator for the quantum harmonic oscillator (Mehler kernel) is given by


The single particle (SP) wavefunctions after the Bragg pulse can then be time evolved by integrating the initial wavefunctions including the cosine phase with the propagator:


where are the groundstate harmonic eigenfunctions


The result in Eq. (51) has been obtained by using the following two identities


.5 Exact momentum distribution at for arbitrary interactions

The one-body density matrix at (after the Bragg pulse) is given by