Quantum quenches to the attractive one-dimensional Bose gas: exact results

We study quantum quenches to the one-dimensional Bose gas with attractive interactions in the case when the initial state is an ideal one-dimensional Bose condensate. We focus on properties of the stationary state reached at late times after the quench. This displays a finite density of multi-particle bound states, whose rapidity distribution is determined exactly by means of the quench action method. We discuss the relevance of the multi-particle bound states for the physical properties of the system, computing in particular the stationary value of the local pair correlation function $g_2$.


Introduction
Strongly correlated many-body quantum systems are often outside the range of applicability of standard perturbative methods. While being at the root of many interesting and sometimes surprising physical effects, this results in huge computational challenges, which are most prominent in the study of the non-equilibrium dynamics of many-body quantum systems. This active field of research has attracted increasing attention over the last decade, also due to the enormous experimental advances in cold atomic physics [1][2][3]. Indeed, highly isolated many-body quantum systems can now be realised in cold atomic laboratories, where the high experimental control allows to directly probe their unitary time evolution [4][5][6][7][8][9][10][11][12][13].
A simple paradigm to study the non-equilibrium dynamics of closed many-body quantum systems is that of a quantum quench [14]: a system is prepared in an initial state (usually the ground state of some Hamiltonian H 0 ) and it is subsequently time evolved with a local Hamiltonian H. In the past years, as a result of a huge theoretical effort (see the reviews [2,[15][16][17][18][19][20][21] and references therein), a robust picture has emerged: at long times after the quench, and in the thermodynamic limit, expectation values of local observables become stationary. For a generic system, these stationary values are those of a thermal Gibbs ensemble with the effective temperature being fixed by the energy density in the initial state [22].
A different behaviour is observed for integrable quantum systems, where an infinite set of local conserved charges constrains the non-equilibrium dynamics. In this case, long times after the quench, local properties of the systems are captured by a generalised Gibbs ensemble (GGE) [23], which is a natural extension of the Gibbs density matrix taking into account a complete set of local or quasi-local conserved charges.
The initial focus was on the role played by (ultra-)local conservation laws in integrable quantum spin chains [24][25][26][27][28][29][30], while more recent works have clarified the role by sets of novel, quasi-local charges [31][32][33][34][35][36][37][38][39][40][41]. It has been shown recently that they have to be taken into account in the GGE construction in order to obtain a correct description of local properties of the steady state [44,45]. Quasi-local conservation laws and their relevance for the GGE have also recently been discussed in the framework of integrable quantum field theories [48,49]. These works have demonstrated that the problem of determining a complete set of local or quasi-local conserved charges is generally non-trivial.
A different approach to calculating expectation values of local correlators in the stationary state was introduced in Ref. [50]. It is the so called quench action method (QAM), a.k.a. representative eigenstate approach and it does not rely on the knowledge of the conserved charges of the system. Within this method, the local properties at large times are effectively described by a single eigenstate of the post-quench Hamiltonian.
One of the most interesting aspects of non-equilibrium dynamics in integrable systems is the possibility of realising non-thermal, stable states of matter by following the unitary time evolution after a quantum quench. Indeed, the steady state often exhibits properties that are qualitatively different from those of thermal states of the post-quench Hamiltonian. The QAM provides a powerful tool to theoretically investigate these properties in experimentally relevant settings.
In this paper we study the quantum quench from an ideal Bose condensate to the Lieb-Liniger model with arbitrary attractive interactions. A brief account of our results was previously given in Ref. [56]. The interest in this quench lies in its experimental feasibility as well as in the intriguing features of the stationary state, which features finite densities of multi-particle bound states. Our treatment, based on the quench action method, allows us to study their dependence on the final interaction strength and discuss their relevance for the physical properties of the system. In particular, as a meaningful example, we consider the local pair correlation function g 2 , which we compute exactly.
The structure of the stationary state is very different from the super Tonks-Girardeau gas, which is obtained by quenching the one-dimensional Bose gas from infinitely repulsive to infinitely attractive interaction [64][65][66][67][68][69][70]. The super Tonks-Girardeau gas features no bound states, even though it is more strongly correlated than the infinitely repulsive Tonks-Girardeau gas, as has been observed experimentally [66]. As we argued in [56], the physical properties of the post-quench stationary state reached in our quench protocol could be probed in ultracold atoms experiments, and the multi-particle bound states observed by the presence of different"lightcones" in the spreading of local correlations following a local quantum quench [71].
In this work we present a detailed derivation of the results previously announced in Ref. [56]. The remainder of this manuscript is organised as follows. In section 2 we introduce the Lieb-Liniger model and the quench protocol that we consider. The quench action method is reviewed in section 3, and its application to our quench problem is detailed. In section 4 the equations describing the post-quench stationary state are derived. Their solution is then obtained in section 5, and a discussion of its properties is presented. In section 6 we address the calculation of expectation values of certain local operators on the post-quench stationary state, and we explicitly compute the local pair correlation function g 2 . Finally, our conclusions are presented in section 7. For the sake of clarity, some technical aspects of our work are consigned to several appendices.

The Hamiltonian and the eigenstates
We consider the Lieb-Liniger model [72], consisting of N interacting bosons on a one-dimensional ring of circumference L. The Hamiltonian reads where m is the mass of the bosons, and c = − 2 /ma 1D is the interaction strength. Here a 1D is the 1D effective scattering length [73] which can be tuned via Feshbach resonances [74]. In the following we fix = 2m = 1. The second quantized form of the Hamiltonian is where Ψ † , Ψ are complex bosonic fields satisfying [Ψ(x), Ψ † (y)] = δ(x − y). The Hamiltonian (1) can be exactly diagonalised for all values of c using the Bethe ansatz method [42,72]. Throughout this work we will consider the attractive regime c < 0 and use notations c = −c > 0. We furthermore define a dimensionless coupling constant by A general N -particle energy eigenstate is parametrized by a set of N complex rapidities {λ j } N j=1 , satisfying the following system of Bethe equations The wave function of the eigenstate corresponding to the set of rapidities {λ j } N j=1 is then where the sum is over all the permutations of the rapidities. Eqns (4) can be rewritten in logarithmic form as where the quantum numbers {I j } N j=1 are integer (half-odd integer) for N odd (even). In the attractive regime the solutions of (6) organize themselves into mutually disjoint patterns in the complex rapidity plane called "strings" [75,76]. For a given N particle state, we indicate with N s the total number of strings and with N j the number of j-strings, i.e. the strings containing j particles (1 ≤ The rapidities within a single j-string are parametrized as [87] λ j,a α = λ j α + where a labels the individual rapidities within the j-string, while α labels different strings of length j. Here λ j α is a real number called the string centre. The structure (8) is common to many integrable systems and within the so called string hypothesis [75,77] the deviations from a perfect string δ j,a α are assumed to be exponentially vanishing with the system size L (see Refs. [78,79] for a numerical study of such deviations in the Lieb-Liniger model). A j-string can be seen to correspond to a bound state of j bosons: indeed, one can show that the Bethe ansatz wave function decays exponentially with respect to the distance between any two particles in the bound state and the j bosons can be thought as clustered together.
Even though some cases are known where states violating the string hypothesis are present [80][81][82][83][84], it is widely believed that their contribution to physically relevant quantities is vanishing in the thermodynamic limit. We will then always assume the deviations δ j,a α to be exponentially small in L and neglect them except when explicitly said otherwise. From (6), (8) a system of equations for the string centres λ j α is obtained [76] jλ j α L − where and where I j α are integer (half-odd integer) for N odd (even). Eqns (9) are called Bethe-Takahashi equations [75,85]. The momentum and the energy of a general eigenstate are then given by

The thermodynamic limit
In the repulsive case the thermodynamic limit was first considered in Ref. [86], and it is well studied in the literature. In the attractive case, the absolute value of the ground state energy in not extensive, but instead grows as N 3 [87,88]. While ground state correlation functions can be studied in the zero density limit, namely N fixed, L → ∞ [76], it was argued that the model does not have a proper thermodynamic limit in thermal equilibrium [75,88]. Crucially, in the quench protocol we are considering, the energy is fixed by the initial state and the limit of an infinite number of particles at fixed density presents no problem.
As the systems size L grows, the centres of the strings associated with an energy eigenstate become a dense set on the real line and in the thermodynamic limit are described by smooth distribution function. In complete analogy with the standard finite-temperature formalism [75] we introduce the distribution function {ρ n (λ)} ∞ n=1 describing the centres of n strings, and the distribution function of holes {ρ h n (λ)} ∞ n=1 . We recall that ρ h n (λ) describes the distribution of unoccupied states for the centres of n-particle strings, and is analogous to the distribution of holes in the case of ideal Fermi gases at finite temperature. Following Takahashi [75] we introduce In the thermodynamic limit the Bethe-Takahashi equations (9) reduce to an infinite set of coupled, non-linear integral equations where a nm (λ) = (1 − δ nm )a |n−m| (λ) + 2a |n−m|+2 (λ) + . . . + 2a n+m−2 (λ) + a n+m (λ) , a n (λ) = 1 2π In the thermodynamic limit the energy and momentum per volume are given by where Finally, it is also useful to define the densities D n and energy densities e n of particles forming n-strings The total density and energy per volume are then additive
As we mentioned before, the energy after the quench is conserved and is most easily computed in the initial state |ψ(0) = |BEC as The expectation value on the r.h.s. can then be easily computed using Wick's theorem. In the thermodynamic limit we have 3 The quench action method

General considerations
Consider the post-quench time evolution of the expectation value of a general operator O. For a generic system it can be written as where {|µ } denotes an orthonormal basis of eigenstates of the post-quench Hamiltonian. In Ref. [50] it was argued that in integrable systems a major simplification occurs if one is interested in the time evolution of the expectation values of local operators O in the thermodynamic limit. In particular, the double sum in the spectral representation (25) can be replaced by a single sum over particle-hole excitations over a representative eigenstate |ρ sp . In particular, we have where we have indicated with lim th the thermodynamic limit N, L → ∞, keeping the density D = N/L fixed.
Here e denotes a generic excitation over the representative state |ρ sp . Finally we have where ω[ρ sp ], ω[ρ sp , e] are the energies of |ρ sp and |ρ sp , e respectively. The representative eigenstate (or "saddle-point state") |ρ sp is described in the thermodynamic limit by two sets of distribution functions {ρ n (λ)} n , {ρ h n (λ)} n . In Ref. [50] it was argued that these are selected by the saddle-point condition where S QA [ρ] is the so-called quench action Here ρ is the set of distribution functions corresponding to a general macro-state, S[ρ] gives the thermodynamically leading part of the logarithm of the overlap and S Y Y is the Yang-Yang entropy. As we will see in section 3.2, we will only have to consider parityinvariant Bethe states, namely eigenstates of the Hamiltonian (1) characterised by sets of rapidities satisfying Restricting to the sector of the Hilbert space of parity invariant Bethe states, the Yang-Yang entropy reads We note the global pre-factor 1/2. From Eq. (26) it follows that the saddle-point state |ρ sp can be seen as the effective stationary state reached by the system at long times. Indeed, if O is a local operator, Eq. (26) gives lim

Overlaps with the BEC state
The main difficulty in applying the quench action method to a generic quantum quench problems is the computation of the overlaps ψ(0)|ρ between the initial state and eigenstates of the post-quench Hamiltonian. At present this problem has been solved only in a small number of special cases [50,[103][104][105][106][107][108][109][110][111].
A conjecture for the overlaps between the BEC state and the Bethe states in the Lieb-Liniger model first appeared in Ref. [51] and it was then rigorously proven, for arbitrary sign of the particle interaction strength, in Ref. [106]. As we have already mentioned, the overlap is non-vanishing only for parity invariant Bethe states, namely eigenstates characterised by sets of rapidities satisfying [105]. The formula reads where The extensive part of the logarithm of the overlap (33) was computed in Ref. [51] in the repulsive regime. A key observation was that the ratio of the determinants is non-extensive, i.e.
In the attractive regime additional technical difficulties arise, because matrix elements of the Gaudin-like matrices G jk , G Q jk can exhibit singularities when the Bethe state contains bound states [111]. This is analogous to the situation encountered for a quench from the Néel state to the gapped XXZ model [57][58][59][60]. In particular, one can see that the kernel K(µ − ν) diverges as 1/(δ n,a α − δ n,a+1 α ) for two "neighboring" rapidities in the same string µ = λ n,a α , ν = λ n,a+1 α , or when rapidities from different strings approach one another in the thermodynamic limit, µ → λ + ic.
These kinds of singularities are present in the determinants of both G Q jk and G jk . It was argued in Refs [57,58,111] that they cancel one another in the expression for the overlap. As was noted in Refs. [57,58,111], no other singularities arise as long as one considers the overlap between the BEC state and a Bethe state without zero-momentum n-strings, (strings centred at λ = 0). Concomitantly the ratio of the determinants in (33) is expected to give a sub-leading contribution in the thermodynamic limit, and can be dropped. The leading term in the logarithm of the overlaps can then be easily computed along the lines of Refs. [57,58] where In the case where zero-momentum n-strings are present, a more careful analysis is required in order to extract the leading term of the overlap (33) [111,112]. This is reported in Appendix A. The upshot of this analysis is that (38) gives the correct leading behaviour of the overlap even in the presence of zero-momentum n-strings.
4 Stationary state

Saddle point equations
As noted before, the stationary state is characterized by two sets of distribution functions {ρ n (λ)} n , {ρ h n (λ)} n , which fulfil two infinite systems of coupled, non-linear integral equations. The first of these is the thermodynamic version of the Bethe-Takahashi equations (16). The second set is derived from the saddle-point condition of the quench action (28), and the resulting equations are sometimes called the overlap thermodynamic Bethe ansatz equations (oTBA equations). Their derivation follows Refs [57][58][59][60]. In order to fix the density D = N/L we add the following term to the quench action (29) As discussed in the previous section, S [ρ] in (29) can be written as where W n (λ) is given in (39). Using (41), (31), and (40) one can straightforwardly extremize the quench action (29) and arrive at the following set of oTBA equations Here a nm are defined in (17), and we have used the notation f * g(λ) to indicate the convolution between two functions Eqns (42) determine the functions η n (λ) and, together with Eqns (16) completely fix the distribution functions {ρ n (λ)} n , {ρ h n (λ)} n characterising the stationary state.

Tri-diagonal form of the oTBA equations
Following standard manipulations of equilibrium TBA equations [75], we may re-cast the oTBA equations (42) in the form Here we have defined η 0 (λ) = 0 and The calculations leading to Eqns (44) are summarized in Appendix B. The thermodynamic form of the Bethe-Takahashi equations (16) can be similarly rewritten. Since we do not make explicit use of them in the following, we relegate their derivation to Appendix B.

Asymptotic relations
Eqns (44) do not fix {η n (λ)} n of Eqns (42), because they do not contain the chemical potential h. In order to recover the (unique) solution of Eqns (42), it is then necessary to combine Eqns (44) with a condition on the asymptotic behaviour of η n (λ) for large n. In our case one can derive from (42) the following relation, which holds asymptotically for n → ∞ Here a 1 (λ) is given in (18) (for n = 1). The derivation of Eqn (47) is reported in Appendix C. The set of equations (44), with the additional constraint given by Eqn (47), is now equivalent to Eqns (42).
5 Rapidity distribution functions for the stationary state

Numerical analysis
Eqns (16), (42) can be truncated to obtain a finite system of integral equations, which are defined on the real line λ ∈ (−∞, ∞). One can then numerically solve this finite system either by introducing a cut-off for large λ, or by mapping the equations onto a finite interval. Following the latter approach, we define where q n (λ) is given by Finally, we have defined h being the Lagrange multiplier appearing in (42). The functions χ n (λ) satisfy the following system of equations where a nm (λ) are defined in (17). We then change variables which maps the interval (0, ∞) onto (−1, 1). Since the distributions χ n (λ) are symmetric w.r.t. 0, they can be described by functions with domain (0, ∞). Using the map (52) they become functions χ n (x) with domain (−1, 1). The set of equations (51) becomes where The thermodynamic Bethe-Takahashi equations (16) can be similarly recast in the form where Θ(x) = ρ t n λ(x) , with λ(x) defined in Eq. (52). The infinite systems (53) and (55), defined on the interval (−1, 1), can then be truncated and solved numerically for the functions χ n (x) and Θ n (x), for example using the Gaussian quadrature method. The functions η n (λ) are recovered from (48) and (52), while the particle and hole distributions ρ n (λ), ρ h n (λ) are obtained from the knowledge of η n (λ) and ρ t n (λ). As γ decreases, we found that an increasing number of equations has to be kept when truncating the infinite systems (53), (55) in order to obtain an accurate numerical solution. As we will see in section 6.2, this is due to the fact that, as γ → 0, bound states with higher number of particles are formed and the corresponding distribution functions ρ n (λ), η n (λ) cannot be neglected in (16), (42). As an example, our numerical solution for γ = 0.25, and γ = 2.5 is shown in Fig. 1, where we also provide a comparison with the analytical solution discussed in section 5.3.
Two non-trivial checks for our numerical solution are available. The first is given by Eq. (24), i.e. the solution has to satisfy the sum rule where ε n (λ) is defined in Eqn (20). The second non-trivial check was suggested in Refs [59,60] (see also Ref. [58]), and is based on the observation that the action (29) has to be equal to zero when evaluated on the saddle point solution, i.e. S QA [ρ sp ] = 0, or where S[ρ] and S Y Y [ρ] are defined respectively in (41) and (31). Both (56) and (57) are satisfied by our numerical solutions within a relative numerical error 10 −4 for all numerically accessible values of h. As a final check we have verified that our numerical solution satisfies, within numerical errors, where τ is defined in (50) and γ =c/D is computed from the distribution functions using (22). Relation (58) is equivalent to that found in the repulsive case [51].

Perturbative expansion
Following Ref. [51] we now turn to a "perturbative" analysis of Eqns (42). This will provide us with another non-trivial check on the validity of the analytical solution presented in section 5.3. Defining ϕ n (λ) = 1/η n (λ) and using (50), we can rewrite (42) in the form where W n (λ) is given in (39). We now expand the functions ϕ n (λ) as power series in τ From (59) one readily sees that ϕ n (λ) = O(τ 2n ), i.e.
Using (62) as a starting point we can now solve Eqns (59) by iteration. The calculations are straightforward but tedious, and are sketched in Appendix D. Using this method we have calculated ϕ 1 (λ) up to fifth order in τ . In terms of the the dimensionless variable x = λ/c we have

Exact solution
In this section we discuss how to solve equations (16), (42) analytically. We first observe that the distribution functions ρ n (λ) can be obtained from the set {η n (λ)} n of functions fulfilling Eqns (42) as where τ is given in (50). This relation is analogous to the one found in the repulsive case in Ref. [51]. To prove (64) one takes the partial derivative ∂ τ of both sides of (59). Combining the resulting equation with the thermodynamic version of the Bethe-Takahashi equations (16), and finally invoking the uniqueness of the solution, we obtain (64). This leaves us with the task of solving (42). In what follows we introduce the dimensionless parameter x = λ/c and throughout this section, with a slight abuse of notation, we will use the same notations for functions of λ and of x. Our starting point is the tri-diagonal form (44) of the coupled integral equations (42). Following Ref. [58] we introduce the corresponding Y -system [113,114] where we define y 0 (x) = 0 and Y n (x) = 1 + y n (x) .
One can prove that the set of functions y n (x) with these properties solve the tri-diagonal form of the integral equations equations (44) [58]. To see this, one has to first take the logarithmic derivative of both sides of (65) and take the Fourier transform, integrating in x ∈ (−∞, ∞). Since the argument of the functions in the l.h.s. is shifted by ±i/2 in the imaginary direction, one has to use complex analysis techniques to perform the integral. In particular, under the assumptions (1), (2), (3) the application of the residue theorem precisely generates, after taking the inverse Fourier transform, the driving term d(λ) in (44) [58].
We conjecture that the exact solution for η 1 (x) is given by Our evidence supporting this conjecture is as follows: 1. We have verified using Mathematica that the functions η n (x) generated by substituting (67) into the Y-system (65) have the properties (1), (2) up to n = 30. We have checked for a substantial number of values of the chemical potential h that they have the third property (3) up to n = 10. (67) agrees with the expansion (63) in powers of τ up to fifth order.

Our expression
3. Eqn (67) agrees perfectly with our numerical solution of the saddle-point equations discussed in section 5.1, as is shown in Fig. 1.
The functions ρ n (x) for n ≥ 3 are always written as rational functions but their expressions get lengthier as n increases.
6 Physical properties of the stationary state 6

.1 Local pair correlation function
The distribution functions ρ n (λ), ρ h n (λ) completely characterize the stationary state. Their knowledge, in principle, allows one to calculate all local correlation functions in the thermodynamic limit. In practice, while formulas exist for the expectation values of simple local operators in the Lieb-Liniger model in the finite volume [115][116][117][118], it is generally difficult to take the thermodynamic limit of these expressions. In contrast to the repulsive case [117,[119][120][121][122][123], much less is known in the attractive regime, where technical complications arise that are associated with the existence of string solutions to the Bethe ansatz equations. Here we focus on the computation of the local pair correlation function g 2 = :ρ 2 (0) : We start by applying the Hellmann-Feynman [119,120,122,124] theorem to the expectation value in a general energy eigenstate |{λ j } with energy E[{λ j }] of the finite system In order to evaluate the expression on the r.h.s., we need to take the derivative of the Bethe-Takahashi equations (9) with respect to c Here a nm is given in Eq. (17) and Taking the thermodynamic limit gives Using the thermodynamic version of the Bethe-Takahashi equations (16) and defining we arrive at The set of equations (80) completely fixes the functions b n (λ), once the functions η n (λ) are known. The right hand side of (75) in the finite volume can be cast in the form Taking the thermodynamic limit, and using (79) we finally arrive at Combining (80) and (82) we can express the local pair correlation function as where the functions b n (x) are determined by In (83), (84) we defined η n (x) = η n (xc) , ρ n (x) = ρ n (xc) , a nm (x) = ca nm (xc).
Using the knowledge of the functions η n (λ) for the stationary state, we can solve Eqns (84) numerically and substitute the results into (83) to obtain g 2 (γ). While (83), (84) cannot be solved in closed form, they can be used to obtain an explicit asymptotic expansion around γ = ∞. To that end we use (19), (20) and (24) to rewrite g 2 (γ) as We then use that large values of γ correspond to small values of τ , cf. (58), and carry out a small-τ expansion of the functions where ϕ n (x) = 1/ η n (x) as in section 5.2. Substituting this expansion into the r.h.s. of (84) and proceeding iteratively, we obtain an expansion for the functions b n (x) in powers of τ . The steps are completely analogous to those discussed in section 5.2 for the functions ϕ n (λ) and will not be repeated here. Finally, we use the series expansions of b n (x) and (1 + η n (x)) −1 in (86) to obtain an asymptotic expansion for g 2 (γ). The result is In Fig. 2 we compare results of a full numerical solution of Eqns (83), (84) to the asymptotic expansion (88). As expected, the latter breaks down for sufficiently small values of γ. In contrast to the large-γ regime, the limit γ → 0 is more difficult to analyze because g 2 (γ) is non-analytic in γ = 0. The limit γ → 0 can be calculated as shown in Appendix E, and is given by As was already noted in Ref. [56], (89) implies that the function g 2 (γ) is discontinuous in γ = 0. Indeed, g 2 (0) can be calculated directly by using Wick's theorem in the initial BEC state This discontinuity, which is absent for quenches to the repulsive regime [51], is ascribed to the presence of multi-particle bound states for all values of γ = 0. The former are also at the origin of the non-vanishing limit of g 2 (γ) for γ → ∞ as it will be discussed in the next section.
Finally, an interesting question is the calculation of the three-body one-point correlation function g 3 (γ) on the post-quench steady state. The latter is relevant for experimental realizations of bosons confined in one dimension, as it is proportional to the three-body recombination rate [125]. For g 3 it is reasonable to expect that three-particle bound states may give non-vanishing contributions in the large coupling limit. While g 3 is known for general states in the repulsive Lieb-Liniger model, its computation in the attractive case is significantly harder and requires further development of existing methods. We hope that our work will motivate theoretical efforts in this direction.

Physical implications of the multi-particle bound states
A particularly interesting feature of our stationary state is the presence of finite densities of n-particle bound states with n ≥ 2. In Fig. 3, their densities and energies per volume are shown for a number of different values of γ. We see that the maximum of D n occurs at a value of n that increases as γ decreases. This result has  Figure 3: Density D n and absolute value of the normalized energies per volume e n /γ of the bosons forming n-particle bound states as defined in (21). The plots correspond to (a) γ = 20, (b) γ = 2, (c) γ = 0.2. The total density is fixed D = 1. The energy densities e n are always negative for n ≥ 2 (i.e. |e n | = −e n for n ≥ 2) while e 1 > 0. a simple physical interpretation. In the attractive regime, the bosons have a tendency to form multi-particle bound states. One might naively expect that increasing the strength γ of the attraction between bosons would lead to the formation of bound states with an ever increasing number of particles, thus leading to phase separation. However, in the quench setup the total energy of the system is fixed by the initial state, cf. (24), while the energy of n-particle bound states scales as n 3 , cf. Eqns (20), (21). As a result, n-particle bound states cannot be formed for large values of γ, and indeed they are found to have very small densities for n ≥ 3. On the contrary, decreasing the interaction strength γ, the absolute value of their energy lowers and these bound states become accessible. The dependence of the peak in Fig. 3 on γ is monotonic but non-trivial and it is the result of the competition between the tendency of attractive bosons to cluster, and the fact that n-particle bound states with n very large cannot be formed as a result of energy conservation.
The presence of multi-particle bound states affects measurable properties of the system, and is the reason for the particular behaviour of the local pair correlation function computed in the previous section. Remarkably, this is true also in the limit γ → ∞. This is in marked contrast to the super Tonks-Girardeau gas, where bound states are absent. To exhibit the important role of bound states in the limit of large γ, we will demonstrate that the limiting value of g 2 (γ) for γ → ∞ is entirely determined by bound pairs. It follows from (83) that g 2 (γ) can be written in the form where g (m) 2 (γ) denotes the contribution of m-particle bound states to the local pair correlation Let us first show that unbound particles give a vanishing contribution In order to prove this, we use that at leading order in 1/γ we have b 1 (x) = x. Using the explicit expressions for η 1 (x), ρ 1 (x) we can then perform the integrations in the r.h.s. of Eq. (92) exactly and take the limit γ → ∞ afterwards. We obtain which establishes (93). Next, we address the bound pair contribution. At leading order in 1/γ we have b 2 (x) = 2x, and using the explicit expression for η 2 (x) we obtain This leaves us with the contribution Although the function ρ 2 (x) is known, cf. Eq. (70), its expression is unwieldy and it is difficult to compute the integral analytically. On the other hand, one cannot expand ρ 2 (x) in 1/γ inside the integral, because the integral of individual terms in this expansion are not convergent (signalling that in this case one cannot exchange the order of the limit γ → ∞ and of the integration). Nevertheless, the numerical computation of the integral in (97) for large values of γ presents no difficulties and one can then compute the limit numerically. We found that the limit in Eq. (97) is equal to 4 within machine precision so that lim γ→∞ g 2 (γ) = 4 = lim γ→∞ g 2 (γ).
Finally, we verified that contributions coming from bound states with higher numbers of particles are vanishing, i.e. g (m)

Conclusions
We have considered quantum quenches from an ideal Bose condensate to the one-dimensional Lieb-Liniger model with arbitrary attractive interactions. We have determined the stationary state, and determined its physical properties. In particular, we revealed that the stationary state is composed of an interesting mixture of multi-particle bound states, and computed the local pair correlation function in this state. Our discussion presents a detailed derivation of results first announced in Ref. [56].
As we have stressed repeatedly, the most intriguing feature of the stationary state for the quench studied in this work is the presence of multi-particle bound states. As was argued in Ref. [56], their properties could in principle be probed in ultra-cold atoms experiments. Multi-particle bound states are also formed in the quench from the Néel state to the gapped XXZ model, as it was recently reported in Refs. [57][58][59][60]. However, in contrast to our case, the bound state densities are always small compared to the density of unbound magnons for all the values of the final anisotropy parameter ∆ ≥ 1 [58].
Our work also provides an interesting physical example of a quantum quench, where different initial conditions lead to stationary states with qualitatively different features. Indeed, a quench in the one-dimensional Bose gas from the infinitely repulsive to the infinitely attractive regime leads to the super Tonks-Girardeau gas, where bound states are absent. On the other hand, as shown in section 6, if the initial state is an ideal Bose condensate, bound states have important consequences on the correlation functions of the system even in the limit of large negative interactions.
An interesting open question is to find a description of our stationary state in terms of a GGE. As the stationary state involves bound states, it is likely that the GGE will involve not yet known quasi-local conserved charges [41,44,45] as well as the known ultra-local ones [126]. In the Lieb-Liniger model technical difficulties arise when addressing such issues, as expectation values of local conserved charges generally exhibit divergences [29,89,126]. In addition, very little is known about quasi-local conserved charges for interacting models defined in the continuum [41,48].
Finally, it would be interesting to investigate the approach to the steady state in the quench considered in this work. While this is in general a very difficult problem, in the repulsive regime the post-quench time evolution from the non-interacting BEC state was considered in [53]. There an efficient numerical evaluation of the representation (26) was performed, based on the knowledge of exact one-point form factors [118]. The attractive regime, however, is significantly more involved due to the presence of bound states and the study of the whole post-quench time evolution remains a theoretical challenge for future investigations.

A Overlaps in the presence of zero-momentum n-strings
In this appendix we argue that Eq. (38) gives the leading term in the thermodynamic limit of the logarithm of the overlap between the BEC state and a parity-invariant Bethe state, even in cases where the latter contains zero-momentum strings.
To see this, consider a parity invariant Bethe state with a single zero-momentum m-string, and K parityrelated pairs of n j -strings. The total number of particles in such a state is then N = 2 j n j + m. In Ref. [111] an explicit expression for the overlap (33) of such states with a BEC state in the zero-density limit (L → ∞ and N fixed) was obtained. Up to an irrelevant (for our purposes) overall minus sign, it reads where λ p is the centre of the p'th string. We see that as a result of having a zero-momentum string, an additional pre-factor L appears. In general, the presence of M zero-momentum strings will lead to an additional prefactor L M [111]. While (99) is derived in the zero density limit, we expect an additional pre-factor to be present also if one considers the thermodynamic limit N, L → ∞, at finite density D = N/L. Importantly such pre-factors will result in sub-leading corrections of order (ln L)/L to the logarithm of the overlaps. This suggests that (38) holds even for states with zero-momentum n-strings.
B Tri-diagonal form of the coupled integral equations

B.1 Tri-diagonal Bethe-Takahashi equations
Our starting point are the thermodynamic Bethe equations (16). For later convenience we introduce the following notations for the Fourier transform of a function We recall that f * g denotes the convolution of two functions, cf. (43). The Fourier transform of a n (λ) defined in (18) is easily computedâ Following Ref. [85], we introduce the symbols We can then perform the Fourier transform of both sides of (16) and obtain where ρ t n (λ) are given in (15). We now definê ρ 0 (k) = 0 .
Iterating the above procedure n times we arrive at − 2h + a n * (f 0 + f 1 ) = a n * ln η n+1 − a n+1 * ln(1 + η n ) − a n+2 * ln(1 + η −1 n+1 ) Fourier transforming and using the definition for f r given in (119) Assuming that η −1 n (λ) is vanishing sufficiently fast as n → ∞ for a generic (and fixed) value of λ, we can drop the infinite sum and the two previous terms, and arrive at Eq. (47).

D Perturbative analysis
In this appendix we sketch the calculations leading to the expansion (63). Throughout this appendix we work with the dimensionless variable x = λ/c. At the lowest order, it follows from Eq. (62) that Since ϕ n (x) ∝ τ 2n , we can neglect ϕ n (x) with n ≥ 2 to compute the third order expansion of ϕ 1 (x). Hence, the infinite sum in (59) for n = 1 can be truncated, at third order in τ , to the first term (m = 1), where we can use the expansion (135) for ϕ 1 (λ). Following Ref. [51] one can then use identity (118) to perform the convolution integral and finally obtain ϕ 1 (x) = τ 2 x 2 x 2 + 1 4 1 − 4τ One can then perform the same steps for higher order corrections, at each stage of the calculation keeping all the relevant terms. For example, already at the fourth order in τ of ϕ 1 (x) one cannot neglect the lowest order contribution coming from ϕ 2 (x) in the r.h.s. of Eq. (59). For higher orders one also has to consider corrections to ϕ n (x) with n ≥ 2.
Our starting point is Eqn (86). Rescaling variables bŷ we have .
The functionsb n (x) satisfy the coupled nonlinear integral equationŝ whereâ Our goal is to determine the limit .
The calculation is non-trivial as we cannot exchange the infinite sum with the limit. However, based on numerical evidence we claim that this limit is finite, and (137) then immediately follows from (139). Note that the numerical computation of g 2 (γ) is increasingly demanding as γ → 0, due to the fact that more and more strings contribute. Accordingly, the infinite systems (83) and (84) have to be truncated to a larger number of equations and the numerical computation takes more time to provide precise results. We were able to numerically compute g 2 (γ) for decreasing values of γ down to γ = 0.025 where g 2 (0.025) 2.11 and approximately 30 strings contributed to the computation. We fitted the numerical data for small γ with G(γ) = α 1 + α 2 √ γ and we correctly found α 1 = 2 within the numerical error.