A coordinate Bethe ansatz approach to the calculation of equilibrium and nonequilibrium correlations of the one-dimensional Bose gas

We use the coordinate Bethe ansatz to exactly calculate matrix elements between eigenstates of the Lieb-Liniger model of one-dimensional bosons interacting via a two-body delta-potential. We investigate the static correlation functions of the zero-temperature ground state and their dependence on interaction strength, and analyze the effects of system size in the crossover from few-body to mesoscopic regimes for up to seven particles. We also obtain time-dependent nonequilibrium correlation functions for five particles following quenches of the interaction strength from two distinct initial states. One quench is from the non-interacting ground state and the other from a correlated ground state near the strongly interacting Tonks-Girardeau regime. The final interaction strength and conserved energy are chosen to be the same for both quenches. The integrability of the model highly constrains its dynamics, and we demonstrate that the time-averaged correlation functions following quenches from these two distinct initial conditions are both nonthermal and moreover distinct from one another.


I. INTRODUCTION
The Lieb-Liniger model of a one-dimensional Bose gas with repulsive delta-function interactions is a paradigmatic example of an exactly solvable continuous, integrable many-body quantum system [1]. In particular, it has served as the context for the development of theoretical tools that have subsequently been widely applied in the study of integrable systems, such as the so-called "thermodynamic Bethe ansatz" functional representation, which provides the exact equation of state, excitation spectrum [1], and bulk parameters [2] of the system in the thermodynamic limit. However, the calculation of correlation functions from the exact solutions provided by the Bethe ansatz is notoriously difficult.
At zero temperature, exact closed-form solutions for some equilibrium correlation functions are known in the Tonks-Girardeau limit of infinite interaction strength [3][4][5][6][7]. This comparatively tractable limit also allows for some strong-coupling expansion results for large but finite interactions [7][8][9][10]. In the opposite weakly interacting quasi-condensate regime, a mean-field approach can be used to describe the system [11] and a Bogoliubov method can be used to determine the low-lying excitation spectrum [12], relying on small density fluctuations. Fewer results are available for intermediate interaction strengths, away from the strongly-interacting and weakly-interacting regimes. The development of the Luttinger liquid description of quantum fluids [13] and the related formalism of conformal field theory [14,15] have lead to the prediction of power-law scaling for first-order * j.zill@uq.edu.au correlations at large distances, with an exponent given in terms of the equation of state that is known exactly from the thermodynamic Bethe ansatz [16]. The algebraic Bethe ansatz provides a determinantal representation of correlations, from which their asymptotic behavior can be extracted [17]. More recently, exact expressions for local second-and third-order correlations [18][19][20], together with exact results for the one-body correlation function at asymptotically short distances [21] in terms of the equation of state have been derived.
Away from the asymptotic short-and long-range regimes, the behavior of correlation functions is less well known. For intermediate interaction strengths and arbitrary length scales one must resort to numerics to determine the correlation functions. Results for the latter have been obtained using numerical methodologies including quantum Monte Carlo [22,23], and density matrix renormalization group approaches [24]. A recently developed, integrability-based approach combines the decomposition of correlation functions into sums over matrix elements (form factors) of certain simple operators between Bethe ansatz eigenstates [25,26]. This approach has generated results, for example, for static and dynamical equilibrium correlations at zero and finite temperature for systems of up to N ≈ 100 particles [27]. Other finite temperature results for correlation functions have been obtained using imaginary time stochastic gauge methods [28,29], taking the non-relativistic limit of a relativistic field theory [30], utilizing Fermi-Bose mapping for the strongly interacting gas [9,31,32], employing perturbative expansions in temperature and interaction strength [33], as well as combining the thermodynamic Bethe ansatz with the Hellmann-Feynman theorem [34].
Experiments with ultracold quantum gases are able arXiv:1601.00434v2 [cond-mat.quant-gas] 14 Apr 2016 to realize effectively one-dimensional systems by tightly confining the gas in two of the three spatial dimensions, either using optical lattice potentials or atom-chip traps [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. These experiments are now probing the predictions of the Lieb-Liniger model. The configurability of quantum-gas experiments allows for so-called quenches of the system, in which Hamiltonian parameters of the system are abruptly changed, and thus for the study of the Lieb-Linger model out of equilibrium, providing even greater challenges for theory.
In this paper we take a different approach, and calculate correlation functions of the Lieb-Linger model, both in and out of equilibrium, by calculating matrix elements between Lieb-Liniger eigenstates directly within the coordinate Bethe ansatz formalism. Given the known expressions for the coordinate-space forms of Lieb-Liniger eigenstates, we generate symbolic expressions for matrix elements of operators between these states in terms of the Bethe rapidities. The numerically obtained values of the rapidities can then be substituted to yield essentially numerically exact values for the matrix elements.
In our previous work we applied this methodology to quenches from the ideal gas ground state to positive γ for up to N = 5 particles [90]. In Sec. II we provide the details of the methodology, and describe how it can be used to calculate the matrix elements of the Lieb-Liniger eigenstates. These symbolic expressions, and thus the computational cost of evaluating them, grow combinatorially with particle number, restricting the method to systems of only a few particles. However for small particle numbers N ≤ 7 we obtain numerically exact results for ground-state correlations, which are described in Sec. III. Our results demonstrate that local correlations in the strongly interacting regime are already close to their thermodynamic-limit values for these few-body to mesoscopic systems.
An additional advantage of our methodology is that it can also calculate overlaps between Lieb-Liniger eigenstates corresponding to any two interaction strengths, which allows us to study the dynamics of quenches of the interaction strength between arbitrary values. In Sec. IV we utilize this property to study the effects of integrability on the relaxation of the Lieb-Liniger model following such a quench. In particular, we compare two nonequilibrium quench scenarios with the same final Hamiltonian and state energy, but beginning from starkly different initial states. Statistical mechanics would predict that the system would relax to the same thermal state in both cases, but due to the integrability of the Lieb-Linger model not only are the time-averaged states following the two quenches non-thermal, they are also distinct. After characterizing and comparing the nonequilibrium dynamics following both quenches, we conclude in Sec. V.

II. COORDINATE BETHE-ANSATZ METHODOLOGY
A. Lieb-Liniger model eigenstates The Lieb-Liniger model [1] describes a system of N indistinguishable bosons subject to a delta-function interaction potential in a periodic one-dimensional (1D) geometry of length L. We work in units such that = 1 and the particle mass m = 1/2, and so the Hamiltonian of this system readŝ where c is the interaction strength. The coordinate Bethe ansatz yields eigenstates |{λ j } of Hamiltonian (1) with spatial representation [17] where the rapidities λ j (or quasimomenta) are solutions of the Bethe equations The quantum numbers m j are any N distinct integers (half-integers) in the case that N is odd (even) [2], and σ denotes a sum over all N ! permutations σ = {σ(j)} of {1, 2, . . . , N }. The normalization constant reads [17] The rapidities determine the total momentum P = N j=1 λ j and energy E = N j=1 λ 2 j of the system in each eigenstate. The ground state of the system corresponds to the set of N rapidities that minimize E and constitute the (pseudo-)Fermi sea of the 1D Bose gas [17]. The Fermi momentum is the magnitude of the largest rapidity occurring in the ground state in the Tonks-Girardeau limit of strong interactions [3]. The only parameter of the Lieb-Liniger model in the thermodynamic limit is the dimensionless coupling γ ≡ c/n, where n ≡ N/L is the 1D density. In finite systems, physical quantities also depend on the particle number N (see, e.g., Sec. III C), whereas the length L of our system, and therefore also the density n, are arbitrary. Consequently, in this article we will specify both N and γ. Unless specified otherwise, we measure time in units of k −2 F , energy in units of k 2 F , and length in units of k −1 F .

B. Calculation of correlation functions and overlaps
As the eigenstates |{λ j } form a complete basis [91] for the state space of the Lieb-Liniger model, the expectation value Ô t = Tr{ρ(t)Ô} of an arbitrary op-eratorÔ in a Schrödinger-picture density matrixρ(t) can be expressed as a sum of matrix elements ofÔ between the states |{λ j } . In particular, in a pure state whereas in a statistical ensemble with density matrix In this article, we focus in particular on the normalized m th -order equal-time correlation functions whereΨ ( †) (x) is the annihilation (creation) operator for the Bose field andn(x) ≡Ψ † (x)Ψ(x). Here and in the following we drop the time index t of the state vectors.
Since the Hamiltonian we consider in this article is translationally invariant along the periodic volume of length L, the mean density n(x) ≡ n is constant in both time and space, and The correlation functions g (m) (x 1 , . . . , x m , x 1 , . . . , x m ; t) can therefore be expressed as the expectation values of the operatorsĝ (m) (x 1 , . . . , x m , x 1 , . . . , x m ) ≡ Ψ † (x 1 ) · · ·Ψ † (x m )Ψ(x 1 ) · · ·Ψ(x m )/n m . We note that for the same reasons as above the matrix elements {λ j }|ĝ (m) (x 1 , . . . , x m , x 1 , . . . , x m )|{λ j } are invariant under global coordinate shifts x → x + d and thus, without loss of generality, we can set one of the spatial variables to zero. For the first-order correlation function, the matrix elements are The evaluation of the integral in Eq. (10) is complicated by the sign function in Eq. (2) and the associated nonanalyticities in ζ {λj } ({x i }) where any two particle coordinates x k and x l coincide. However, we can use the Bose symmetry of the wave function ζ {λj } ({x i }) to reexpress this matrix element as a sum of integrals over the ordered domains [10] R M,j (x) : 0 ≤ x 1 < · · · < x j < x < x j+1 < · · · < x M ≤ L.
where (for the repulsive interactions c > 0 considered in this article) the κ m are real numbers. A single closed form for this integral does not exist, as in general one or more κ m may vanish, and this must be handled separately from the case of κ m = 0. However, given knowledge of the particular sets of rapidities {λ j } and {λ j } (and permutations σ and σ ), and thus of the locations of zero exponents κ m = 0 in Eq. (16), each individual integral of this form can be reduced to an algebraic expression in terms of {κ m }. More specifically, each successive integration dx m yields a term (involving, in general, x m+1 ) arising from the primitive integral [92] dx x p e ikx = −(i/k) p+1 Γ(p + 1, −ikx) in the case that κ m is nonzero, or from dx x p otherwise. In our calculations, the construction of algebraic expressions for the integrals occurring in Eqs. (13)- (15) in terms of the rapidities λ j is efficiently performed by a simple computer algorithm that accounts for and combines the symbolic terms that arise from these succes-sive reductions. We note that, e.g., each matrix element {λ j }|ĝ (1) (0, x)|{λ j } is a sum of N integrals over (N −1)dimensional domains and that the integrand in each case comprises (N !) 2 terms [10], illustrating the dramatically increasing computational cost of evaluating correlation functions with increasing N . Nevertheless, the explicit closed-form expression for the integral produced by our algorithm can be evaluated to obtain a numerically exact result by substituting in the values of the rapidities. The latter are obtained by solving Eq. (3) numerically using Newton's method, starting in the Tonks-Girardeau regime of strong interactions γ 1 and iteratively pro-gressing to smaller values of γ using initial guesses given by linear extrapolation of the solutions at stronger interaction strengths. We note that this algorithmic approach also provides for the efficient and accurate calculation of the overlaps {λ j }|{µ j } between eigenstates of Hamiltonian (1) corresponding to different values of γ, which we make use of in our analysis of nonequilibrium dynamics in Sec. IV. In particular, the overlap between an arbitrary eigenstate |{λ j } ofĤ at a finite interaction strength γ > 0 and the noninteracting ground state |0 , with constant spatial representation {x i }|0 = L −N/2 , is simply given by which can easily be evaluated semi-analytically using our algorithm. In practice we find that the results we obtain for the overlaps from our evaluation of Eq. (18) agree with the recently derived closed-form expressions for these quantities [84,[93][94][95], which imply in particular that {λ j }|0 ∝ 1/λ 2 j as any λ j → ∞.

III. GROUND-STATE CORRELATION FUNCTIONS
As a first application of our methodology we calculate the correlation functions of the Lieb-Liniger model in the ground state for up to N = 7 particles. In this case, we need to evaluate only the diagonal elements of Eqs. (13)- (15) in the ground-state wave function, thereby obtaining exact algebraic expressions for correlation functions in terms of the ground-state rapidities, which are themselves determined to machine precision (Sec. II B). The ground-state correlations of the Lieb-Liniger model have been considered extensively in previous works (see Refs. [96,97] and references therein), and we compare our exact mesoscopic results to those obtained with various other methods and approximations, for finite system sizes as well as in the thermodynamic limit. This allows us to clarify the utility and limitations of calculations, such as ours here and in Ref. [90], that involve only small particle numbers.

A. First-order correlations
We begin by considering the first-order correlation function g (1) (x) ≡ g (1) (0, x) in the ground state of the Lieb-Liniger model. In Fig. 1(a) we plot g (1) (x) for N = 7 particles for a range of interaction strengths γ, which exhibits the expected decrease in spatial phase co-herence with increasing γ [16]. As is well known, true long-range order, lim x→∞ g (1) (x) = n 0 > 0 [98,99], is prohibited in an interacting homogeneous 1D Bose gas in the thermodynamic limit, even at zero temperature (see Ref. [97] and references therein). Indeed the Lieb-Liniger system is quantum critical at zero temperature, and the asymptotic long-range behavior of g (1) (x) is a power-law decay (so-called quasi -long-range order) [17].
This power-law scaling of g (1) (x) is only expected to be realized at separations x large compared to the healing length ξ = 1/ √ γ and, in a finite periodic geometry such as we consider here, is curtailed by the finite extent L of the system (see, e.g., Ref. [16]). Indeed, for γ = 0.1, the power-law decay is not visible in our finite-sized calculation, although as the interaction strength γ increases g (1) (x) exhibits behavior consistent with power-law decay over an increasingly large range of x, see Fig. 1(a). In particular, for γ 10, our results for g (1) (x) seem to converge toward the asymptotic scaling of the Tonks-Girardeau limit (black dot-dashed line) with increasing γ.
Due to the translational invariance of our system, the first-order correlations of the Lieb-Liniger ground state are encoded in the momentum distribution which, in our finite periodic geometry, is only defined for discrete momenta k j = 2πj/L, with j an integer. In Fig. 1(b) we plot the momentum distributions n(k j ) corresponding to the first-order correlation functions g (1) (x) shown in Fig. 1(a). The first feature that we note in Fig. 1(b) is that for all interaction strengths, n(k) exhibits a power-law decay n(k) ∝ k −4 (dot-dashed black line) at high momenta. This is a universal result for delta-function interactions in one dimension [21,89,100] (and indeed also in higher dimensions [101]). The effects (Color online) One-and two-body correlations in the Lieb-Liniger ground state, for N = 7 particles. (a) Non-local first-order coherence g (1) (x). The black dot-dashed line indicates the asymptotic long-range behavior g (1) (x) ∝ |x| −1/2 of a Tonks-Girardeau gas in the thermodynamic limit. (b) Corresponding zero-temperature momentum distribution n(kj). The black dot-dashed line indicates the universal highmomentum power-law scaling n(k) ∝ k −4 common to all positive interaction strengths [21]. (c) Non-local second-order coherence g (2) (x). (d) Corresponding static structure factor S(k).
of the finite extent L of the system on the first-order correlations are again evident in this momentum-space representation. For γ = 0.1, no deviation from the ∝ k −4 scaling is observed for the smallest (nonzero) momenta k j that can be resolved in the periodic geometry. For larger values of the interaction strength, n(k) departs from the ∝ k −4 scaling at increasingly large values of k with increasing γ, and develops a hump at momenta near k F for γ 10 [89]. We note that although the small-k behavior of n(k) tends towards the ∝ k −1/2 scaling exhibited by the Tonks-Girardeau gas in the thermodynamic limit, the rounding off of the power-law decay of g (1) (x) as x → L/2 precludes n(k) from reaching the known asymptotic k → 0 behavior in our finite geometry.
B. Second-, third-, and fourth-order correlations In Fig. 1(c), we present the nonlocal second-order coherence g (2) (x) ≡ g (2) (0, x, x, 0), which provides a measure of density-density correlations, for N = 7 particles at a range of interaction strengths γ. In the limiting case of an ideal gas (γ = 0), the ground state of the system is a Fock state of N particles in the zeromomentum single-particle mode, and the second-order coherence g (2) γ=0 (x) = 1 − N −1 (horizontal dashed line) is therefore independent of x. As the interaction strength γ is increased, the second-order coherence is increasingly suppressed at zero spatial separation and correspondingly enhanced at separations x 2k −1 F . Oscillations in g (2) (x) develop at finite x as the system enters the strongly interacting regime γ 1 [9,17] and, in particular, for γ = 100 (dashed cyan line), our numerical results are practically indistinguishable from the exact Tonks-Girardeau limit result (solid black line) [3].
An alternative representation of the second-order correlations of the ground state is given by the static structure factor S(k), which is related to g (2) (x) by [11] In Fig. 1(d) we present the structure factors S(k) corresponding to the correlation functions g (2) (x) shown in Fig. 1(c). For all values of γ, S(0) = 0 due to particlenumber conservation and translational invariance. In the ideal-gas limit (red circles) S(k j ) = 1 for all nonzero k j .
In the opposite limit of a Tonks-Girardeau gas which tends, in the thermodynamic limit, to the wellknown result (see, e.g., Ref. [9]) S(k) = |k|/2k F for |k| ≤ 2k F , and S(k) = 1 for |k| > 2k F . Just as for g (2) (x), we observe that for γ = 100 (cyan plus symbols), our numerical results for S(k) are almost identical to the known exact expression [Eq. (21)] for the Tonks-Girardeau limit (black crosses). For smaller values of γ our mesoscopic results for S(k) appear consistent with those of Refs. [22,25], obtained using quantum Monte Carlo and algebraic-Bethe ansatz techniques, respectively.
We now focus in more detail on local correlation functions. We note that the local second-order coherence has recently been proposed as a measure of quantum criticality in the 1D boson system [102], while the local third-order correlations have received increasing attention both in theory [103] and experiment [47,[104][105][106]. The local fourth-order correlations for the Lieb-Liniger model have also been investigated [107]. In Fig. 2, we plot the local second-order coherence g (2) (0) (solid red line), together with the local third-order coherence g (3) 3 /n 3 (dotted green line), and the local fourth-order coherence g (4) 4 [Ψ(0)] 4 /n 4 (dashed blue line) for N = 7 particles and a broad range of interaction strengths γ. For comparison, we also plot the asymptotic results obtained in the Bogoliubov limit of weak interactions (γ → 0) in the thermodynamic limit [12,18] (left-hand dot-dashed lines). The numerical results for small γ are broadly comparable to these thermodynamic-limit results. However, for the small particle numbers considered here, the suppression of g (2) (0), 10 −16 (Color online) Interaction-strength dependence of the local second-, third-and fourth-order coherence in the Lieb-Liniger ground state, for N = 7 particles. To aid visibility, we plot g (2) (0) scaled by a factor of 10 1 , and g (4) (0) scaled by a factor of 10 −1 . Dot-dashed lines indicate asymptotic weak-(γ 1) and strong-coupling (γ 1) expressions for g (2) (0), g (3) (0) and g (4) (0) in the thermodynamic limit (see text). g (3) (0), and g (4) (0) due to interactions in the limit of small γ is overshadowed by the suppression due to the finite population of the system [20]. At larger γ, the effects of interactions dominate, and the numerical results converge closely to the appropriate strong-coupling expressions [18] (right-hand dot-dashed lines). We note, therefore, that the local correlations of the Lieb-Liniger ground state, and particularly their scaling with γ, appear to be quite insensitive to the infrared cutoff imposed by the finite extent of our system in the strongly interacting regime γ 1.

C. System-size dependence
The results we have obtained so far indicate that, as expected, the small size of our system leads to corrections to correlation functions as compared to their known asymptotic forms in the thermodynamic limit. However, our results also suggest that the effects of finite system size are comparatively less important for local correlations, particularly in the limit of large interaction strengths γ 1. To further elucidate the potential significance of finite-size effects in our calculations of nonequilibrium dynamics [90], here we give a brief characterization of the dependence of correlation functions of the Lieb-Liniger ground state on the particle number N at a fixed value of the interaction strength γ.
Specifically we consider the case for γ = 10, as this value places the system in the strongly interacting regime γ 1 (which appears less sensitive to finite-size effects than the weakly interacting regime γ 1), while still exhibiting significant deviations from the Tonks-Girardeau limit (see, e.g., Ref. [9]). Whereas elsewhere in this paper we quote momenta (lengths) in units of k F (k −1 F ), in comparing results between systems with different particle numbers N we quote momenta (lengths) in units In Fig. 3(a) we plot g (1) (x) for particle numbers N = 3, 4, 5, 6, and 7. For small x, the curves fall nearly perfectly on one line. The same behavior can be observed for the large-k tail of the corresponding momentum distribution n(k), which we plot in Fig. 3(b). Indeed, at larger momenta k 2πn, n(k) appears to exhibit a rapid collapse to a single curve with increasing N [21,109]. However, the differences in n(k) are so small that they can not be seen in Fig. 3(b). For small momenta, our choice of units implies an increasing resolution with increasing particle number, specifically k 1 = 2π/L × (πn) −1 = 2/N . However, this lowest resolvable momentum seems to fall on one line for increasing particle number, indicating that the infrared behavior of large systems can be at least partly accessed by our mesoscopic system sizes.
Luttinger-liquid theory predicts a long-range powerlaw decay g (1) (x) ∝ |x| −1/2K , where the Luttinger parameter K can be calculated from the thermodynamic limit of the Bethe ansatz solution (see, e.g., Refs. [16,17] and references therein). For our parameters we have K = 1.40, implying an asymptotic scaling g (1) (x) ∝ |x| −0.357 [black dot-dashed line in Fig. 3(a)]. This corresponds to a power-law behavior n(k) ∝ |k| −1+1/2K = |k| −0.643 [16] [dot-dashed line in Fig. 3(b)] for small momenta. We note that this infrared scaling is a true many-body effect and as such does not show up for N = 2 particles. Indeed, one can show analytically that, for N = 2, the momentum distribution n(k) ∝ (λ 2 1 −k 2 ) −2 and thus k −4 is the highest power in the series expansion of n(k).
In Fig. 3(c) we plot the nonlocal second-order coherence g (2) (x) for γ = 10 and N = 3, 4, 5, 6, and 7. The corresponding static structure factor S(k) is shown in Fig. 3(d). In Fig. 3(d) we also plot (black dot-dashed line) the form of S(k) resulting from the phenomenological expression proposed in Ref. [108] (see also Ref. [110]). This expression involves the limiting dispersions and edge exponents of the Lieb-Liniger model, which we obtain by numerically solving the appropriate integral equations [1,111]. We also plot the corresponding prediction for g (2) (x) (black dot-dashed line) in Fig. 3(c). We note that the numerical results for our mesoscopic systems are, in general, rather close to the phenomenological thermodynamic-limit expressions even for the relatively small particle numbers considered here.

IV. APPLICATION TO NONEQUILIBRIUM DYNAMICS
We now apply our methodology to the nonequilibrium dynamics of the Lieb-Liniger model. Specifically, we consider the evolution of a system, initially prepared in the ground state of Hamiltonian (1) with interaction strength γ 0 , following an abrupt change, at time t = 0, of the interaction strength to a distinct value γ = γ 0 -a socalled "interaction quench". The evolution of the system following such a quench is generated by Hamiltonian (1) with interaction strength γ, which we denote byĤ(γ) hereafter. The time-evolving state is given at all times t > 0 by where |{λ j } are the eigenstates ofĤ(γ) with energies E {λj } , and C {λj } ≡ {λ j }|ψ 0 are the overlaps of the |{λ j } with the initial state |ψ 0 . The expectation value of an arbitrary operatorÔ in the state |ψ(t) is given by We use the methodology described in Sec. II to evaluate both the overlaps C {λj } and the matrix elements {λ j }|Ô|{λ j } that appear in Eq. (23). One of the features of our methodology is that it allows us to describe quenches between arbitrary interaction strengths. In this paper we consider two interaction-strength quenches, from different initial interaction strengths γ 0 , to a common final value of the coupling γ. Specifically, we consider a quench from the non-interacting limit γ 0 = 0 (similar to those previously studied in Refs. [63,78,84,90,[112][113][114][115]) and a quench from the correlated ground state obtained for a strong interaction strength γ 0 = 100. AsĤ(γ) is time independent following the quench, energy is conserved during the dynamics. We choose the final interaction strength after the two quenches such that the postquench energy is the same in both cases.
The statistical description of the dynamics of sufficiently ergodic systems is usually based on the assumption that the energy is the sole integral of motion, such that the equilibrium system is entirely determined by its energy. If this would be the case for our system, the two quenches would lead to the same equilibrium state. However, the dynamics according to the integrable Lieb-Liniger Hamiltonian are strongly constrained by the conserved quantities other than the total energy. By performing two different quenches to the same final Hamiltonian and energy, we investigate the effects of integrability on the postquench evolution of the Lieb-Liniger system.
Here, we consider quenches of N = 5 particles, and determine this final interaction strength to machine precision, inferring a value γ * = 3.7660 . . . from numerical solutions for the energy and local second-order coherence of the ground state at finite γ (Sec. III B). We note that although the overlaps C {λj } of the initial state |ψ 0 with the eigenstates ofĤ(γ * ) can be calculated analytically in the case of the quench from γ 0 = 0 [93][94][95], for the quench from γ 0 = 100 no closed-form expressions for these quantities are known, and thus their numerical values must be determined using the semi-analytical Neglecting degeneracies in the spectrum ofĤ(γ * ) (see discussion in Appendix B), such averages are given by the expectation values Ô DE = Tr{ρ DEÔ } of operatorŝ O in the density matrix of the diagonal ensemble [116,117]. Formally, the sums in Eq. (22), (23), and (26) range over an infinite number of eigenstates |{λ j } , and thus the basis over which |ψ(t) is expanded must be truncated in our numerical calculations. By only including eigenstates with an absolute initial-state overlap |C {λj } | larger than some threshold, we consistently neglect small contributions to correlation functions from weakly occupied eigenstates and minimize the truncation error for a given basis size. We quantify this truncation error by the violations of the normalization and energy sum rules, as we discuss in Appendix A.

A. Evolution of two-body correlations
In Fig. 4 we plot the time evolution of the local secondorder coherence g (2) (0, t) for N = 5 particles following quenches of the interaction strength from initial values γ 0 = 0 (red dotted line) and γ 0 = 100 (blue dashed line) to the common final value γ * . For the quench from the noninteracting initial state (γ 0 = 0), as time evolves the local second-order coherence decays from its initial value g (2) (0, t = 0) = 1 − N −1 before settling down to fluctuate about the diagonal-ensemble expectation value g (2) DE (0) (horizontal dot-dashed line). This behavior is consistent with results obtained for similar quenches of the interaction strength from zero to a positive value in Ref. [90]. For the quench from γ 0 = 100, the value of g (2) (0) in the initial "fermionized" state is g (2) (0) ≈ 10 −3 . In this case g (2) (0, t) rises as time progresses, and then exhibits somewhat irregular oscillations about g (2) DE (0) (horizontal solid line). We observe that the decay (growth) of g (2) (0, t) to its diagonal-ensemble value and the onset of irregular oscillations about this value occur on comparable time scales in the two quenches.
We note that the predictions of the diagonal ensemble for the local second-order coherence g (2) DE (0) are very similar for the two quenches, despite the significant difference between the values of g (2) (0) in the two initial states. However, they are clearly distinct -g (2) DE (0) for the quench from the noninteracting state is in fact larger than that for the quench from the correlated state by an amount ≈ 0.0125, demonstrating that the system retains some memory of its initial state in the long time limit as is expected for an integrable system. We analyze this difference in more detail in Sec. IV C.
We now turn our attention to the time evolution of the full non-local second-order correlation function g (2) (x, t). In Fig. 5(a) we show the dependence of g (2) (x, t) on separation x for the quench from the noninteracting initial state at four representative times. [Note that the upper limit x = 2πk −1 F of the x axis in Fig. 5(a) corresponds to x = L/2 in the present case of N = 5 particles.] At t = 0 (horizontal solid line), the second-order coherence has the constant form of the noninteracting ground state. At short times (e.g., t = 0.01 k −2 F , red dashed line) a minimum in g (2) (x) develops at zero separation, together with the corresponding maximum required by the conservation of L 0 dx g (2) (x, t) [66]. As time progresses a wave pattern of maxima and minima develops and propagates away from the origin (e.g., t = 0.1 k −2 F , green dotted line). By time t = 1 k −2 F (blue dot-dashed line), the distinct maxima and minima of g (2) (x, t) have broadened in such a way that they are no longer clearly distinguishable and the correlation function agrees reasonably well with its diagonal-ensemble form (black dotdashed line) for small separations x 0.25 × 2πk −1 F . In Fig. 5(b) we show the full space and time dependence of g (2) (x, t) following a quench from γ 0 = 0, which gives a more complete picture of the development of a correlation wave at short length scales and its propagation to larger values of x as time progresses. The correlation wave we observe here is consistent with the results of previous investigations of the dynamics following the sudden introduction of repulsive interactions in an initially noninteracting gas [63,64,66,78,118].
In Fig. 5(d) we plot the spatial form of g (2) (x, t) for the quench from γ 0 = 100 at the same four representative times considered in Fig. 5(a). Despite the obvious distinction that the initial (t = 0, solid grey line) correlation function is in the fermionized regime with g (2) (0) 1, the behavior of g (2) (x, t) for this quench is qualitatively similar to that observed for the quench from γ 0 = 0, in that at early times (e.g., t = 0.01k −2 F , red dashed line), deviations from g (2) (x, t = 0) occur only at small separations x 2πk −1 F . Moreover, as time evolves and g (2) (0, t) increases towards g (2) DE (0), larger modulations of g (2) (x, t) about its initial functional form develop (e.g., t = 0.1k −2 F , green dotted line). At later times (e.g., t = 1k −2 F , blue dot-dashed line), g (2) (x, t) is close to g (2) DE (x) at small separations x 0.25 × 2πk −1 F , but exhibits large excursions away from it at larger x. In Fig. 5(e) we plot the full space and time dependence of g (2) (x, t) following the quench from γ 0 = 100. Although the behavior of g (2) (x, t) here obviously differs from that following a quench from the noninteracting initial state [ Fig. 5(b)], with the "fermionic" depression around x = 0 lessening rather than growing in magnitude, a similar pattern of propagating correlation waves in g (2) (x, t) can again be seen.
The correlation-wave pattern common to both quenches is more clearly exhibited by the change g (2) (x, t) − g (2) (x, 0) in the correlation function following the quench, which we plot in Figs. 5(c) and 5(f). This representation of the postquench second-order coherence of the system reveals a remarkably similar pattern of propagating waves in both cases, although the maxima and minima of the two wave patterns are inverted relative to one another. Fitting a power law to the position x(t) of the first propagating extremum of each of the two correlation waves, we find x ∝ t 0.516±0.012 for the quench from γ 0 = 0 and x ∝ t 0.496±0.005 for the quench from γ 0 = 100, which we indicate by the solid black lines in Figs. 5(c) and 5(f). These power-law trajectories are consistent with the "telescoping" x ∝ t 1/2 behavior obtained for a quench γ = 0 → ∞ in Ref. [63], and for quenches x (units of 2πk −1 F ) from finite repulsive interactions to the noninteracting limit in Ref. [119] (see also Ref. [120]). The small scale features on top of the main propagating extrema differ for the two quenches, with fast oscillations appearing more pronounced for the quench γ = 0 → γ * in Fig. 5(c). Even though hardly visible in Fig. 5(f), they are still present for the quench from γ = 100 → γ * , but due to the different distribution of overlaps in the final basis compared to the quench from γ 0 = 0 (cf. Sec. IV C), they contain more high-frequency components and therefore the fine structure differs.

B. Time-averaged correlations
We now compare the time-averaged second-order correlation functions following the two quenches with the form of this function that would be obtained if, following the quench, the system relaxed to thermal equilibrium. As in Ref. [90] we make use of the canonical ensemble, for which the density matrix is given bŷ where the partition function The inverse temperature β is determined implicitly by fixing the mean energy in the stateρ CE to the common postquench energy, i.e., Tr{ρ CEĤ (γ * )} = E 0→γ * . The sum in Eq. (27), like that in Eq. (26), formally ranges over an infinite number of eigenstates. We therefore truncate this sum by applying a cutoff in energy, as described in Appendix A.
In Fig. 6(a) we plot the second-order correlation function g (2) CE (x) = Tr{ρ CEĝ (2) (0, x)} in the canonical ensemble (black dot-dashed line), along with the diagonalensemble predictions g (2) DE (x) for the quenches from γ 0 = 0 (red solid line) and from γ 0 = 100 (blue dotted line). For comparison we also plot the correlation functions in the initial states with γ 0 = 0 (horizontal line), γ 0 = 100 (grey dashed line), as well as the ground state for γ = γ * (solid black line). For the quench from γ 0 = 0, the timeaveraged value g (2) DE (0) is smaller than the corresponding thermal value g (2) CE (0), consistent with the results of Refs. [84,90,112]. In fact g (2) DE (x) is suppressed below g DE (x) > g (2) CE (x) at larger separations x due to particle number and momentum conservation. For the quench γ = 100 → γ * , the diagonal-ensemble coherence function g (2) DE (x) is similar in shape to that of the quench from γ 0 = 0. However, it is somewhat smaller at x = 0, and correspondingly larger at large x. This indicates some memory of the initial state preserved by the dynamics of the integrable Lieb-Liniger system [58,85]. Despite these differences, on the whole both functions g (2) DE (x) are comparable to g (2) CE (x) (cf. also Ref. [66]). We note, however, that they are also both reasonably close to the ground state result for g (2) (x) at interaction strength γ * (solid black line), although the local value g (2) DE (0) for both quenches is much closer to the thermal value than the ground state value.
Since the system is in its ground state before the quench for both γ 0 = 0 and γ 0 = 100, and the total momentum operatorP commutes with the Hamiltonian, the postquench states at γ * only have support on eigenstates with total momentum P = 0. Furthermore, the spatially structureless initial state at γ 0 = 0 implies additional parity-invariance ({λ j } = {−λ j }) in Bethe rapidity space for the postquench eigenstates [93][94][95]. Thus an interesting question to ask is if we constructed a canonical density matrix (27) restricted to P = 0 states, or one further restricted to parity-invariant states (which are a subset of the P = 0 states), would these yield better agreement with the diagonal ensemble predictions for the quenches? We have performed these constructions with the temperature in both cases fixed via the postquench energy in the same way as for the canonical ensemble, cf. Eq. (27) and the following text.
In Fig. 6(b), we plot the resulting second-order correlation function g (2) CE (x) = Tr{ρ CEĝ (2) (0, x)} for the standard canonical ensemble (black dot-dashed line), as well as in the restricted P = 0 ensemble (solid black line), and the parity-invariant ensemble (solid grey line). We also include the diagonal-ensemble predictions g (2) DE (x) for the quenches from γ 0 = 0 (red solid line) and from γ 0 = 100 (blue dotted line). It can be seen that the restricted ensembles give results for the correlation function that are quite close to the standard canonical ensemble, and are no closer to the diagonal ensemble results.

C. Contributions to relaxed correlation functions
The relaxation of the nonlocal correlations g (2) (x, t) takes place on a similar time scale to that of the local coherence g (2) (0, t) for both of the quenches considered here. This should be contrasted with, e.g., the behavior following a quench from the noninteracting limit to γ = 100 reported in Ref. [90], in which g (2) (0, t) decays rapidly and the development and propagation of correlation waves occurs over a significantly longer time scale. We identify the absence of a significant separation of the time scales of local and nonlocal evolution here as a consequence of the fact that only a small number of eigenstates contribute significantly to the postquench dynamics (cf. Ref. [90] and references therein). Indeed, we find that the purity Γ DE ≡ Tr{(ρ DE ) 2 } of the diagonal-ensemble density matrix takes values ≈ 0.52 for the quench γ = 0 → γ * and ≈ 0.63 for the quench γ = 100 → γ * , indicating rather weak participation of the eigenstates |{λ j } in the dynamics. The difference in the purities can largely be attributed to the somewhat greater occupation of the ground state ofĤ(γ * ) following the quench from the γ = 100 initial state.
To further illustrate the difference in the final states, in Fig. 7 we plot the occupations of eigenstates with energy E {λj } for the quenches from γ 0 = 0 (red crosses) and γ 0 = 100 (blue squares). For the quench from γ 0 = 100, significantly more eigenstates have occupations above a given threshold than in the case of γ 0 = 0, resulting in a much larger basis size in this case. However, the occupation of the ground state ofĤ(γ * ) is somewhat larger for the quench from γ 0 = 100 than for γ 0 = 0, and the low-lying excited states are comparatively weakly occupied for γ 0 = 100, cf. Fig. 7(b). This result is reasonably intuitive, as the ground state for γ = γ * is moderately correlated, and will be more similar to the γ = 100 than  the γ = 0 ground state. The distribution of normalization over eigenstates |{λ j } is thus more sharply "localized" on the ground state in this case, resulting in the somewhat larger value of the purity Γ DE following this quench.
For comparison, we also plot the occupations of the three ensembles introduced in Sec. IV B in Fig. 7. The restrictions lead to a reduction in available eigenstates for any given energy-window, and correspondingly the temperature of the canonical ensemble is smaller than that of the P = 0 ensemble, which is in turn smaller than that of the parity-invariant ensemble. The occupations of eigenstates for the quench from γ 0 = 0 (red crosses) and from γ 0 = 100 (blue squares) are suggestive of power-law decay at high energies. For small energies on the other hand, Fig. 7(b) shows that the functional form is not incompatible with exponential decay.

V. CONCLUSIONS
We have described a method to calculate matrix elements between eigenstates of the Lieb-Liniger model of one-dimensional delta-interacting bosons. This method is based on the coordinate Bethe ansatz, which generates a complete set of energy eigenfunctions for any fixed coupling strength. This allows us to obtain overlaps between eigenstates of different Hamiltonians, as well as expressions for correlation functions. By introducing periodic boundary conditions, we obtained expressions amenable to numerical evaluation. We applied our methodology to the evaluation of first-, second-, third-, and fourth-order correlation functions in the ground state of the Lieb-Liniger model for various values of the interparticle interaction strength. Our results indicate that although the correlations of the system are in general distorted by the small system size, finite-size effects become increasingly less significant with increasing interaction strength and decreasing spatial separation.
Out of equilibrium, we investigated the dynamics of relaxation after a quantum quench of the interparticle interaction strength towards a non-thermal steady state. Starting from two different initial states, we quenched to a common final interaction strength γ * chosen in such a way that both postquench energies were the same. Our calculations reveal a similar relaxation process for the second-order coherence g (2) (x, t) for both initial states: the build-up of correlations on short interparticle distances and their propagation through the system as time progresses. The time-averaged second-order correlation functions in both cases disagreed with the prediction for thermal equilibrium and were biased, relative to one another, towards their pre-quench forms -an intuitive result given the integrability of the system. In the future it would be interesting to study quenches from other initial states with the same final energy to explore how the memory of the initial state is manifest in different situations.
Although our method is restricted to small system sizes due to computational complexity and here only applied to five particles out of equilibrium, we were able to obtain the dynamical evolution as well as time-averaged correlation functions to high precision. Finally we note that the evaluation of matrix elements of the Lieb-Liniger model with this method is not restricted to real-valued Bethe rapidities, opening the door to investigating the nonequilibrium dynamics of attractively interacting systems (where the rapidities become complex-valued) and that following quenches from more complex initial states. The Hilbert space of the Lieb-Liniger model is infinite dimensional, and therefore the sums in Eqs. (22), (23), and (26) must be truncated for numerical purposes. Here, we provide details of the truncation scheme for the two different initial states we considered in Sec. IV, and explain how we quantify the error resulting from this truncation.
For the quench from γ 0 = 0, the initial state |ψ 0 only has nonzero overlap with eigenstates |{λ j } ofĤ(γ * ) that are parity invariant (i.e., eigenstates for which {λ j } = {−λ j }) and, a fortiori , have zero total momentum P [114]. The strongly-correlated initial state of the quench from γ 0 = 100 similarly has zero overlap with eigenstates |{λ j } with nonzero total momentum, but in this case states contributing to |ψ(t) , and thusρ DE , need not be parity-invariant in general. For γ 0 = 0 our results for the overlaps agree with recently obtained analytical expressions [94,95], which predict real positive overlaps, given the phase convention implicit in Eq. (2), for quenches to γ > 0. For γ 0 = 100, we find that the overlaps are still real, but are no longer restricted to positive values.
We briefly summarize our procedure to determine the cutoff here -see Appendix A of Ref. [90] for an extended discussion for the case of parity-invariant states. It can be shown [2] that the solutions {λ j } of the Bethe equations (3) are in one-to-one correspondence with the numbers m j that appear in Eq. (3). This allows us to uniquely label states by the set {m j }. Without loss of generality, we order the numbers m j such that m 1 > m 2 > · · · > m N −1 > m N , and we only need consider states for which j m j = 0, corresponding to zero total momentum P . We specialize hereafter to the case N = 5, which is the largest N for which we consider the dynamics in this article. The states can be grouped into families, labelled by m 1 . We have found empirically that within each such family, the eigenstate (m 1 , 1, 0, −1, −m 1 ) has the largest absolute overlap | {λ j }|ψ 0 | with the initial state, for both initial states we consider (γ 0 = 0 and γ 0 = 100). Furthermore, this overlap is larger than that of the most significantly contributing eigenstate (m 1 +1, 1, 0, −1, −m 1 −1) of the following family (m 1 +1). We therefore construct the basis by considering in turn each family m 1 and including all states within that family for which the overlap with the initial state exceeds our chosen threshold value C min . Eventually, for some value of m 1 , even the eigenstate (m 1 , 1, 0, −1, −m 1 ) has overlap with |ψ 0 smaller than the threshold, at which point all states that meet the threshold have been accounted for.
We note that the Lieb-Liniger model has an infinite number of conserved charges [Q (m) ,Ĥ(γ)] = 0; m = 0, 1, 2, . . . , with eigenvalues given byQ (m) |{λ j } = N l=1 λ m l |{λ j } . However, for a quench from γ 0 = 0 their expectation values in the diagonal ensemble Q (m) DE diverge for all even m ≥ 4 [94,95]. Our numerical results suggest that this is also the case for quenches from γ 0 > 0 (indeed, they diverge for almost all states but eigenstates [121,122]). For all odd values m, the expectation values of the corresponding conserved chargeŝ Q (m) are identically zero for our initial states and quench protocol. Thus, the only nontrivial and regular conserved quantities are the particle number (m = 0) and energy (m = 2). As in Ref. [90], we quantify the saturation of the normalization and energy sum rules by the sum-rule violations respectively, where E γ0→γ is the exact postquench energy [Eq. (24)]. We note that the calculation of timedependent observables involves a double sum over {λ j }, and is therefore more numerically demanding than the calculation of expectation values in the DE. Moreover, the calculation of the local coherence g (2) (0, t) is much less demanding than that of the full nonlocal g (2) (x, t).
We therefore use different thresholds C min , resulting in different basis sizes and sum-rule violations, in the calculation of g (2) (0, t), g (2) (x, t), and g DE (x), as indicated in Table I. We note that the energy sum rule is in general less well satisfied than the normalization sum rule, due to the ∝ λ −4 tail of the diagonal-ensemble distribution of eigenstates [90]. We find also that both sum rules are less well satisfied for the quench γ = 100 → γ * , despite the truncation procedure described above resulting in more than five times as many basis states being employed in its solution than are used in the quench γ = 0 → γ * .
For expectation values in the CE [Eq. (27)], we truncate the basis by retaining all states with energies below some cutoff E cut . The inverse temperature β is then chosen to minimize the energy sum-rule violation ∆E. The normalization sum rule is fulfilled by construction. Since all states (not only those with zero momentum) contribute to this sum, the number of eigenstates involved in canonical-ensemble calculations is much larger than that in diagonal-ensemble calculations. For the canonical-ensemble correlation function plotted in Fig. 6 we used an energy cutoff of 3.2 × 10 2 k 2 F , which yields a basis of 2.1 × 10 6 eigenstates |{λ j } . We checked that this cutoff is sufficiently large to ensure saturation of g (2) CE (x) (Fig. 6). For the ensemble restricted to P = 0 eigenstates [ Fig. 6(b)], we used an energy cut-off of 6.4 × 10 5 k 2 F , corresponding to 44530 eigenstates, while for the parity-invariant ensemble we used an energy cutoff of 8.5 × 10 6 k 2 F , corresponding to 64204 eigenstates.  I. Basis-set sizes and sum-rule violations for full nonlocal, time-evolving second-order coherence g (2) (x, t), for local, time-evolving second-order coherence g (2) (x = 0, t), and for time-averaged second-order coherence g (2) DE (x) following quenches from γ0 = 0, and γ0 = 100 to γ * = 3.7660 . . . .
The first term in Eq. (B1) is simply the diagonalensemble density matrixρ DE [Eq. (26)], to whichρ re-duces in the absence of degeneracies in the spectrum of H(γ). This is the case for the quench from γ 0 = 0, as the only eigenstates ofĤ(γ * ) with nonvanishing overlaps with |ψ 0 in this case are the parity-invariant states |{λ j } with {λ j }={−λ j }, which are nondegenerate (see Ref. [90] and references therein). By contrast, in a quench from γ 0 > 0, |ψ(t) has support on non-parity-invariant states |{λ j } , which are degenerate with their parity conjugates |{−λ j } . In general such degeneracies can have observable consequences for time-averaged expectation values [117]. However, as can be seen from Fig. 8, the correction to g (2) DE (x) due to the contributions of degenerate eigenstates in the case of the quench from γ 0 = 100 is small. It is straightforward to show that the elements {−λ j }|ĝ (2) (0)|{λ j } of the local second-order coherence between parity-conjugate states must vanish due to symmetry considerations. At larger separations x, the matrix elements between these pairs of states are nonzero, as illustrated in Fig. 8(a). However, these contributions are small compared to the diagonal-ensemble result g (2) DE (x), and indeed the total contribution of all parity-conjugate states in our finite-basis description [ Fig. 8(b)] would yield a barely visible correction to the function g (2) DE (x) plotted in Fig. 6. We note also that the substitution of ρ DE for the time-averaged density matrixρ introduces negligible error in the calculation of the purity of this matrix (Sec. IV C).