Quench dynamics of a TonksGirardeau gas released from a harmonic trap
Abstract
We consider the nonequilibrium dynamics of a gas of impenetrable bosons released from a harmonic trapping potential to a circle. The many body dynamics is solved analytically and the time dependence of all the physically relevant correlations is described. We prove that, for large times and in the thermodynamic limit, the reduced density matrix of any subsystem converges to a generalized Gibbs ensemble as a consequence of the integrability of the model. We discuss the approach to the stationary behavior at late times. We also describe the timedependence of the entanglement entropy which attains a very simple form in the stationary state.
1 Introduction
Recent experiments on trapped ultracold atomic gases have shown that it is possible to follow and measure the unitary nonequilibrium evolution of an isolated quantum system [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. A particular class of these nonequilibrium problems which is experiencing an enormous theoretical activity is that of a sudden quench of a Hamiltonian parameter. In a global quantum quench, the initial condition is the ground state of a translationally invariant Hamiltonian which differs from the one governing the evolution by an experimentally tunable parameter such as a magnetic field. In these experiments the two key questions are: i) how the correlations and entanglement spread into the system with time [12, 13, 14, 15, 16, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29] and ii) whether the system relaxes (in some sense) to a stationary state, and if it does, how to characterize from first principles its physical properties at late times. For the latter question, it is widely believed that, depending on the integrability of the Hamiltonian governing the time evolution, the behavior of local observables can be described either by an effective thermal distribution or by a generalized Gibbs ensemble (GGE), for nonintegrable and integrable systems respectively (see e.g. [12] for a review). This scenario is corroborated by many investigations [30, 31, 14, 15, 32, 33, 34, 35, 36, 37, 24, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58], but still a few studies [59, 60, 61, 62, 63, 64, 65, 66] suggest that the behavior could be more complicated. Indeed, it has been argued that the initial state can affect this scenario, in particular if it breaks some symmetries of the Hamiltonian governing the subsequent evolution which tend to be recovered in a statistical ensemble such as thermal or GGE. The case that has been most largely studied is that of a nontranslationally invariant initial state generically referred to as inhomogeneous quenches [67, 68, 69, 70, 71, 72, 73, 74].
Among the inhomogeneous quenches, a particularly relevant one which has been already experimentally realized (also in one dimension [8, 9]) is the nonequilibrium dynamics of a gas released from a parabolic trapping potential. A very interesting experimental finding is that the spreading of correlations is ballistic for an integrable system and diffusive for a nonintegrable one [9]. Both experimental [8, 9] and theoretical [64, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86] analyses concentrated on the expansion in the full onedimensional space, which has the advantage to avoid unwanted finitesize effects. However, if a repulsive gas expands on the full line, its density will decrease as time passes and for infinite time it goes to zero, making senseless to distinguish thermal and GGE states.
An alternative proposal by J.S Caux and R. Konik [87] is that of considering the release of a gas from a parabolic trap not in free space but on a closed circle of length (as sketched in Fig. 1), so that the gas has finite density even for infinite time. However, in a gas with a finite number of particles (or more generally in a system with a finite number of degrees of freedom) a stationary state cannot be approached because of revival and recurrence effects (i.e. the system is always quasiperiodic). To circumvent this, the thermodynamic (TD) limit should be defined properly: for fixed final density , when , at the same time the frequency of the initial confining potential should also vanish, i.e. , but with constant (i.e. fixed ‘initial density’) [87, 88]. In order to tackle quite generally this problem, a new approach based on integrability has been developed and applied explicitly to the LiebLiniger Bose gas [89] (a preliminary analysis for nonintegrable models has also been presented [90]). Using this approach, in Ref. [87], it has been shown numerically that for a LiebLiniger gas, the time averaged correlation functions are well described by a GGE, apart from finite size effects (the maximum number of particles considered in Ref. [87] is ). In principle this new approach allows also the study of the time evolution, but it is much more computationally demanding and it has not yet been done. As a consequence it has not yet been established whether (and in which sense) an infinite time limit exists and, if yes, how it is approached.
In order to overcome these limitations, in a recent letter [88] we presented a full analytic solution of the nonequilibrium trap release dynamics in the limit of strong coupling, i.e. in the TonksGirardeau regime [91]. This allowed us to understand that also the infinite time limit should be handled with care: in the trap release dynamics, a stationary behavior is possible because of the interference of the particles going around the circle many times (see Fig. 1), i.e. to observe a stationary value we must require (with the expansion velocity of the gas). This is very different from equilibration in standard global quenches where, in order to avoid revival effects, the time should be such that the boundaries are never reached (i.e. one first considers the TD limit and only after the infinite time limit , which, in finite systems, corresponds to the condition , see e.g. [39]). In the trap release problem the revival scale is (see also [92]) and so the infinite time limit in which a stationary behavior can be achieved is provided . In Ref. [88], we have showed that, in the TD limit, the reduced density matrix of any finite subsystem converges for long times (in the sense just explained) to the GGE one. This implies that any measurable local observable will converge to the GGE predictions. In this manuscript, we extend the previous letter [88] in several aspects. First of all, we give complete derivation of all results in the GGE previously presented. Secondly, for many observables we will characterize the full asymptotic time dependence and not restrict to the stationary results. As particularly important new aspects absent in Ref. [88], we study the time evolution of the entanglement entropies and we construct the GGE in terms of local integrals of motion.
The manuscript is organized as follows. In Sec. 2 we introduce the model under investigation and the quench protocol. In Sec. 3 we calculate the time evolution of the twopoint correlation function and prove that for infinite time a stationary value is approached. We also discuss the approach to the stationary value. In Sec. 4 we show that the stationary values of all local observables are described by a GGE both in fermionic momentum occupation numbers and in the local integrals of motion. In Sec. 5 we compute the densitydensity correlation and in Sec. 6 the bosonic oneparticle density matrix (Fourier transform of the momentum distribution function). In Sec. 7 we move our attention to the entanglement entropies. The trap release dynamics from a trapped gas to a larger trap is addressed in Sec. 8. Finally in Sec. 9 we draw our conclusions.
2 The Model and quench protocol
The LiebLiniger model describes a system of identical bosons in one dimension (1D) interacting via a pairwise Diracdelta potential. In first quantization language, the Hamiltonian is given by [89]
(1) 
where is the coupling constant and we set . For definiteness, we consider a system of length with periodic boundary conditions (PBC). In the repulsive regime, , and in the TD limit, the equilibrium physics of the model depends on the single parameter where is the particle density. Then in 1D, in stark contrast to higher dimensions, low densities lead one to the strongcoupling regime of impenetrable bosons , known as the TonksGirardeau limit [91]. In the attractive regime, , the physics of the model is completely different (see e.g. [93, 94]) and will not be considered here.
In second quantization language, the Hamiltonian (1) can be rewritten as a quantum nonlinear Schrodinger equation
(2) 
where and are the bosonic annihilation and creation field operators respectively.
The LiebLiniger model is Bethe ansatz integrable [89], but the analytic calculation of the nonequilibrium dynamics in the TD limit is still a formidable task, despite the numerous attempts in the literature [77, 87, 79, 92, 95, 96]. For this reason, as already anticipated, we concentrate here in the impenetrable limit in which the Hamiltonian (2) can be simply written as
(3) 
and where the infinite coupling is encoded in the hardcore constraint , i.e. the condition that two bosons cannot occupy the same position. At operator level, this constraint can be imposed by requiring that and commute at different spatial points and they anticommute when evaluated at the same point. In other words, they are similar to fermionic operators but they commute on different space positions. In order to restore a genuine Fermi algebra, fermionic field operators and are built through a JordanWigner transformation
(4)  
which by construction satisfy and , with the density operator of both fermions and bosons. This is the standard mapping between impenetrable bosons and free fermions [91] which ensures that all spectral and thermodynamical properties of the bosons can be simply obtained from free fermions. However, being the transformation (4) nonlocal, bosonic correlation functions are different from fermionic ones and they should be reconstructed with the help of Wick theorem, as explicitly done in the following.
The Hamiltonian (3) of impenetrable bosons is the one governing the time evolution in our problem and we have now to fix the initial manybody state.
2.1 The initial state
The initial state we consider is the ground state of the TonksGiradeau gas in a harmonic confining potential, i.e. the ground state of the Hamiltonian
(5) 
with and for . The translationally invariant LiebLiniger model is recovered for . In the TonksGiradeau limit, the corresponding free fermionic Hamiltonian is
(6) 
The many bodyground state is the Slater determinant built with the lowest energy oneparticle eigenfunctions. This is easily worked out from the diagonalization of the singleparticle Hamiltonian
(7) 
Let us assume for the moment that so that the eigenfunctions of , for the parabolic potential , are the ones of the onedimensional harmonic oscillator
(8)  
with the Hermite polynomials and a nonnegative integer number. Introducing now the fermionic operators as
(9) 
satisfying the canonical anticommutation relations , the manybody Hamiltonian is diagonal in the representation
(10) 
Clearly, all previous results remain valid for any external potential as long as one uses the corresponding eigenfunctions of the onebody Hamiltonian (7).
In Fock space, the manybody ground state of impenetrable bosons in a parabolic trap is
(11) 
where is the vacuum state annihilated by for all .
Let us now consider the initial density profile
(12) 
The sum over can be analytically carried out using the ChristoffelDarboux formula for the Hermite polynomials
(13) 
which in the limit leads to
(14) 
Thus the TD initial density profile is
(15) 
which is the wellknown ThomasFermi profile (straightforwardly obtained for free fermions also by local density approximation). Notice that for larger than the ThomasFermi radius the gas density is exactly zero in the TD limit. Also the trapped fermionic twopoint correlation function is straightforwardly obtained from ChristoffelDarboux formula
(16) 
with the single particle wavefunction in Eq. (8). When the two points are very close to the center of the trap, i.e. the above formula simplifies to
(17) 
which is the translationally invariant result with .
The vanishing of the density and of the manybody wavefunction, in the TD limit, for is the fundamental property allowing us to treat analytically also the time evolution in a finite circle of length . Indeed, we now make the only crucial physical assumption of our treatment: we impose that the space initially occupied by the trapped gas as a whole is within the external box of length , i.e. the PBC are irrelevant for the gas in the initial state which only “sees” the parabolic trap. This means that the extension of the gas in the trap, in Eq. (15), must be smaller than the box size :
(18) 
In terms of the number of particles , this condition means that must be smaller than the first level of the parabolic potential that is affected by the PBC. Furthermore, this is the hypothesis which allows us to talk about release of the gas, because if the gas would feel the PBC before the quench, we would not have a trapped gas, but something more complicated. Clearly, under the assumption (18), the manybody ground state in infinite space (11) is also the ground state for a finite circular geometry in the presence of the trap.
The traprelease condition can be also written in a maybe more transparent way in terms of the initial average density and final density . Indeed, by definition we have
(19) 
and the traprelease condition becomes
(20) 
i.e. that the initial average density is larger than the final one signaling that the gas expands.
2.2 The quench protocol
In this section we describe the nonequilibrium dynamics which is the focus of the paper. The initial state is in Eq. (11) and the Hamiltonian governing the evolution for is the TonksGirardeau in Eq. (3) clearly with periodic boundary conditions. In practical terms, this protocol is a quench of the trapping potential from a given to , i.e. a trap release at .
The Hamiltonian (3) in terms of the fermionic field operators is
(21) 
which is diagonalized by Fourier transform in terms of the free fermionic operators and (with and integer)
(22) 
The time evolution of an observable is obtained in a standard way:

Write the desired observable in terms of the postquench mode operators , whose timeevolution, in Heisenberg representation, is
(23) 
Write the postquench mode operators as a function of the prequench operators , whose action on the initial state is trivial.
The relation between prequench and postquench modeoperators can be written as
(24) 
where we introduced the overlap between the prequench and the postquench oneparticle eigenfunctions
(25) 
The inverse relation is
(26) 
where the term is due to the different domain of integration between the prequench and the postquench Hamiltonians.
All this derivation is completely general and it is the practical way we construct the exact time evolution for finite number of particles , finite and . However, the results greatly simplify in the TD limit if this is properly defined as follows. We should consider at fixed density and, at the same time, with constant (i.e. fixed ‘initial density’). This is exactly the same TD limit defined in Ref. [87]. In terms of these TD quantities the trap release condition reads . Under this condition, the functions entering in the definition of (i.e. with ) are exponentially small outside the interval and the mapping between the operators and in Eq. (26) is exact also with the integration domain in Eq. (25) extended to . Thus, if the trap release condition is satisfied, the overlaps are simply the Fourier transforms of the eigenfunctions of the onedimensional harmonic oscillator, i.e.
(27) 
A very important quantity for the nonequilibrium dynamics is the expansion velocity of the gas in full space. This is obtained straightforwardly from the analytic solution of the dynamics [75] which we will discuss later, but can be also simply written down from elementary arguments. Indeed this velocity is determined by the maximum energy singleparticle occupied level in the initial state with energy . In terms of the postquench Hamiltonian with single particle spectrum , corresponds to an initial Fermimomentum . Since , we have for the Fermi velocity
(28) 
Notice how the expansion velocity remains finite in the proper TD limit with constant.
3 The twopoint fermionic correlation function
The easiest observable that we can calculate is the twopoint fermionic correlator
(29) 
Indeed, since the Hamiltonian is quadratic in the fermionic operators, the evolved state is a Slater determinant and Wick’s theorem applies allowing to obtain (with some work as we shall see) all other observables. In terms of oneparticle wave functions the fermionic correlator is
(30)  
which is a wellknown result for Slater determinants.
The time evolved oneparticle wave functions are the solutions to the Schrödinger equation
(31) 
where is the single particle Hamiltonian with PBC. In terms of the overlaps in Eq. (27) the solutions to this equation read
(32) 
3.1 The time average
Let us first compute the time average of the fermionic correlation function since if a large time limit of Eq. (30) exists, it should be equal to its time average. To this aim, it is convenient to split the double momentum sum in Eq. (30) in a term with and another with . After timeaveraging only the latter terms survive to give
(33) 
where we defined
(34) 
in terms of the overlaps in Eq. (27).
We now calculate the TD limit of the three pieces , and separately. The writing is simplified by the use of the Dirac notation for the oneparticle states
(35) 
(we use instead of to avoid confusion with the state ). Notice that, because of our normalization, the momentum operator acts on freewaves as .
Let us first compute :
(36)  
These matrix elements can be calculated using (with the harmonic oscillator ladder operators) and the BakerCampbellHausdorff formula:
(37)  
where is the hypergeometric function. To evaluate in the TD limit we need to calculate the coefficients of the powers of with , and
(38)  
The large limit of the expression in the last line in the round parenthesis is just and therefore we obtain
(39) 
with the Bessel function. Notice , the initial Fermi momentum.
In order to evaluate we notice from Eq. (27) that . Thus the only differences compared to are (i) an extra sign and (ii) the replacement of with . Therefore the calculation of follows the same steps as for up to the middle line in Eq. (38), which now becomes
(40) 
but unlike , the piece in the round parenthesis does not grow like for large , but it has a finite limit because . Therefore, for , we have
(41) 
Thus, unlike which is finite in the TD limit, decays to zero as giving just a finitesize correction to the correlation function.
The calculation of the last term is straightforward and in the TD limit we have
(42) 
which also vanishes for large .
Summing up, the time average of the fermionic correlation function in the thermodynamic limit gets a non vanishing contribution only from and so it is
(43) 
3.2 The time dependent oneparticle problem
As detailed in A, Fourier analysis allows us to rewrite the oneparticle evolution in Eq. (32) in terms of the time evolved wave function in the infinitespace
(44) 
This formula is valid for any oneparticle timedependent problem with PBC and its physical meaning is very simple: in a circle of length with PBC, the time evolution is equivalent to the sum (superposition) of the time evolution of infinite copies (replicas labeled by in Eq. (44)) of the initial state in infinite space but shifted by integer multiples of (see Fig. 2 for a pictorial representation). For the particular case at hand we have
and using the following property of the Hermite polynomials
we get
and therefore
(47)  
which coincides with the result in Ref. [75]. The full time dependence in the ring is obtained by plugging the above equation (47) in Eq. (44).
3.3 The time evolution of the density profile
We start the timedependent analysis of the manybody problem from the diagonal part of the fermionic correlation function, i.e. the particle density (both for fermions and bosons). Plugging the oneparticle wavefunctions (47) and (44) into Eq. (30) we obtain
(48)  
which is an exact formula for any .
The physical interpretation of Eq. (48) is again the one suggested by Fig. 2. The wavefunctions of periodically placed replicas expand in infinite space and eventually overlap with each other when they reach a boundary between two replicas, that happens (approximately) at times which are integer multiples of the characteristic time . Therefore, if the boundaries are not reached and the system does not feel the PBC. When the time becomes much larger than , the overlap between replicas (or the turning of the particles around the circle) leads, as we shall prove, to equilibration manifested as a uniform distribution. However, as already stressed in the introduction, the infinite time limit should be handled with care because revivals will take place for larger time scales of the order of (because the fundamental frequency of the double momentum sum is , i.e. the revival period is , see also [92] for a general treatment). Therefore in the TD limit, the physically relevant scaling regime is obtained by taking and and consequently , so that the revivals are eliminated. In this regime the leading behavior of Eq. (48) is extracted by stationary phase arguments and it comes from the ‘diagonal replicas’ , i.e. the terms with . Indeed, in the TD limit with , the phase of the exponent in (48) is stationary only for , but the terms with give a finitesize correction going like . Thus, in the TD limit, the leading behavior of the timedependent density profile is given by
(49) 
To perform the sum over , we use the ChristoffelDarboux formula for the Hermite polynomials in Eq. (13) which in the limit leads to Eq. (14). Thus Eq. (49) can be written in terms of the particle density at initial time in Eq. (15) as
(50) 
showing that the density profile in the TD limit is simply given by the sum of the replicated densities and all the interference effects are subleading in , as probably expected.
In Figs. 3, 4 and 5 we show the numerically calculated exact time dependent density for finite but large . For large enough systems, the numerical data perfectly agree with the above TD prediction for any time. The infinitetime limit of Eq. (50) is straightforward and gives the expected result .
3.4 The time evolution of the twopoint fermionic correlation and its largetime limit
The calculation of the time evolution of the twopoint fermionic correlator is similar to the one just reported for the density. For finite , plugging Eq. (44) into Eq. (30) we have an exact starting point
(51)  
which once again can be easily interpreted in terms of replicas.
As for the density, a stationary phase argument allows us to conclude that only diagonal terms contribute to the TD limit. (To quantitatively support this statement, in Fig. 6 we compare the full sum with the one restricted over the diagonal terms: for small differences are visible, but they are negligible already for .) This leads to