Exact time evolution of space- and time-dependent correlation functions after an interaction quench in the 1D Bose gas

We consider the non-equilibrium dynamics of the interacting Lieb-Liniger gas after instantaneously switching the interactions off. The subsequent time evolution of the space- and time-dependent correlation functions is computed exactly. Different relaxation behavior is observed for different correlation functions. The long time average is compared with the predictions of several statistical ensembles. The generalized Gibbs ensemble restricted to a fixed number of particles is shown to give correct results at large times for all length scales.

With the recent advances in ultracold atomic gases it is now possible to realize isolated quantum systems with long coherence times.These experiments are ideal to study non-equilibrium physics [1][2][3].The outcome of these experiments initiated the ongoing debate concerning the circumstances under which an isolated quantum system will thermalize, which has led to many theoretical studies on non-equilibrium dynamics of quantum systems.Most of these studies considered a so-called quantum quench: an isolated system is prepared in the ground state of a Hamiltonian H(c), where c is some controllable parameter, like an interaction strength.At t 0 = 0 the parameter c is instantaneously switched to a different value resulting in unitary time evolution governed by the new Hamiltonian.Considerable progress has been made (for a recent review see [4] ), however there are still many open questions.Under what conditions can a system equilibrate and what is the relevant mechanism?What is the role played by integrability or its absence?Do different kinds of correlation functions show different relaxation behavior?And what is the importance of the initial state?
The difficulty with answering these questions is that one must rely on case by case studies, thus, one is faced with the problem of distinguishing universal behavior from case-specific results.Moreover, most methods used so far suffer from the fact that they are either valid for short times or for large ones.What we propose here is to study a specific example where we can study the behavior of correlation functions at all time scales, using an exact method.We can, therefore, not only determine what the equilibrium state is, if there is one, but also how it is reached.Furthermore, the method we use allows us to consider various types of correlation functions, thereby obtaining a more complete picture for this case.
The model we consider is the Lieb-Liniger model, which has experimentally been realized in various circumstances [2,3,5,6].We bring the model out of equilibrium by instantaneously turning off the interactions; this is a special case of an interaction quench.For short times one can think of this as a simulation of experiments where an ultracold Bose gas is released from a tight transverse confinement [7,8].
Several studies on the non-equilibrium dynamics of the Lieb-Liniger model have appeared before.The effect of instantaneously turning on the interactions was investigated in [9] using the Bethe Ansatz in combination with a Monte-Carlo sampling technique, and in [10] using the numerical time-evolving block decimation algorithm.The expansion starting from a regular array was studied in [11], where the properties of this specific initial state were exploited using the coordinate Bethe Ansatz.In [12] the possibility of studying quenches in the Lieb-Liniger model using an integrable field theory is discussed.
The outline of this article is as follows.In section I we describe the setup and explain the methods to compute correlation functions after the quench.This is followed by section II where we give the results for the time evolution of the non-local pair correlation and the autocorrelation.We also introduce a new type of correlation function that is a measure for the correlation between the pre-and post-quench states, which we will call quenchstraddling correlations.In section III we compare the long time average of the correlation functions with the predictions of various statistical ensembles.First we discuss why the Generalized Gibbs (GGE) ensemble fails in this case.Considering the generalized canonical ensemble (GCE), that is the GGE but with the total number of particles explicitly fixed, we show that it yields the correct results at all length scales.We end with a discussion and conclusions.

I. SETUP AND METHOD
The model we discuss in this paper is the onedimensional Bose gas which is described by the Lieb-Liniger Hamiltonian [13] with L the length of the system and c the interaction strength.For simplicity we will work with periodic boundary conditions.We will only consider the repulsive case c > 0. Two limiting regimes of the Lieb-Liniger arXiv:1201.1885v2[cond-mat.stat-mech]24 May 2012 model are: noninteracting Bosons (c = 0), and the socalled Tonks-Girardeau regime (c = ∞) [14] where the bosons effectively behave as free fermions.In order to study the non-equilibrium dynamics, we quantum quench the system [15,16].We do this by first preparing the system in the ground state for c > 0 and then instantaneously turning off the interactions: H(c) → H(0) at t 0 = 0.
The equal-time 2-point function Ψ † (x)Ψ(0) will not evolve in time for this specific quench.The Fouriertransform of this 2-point function consists of computing Ψ † k Ψ q .Since the states before and after the quench are translationally invariant, one only needs to compute the diagonal components (q = k).However, when c = 0, the momentum occupation operators Ψ † k Ψ k are constants of motion, hence the 2-point function will not be affected by this quench.For higher-point correlations or dynamical ones this will not be the case.As an illustration we will consider the measurable non-local pair correlation The non-local pair correlation for finite size systems has been computed in [17] using exact methods for the ground state.The case of finite temperature was studied in [18] using perturbative techniques.Without loss of generality, we will fix in the following the density at unity Ψ † (0)Ψ(0) = N/L = 1.In order to study the time evolution after the interactions are turned off, it is useful to write the final Hamiltonian as The time evolution of the field operators readily follows: e iH(0)t Ψ † k e −iH(0)t = Ψ † k e ik 2 t .So instead of acting with the time evolution operator on the initial state, we can act with it on the operators.The upshot of this is that we can avoid decomposing the initial state in terms of eigenstates of the final Hamiltonian.This decomposition typically involves a number of states that grows exponentially with the total number of particles.
Using the Fourier-transform of the field operators k the Heisenberg picture allows us to express the time evolution of the dynamical correlation function after the quench as where k ij = k i − k j and we used k 1 + k 2 = k 3 + k 4 because of translational invariance.The initial state |φ is the ground state of the fully interacting Lieb-Liniger model (c > 0).The quench takes place at t 0 = 0 and we evaluate the correlation function at times 0 ≤ t 1 ≤ t 2 .
In equilibrium (at t 2 = t 1 = t 0 ), this is the normal ordered 2-point function of current Ψ † (x)Ψ(x) operators and can efficiently be obtained using the techniques described in [17].In order to study the time evolution, we need to compute a four-point function of field operators Ψ † k .This four-point function of the initial state |φ can be computed in the framework of the Algebraic Bethe Ansatz [19] by generalizing the methods of [20] for their computation of two-point functions.We will briefly discuss how the computations are done without going in too much detail.First, we express the correlation function in terms of matrix elements.By inserting a resolution of the identity operator between every pair of adjacent field operators, one can write the four-point function in (3) as in terms of the matrix elements n|Ψ † k |m of the field operator.Here |φ is the initial state for N particles.The states |n 1 , |n 2 and |n 3 are intermediate states with N − 1, N − 2 and N − 1 particles respectively.In order to compute the matrix element n|Ψ † k |m one first solves the Bethe equations for the states |n and |m which results in a set of so-called rapidities for both states.The matrix elements can then be evaluated by computing the determinant of a matrix whose entries are rational functions of the rapidities of the two eigenstates involved.The explicit expression for the matrix element can be found in [21] as a sum of determinants and in [20] (see eq. 41) in terms of a single determinant.What remains to be performed are the actual summations over n 1 , n 2 and n 3 .Starting from the initial state |φ the Fock space of intermediate states |n 1 is scanned by navigating through choices of sets of quantum numbers.For each individual intermediate state, the Bethe equations are solved, and the matrix element is computed.The search for important states is done via the ABACUS method [22], and is close to optimal.In order to verify the accuracy of our results, we keep track of the sum-rule The sum over n 1 is truncated after a desired precision is achieved.Next, for every state |n 1 , we repeat the process by performing a summation over n 2 , for which we can also compute a sum-rule.We do not need to perform the summation over n 3 explicitly because we can use the hermitian conjugate of the matrix element φ|Ψ † k1 Ψ † k2 |n 2 , thereby constructing the full four-point function.For a system of N = 20 particles and large interaction c = 1000, which is the most difficult case, a number of intermediate states |n 2 of the order of 10 8 are needed in order to saturate the sum-rule at 99.6%.Once the four-point function is obtained, dynamical correlation functions (3) can be computed straightforwardly.

II. RESULTS
In this section, we study the dynamical correlation function g 2 (x, t 1 , t 2 ) (3) after the quench, for all times.We present the results for a system of size L = N = 20 starting from the ground state for various c > 0. For clarity, we specialize to the cases of the non-local pair correlation g 2 (x, t) = g 2 (x, t, t) and the auto-correlation g 2 (t 1 , t 2 ) = g 2 (0, t 1 , t 2 ).We conclude this section by considering what we will call the quench-straddling correlations.

Short Times
First, we consider the non-local pair correlation.We plot g 2 (x, t) for various times t after the quench as a function of x in fig. 1.For clarity we also plot how g 2 (x, t) for various fixed values of x, develops in time in fig. 2. We display here only the results for c = 1000 because for smaller c the behavior is roughly the same but less pronounced.In equilibrium (t = 0) we see that g 2 (x) vanishes for small x, as a result of destructive interference due to the fermionic character of the Tonks-Girardeau gas (this is also called anti-bunching).For large x, correlations decay, no destructive interference takes place and g 2 (x) → 1.If we now focus on the behavior for x = 0 we see that because of loss of coherence the correlation builds up.Eventually, for large times there will be constructive interference ( g 2 (x) > 1), which is called bunching, and is typical for free bosons (see fig. 2).For small x and small t we also expect that correlations will grow because of loss of coherence, however we have to take into account that L 0 g 2 (x, t)dx = 1 − 1/N is a conserved quantity which results in non-monotonic behavior for x > 0 at small times as can be clearly seen in fig. 2  In order to study how the correlations propagate through the system we plot |g 2 (x, t) − g 2 (x, 0)| in fig. 3.There is no light-cone effect because of the absence of a maximal velocity.As reference we have plotted x = v g t (dashed line) where the group velocity is computed as We see that the dominant changes take place after t ∼ x/v g .

Finite size effects
From a finite size study we can conclude that the behavior for short times, discussed in the previous section, is independent of the system size.If we now turn to large times finite size effects become more and more pronounced.Because of the simple dispersion relation ω(k n ) = (2πn/L) 2 exact revival occurs at t rev = L 2 /(4π).Partial revivals of decreasing strengths also occur at higher harmonics of the revival time: t rev /2, t rev /3, t rev /4, . ... Due to the revival, g 2 (0, t) gets suppressed for large times, as can be seen in fig. 4. In the thermodynamic limit there will be no suppression and one can easily show that g 2 (0) = 2 for free bosons in equilibrium.

Asymptotic behavior
In order to understand the intermediate time behavior, which is after the non-monotonic behavior and before finite size effects start to dominate, we fit the results for small x with the following test function We motivate this function as follows.The first part a(x; c) is the stationary part.For the time-dependent part we expect algebraic decay, since the initial state is quantum critical.The exact revival, due to the finite size, is taken into account by means of the sine-function.
The fits for small x and t are compared with the exact result in fig. 2. The data suggests that a(x, c) is a decreasing function of x which agrees with the result for the long time average which we will discuss in section III.We find that α(c) is a decreasing function of c and α ∼ 1/2 for c = 1000.However, since we have only a limited range where we can fit the data and α is small, no firm predictions can be made.After comparing fits for various system sizes, our data seem to be consistent with the dependence on L entering only via t rev .This agrees with the intuition that for small x and small times t after the quench finite size effects are irrelevant.Assuming this is indeed the case, one can from the finite size data make predictions for the thermodynamic limit by sending t rev → ∞ in (5) while keeping the other parameters fixed.The stationary part a(x; c) would then be the large time limit in the thermodynamic limit, since the time-dependent part now vanishes for large t, which is not the case for finite size systems.For free bosons in equilibrium one can easily show that g 2 (0) = 2, which is compatible with what we find from the fitting data: a(0, c) ∼ 2. In fig. 4 we compare the exact results for g 2 (0, t) with the fit (5) both for finite size and our predictions for the thermodynamic limit.
Since the post-quench state is completely described by the correlations on the initial ground state, one can try to predict asymptotic behavior, such as the exponent in (5), using low-energy effective theories like bosonization.If we write (3) as one can see that because of the integral over z one needs the correlation function at all length scales.Hence whether low-energy descriptions can be used or not remains an open question.

B. Auto-correlation
The auto-correlation g 2 (t 1 , t 2 ) is plotted in fig. 5 as a function of t 2 after various times t 1 after the quench.As for the non-local pair correlation, there is an antibunching bunching transition.However, in contrast to the evolution of the non-local pair correlation, the autocorrelation does increase monotonically in time after the quench, as can be seen in fig. 5.This is consistent with the fact that the integrated auto-correlation is not conserved in time, in contrast with the non-local pair correlations.

C. Quench-straddling correlations
The versatility of our method allows the computation of various different dynamical correlation functions.As an example we will study a correlation function which is a measure of how the pre-and post-quench states are correlated, what we will call a quench-straddling correlation function.We consider the situation where the quench takes place at time t 0 = 0. We then compute a two-point function where we evaluate Ψ(0) before the quench at time t − < t 0 and Ψ † (x) after the quench, at time t + > t 0 g straddle (x, t − , t + ) = φ|e iH(0)t+ Ψ † (x)e −iH(0)t+ e −iH(c)t− Ψ(0)e iH(c)t− |φ .(7) To simplify this expression we first use the Fouriertransform of the field operators Ψ † (x) =  5.The auto-correlation g2(0, t2; t1) for various times t1 = 0.00, 0.05, 0.10, 0.14, 0.24, 0.53 as a function of t2 for a system with N, L = 20 and c = 1000.
As before, we can then handle the time-evolution of the post-quench Hamiltonian H(0) using the Heisenbergpicture.For the time evolution of the pre-quench Hamiltonian H(c) we use the Schrödinger picture; this is achieved by inserting a resolution of the identity n |n n|, in terms of eigenstates |n of the pre-quench Hamiltonian between e −iH(0)t+ and e −iH(c)t− .By pulling out all the different phases and using translation invariance we arrive at Here E 0 is the energy of the initial state |φ and E n , k n are the energy and momentum of the intermediate state |n respectively.In fig.6 we plot Re{g straddle (0, t − , t + )} for various t − as a function of t = t + .A first observation is that in the case of no quench, t + = t 0 , the correlation decays much faster than the case after the quench t − = t 0 .This can be understood from that fact that, for a sufficiently large interaction strength c, the dispersion E(k) is linear for small k.Therefore, the majority of the states |n in ( 7) have an energy E(k) > k 2 .When considering the general quench times t 0 − t − , a striking observation we can make is that just after the quench the correlation function suddenly seems completely relaxed.The explanation is that for a quench-straddling correlation dephasing takes place both via the energy and the momenta, which can be considered as orthogonal directions.This speeds up the relaxation significantly for a short period right after the quench.The only coherence that is left is for very small k, resulting in an extremely slow relaxation compared to the case if there were no quench.

III. STATISTICAL ENSEMBLES
In this section we discuss whether the large time behavior of correlation functions after the quench can be described by a statistical ensemble.For simplicity we restrict ourselves in this section to the analyses of g 2 (x, t).Due to the finiteness of the system, no actual relaxation occurs.However, for most of the time, g 2 (x, t) oscillates around the same mean value, as is seen in fig. 4. Therefore, it is still useful to consider the long time average (LTA) of correlation functions Note that because of the exact revival we can take T LT A = t rev .The result of the time average is presented in fig. 7. It is well-known that quantum integrable models do not always relax to a state of thermal equilibrium.
As an alternative the Generalized Gibbs ensemble (GGE) [23] has been proposed.The GGE is based on the idea of Jaynes [24] to construct a statistical ensemble with maximal entropy subject to constraints for the expectation values of the conserved charges.Expectation values of observables in the GGE are then computed via O(x) GGE = Tr O(x)e − n βnQn /Z GGE (10) where Z GGE = Tr e − n βnQn is the generalized partition function.The Langrange multipliers β n corresponding to the conserved charges Q n are fixed via the initial conditions φ|Q n |φ = Q n GGE .So far, the GGE has only been used for models which are effectively free, namely the Hamiltonian can be written as

with the momentum occupation numbers Ψ †
k Ψ k as the obvious choice for the conserved charges.Many successful examples exist [16,23,[25][26][27][28][29][30][31].Usually, the GGE is implemented as a grand canonical ensemble, which simplifies the results because the partition function can be written as a product of single states.For example, the partition function can be written as where the + and − signs corresponds to fermions and bosons respectively.The Lagrange multipliers are determined via φ|Ψ † k Ψ k |φ = 1/(e β k ± 1), (+ fermions and − bosons).However, in this procedure the correlation between different conserved charges are not kept; that is, GGE for k = q.Despite these shortcomings, the GGE has been successful in many cases, although exceptions are known, for example when translation invariance is broken [32,33].Let us now turn to the case of the construction of the GGE for the LTA of g 2 (x).First we explicitly write down the LTA of g 2 (x) using ( 3) and ( 9) One can see that the LTA for the g 2 (x) is explicitly written as a sum over expectation values of products of the conserved charges.From this we can see that the GGE is not applicable in this case, since the initial state is the ground state of an interacting system and Wick's theorem does not apply here.Instead we considered the generalized canonical ensemble (GCE) by keeping the total number of particles fixed.
A. The generalized canonical ensemble Consider a system where the Hamiltonian can be diagonalized in terms of free particles H = k ω k Ψ † k Ψ k .We impose no restrictions on the dispersion relation ω k and the operators Ψ † k can have either fermionic or bosonic commutation relations.The partition function for such a system in the canonical ensemble (CE) can be written as where n max = 1, N for fermions and bosons respectively.The presence of δ j nj ,N makes it very complicated to evaluate the sums directly.Fortunately, one can compute the partition function Z N for a system of N particles via a well known recursion relation (see for instance [34][35][36]) with Z n = 0 for n < 0, Z 0 = 1 and The minus and plus signs correspond to fermions and bosons respectively.For a given β and N the partition function can be evaluated numerically.The computational complexity is of the order of n cutoff N 2 , where n cutoff is the number of one-particle states considered in the computation of z n .One can easily generalize the computation of the canonical ensemble to what we will call the generalized canonical ensemble (GCE) by making the replacement βω k → β k .Since the partition function now does not factorize, as is the case for the GGE, the expectation values of products of conserved charges are not necessarily uncorrelated.From the partition function one can derive the probability P ≥ k (n) of having at least n particles in the state k.This is done by starting the summation of n k at n in (13) and then using the recursion relation (14) to obtain The probability of having exactly n particles in the state k is therefore Similarly, one can derive the probability of having at least n k states in k and n q in q with k = q: from which an expression for p k,q (n k , n q ) can be obtained.Expectation values are now computed as follows: One can numerically solve these coupled equations in order to obtain values for β k .Once all β k are determined we can compute the other expectation values, for instance  8.The long time average (LTA) of the g2(x) compared with the predictions of the generalized canonical ensemble (GCE), the generalized Gibbs ensemble (GGE) and the canonical ensemble (CE), for a system with N, L = 20 and c = 1000.

B. Results
In fig.8 we compare the LTA, GGE, GCE and CE for g 2 (x).As expected, the predictions from CE and GGE are completely off.The predictions of GCE agree extremely well with the LTA see (fig.7), the relative error being less than 0.5% for all x.This may be ascribed to the fact that the sum-rule for g 2 (x) is only saturated to 99.6%.The reason why the GCE works is not completely obvious.We would like to stress that both for the GGE and the GCE only the expectation values of Ψ † k Ψ k are fixed via corresponding Lagrange multipliers β k .In principle one could expand the GGE or GCE by including products of conserved charges Ψ † k1 Ψ k1 Ψ † k2 Ψ k2 with their own Lagrange multiplier β k1,k2 , but that is not what is done here.Apparently, by demanding that N 2 = N 2 one almost completely fixes higher order expectation values.To gain more insight, we compare the expectation values Ψ † k Ψ k Ψ † k+q Ψ k+q for the LTA, GCE and the GGE plotted in fig. 9 for k = 0, 2π/L.If we first focus on the top graph (k = 0), we see that all q the LTA and GCE agree extremely well.The GGE is only valid for large q > 10, from which we conclude that at this point the conserved charges become uncorrelated.Furthermore, from the inset we see that the fluctuations of Ψ † 0 Ψ 0 differ significantly for the GGE.In the bottom graph (k = 2π/L), the correlations are much weaker compared to the case k = 0. We see that the LTA and GCE differ slightly now, but the GCE is still better than the GGE.
The results here are presented for finite size.In the thermodynamic limit one might expect that the prediction of the GCE and GGE coincide for local observables.On the other hand, correlations like Ψ † 0 Ψ 0 Ψ † q Ψ q for small q cannot be approximated using Wick's theorem; for these correlations the GGE remains invalid.To study how correlation functions like g 2 (x) behave in the thermodynamic limit, we compared the predictions of GCE and GGE for a system with N, L = 150 in fig.10.We see that the difference between GCE and GGE is still present for small x, but it is considerably smaller than for the case of N, L = 20.This leads us to believe that in the thermodynamic limit the GCE and GGE yield equivalent predictions for g 2 (x).

IV. CONCLUSIONS AND DISCUSSION
In this paper we have considered an interaction quench by turning off the interactions in the Lieb-Liniger model.By using the Heisenberg picture explicitly, the time evolution was computed via the correlation functions on the fully interacting initial state.The exact method we have used allowed us to study various correlation functions after the quench at all time and length scales.Although this paper focused on the results for the Lieb-Liniger model, we would like to stress that the methods presented in this paper can be applied to any model for which the Heisenberg picture can be used to efficiently obtain the time evolution and for which the correlation functions on the initial state can be computed exactly.For example, one could prepare the system in the ground state of the anisotropic Heisenberg spin chain.By switching off the anisotropy term, the time evolution is accessible from that of free fermions.
In case of the Lieb-Liniger model, we have studied several types of correlation functions after the quench.As initial state, the ground state of the Lieb-Liniger model was used for various interaction strengths c > 0. The time behavior of the correlation functions was studied both in the short and long time regimes.As expected, the results are most pronounced for large c, although the qualitative features do not strongly depend on the initial interaction strength.As a byproduct we have computed the four-point function of the Lieb-Liniger model for the first time.
The long time average has been compared with various statistical ensembles.As expected, the canonical ensemble fails to make correct predictions, as is usually the case for integrable models.The GGE gives better results.However, for an intermediate number of particles it still fails to capture all the features of the LTA.This discrepancy can be understood by the fact that the fluctuation of the total number of particles is large, while for the quench under consideration the total number of particles remains constant.By introducing the GCE, which keeps the total number of particles explicitly fixed, correct predictions for the LTA are obtained.The failure of the GGE can also be explained by the fact that correlations between the conserved charges are not kept.The GCE, which is essentially the GGE with one additional constraint, seems to give the correct correlations between the conserved charges; why this is the case remains unclear on a more formal level.It is worth to mention that the GCE can be applied in the same cases as the GGE.A firmer test of the GCE would be to study even higher order correlation functions such as g n (x) = (Ψ † (x)) n (Ψ(0)) n .In the thermodynamic limit it is expected that the particle number fluctuations of the GGE become irrelevant for local observables such as the g 2 (x); a comparison of GGE and GCE for a large number of particles confirm this.
It would also be interesting to see what the effect is of considering a different basis of conserved charges.For example in [28] the validity of the GGE was studied using two different bases for the conserved charges, of which only one agreed with the long time average.To make a connection with the interacting Lieb-Liniger model with c > 0, an appropriate basis would consist of the conserved charges in terms of derivatives of the transfer matrix, while sending c → 0. This leads to the following set of conserved charges: However, for this quench problem the expectation values φ|Q n |φ for n ≥ 3 are not well-defined because of the 1/k 4 behavior of φ|Ψ † k Ψ k |φ for large k [37] which reflects the ultralocality of the Lieb-Liniger model.This does not imply that the GGE fails in this case, as one could regularize the results by introducing a lattice spacing a while keeping the system integrable.One can then write down the expectation values as function of a and send a → 0 at the very end of the calculation.This will be investigated in future work.
A way to generalize the results of this paper is by considering quenches starting from arbitrary interaction strength c 1 and ending in a different arbitrary c 2 .In this case, one cannot use the Heisenberg picture to compute the time evolution of observables as was done in this paper.One could try to solve the time evolution in the Schrödinger picture by making a spectral decomposition of the initial state in terms of the eigenstates of the Hamiltonian after the quench.The first problem in this approach is the need for the overlap coefficients between the initial state and final states.For the Lieb-Liniger model, these overlap coefficients are only known in one particular case [9]; the general problem is still unsolved.Secondly, the spectral decomposition typically involves an exponential number of states as function of the system size.Unless there are huge degeneracies present in the spectrum, as is only the case for the special points c = 0 and c = ∞, performing the spectral sum seems intractable.We will return to these and further issues in future publications.

FIG. 4 .
FIG.4.Complete and partial revival of g2(x = 0, t) compared with the fits for a system of finite size and infinite size, while keeping the density fixed (N/L = 1).

FIG. 7 .
FIG. 7. The long time average g2(x) for various c for a system with N, L = 20.