Quench dynamics of a Tonks-Girardeau gas released from a harmonic trap

Quench dynamics of a Tonks-Girardeau gas released from a harmonic trap

Mario Collura, Spyros Sotiriadis, and Pasquale Calabrese Dipartimento di Fisica dell’Università di Pisa and INFN, 56127 Pisa, Italy
July 9, 2019

We consider the non-equilibrium 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 time-dependence of the entanglement entropy which attains a very simple form in the stationary state.

1 Introduction

Recent experiments on trapped ultra-cold 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 non-integrable 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 non-translationally 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 non-equilibrium 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 non-integrable 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 one-dimensional space, which has the advantage to avoid unwanted finite-size 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.

Figure 1: Sketch of the trap release dynamic in a circle.

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 quasi-periodic). 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 Lieb-Liniger Bose gas [89] (a preliminary analysis for non-integrable models has also been presented [90]). Using this approach, in Ref. [87], it has been shown numerically that for a Lieb-Liniger 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 Tonks-Girardeau 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 two-point 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 density-density correlation and in Sec. 6 the bosonic one-particle 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 Lieb-Liniger model describes a system of identical bosons in one dimension (1D) interacting via a pairwise Dirac-delta potential. In first quantization language, the Hamiltonian is given by [89]


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 strong-coupling regime of impenetrable bosons , known as the Tonks-Girardeau 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 non-linear Schrodinger equation


where and are the bosonic annihilation and creation field operators respectively.

The Lieb-Liniger model is Bethe ansatz integrable [89], but the analytic calculation of the non-equilibrium 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


and where the infinite coupling is encoded in the hard-core 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 anti-commute 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 Jordan-Wigner transformation


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) non-local, 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 many-body state.

2.1 The initial state

The initial state we consider is the ground state of the Tonks-Giradeau gas in a harmonic confining potential, i.e. the ground state of the Hamiltonian


with and for . The translationally invariant Lieb-Liniger model is recovered for . In the Tonks-Giradeau limit, the corresponding free fermionic Hamiltonian is


The many body-ground state is the Slater determinant built with the lowest energy one-particle eigenfunctions. This is easily worked out from the diagonalization of the single-particle Hamiltonian


Let us assume for the moment that so that the eigenfunctions of , for the parabolic potential , are the ones of the one-dimensional harmonic oscillator


with the Hermite polynomials and a non-negative integer number. Introducing now the fermionic operators as


satisfying the canonical anti-commutation relations , the many-body Hamiltonian is diagonal in the representation


Clearly, all previous results remain valid for any external potential as long as one uses the corresponding eigenfunctions of the one-body Hamiltonian (7).

In Fock space, the many-body ground state of impenetrable bosons in a parabolic trap is


where is the vacuum state annihilated by for all .

Let us now consider the initial density profile


The sum over can be analytically carried out using the Christoffel-Darboux formula for the Hermite polynomials


which in the limit leads to


Thus the TD initial density profile is


which is the well-known Thomas-Fermi profile (straightforwardly obtained for free fermions also by local density approximation). Notice that for larger than the Thomas-Fermi radius the gas density is exactly zero in the TD limit. Also the trapped fermionic two-point correlation function is straightforwardly obtained from Christoffel-Darboux formula


with the single particle wave-function in Eq. (8). When the two points are very close to the center of the trap, i.e. the above formula simplifies to


which is the translationally invariant result with .

The vanishing of the density and of the many-body wave-function, 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 :


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 many-body ground state in infinite space (11) is also the ground state for a finite circular geometry in the presence of the trap.

The trap-release 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


and the trap-release condition becomes


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 non-equilibrium dynamics which is the focus of the paper. The initial state is in Eq. (11) and the Hamiltonian governing the evolution for is the Tonks-Girardeau 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


which is diagonalized by Fourier transform in terms of the free fermionic operators and (with and integer)


The time evolution of an observable is obtained in a standard way:

  1. Write the desired observable in terms of the post-quench mode operators , whose time-evolution, in Heisenberg representation, is

  2. Write the post-quench mode operators as a function of the pre-quench operators , whose action on the initial state is trivial.

The relation between pre-quench and post-quench mode-operators can be written as


where we introduced the overlap between the pre-quench and the post-quench one-particle eigenfunctions


The inverse relation is


where the term is due to the different domain of integration between the pre-quench and the post-quench 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 one-dimensional harmonic oscillator, i.e.


A very important quantity for the non-equilibrium 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 single-particle occupied level in the initial state with energy . In terms of the post-quench Hamiltonian with single particle spectrum , corresponds to an initial Fermi-momentum . Since , we have for the Fermi velocity


Notice how the expansion velocity remains finite in the proper TD limit with constant.

3 The two-point fermionic correlation function

The easiest observable that we can calculate is the two-point fermionic correlator


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 one-particle wave functions the fermionic correlator is


which is a well-known result for Slater determinants.

The time evolved one-particle wave functions are the solutions to the Schrödinger equation


where is the single particle Hamiltonian with PBC. In terms of the overlaps in Eq. (27) the solutions to this equation read


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 time-averaging only the latter terms survive to give


where we defined


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 one-particle states


(we use instead of to avoid confusion with the state ). Notice that, because of our normalization, the momentum operator acts on free-waves as .

Figure 2: Pictorial representation of the solution of the trap release dynamics in a ring as a superposition of replicas of the infinite-space time-evolved function periodically shifted in space. See A for the mathematical derivation.

Let us first compute :


These matrix elements can be calculated using (with the harmonic oscillator ladder operators) and the Baker-Campbell-Hausdorff formula:


where is the hypergeometric function. To evaluate in the TD limit we need to calculate the coefficients of the powers of with , and


The large limit of the expression in the last line in the round parenthesis is just and therefore we obtain


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


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


Thus, unlike which is finite in the TD limit, decays to zero as giving just a finite-size correction to the correlation function.

The calculation of the last term is straightforward and in the TD limit we have


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


3.2 The time dependent one-particle problem

As detailed in A, Fourier analysis allows us to rewrite the one-particle evolution in Eq. (32) in terms of the time evolved wave function in the infinite-space


This formula is valid for any one-particle time-dependent 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


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).

Figure 3: Color plot of the numerically calculated density evolution for (from left to right) at and as a function of the rescaled space variable and for rescaled times
Figure 4: In each panel we plot the time evolution of the density as function of the rescaled time at fixed : for different (large enough) sizes the curves collapse on top of each other. Dashed red lines indicate the equilibration value reached at infinite time. The symbols are the exact dynamics [cf. Eq. (48)] for finite , while full black lines are the TD limit in Eq. (50).
Figure 5: In each panel we report the density profile as function of for fixed : again the results for various sizes collapse on a single curve when increases. Symbols are the exact dynamics for finite [cf. Eq. (48)], while full black lines are the TD limit in Eq. (50). As the time increases, the profiles tighten close to the equilibration value (note the vey different vertical scale in the four panels).

3.3 The time evolution of the density profile

We start the time-dependent analysis of the many-body problem from the diagonal part of the fermionic correlation function, i.e. the particle density (both for fermions and bosons). Plugging the one-particle wave-functions (47) and (44) into Eq. (30) we obtain


which is an exact formula for any .

The physical interpretation of Eq. (48) is again the one suggested by Fig. 2. The wave-functions 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 finite-size correction going like . Thus, in the TD limit, the leading behavior of the time-dependent density profile is given by


To perform the sum over , we use the Christoffel-Darboux 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


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 infinite-time limit of Eq. (50) is straightforward and gives the expected result .

Figure 6: Top: Color plot of the numerically calculated fermionic correlation function for considering both the full replica sum on (Full) in Eq. (51) and only the diagonal part (Diag.) in Eq. (52). We fix and . Notice how the differences between the full and the diagonal calculation, visible for , disappear already for . Bottom: For we report the correlation function (calculated as sum over only the diagonal term in Eq. (52)) for several values of TD parameters. From left to right: (a) ; (b) ; (c) ; . Notice (i) the slope of the signal lines (yellow) is equal to integer multiples of (ii) the difference in the time-scale (at the bottom of the figure) which depends on .
Figure 7: Snapshots of the correlation at different rescaled times and sizes. The data for different large enough sizes nicely collapse on top of each other. In panel (a) for , the full line is the initial correlation in the TD limit (cf. Eq. (17)). Panels (b) and (c) show that, as time increases, two symmetric peaks are expelled from the central region. In the panel (c) for , the full line is the stationary value plus the first-order correction in Eq. (59) which correctly describes the position of the two moving peaks, but not their amplitudes. For comparison, in panel (c) we report also the leading contribution for infinite time (red dashed line). In panel (d) with , in the considered spatial region , the time evolved data are almost indistinguishable from the stationary values (full line).

3.4 The time evolution of the two-point fermionic correlation and its large-time limit

The calculation of the time evolution of the two-point 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


which once again can be easily interpreted in terms of replicas.

Figure 8: Time dependence of the correlation for (a,c) and (b,d). The dynamics is rather irregular before the propagating peak travels the distance (in a time ). For later times instead there is a simple damped oscillatory behavior around the stationary value (dashed line in top panels) which is zoomed in the insets for clarity. In (a,b) the points represent the full correlation function in Eq. (51) while the full lines are the diagonal sum in Eq. (52). In the two bottom panels (c,d), the points represent from the diagonal sum in Eq. (52) for larger time scales. They are compared with the asymptotic expansion in Eq. (60) (full lines) which is almost indistinguishable from the data soon after the moving peak passed through. Also the asymptotic behavior for the envelopes of maxima and minima is reported (dashed lines cf. Eq. (61)).

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