Abstract

Quantum quench dynamics of the attractive one-dimensional Bose gas via the coordinate Bethe ansatz

J. C. Zill1, T. M. Wright1, K. V. Kheruntsyan1, T. Gasenzer2,3, M. J. Davis4,5*

1 School of Mathematics and Physics, The University of Queensland, Brisbane QLD 4072, Australia

2 Kirchhoff-Institut für Physik, Universität Heidelberg, Im Neuenheimer Feld 227, 69120 Heidelberg, Germany

3 ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung, 64291 Darmstadt, Germany

4 ARC Centre of Excellence in Future Low-Energy Electronics Technologies, School of Mathematics and Physics, The University of Queensland, Brisbane QLD 4072, Australia

* mdavis@physics.uq.edu.au

July 5, 2019

## Abstract

We use the coordinate Bethe ansatz to study the Lieb–Liniger model of a one-dimensional gas of bosons on a finite-sized ring interacting via an attractive delta-function potential. We calculate zero-temperature correlation functions for seven particles in the vicinity of the crossover to a localized solitonic state and study the dynamics of a system of four particles quenched to attractive interactions from the ideal-gas ground state. We determine the time evolution of correlation functions, as well as their temporal averages, and discuss the role of bound states in shaping the postquench correlations and relaxation dynamics.

Contents\@mkbothCONTENTSCONTENTS

toc

## 1 Introduction

The near-perfect isolation and exquisite control possible for many experimental parameters in ultra-cold atomic gases has enabled the study of nonequilibrium dynamics of closed many-body quantum systems [1]. A number of different trapping geometries have led to the realization of quasi-one-dimensional systems [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17] that are well described by the paradigmatic exactly solvable Lieb–Liniger model of pointlike interacting bosons [18, 19, 20]. As this model is integrable, the various forms of the Bethe ansatz provide powerful methodologies with which to investigate the physics it describes [18, 19, 21, 22, 23, 24].

One of the simplest methods of taking a quantum system out of equilibrium is to effect an instantaneous change of a parameter in its Hamiltonian — a so-called quantum quench. Several authors have considered the nonequilibrium dynamics of repulsively interacting systems, where one particularly well-studied scenario is an interaction quench starting from the zero-temperature ideal gas [25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Here we study quantum quenches in which a one-dimensional Bose gas, initially prepared in its noninteracting ground state, is subjected to the abrupt introduction of attractive interparticle interactions [35, 36].

The ground-state wave function for the attractive one-dimensional (1D) Bose gas on the infinite line with finite particle number was constructed by McGuire [37] and consists of a single bound state of all the particles. For systems with finite spatial extent, the coordinate Bethe ansatz provides solutions in terms of quasi-momenta (or rapidities), which for attractive interactions are in general complex-valued. Ground-state solutions on a finite ring were found numerically in Refs. [38, 39].

Since the energy of the ground state is proportional to , where is the particle number, a proper thermodynamic limit with and fixed density does not exist [18, 40, 23]. However, the limit with is well defined, and was recently analysed in Ref. [41]. The zero-density limit , is also well defined and nontrivial for attractive interactions. In this limit, some correlation functions are accessible with the algebraic Bethe ansatz [42, 43].

An alternative large-system limit is given by in a finite ring of circumference . In particular, in the Bogoliubov limit , , , where is the interaction strength, a mean-field Gross–Pitaevskii description of the finite-circumference system predicts the appearance of a localized bright-soliton state beyond some threshold interaction strength [44, 45]. This has been interpreted as evidence for spontaneous breaking of translational symmetry in the infinite-, finite- limit [44, 46, 47]. However, Bogoliubov theory predicts a diverging quantum depletion in the vicinity of the threshold interaction strength, invalidating the mean-field description in this regime [44].

A many-body analysis for finite reveals a smooth crossover between a uniform condensate and a state with solitonic correlations, as expected in a finite system [44, 48, 49, 46]. Such an analysis also indicates that the gap at the crossover point vanishes as  [44]. The Bogoliubov-theory prediction of a vanishing gap at the crossover point in the semiclassical limit is thus regained. The crossover to the correlated state has therefore been interpreted [44] as a kind of effective quantum phase transition in the finite- system, though it should be stressed that the crossover in a system of finite particle number cannot be considered a finite-size precursor of a true quantum phase transition, as no proper thermodynamic limit exists.

In a full many-body quantum-mechanical treatment, energy eigenstates on the localized side of the crossover respect the symmetry of the Hamiltonian, but may contain solitonic structure in (pair) correlations. Localized bright solitons can thus be constructed from superpositions of certain exact many-body wave functions [50, 51, 52], which are given by the Bethe ansatz [18, 19, 37]. An integral equation for the density of Bethe rapidities of the ground state for particle number , valid across the crossover, has recently been derived and signatures of the crossover were observed in this density [47]. Bright-soliton-like structures have also been observed experimentally in elongated quantum-gas samples [53, 54, 55, 56, 57, 58, 59].

A particular nonequilibrium scenario for the attractive 1D Bose gas was proposed in Refs. [60, 61] and subsequently realized experimentally in Ref. [7]. In the latter work the system was prepared near the ground state at strong repulsive interactions, before the interactions were suddenly switched to strongly attractive using a confinement-induced resonance [20]. In doing so a metastable state was created: the so-called super-Tonks gas [60, 61, 62, 63, 64]. This highly excited state of the attractive gas has a “fermionized” character [62] that both stabilizes it against decay via recombination losses and implies a large overlap with the Tonks–Girardeau-like prequench state, leading to efficient state preparation via the interaction quench [63, 64]. This comparatively tractable regime also allows for a Luttinger-liquid description [65], as well as numerical studies with algebraic Bethe-ansatz [65] and tensor-network methods [66]. Local correlations in the super-Tonks regime can be obtained via an identification of the Lieb–Liniger gas with a particular nonrelativistic limit of the sinh-Gordon model [67], as well as by combining the equation of state of the super-Tonks gas with the Hellmann–Feynman theorem [63].

There are fewer results available for more general quench scenarios of the one-dimensional Bose gas involving attractive interparticle interactions. References [68, 69] introduced a Bethe-ansatz method, based on the Yudson contour-integral representation [70], for calculations of nonequilibrium correlation functions in systems of a few particles in the infinite-volume limit. Recently, the local second-order correlation function in the relaxed state following a quench from the ideal-gas ground state to attractive interactions was determined in the thermodynamic limit111The quench from the ideal gas to attractive interactions leaves the system with a finite energy per unit length and the thermodynamic limit is therefore well defined in this case [35, 36]. [35, 36] using the quench-action method [71, 72].

In Refs. [32, 33] we developed a methodology for the calculation of equilibrium and nonequilibrium correlation functions of the repulsively interacting Lieb–Liniger gas based on the semi-analytical evaluation of matrix elements between the eigenstates of the Lieb–Liniger Hamiltonian given by the coordinate Bethe ansatz. Here we extend this approach to the attractively interacting gas, for which the Bethe rapidities that characterize the eigenstates are in general complex-valued, indicating the presence of multiparticle bound states. We apply our method to calculate results for the time evolution of correlation functions following a quench to attractive interactions from the ideal-gas ground state, for a system of four particles. As in our previous studies of quenches to repulsive interactions [32, 33], we find that finite-size effects are significant for quenches to weak final interaction strengths. For strong final interaction strengths our results for the time-averaged local second-order correlation function are consistent with the stationary values in the thermodynamic limit calculated in Refs. [71, 72]. In contrast to that work, however, our approach allows us to also calculate the time-averaged value of the postquench third-order correlation function, which we find to be dramatically enhanced over the ideal-gas value, implying that three-body recombination losses would be significant in experimental realizations of the quench. Our approach also allows us to calculate the dynamical evolution of correlation functions following the quench, and for a quench to strong attractive interactions we observe behaviour similar to that following a quench to repulsive interactions of the same magnitude, superposed with characteristic contributions of bound states at small interparticle separations.

This paper is organised as follows. We provide a brief summary of the Lieb-Liniger model in Sec. 2. We also discuss the complications that arise in numerically solving the Bethe equations due to the appearance of complex Bethe rapidities, and explain how we manage these. In Sec. 3, we calculate ground-state correlation functions for up to seven particles in the vicinity of the mean-field crossover point where solitonic correlations emerge. We also present results for the ground state of four particles subject to strongly attractive interactions. In Sec. 4, we compute representative nonequilibrium correlation functions following quenches of the interaction strength from zero to attractive values for up to four particles. We discuss quenches to the weakly interacting regime in the vicinity of the mean-field crossover, as well as those to the more strongly interacting regime. We also compare the nonequilibrium dynamics to that following an interaction quench to repulsive interactions of the same magnitude. In Sec. 5 we present results for time-averaged correlation functions, before concluding in Sec. 6.

## 2 Methodology

### 2.1 Lieb–Liniger model

The Lieb–Liniger model [18, 19] describes a system of indistinguishable bosons subject to a delta-function interaction potential in a one-dimensional geometry. The Hamiltonian is

 ^H=−N∑i=1∂2∂x2i+2cN∑i

where is the interaction strength, and we have set and the particle mass . The interactions are attractive for , and repulsive for . The eigenstates of Hamiltonian (1) in the ordered spatial permutation sector () are given by the coordinate Bethe ansatz in the form [22]

 ζ{λj}({xi}) ≡⟨{xi}|{λj}⟩ =A{λj}∑σ(−1)[σ]a(σ)exp[iN∑m=1xmλσ(m)], (2)

where the sum runs over all permutations of , denotes the sign of the permutation, and the scattering factors are

 a(σ)=∏k>l(λσ(k)−λσ(l)−ic). (3)

The quantities are termed the rapidities, or quasimomenta of the Bethe-ansatz wave function. The normalization constant is given by [22]

 A{λj}=[N!det{M{λj}}∏k>l[(λk−λl)2+c2]]−1/2, (4)

where is the matrix with elements

 [M{λj}]kl =δkl(L+N∑m=12cc2+(λk−λm)2) −2cc2+(λk−λl)2. (5)

Imposing periodic boundary conditions leads to a set of equations for the rapidities, the so-called Bethe equations

 eiLλj=∏l≠j(λj−λl)+ic(λj−λl)−ic, (6)

where is the length of the periodic geometry. The rapidities determine the total momentum and energy of the system in each eigenstate. The ground state of the system for attractive interactions is an –body bound state (the finite-system analogue of the McGuire cluster state [37]) and has purely imaginary rapidities [38, 39]. All eigenstates corresponding to bound states have some Bethe rapidities with imaginary components. This is in contrast to the repulsively interacting system (), for which the solutions to the Bethe equations (6) are purely real. These are usually parameterized by a set of quantum numbers , which for are proportional to , see e.g. Ref. [22]. For the attractively interacting gas, it is more convenient to enumerate the solutions of the Bethe equations (6) by their corresponding ideal-gas (i.e., ) quantum numbers , where are the quantized free single-particle momenta in the finite ring and is an integer222The energy of an eigenstate with for connects to the energy of the eigenstate with for . Here, are the quantum numbers of the “Fermi-sea” ground state for . In the remainder of this article, we will label states of the repulsive gas by their reduced quantum numbers . [39]. In this paper, in which we consider ground-state correlations and quenches from the ideal-gas ground state, we only need to consider eigenstates that are parity invariant, i.e., those for which we can order the such that for . Thus, we can label all eigenstates by quantum numbers , where is the floor function. By convention we choose these numbers to be the nonnegative values , which we regard as sorted in descending order (for odd , ).

Our results depend explicitly on the number of particles in our system, though the extent of our periodic geometry, and consequently the density of the gas, is arbitrary. We follow Refs. [18, 19] in absorbing the density into the dimensionless interaction-strength parameter . Our finite-sized system is then identified by the specification of both and . The Fermi momentum , which is the magnitude of the largest rapidity in the ground state in the Tonks–Girardeau limit of infinitely strong repulsive interactions [22], is a convenient unit of inverse length and so we specify lengths in units of , energies in units of , and times in units of .

### 2.2 Correlation functions

The static and dynamic behaviour of the Lieb-Liniger gas can be characterized by the normalized -order correlation functions

 g(m) (7)

where is the annihilation (creation) operator for the Bose field, is the particle-density operator, and denotes an expectation value with respect to a Schrödinger-picture density matrix . Due to the translational invariance of the system the density is constant [i.e., ] and the correlation functions are invariant under global coordinate shifts . Without loss of generality, we therefore set one of the spatial coordinates to zero and focus on the first-order correlation function , the second-order correlation function , and the local third-order correlation . We also consider the momentum distribution

 ˜n(k)=n∫L0dxe−ikxg(1)(x), (8)

which we evaluate at the discrete momenta .

For a system in a pure state , Eq. (7) reads

 g(m)(x1,…,xm,x′1,…,x′m;t)=1nm⟨ψ(t)|^Ψ†(x1)⋯^Ψ†(xm)^Ψ(x′1)⋯^Ψ(x′m)|ψ(t)⟩, =N!∫L0dxm+1⋯dxNnm(N−m)!ψ∗(x1,…,xm,xm+1,…,xN,t)ψ(x′1,…,x′m,xm+1,…,xN,t). (9)

By expressing the wave function in terms of Lieb–Liniger eigenstates [Eq. (2.1)], we can calculate the integrals in Eq. (2.2) semi-analytically with the methodology of Ref. [33]. This approach also allows for the evaluation of the overlaps of the initial state with Lieb–Liniger eigenstates necessary for our nonequilibrium calculations in Sec. 4333We note that direct evaluation of the normalization constant via Eq. (4) is susceptible to catastrophic cancellations similar to those discussed in Appendix B. In practice, we therefore obtain the constants by evaluating the self-overlaps of unnormalized Bethe eigenfunctions using the methodology of Ref. [33].. In Sec. 5, we consider the relaxed state of the system, as described by the diagonal-ensemble [73] density matrix , for which Eq. (7) reads

 g(m)DE(x1,…,xm,x′1,…,x′m)=1nmTr{^ρDE^Ψ†(x1)⋯^Ψ†(xm)^Ψ(x′1)⋯^Ψ(x′m)}, =1nm∑{λj}ρDE{λj}⟨{λj}|^Ψ†(x1)⋯^Ψ†(xm)^Ψ(x′1)⋯^Ψ(x′m)|{λj}⟩, =N!∑{λj}ρDE{λj}∫L0dxm+1⋯dxNnm(N−m)!ζ∗{λj}(x1,…,xm,xm+1,…,xN) ×ζ{λj}(x′1,…,x′m,xm+1,…,xN). (10)

### 2.3 Numerical considerations

For repulsive interactions the solutions to the Bethe equations (6) are characterized by purely real rapidities , and finding these numerically is relatively straightforward — see, e.g., Ref. [32]. However, for attractive interactions solutions with complex rapidities are possible, and the associated Yang-Yang action [21] of the problem is nonconvex (see, e.g., Ref. [43]), which significantly complicates the root-finding procedure.

To find the rapidities for attractive interactions, we start our root-finding routine close to . Here the rapidities are close to the free-particle momenta corresponding to , and these can be used as an initial guess for a Newton-method root finder. We then decrease in small steps, using linear extrapolation of the previous solutions to form initial guesses for the rapidities at each new value of . We have found that this procedure gives good convergence of the rapidities to machine precision.

Eigenstates with complex rapidities arrange themselves in so-called string patterns in the complex plane for large values of , with deviations from these strings exponentially small in the system length  [74, 23, 39, 42, 43]. For these states, some of the scattering factors in Eq. (3) become increasingly smaller with increasing , cancelling the extremely large exponential factor to give a finite result. Naïve evaluation of the wave function would therefore lead to numerical inaccuracies due to catastrophic cancellations as soon as the string deviations shrink to the order of machine precision. This problem can be overcome by using the Bethe equations (6) to rewrite the problematic factors in in terms of exponentials, thereby rendering the expressions more amenable to numerical calculation, as we discuss in Appendix B. For , this enables us to calculate correlation functions for attractive interaction-strength values using standard double-precision floating-point arithmetic, with the exception of a single eigenstate that we treat with high-precision arithmetic, as we discuss in Appendix B.3. For larger values of , the bound states become increasingly localized, leading to factors in Eq. (2.1) that are too large to be represented with double-precision floating-point arithmetic. We could in principle treat systems with through extensive use of high-precision arithmetic, but find that the regime to which we restrict our analysis reveals many important features of the physics of the attractively interacting system.

## 3 Ground-state correlation functions

The ground-state correlation functions of the one-dimensional Bose gas with attractive interactions have so far been investigated both in the mean-field regime [75, 45, 44] and with beyond-mean-field methodologies [44, 76, 48, 49, 46]. The corresponding Bose-Hubbard lattice approximation was considered in Ref. [77]. Systems in the limit were studied in Refs. [40, 51, 78, 79, 80], while in Ref. [81] correlation functions for up to particles under hard-wall boundary conditions were obtained via the coordinate Bethe ansatz. References [42, 43] used the algebraic Bethe ansatz to calculate the dynamic structure factor to first order in the string deviations under periodic boundary conditions. Piroli and Calabrese recently computed the local two- and three-body correlations in the limit where the interaction strength goes to zero as the system size increases at fixed particle density [41].

Here we compute exact correlation functions for a finite system of length with periodic boundary conditions and compare them with the predictions of mean-field theory, first for particles in the vicinity of the uniform-density to bright-soliton crossover , before considering more strongly attractive systems of particles with .

### 3.1 Correlations near the crossover

In Fig. 1 we plot the first- and second-order correlation functions of the ground state for particles for a range of . Figure 1(a) shows the first-order correlation in the spatial domain. For (red dashed line), the proximity to the noninteracting gas results in a nearly constant . For more attractive values of , begins to decay towards zero at larger separations . For (pink dot-dashed line), comes close to zero for , which corresponds to for . [Due to the periodic nature of our geometry, is symmetric around , and we therefore only show up to this point.]

Mean-field theory predicts a crossover from a uniform mean-field wave function to a localized bright-soliton state at an interaction strength  [44, 45, 46, 47]. In our exact quantum-mechanical treatment of the translationally invariant (and particle-number conserving) system, the density is necessarily constant. However, a signature of the bright-soliton-like state can be found in the first-order correlation function. In the finite-sized system the crossover is broad, but there is clearly a significant change in between [red dashed line in Fig. 1(a)] and (blue dot-dashed line). In the mean-field description, the many-body wave function is approximated by a translationally symmetrized Hartree-Fock product of single-particle wave functions [79, 80]. In this approximation correlation functions for the small system sizes we consider here are comparatively straightforward to compute numerically, see Appendix A for details.

Whereas the mean-field analysis predicts a sharp transition to the localized regime at the threshold interaction strength, the inclusion of quantum fluctuations leads to a smooth crossover between the delocalized and localized regimes in a system of finite  [44, 39]. To characterize the breadth of the crossover in our system, we calculate the single-particle entanglement entropy; i.e., the von Neumann entropy of the single-particle density matrix  [82]. In translationally invariant systems , where the are the momentum-mode populations.

In the (symmetrized) mean-field description, the ground state for is a pure product state, and hence . For , the ground state is a superposition of bright solitons, and  [46]. This can indeed be seen in the inset of Fig. 1(d), where we plot the single-particle entanglement entropy of the exact solution (black line) and of the mean-field solution (grey line) for particles. The mean-field entropy exhibits a slope discontinuity at the crossover point, whereas the von Neumann entropy of the exact ground state (black line) varies smoothly.

For the mean-field wave function is uniform, leading to a constant . In Fig. 1(a) we compare our exact results to the mean-field solution just on the localized side of the crossover at (green crosses), and find that the exact many-body solution (green dotted line) is slightly more localized. By contrast, for , i.e., further from the crossover point, the mean-field solution (blue diamonds) is more localized than the exact solution (blue dot-dashed line). For the mean-field solution (pink triangles) and the exact (pink dot-dot-dashed line) are reasonably similar, though the mean-field solution is again somewhat more localized than the exact solution. We note that this behaviour is consistent with that of the entanglement entropy [inset to Fig. 1(d)], which is smaller for the exact solution than for the mean-field approximation for . By contrast, at weaker interaction strengths finite-size rounding of the crossover yields an entropy for the exact system larger than the mean-field value.

In Fig. 1(c), we plot the momentum distribution corresponding to the first-order correlations shown in Fig. 1(a). [For our system and hence we only plot positive momenta.] We note that for all interaction strengths we consider here, the exact momentum distributions exhibit a power-law decay at high momenta — the universal large-momentum behaviour for systems with short-range interactions [83, 84, 85]. For the case of (red empty circles), interactions are sufficiently weak that no visible deviation from this scaling is visible at the smallest nonzero momenta resolvable in our finite geometry. By contrast, for (green triangles), less trivial behaviour of the momentum distribution can be seen, with the lowest nonzero momentum modes deviating visibly from the scaling. As increases, the deviations from this scaling extend to higher momenta, and a broad hump in the momentum distribution develops. This broadening can be more clearly seen in Fig. 1(d), where we plot the momentum distribution for low momenta on a linear scale. For (red empty circles), the zero-momentum occupancy is close to its ideal-gas value of . The zero-momentum mode occupation decreases with increasing and much of this population is redistributed to the first few nonzero momentum modes, resulting in, e.g., a broad distribution for (pink empty squares).

The ground-state mean-field momentum distributions in Fig. 1(c) do not show the scaling for large — this feature appears with a first-order Bogoliubov analysis [86]. For an interaction strength , i.e., close to the crossover point, the exact (green dotted line) and the mean-field solution (green crosses) are clearly different away from . For larger attractive values of , however, the two momentum distributions start to agree more closely. For example, from Figs. 1(c) and 1(d) we observe reasonable agreement between the exact and mean-field results for the lowest three modes at (blue diamonds for mean-field solution, blue dot-dashed line for exact solution). Even closer agreement is observed for , where the lowest six modes of the exact solution (pink dot-dot-dashed line) agree well with the mean-field solution (pink triangles), before the tail of the exact momentum distribution takes over.

In Fig. 1(b), we plot the second-order correlation for the same values of as before. For (red dashed line), is close to the ideal-gas value (horizontal grey line). For (green dotted line), is increased over the ideal-gas value at distances and correspondingly decreased at larger distances. This behaviour is even more pronounced for (blue dot-dashed line), and the trend continues for larger attractive values of , for which there is significant bunching of particles. Comparing the exact results to the mean-field solutions, we again observe a clear difference at , where the exact solution (green dotted line) is more localized than the mean-field solution (green crosses). For , the exact solution (blue dot-dashed line) has a slightly increased value at zero separation compared to the mean-field solution (blue diamonds), but at intermediate separations the latter is marginally broader. For , the local value of the exact solution (pink dot-dot-dashed line) is again slightly larger than the mean-field value (pink triangles). At separations , the mean-field and exact distributions show good agreement.

### 3.2 Correlations for strongly interacting systems

In Fig. 2, we plot the first- and second-order correlation functions of the ground state for particles and for a larger range of values of the interaction strength . For , the mean-field critical interaction strength is , and all ground states we consider here are therefore well in the localized regime. Figure 2(a) indicates the first-order correlation function , which shows that the soliton-like state becomes increasingly tightly localized with increasing . This can also be observed in momentum space, Fig. 2(c), where the corresponding momentum distributions become broader with increasing . We note that the momentum distributions for the most strongly interacting systems considered here are much broader than the “hump” that forms in the ground-state momentum distribution of the repulsive gas in the strongly interacting Tonks limit, which extends to  [87, 88, 33]. For comparison, we also plot the mean-field correlation functions for in Figs. 2(a) and (c) (grey diamonds). The mean-field first-order correlation function is similar to that of the exact solution but slightly more localized, and its momentum distribution is correspondingly somewhat broader than the exact distribution for small values of . Nevertheless, the two momentum distributions agree well over a wide range of momenta up to , where the universal scaling of the exact momentum distribution begins.

Figure 2(b) shows the second-order correlation for separations up to (which corresponds to for ). We again observe that the system becomes more tightly bound with increasingly attractive interactions. In order to ensure that the form of the correlation function at moderate separations is visible in this figure, we have limited the extent of the  axis. The maximum value of the second-order correlation function for (solid black line), , is therefore not shown. The mean-field correlation function for (grey diamonds) shows good agreement with the exact solution, though its value at zero separation (not shown) is reduced compared to that of the exact solution.

Figure 2(d) shows the local second- and third-order correlations for a wide range of interaction strengths. For small values of , these correlations are close to their respective ideal-gas values, and  [89]. In the vicinity of the mean-field crossover point (indicated by the vertical grey line), both and begin to increase significantly with increasing . For larger values of , we observe a linear scaling of the second-order correlation and a quadratic scaling of the third-order correlation , both of which we indicate by black dot-dashed lines in Fig. 2(d). The former scaling can be understood by noting that the McGuire cluster energy scales as  [37], and that  [90].

In summary, the exact finite-system correlation functions show behaviour consistent with a broad crossover around the mean-field critical value. At stronger interactions, our exact results for small atom numbers are in close agreement with the predictions of mean-field theory.

## 4 Dynamics following an interaction quench

In this section we investigate the nonequilibrium evolution of the attractively interacting Lieb–Liniger gas following an interaction quench for particles at time . Initially the system is prepared in the ideal-gas ground state, for which the wave function is constant in space, . Formally, the state of the system at time is given by

 |ψ(t)⟩=∑{λj}C{λj}e−iE{λj}t|{λj}⟩, (11)

where the are the overlaps of the initial state with the Lieb–Liniger eigenstates at the postquench interaction strength , and the are the corresponding energies. The evolution of equal-time correlation functions (Sec. 2.2) is calculated by noting that the time evolution of the expectation value of an arbitrary operator in the time-dependent state is given by

 ⟨^O(t)⟩ ≡⟨ψ(t)|^O|ψ(t)⟩ (12) =∑{λj}∑{λ′j}C∗{λ′j}C{λj}ei(E{λ′j}−E{λj})t⟨{λ′j}|^O|{λj}⟩.

The matrix elements and overlaps are calculated with the method described in Ref. [33].

Numerically it is necessary to truncate the infinite sum in Eq. (12), and our truncation procedure is analogous to that described in Appendix A of Ref. [32]: we include all eigenstates for which the populations are larger than some threshold value, thereby minimizing the normalization sum-rule violation for the corresponding basis size. For calculations of and for interaction-strength quenches to we use a cutoff , leading to a sum-rule violation of . All other correlation functions are calculated with a more stringent cutoff , and the sum-rule violations are correspondingly smaller. We have checked that increasing the cutoff does not visibly alter any of our results.

### 4.1 Influence of bound states following a quench

Before investigating the detailed nonequilibrium dynamics of the Lieb–Liniger gas following a quench to attractive interactions, we first consider the populations of the eigenstates of the postquench Hamiltonian, which are constant at all times [cf. Eq. (11)]. Comparing these populations to those resulting from quenches to repulsive interactions helps provide an understanding of the contribution of bound states to the nonequilibrium dynamics in the attractive case.

In Fig. 3 we plot the populations of several representative Lieb-Liniger eigenstates following quenches of the interaction strength from zero to a wide range of final interaction strengths . [Recall from Sec. 2.1 that for there are two independent to be specified, which we indicate by the legend in Fig. 3(b)444Note that for repulsive interactions the quantum-number pairs quoted here refer to the “reduced” quantum numbers, i.e., the excitation numbers relative to the Fermi-sea ground state (cf. Sec. 2.1)..] For attractive interactions [Fig. 3(a)] several eigenstates containing bound states have significant populations for small values of . (Note that the number of particles in the bound state can be inferred from the distribution of the rapidities in the complex plane.) The populations of the ground state (red solid line), which is a four-particle bound state, and the three-particle bound state (green dotted line) are dominant for quenches to . However, their populations decrease rapidly with increasing absolute interaction strength beyond .

At intermediate interaction strengths , two-body bound states start to dominate the populations [e.g.,  the states with (blue dot-dashed line) and (pink dot-dot-dashed line)]. For increasingly attractive values of , the populations of gas-like states with no bound-state component grow [e.g.,  (black solid line) and (pink dotted line)]. Indeed, at , the population of the super-Tonks state — the lowest-energy gas-like state at strong interactions — begins to dominate. However, the two-body bound state with (blue dot-dashed line) still has a significant population in the strongly interacting regime555We note that at this state has an energy of , which is close to the energy of the two-particle McGuire cluster state with  [37].. Consequently, we expect bound states to influence the dynamical evolution of correlation functions following a quench from the ideal gas to all attractive interaction strengths that we consider. Comparing the populations of eigenstates for attractive postquench interactions to those for repulsive interactions, Fig. 3(b), we can see that there is significantly less structure in the latter, which are all gas-like. We observe that the populations of excited gas-like eigenstates increase monotonically with increasing for both repulsive and attractive interactions, whereas the results of Fig. 3(a) suggest that the populations of the eigenstates containing bound states all eventually decrease as . We note that although scattering states of the attractive gas connect adiabatically to states of the repulsive gas in the limit  [66], the quantum-number labels of the states differ on either side of the infinite-interaction-strength limit. For example, for particles, the super-Tonks state with connects on to the ground state for repulsive interactions, .

To better understand the eigenstate contributions to the nonequilibrium dynamics following a quench to attractive interactions, we focus on quenches of particles from the ideal-gas ground state to attractive and repulsive interactions with , and plot in Fig. 4 the populations of the contributing eigenstates against their energies . We see that there are additional families of populated states for the attractive gas (sequences of blue crosses that extend to negative energies) that are not present for the repulsive gas (red circles). These are due to four different types of contributing bound states, which we now describe.

The first two types of bound states are four-body and three-body bound states, and each of these types contains only a single populated state. These are, respectively, the ground state at with and the first parity-invariant excited state at with . We note that the parity invariance of eigenstates for quenches from the initial ideal gas [28] restricts the appearance of bound states with more than two bound particles to only these two states.

The third type is represented by the eigenstate with , which has two bound particles and two free particles, and is the first in a family of similar states ( a nonnegative integer) whose populations decrease gradually with increasing . The fourth type is represented by the eigenstate with , which contains two two-particle bound states, and is the first in a family with decreasing populations for higher excitations which alternate between the quantum numbers and , with a positive integer. For larger , the two two-body bound states have higher “centre-of-mass” momenta with opposite sign (recall that only eigenstates with total momentum have nonzero occupations following the quench), and for the corresponding positive centre-of-mass energy of the pairs exceeds their binding energy.

We can see from Fig. 4 that the distributions of populations over gas-like eigenstates are similar for quenches to , aside from a shift in energy and a small decrease in populations for the attractive gas due to the appearance of the additional bound states. In particular, the number of eigenstates with populations is () for the attractive (repulsive) gas. The shift in energy can be explained by noting that for , the system is in the strongly interacting regime and the Bethe rapidities of scattering states (i.e. states with no bound particles) can be obtained by a strong-coupling expansion around the Tonks–Girardeau limit of infinitely strong interactions (see, e.g., Ref. [91]). This yields , where the are the Tonks–Girardeau values, implying opposite energy shifts in the attractive and repulsive cases.

### 4.2 Dynamics of local correlations

We now consider the nonequilibrium dynamics following the quench. In Fig. 5(a) we plot the local second-order correlation for particles following a quench from to four representative final interaction strengths. Initially, (cf. Sec. 3.1). For a quench to (pink dot-dashed line), shows nearly monochromatic oscillatory behaviour. This is similar to the behaviour following quenches to small repulsive interaction strengths analyzed in Ref. [32]. Because the difference between the postquench energy  [92, 32] and the ground-state energy of the system is small compared to the finite-size energy gap to the first (parity-invariant) excited state, the ensuing dynamics are dominated by these two states, and the energy difference between them determines the dominant frequency of the oscillations.

Quenches to more attractive values of show the generic behaviour of an initially rising that eventually fluctuates about a seemingly well-defined average value. The frequencies of the oscillations are determined by the energy differences between the Lieb-Liniger eigenstates with the largest populations. For example, for (solid red line), the postquench wave function is dominated by the super-Tonks state and the first two-body bound state, cf. Fig. 3, and the dominant frequency in the oscillations at early times matches the energy difference between these two eigenstates. At later times, the shape of is more irregular, but the large oscillations due to the two dominant eigenstates persist.

In Fig. 5(b) we plot the local third-order correlation for particles following a quench from to the same four final interaction strengths as before. Initially, (see Sec. 3.2). For small postquench interaction strengths, (pink dot-dashed line) and (blue dashed line), the evolution is similar to that of for the same interaction strengths. For larger attractive values of the postquench interaction strength, on the other hand, the shape of is more regular compared to , reflecting the fact that only one three-body bound state contributes to the postquench wavefunction, whereas multiple states containing bound pairs are present. Indeed for (green dotted line) and (solid red line), is dominated by a single frequency, given by the energy difference between the three-body bound state and the predominant two-body bound state . The initial rise of both and terminates on an increasingly shorter time scale with increasingly attractive postquench interaction strength. This time scale corresponds to about half the period of the ensuing oscillations and is proportional to , corresponding to the scaling of the energy of eigenstates containing bound states [37].

For quenches from the ideal-gas initial state, we find that the population of the bound states leads to significantly increased values of both and — in stark contrast to the decay of the same quantities following quenches to repulsive interactions [32] due to the “fermionization” of the system. Such large values of these local correlation functions would lead to strong particle losses in experiments [7, 93, 94]. This is in contrast to the observations in the quench experiments performed in Ref. [7], where the quasi-one-dimensional gas was quenched from strongly repulsive interactions to strongly attractive interactions, and no significant losses were observed. In such a scenario the overlap of the initial strongly repulsive ground state with the super-Tonks state is dominant, and the bound states thus acquire only small populations in the course of the quench [7, 66, 64, 63].

To investigate the influence of the initial state on the populations of the two most dominant postquench eigenstates (cf. Fig. 3), we find the (correlated) ground state of the system at and then compute the populations of the eigenstates following a quench to . In Fig. 6, we plot the populations and of the aforementioned two-body bound state and the super-Tonks state, respectively, for a wide range of initial values . Starting in the strongly interacting regime , the overlap between the initial (Tonks–Girardeau) state and the super-Tonks state is close to unity. As is decreased, the population of the super-Tonks gas decreases, while the population of the bound state increases. At , the two populations are already near their respective values following a quench from the ideal-gas initial state (indicated by black arrows on the left-hand side). The results of Fig. 6 suggest that the postquench values of and would be much smaller for quenches from initial values of compared to those from the noninteracting initial state.

### 4.3 Dynamics of the momentum distribution

We now turn our attention to the postquench dynamics of the momentum distribution. Quenches from the ideal-gas ground state with particles to three different values of are compared in Fig. 7. In each case we plot the time evolution of the momentum-mode occupations [cf. Eq. (8)] for the first six nonnegative momentum modes (). Initially, all particles occupy the zero-momentum single-particle orbital, . At times , the interaction quench leads to a redistribution of this population over other single-particle modes. At early times, all nonzero modes rise with the same rate, independent of , due to the local nature of the interaction potential, which corresponds to a momentum-independent coupling [95]. This applies to all postquench interaction strengths , but the time at which deviations from this behaviour first appear depends on .

All quenches show the same generic behaviour — the momentum-mode populations eventually level off and fluctuate about a well-defined value. These populations undergo oscillations with frequencies determined by the energy differences between the dominant Lieb-Liniger eigenstates. For example, for the case of Fig. 7(c) each mode exhibits fast oscillations at a single frequency given by the energy difference between the super-Tonks state and the two-body bound state , superposed with some irregular envelope function.

In Fig. 8, we compare for quenches from the ideal gas to repulsive and attractive interaction strengths of the same magnitude. In Fig. 8(a), we plot the time evolution of the zero-momentum mode occupation for quenches from to (solid red line) and (blue dashed line). The envelope of for attractive interactions is similar to the shape of for repulsive interactions. On top of this envelope for quenches to attractive interactions, shows large regular oscillations. This also applies for quenches to , Fig. 8(b), but the oscillations for quenches to (solid red line) are faster than for quenches to . The correspondence between following a quench to strong attractive interactions and that following a quench to equally strong repulsive interactions reflects the fact that the two postquench wave functions are similar in their composition, aside from the additional presence of two-body bound states for attractive interactions, as illustrated in Fig. 4.

We also observe a partial revival in for quenches to . This revival is due to the proximity of the system at to the Tonks–Girardeau limit of infinitely strong interactions, where the spectrum of the repulsive Lieb–Liniger model is identical to that of free fermions [96]. This also applies to the scattering states of the attractive system. For , this would lead to recurrences at integer multiples of  [32] due to the commensurability of eigenstate energies [97]. However, for the finite interaction strengths considered here, the revival time is shifted to a later time for repulsive interactions [32] and to an earlier time for attractive interactions, due to the finite-coupling corrections to the Bethe rapidities discussed in Sec. 4.1.

### 4.4 Dynamics of nonlocal pair correlations

We now consider the evolution of the full nonlocal second-order correlation . In Fig. 9 we plot the behaviour of this quantity for an interaction quench from zero to for particles. Figure 9(a) shows at four representative times . Initially, (horizontal line). At (red dashed line), the local value is already greatly enhanced, , cf. Fig. 5(a). [The scale of the  axis is chosen so that the long-range features of are visible, and the large values for are therefore cut off.] In addition to the central peak, at separations a secondary peak emerges, while at larger distances exhibits a decaying oscillatory structure. As time progresses, this secondary peak propagates away from the origin and broadens as can be seen at, e.g. (green dotted line) and (blue dot-dashed line).

The build-up of this secondary correlation peak and its propagation through the system can be more clearly seen in Fig. 9(b), where we plot the time-evolution of up to . The propagation of this peak is consistent with , which was also observed for quenches from the same initial state to strongly repulsive interactions [32, 29]. (Note that the color scale is chosen so that the long-range behaviour is visible, and the local second-order correlation is again not resolved.) Figure 9(c) shows for longer times up to . The overall structure on this longer time scale is more complicated, with several soliton-like correlation dips propagating through the system [32] and a partial revival of at [cf. Figs. 7(c) and 8(b)]. Besides the largely increased value at small distances, the behaviour of is strikingly similar to the results obtained in Ref. [32] for quenches from the same noninteracting ground state to repulsive final interaction strengths.

In summary, quenches from the ideal-gas ground state to attractive values of result in the occupation of energy eigenstates containing bound states in addition to the gas-like scattering states of the attractively interacting model, which are analogous to the eigenstates of the repulsively interacting Lieb–Liniger gas. As the magnitude of the final interaction strength is increased, the postquench occupations of the gas-like excited states approach those of their counterparts following a quench to the corresponding repulsive interaction strength, and the occupations of bound states eventually decrease. However, these bound states significantly influence the dynamics of postquench correlation functions for all final interaction strengths we have considered, causing large oscillations in local correlations and in the occupation of the zero-momentum mode. For large attractive values of , bound states are highly localized and thus influence the second-order correlation function only at small separations, whereas at larger separations this function exhibits postquench dynamics similar to those observed following quenches to repulsive interactions [32].

## 5 Time-averaged correlations

A closed quantum-mechanical system prepared in a pure state will remain in a pure state for all time. However, for a nondegenerate postquench energy spectrum, as is the case here (cf. Refs. [32, 33]), the energy eigenstates will dephase, and the time-averaged expectation value of any operator can be expressed in terms of its diagonal matrix elements between energy eigenstates

 ⟨^O⟩DE =limτ→∞1τ∫τ0dt⟨ψ(t)|^O|ψ(t)⟩, =∑{λj}|C{λj}|2⟨{λj}|^O|{λj}⟩. (13)

This quantity can be viewed as the expectation value of in the diagonal-ensemble density matrix [73]

 ^ρDE =∑{λj}|C{λj}|2|{λj}⟩⟨{λj}|. (14)

We note that in practice the sum in Eq. (14) runs over a finite set of energy eigenstates with populations exceeding some threshold value. If the expectation value of an operator relaxes at all, it must relax to the corresponding diagonal-ensemble value [98]. Although expectation values may exhibit rather large fluctuations around their time-averaged values for system sizes as small as those considered here, in general the relative magnitude of these fluctuations should decrease with increasing system size and vanish in the thermodynamic limit. However, establishing this behaviour is beyond the scope of the current work and we will simply regard the diagonal ensemble defined by Eq. (14) as the ensemble appropriate to describe the relaxed state of our finite-sized system. In the following we consider the time-averaged properties of the quenched system.

### 5.1 Local correlations

In Fig. 10(a), we plot the enhancement of the diagonal-ensemble value of the local second-order correlation over the initial noninteracting value of this function following an interaction quench from zero to for particle numbers , , and . For all particle numbers considered, as is increased from the ideal-gas limit, initially increases rapidly before reaching a local maximum, which occurs at smaller values of for larger particle numbers . For particles (solid blue line) this local maximum in occurs at and coincides with the crossing of the population of the three-particle bound state and that of the ground state [see Fig. 3(a)]. The local minimum of at coincides with the maximum population of this three-particle bound state, and as soon as the population of this state starts to decrease, begins to increase monotonically with increasing .

For large attractive values of , the local second-order correlation tends to a constant value , which is much larger than the ideal gas and super-Tonks values [67]. The decrease of with increasing particle number at fixed large appears consistent with an approach toward the quench-action thermodynamic-limit strong-coupling value obtained to third order in in Refs. [35, 36], indicated by the solid grey line, as .

Using the quench-action approach [71] in the thermodynamic limit, Refs. [35, 36] found that for . Our methodology does not recover this result for small values of , as our small system sizes lead to a finite-size gap for excitations and therefore the energy added by the quench is small in this case. Additionally, eigenstates with more than four bound particles are trivially absent in our calculations, whereas for small postquench values of they contribute significantly in the analysis of Refs. [35, 36]. For larger values of , however, states with more than two bound particles are strongly suppressed and we expect our results to be less influenced by finite-size effects [33].

In Fig. 10(b), we plot the enhancement of the diagonal-ensemble value of the local third-order correlation over its noninteracting initial value following an interaction quench from zero to for particle numbers and . The qualitative behaviour is similar to that of . For strong interactions, also tends to a constant value that is much larger than the initial value. Whether this result persists for larger atom numbers is an important open question, given that large values of lead to strong recombination losses in experiments with ultracold gases [93, 94].

### 5.2 Nonlocal correlations

In Fig. 11(a) we plot the momentum distribution in the diagonal ensemble for particles and for several postquench interaction strengths . At high momenta and for all interaction strengths , exhibits a scaling of . This behaviour is due to the universal character of short-range two-body interactions [83, 84, 85]. For (pink squares), the functional form of is nearly perfectly given by this scaling, and only the three lowest resolvable nonzero momentum modes in our finite periodic system deviate slightly from it.

For a quench to (blue filled circles), the low-momentum part of