Fermionisation dynamics of a strongly interacting 1D Bose gas after an interaction quench

We study the dynamics of a one-dimensional Bose gas after a sudden change of the interaction strength from zero to a finite value using the numerical time-evolving block decimation (TEBD) algorithm. It is shown that despite the integrability of the system, local quantities such as the two-particle correlation $g^{(2)}(x,x)$ attain steady state values in a short characteristic time inversely proportional to the Tonks parameter $\gamma$ and the square of the density. The asymptotic values are very close to those of a finite temperature grand canonical ensemble with a local temperature corresponding to initial energy and density. Non-local density-density correlations on the other hand approach a steady state on a much larger time scale determined by the finite propagation velocity of oscillatory correlation waves.


Introduction
The dynamics of interacting quantum systems from an initial non-equilibrium state constitutes a major challenge for many-body theory. In particular the question of thermalisation of integrable models regained attention recently due to the experimental progress in ultra-cold gases. As demonstrated in a beautiful experiment by Kinoshita et al. [1] for the example of a 1D ultra-cold Bose gas, integrable systems do not thermalise in the usual sense, i.e., reduced density-matrices relax on considerably different time scales than in the absence of integrability.
The speciality of the relaxation dynamics of integrable systems has been attributed to the presence of an infinite set of constants of motion with local character, i.e., which can be written as sums of operators acting only over a finite spatial range. Although thermalisation has been studied in a large body of theoretical papers it remains a largely unsolved problem. Most studies of specific models have been done either for non-interacting particles [2] or systems that can directly be mapped to free systems such as hard-core bosons [3], the Luttinger model [4], or the 1/r fermionic Hubbard model [5].
In the present letter we analyse the dynamics of a 1D Bose gas with s-wave scattering interactions, described by the Lieb-Liniger (LL) model, after a sudden quench of the interaction strength from zero to a finite value, covering the full range from weak to strong interactions. Performing numerical simulations using the time evolving block decimation algorithm (TEBD) [6,7], we show that local quantities, in particular the local two-particle correlation g (2) (0, 0; t), do approach steady-state values on a short time scale determined only by the Tonks parameter γ and the particle density ̺. This shows that although non-local quantities such as the momentum distribution do not approach a steady state over long times [8], there is an equilibration in a local sense. Furthermore the asymptotic values of g (2) (x, x) are very close to those obtained from a thermal Gibbs ensemble [9], with temperature and chemical potential determined by the initial conditions and the amplitude of the interaction quench. Thus it is possible to define local temperature and chemical potential and the influence of constants of motion other than total energy and particle number is very small, if present at all. Non-local quantities such as the density-density correlation approach a steadystate distribution on a larger time scale by way of correlation waves propagating out of the sample.

LL model and lattice approximation
A Bose gas in one spatial dimension is described by the Hamiltonian in units were = m = k B = 1. HereΨ(x) is the field operator of the Bose gas in second quantisation, V (x) some possible trap potential, and g the strength of the local particle-particle interaction. The latter is characterised by the dimensionless Tonks parameter γ = g/̺, where is the 1D density of particles. Specifically we consider here a system initially prepared in the non-interacting ground state. The initial canonical state has locally only diagonal elements. In the course of interactions non-diagonal elements are not created. Thus the reduced local density matrix is entirely determined by the number distribution, and the quantities of interest are the density ρ and the local two-particle correlation g (2) In principle also the corresponding higher-order moments are nonzero. They will, however, not be considered here.
We study the dynamics after the interaction quench numerically by means of the the time-evolving block decimation algorithm [6,7]. As for any numerical methods this requires the use of a discretised version of the Lib-Liniger model. The discretization of the LL model leads to the (non-integrable) Bose-Hubbard model [10] Here J = 1/(2∆x 2 ), U = g/∆x, and D = 1/∆x 2 , with ∆x being the lattice constant of the discretization grid. The appropriateness of discretised lattice models to describe continuous interacting Bose or Fermi gases in the limit ∆x → 0 has been discussed and verified in a number of earlier papers [10,11]. Note furthermore that using the boson-fermion duality in 1D the LL model can be mapped to an integrable lattice model, the spin 1/2 XXZ model [11].
For some data sets we used both, Equation (3), and the XXZ lattice model to verify that in the considered limit the non-integrability of (3) has no influence on the results. It turned out that for the simulation of dynamics the necessary grid sizes are much smaller than for equilibrium simulations [10,11]. Empirically we found that in order to minimise lattice artifacts resulting into numerical errors, the average number of particles per lattice site n should be small compared to 1/γ, where γ = g/̺ corresponds to the density at the centre of the cloud, i.e. ̺ → ̺(x = 0). This can be explained as follows: The interaction energy of a two-particle collision in the lattice, i.e., g/∆x should be smaller than the bandwidth of the lowest Bloch band, ∼ 1/∆x 2 , in order not to see lattice artifacts. At the edges of the cloud where collisions become less probable the condition is less strong. To accommodate the requirement of a very small n at the centre of the cloud we used space dependent grid sizes such that the average boson number per site was constant for the centre part of the particle distribution. Nevertheless to approximate the continuous model sufficiently well, very fine grids were needed leading to rather large lattice sizes on the order of up to 2880, which is rather challenging. In order to illustrate the effects of discretization we have plotted in Figure  1a the local two-particle correlation (see following section) for γ = 200/9 and increasing lattice sizes L, corresponding to finer grids. One clearly recognises oscillation artifacts which only slowly disappear with increasing L.
The convergence of the TEBD scheme was checked by varying the bond dimension χ of the matrix product state (MPS) and calculating the truncation error in the state norm a) accumulated during the time evolution. In Figure 1b the accumulated truncation error is plotted for γ = 200/9 and increasing values of χ from 25 to 200. One recognises that for the maximum value of χ = 200 which we were able to use, the truncation error is below the level of 10 −3 for the time scale of interest. This value is larger then the accuracy typically reached in ground state calculations. However we are not at the point where the cut-off explodes, which typically happens in dynamical simulations at some point. Finally the matrix dimension required to achieve a given accuracy does not depend on the discretization length, i.e. the number of lattice sites used. It is rather the number of particles which determines the complexity of the calculations. Thus the restriction to a moderate particle number allows us to work on lattice large compared to other applications of the algorithm.

Local Dynamics
In order to be able to perform numerical simulations with a fixed number of particles (up to 18, which corresponds to the experiments in references [12,13]) we have to work with a finite size system. Therefore we assumed an initial weak harmonic trapping potential V (x) = 1 2 ω 2 x 2 . Open or periodic boundary conditions, which can also be dealt with by the TEBD algorithm, could be used alternatively. Initially the Bose gas is in the canonical ground state (T = 0) of non-interacting bosons in the trap, for which the matrix product representation is analytically known, since it is a product of single particle states. At t = 0 we suddenly switch the interaction strength from zero to a finite value. At the same time the trap, the only purpose of which is was the preparation of an appropriate initial state, is switched off. On the time scales we are interested in, the density distribution does not change, so that the presence of a trap would be of no relevance. This also allows to apply the results of the present analysis to a homogeneous gas.
The initial state has a Gaussian density distribution ̺(x) = N ω π e −ωx 2 with l osc = 1 √ ω being the oscillator length. Figure 2 shows the time evolution of g (2) (0, 0; t) with time normalised to a characteristic time scale t ia for different values of γ. One recognises after an initial phase a power-law decay with an exponent that is monotonous in the interaction parameter. At times close to t ia a steady state value is attained indicating that a local equilibrium is reached. I.e., although globally a LL gas does not thermalise [1], local quantities do. The time scale t ia of the local dynamics can be estimated from Equation (3). The repulsive interaction Un(n−1) causes particle number fluctuations to be driven out of a given lattice site. This happens in the following way: Initially all components of the state vector have the same phase and tunnelling has no effect. However, due to the interaction, components with different particle number attain a differential phase shift and are subsequently coupled to states in adjacent lattice sites by tunnelling with rate J. Since in the limit ∆x → 0 we have J ≫ U, the maximum rate of this process is limited by the average interaction energy per particle U n . Thus we have Note that already for moderate interaction strength, γ ≫ 1/N 2 , this time is much shorter as for example the oscillation time t osc in the trap.
Note furthermore that although the characteristic time of the expansion of the gas after switching off the trap becomes much shorter than t osc for larger interactions, we found that the density profile did not change on the scale of t ia even for the largest values of γ used, since also t ia ∼ 1/γ.Note furthermore that although the characteristic time of the expansion of the gas after switching off the trap becomes much shorter than t osc for larger interactions, it will be large compared to t ia . This is because the kinetic energy transferred to the particles will be of the order γ and therefor their characteristic speed will be of the order √ γ only. Accordingly we found numerically that the density profile did not change on the timescale t ia even for the (0, 0, t) γ = 1/9 γ = 2/9 γ = 10/9 γ = 20/9 γ = 100/9 γ = 200/9 γ = 500/9 γ = 1000/9 largest values of γ used. Whether or not the thermalised local correlation will adiabatically follow the density evolution after longer times, i.e. when the expansion of the cloud sets in, cannot be concluded from our simulations. We would however expect such a behaviour. The fluctuations in the plots are artifacts of the discretization, which leads to an oscillatory behaviour of g (2) on top of the continuous-system time evolution. These artifacts, which are most pronounced for larger interactions, could not be eliminated completely even for the smallest grid sizes used. As a result the asymptotic values of g (2) (0, 0, t) can only be given with a certain error.
In figure 3 we have plotted the exponents obtained from a fit to the curves in Figure 2 which for intermediate times follows a power law The exponent is a monotonous function of the interaction strength and slowly approaches the limit −1 for γ → ∞, i.e., for a Tonks-Girardeau gas [14]. We now want to analyse the local state of the system in the stationary limit. In particular we will show that the local steady-state can be well described by the usual finite-temperature Gibbs state. To this end we calculate the expected asymptotic value g (2) YY (0, 0) from the thermodynamic Bethe Ansatz [9]. The system is initially prepared in its non-interacting ground state, so we have g (2) ini (0, 0) = 1 − 1/N, which in the thermodynamic limit N → ∞ approaches unity. The energy of this state with respect to the non-interacting Hamiltonian is 0. At time t = 0 the interaction is switched to a finite strength g > 0 and the expectation value of the interaction energy immediately after the quench is given by Since in a homogeneous system there is no x-dependence the energy per particle is Here we have introduced the critical temperature T c in one dimension T c = ̺ 2 /2. One recognises that the system is in a highly excited non-equilibrium state after the quench if γ 1. Using the energy per particle, the density ̺ and the Tonks parameter γ as input parameter, we can extract a temperature T of a corresponding thermal Gibbs state by inverting the Yang-Yang equations of the thermodynamic Bethe Ansatz [9] for the excitation-energy and particle densities ǫ(q), n(q) in momentum space Here K(λ, ξ) = 2g g 2 +(λ−ξ) 2 and ̺ = dλ n(λ). With the help of the Hellmann-Feynman theorem we can then obtain the value g (2) YY (0, 0) corresponding to the Gibbs state at temperature T [15]: with f (γ, T ) being the free energy per particle In the limit 1 ≪ τ ≪ γ 2 , where τ = T /T c , (9) attains the simple form [15] g YY (0, 0) = 2τ γ 2 .  In Figure 4 we have plotted the values of g (2) YY (0, 0) from the thermal Gibbs state in the thermodynamic limit as function of the interaction strength γ (solid line). Also shown are the steady-state values obtained from the numerical simulation in Figure 2. The error bar indicates uncertainties which are here due to discretization artifacts and error estimates obtained from comparing simulations with MPS bond dimensions χ = 100 and 200. It is available only for one parameter set, since the variation of the discretization length and the bond dimension is numerically expensive. However we expect it to be of about the same relative size for all data points. One recognises that g (2) (0, 0; t) attains in the long-time limit values which are close to that of the thermal Gibbs state. One should note that the steady-state values for the largest values of γ (γ = 500/9 and 18 bosons as well as γ = 1000/9 and 9 bosons) are slightly overestimated in the simulation due to the remaining grid artifacts since here n γ ≈ 0.73 is no longer small compared to unity. Also shown is the asymptotic local temperature of the gas in units of the degeneracy temperature. For large values of γ, T ≫ T c , i.e., after relaxation the gas is in a state with large local temperature.

Non-local Relaxation
We now discuss the dynamics of non-local quantities. Specifically we consider the non-local two-particle correlation g (2) (0, x; t). In Figure 5 we have plotted g (2) (0, x; t) for different times after the interaction quench. One recognises that while the local correlations attain a steady-state value on a short time scale, the non-local evolution happens much slower. Switching on the particle-particle repulsion leads to a fast reduction of the probability to find two particles at the same position. Associated with this is a correlation flow to larger distances (Colour online) Time evolution of non-local density-density correlations g (2) (0, x; t) for γ = 200/9. x = 0 denotes the centre of the cloud. One recognises the formation of expanding correlation waves. The dashed blue line shows the approximation (13) to the non-local correlation in a thermal Gibbs state from [16] multiplied by g (2) (0, 0, t = 0) = 8/9 to account for the final particle number (N = 9) used in the simulation.
leading to expanding correlation waves. For very short times the propagation velocity of correlation waves is faster than the Fermi velocity v F = π̺. But at the largest time shown in Figure 5 corresponding to t = 0.01t osc , the maximum of the correlation wave has travelled a distance of approximately ∆x = 0.12l osc which is consistent with the speed of sound which for large values of γ approaches The buildup of a maximum that behaves like a wavefront can be understood as follows: The integral over space of g (2) (0, x; t) is (in a homogeneous system) a constant with respect to time due to particle number conservation. The numerator of g (2) (0, x; t) is proportional to the joint probability distribution to find a particle at position x given that there is one particle at the origin.
So as the quench can not change g (2) (0, x; t) significantly outside the light cone given by the Fermi velocity, the reduction of the probability to find two particles close together must be accompanied by an increase at finite distance.
As one can see from Figure 5 also the non-local correlation function approaches at least for smaller distances in the large-time limit that of the thermal Gibbs state with temperature and density given by the initial conditions and the Tonks parameter γ. For comparison we have plotted an approximation to the finite-temperature non-local g (2) from reference [16] which holds in the regime 1 ≪ τ ≪ γ 2 g (2) Here λ T = 4π/τ ̺ 2 is the thermal de Broglie length.

Experimental Observation
In the following we discuss the possibility to test the local relaxation in an experiment. For this we make use of the fact that by energy conservation the interaction energy lost by the decrease of g (2) (0, 0) must be gained as kinetic energy, E kin (t) = E ia (t = 0) − E ia (t) and the kinetic energy therefore directly gives the value of g (2) (0, 0, t) in the homogeneous case: If the interaction is turned on at t = 0 and turned off abruptly at some time t = t 1 the kinetic energy is the only remaining in the system and can be used to measure g (2) (0, 0; t 1 ).
In an the experimental setup, the gas must be confined e.g. by an harmonic trapping potential. So the initial non-interacting state has a Gaussian density distribution. It is also a good assumption, that the correlations decay locally as in the homogeneous system corresponding to the local density provided the density ̺(x) remains constant over the time scale of interest. This is indeed the case, if t int = 1/(γρ 2 ) ≪ l osc /v s ∼ l osc /ρ This means, that the Tonks parameter must be large compared to 1 N , which is of course the case we are interested in. We note that the region in the wings of the density distribution which does not fulfil this constraint gives a negligible contribution to the total interaction energy. Of course measuring the kinetic energy in the trap gives only an average of g (2) (x,x) ̺ 2 (x) over the trap.

Summary
Using the time-evolving block decimation scheme we have numerically analysed the dynamics of a 1D Bose gas (LL-model), after an interaction quench from zero to a finite value. Although globally the 1D Bose gas does not thermalise, we have shown that local quantities attain a steady-state value on a time scale t ia = (γ̺ 2 ) −1 . Within the achievable accuracy these values are consistent with the assumption that local quantities relax to a thermal Gibbs state with local temperature determined by the initial energy and chemical potential. Non-local quantities such as the density-density correlation relax on a much longer time scale set by the velocity of sound by means of correlation waves propagating out of the sample.