Stationary State After a Quench to the Lieb–Liniger from Rotating BECs

We study long-time dynamics of a bosonic system after suddenly switching on repulsive delta-like interactions. As initial states, we consider two experimentally relevant configurations: a rotating BEC and two counter-propagating BECs with opposite momentum, both on a ring. In the first case, the rapidity distribution function for the stationary state is derived analytically and it is given by the distribution obtained for the same quench starting from a BEC, shifted by the momentum of each boson. In the second case, the rapidity distribution function is obtained numerically for generic values of repulsive interaction and initial momentum. The significant differences for the case of large versus small quenches are discussed.


Introduction
In the last ten years, out-of-equilibrium many-body quantum physics has become a very active area of research, especially fuelled by the experimental realizations developed in cold-atomic setups [1][2][3][4]. As opposed to usual solid state devices, with ultracold atoms it is possible to experimentally realize quantum systems in which the effect of decoherence and dissipation is negligible on long time scales [5][6][7][8][9][10]. Essentially unitary time evolution of an initially isolated system can be studied, enough for allowing the observation of the stationary configuration eventually reached [11][12][13], if any. This has led to the possibility of addressing long standing theoretical questions concerning equilibration and thermalisation in many body quantum systems, the role of integrability and universality and the study of the time-dependence of correlation functions [14][15][16][17][18] One of the most studied out-of-equilibrium protocol is the sudden quantum quench [19][20][21][22][23][24][25][26]. After having prepared a system in an initial configuration, eigenstate of a pre-quenched Hamiltonian, from a certain instant of time on it is let evolved isolated from the rest according to a different Hamiltonian (the post-quenched one). In this kind of settings, the main focus is the investigation of the properties of the stationary state reached long after the quench. To make the time evolution non trivial, the initial state should not be chosen as an eigenstate of the post-quenched Hamiltonian. When considering a thermodynamically large system, this initial state has a non vanishing overlap with a very large number of eigenstates of the post-quenched Hamiltonian. If the thermodynamic limit is taken before the long time limit, a dephazing mechanism takes place: as a consequence, any observable relative to a finite subsystem will relax to time-independent values, even if the entire system is isolated [13].
While for non integrable systems the stationary distribution is generally given by the thermal Gibbs one [27,28], it has been conjectured that integrable systems relax to a Generalized Gibbs Ensemble (GGE) [29], obtained imposing the conservation of all local charges throughout the unitary evolution. For free, or mappable-to-free systems, the conserved charges can be identified with the momentum occupation numbers and thus the GGE can be explicitly constructed [30][31][32][33][34][35][36][37][38], [39,40]. In the case of a truly interacting (yet integrable) model, i.e. for those having a dressed two-body scattering matrix, the construction of the GGE is more involved, due to the difficulty in writing the conserved charges in a closed operatorial form [41][42][43][44][45]. As pointed out in [46], the GGE is hard to construct for the Lieb Liniger model because the conserved charges are divergent. Moreover, it has been recently pointed out that a GGE which takes into account only the local conserved charges is not able to capture the stationary state in some cases [47][48][49][50] and one needs to include also the quasi local charges [51][52][53].
It has been shown that there is another way of computing expectation values of local quantities at any time after a quench in integrable models, the so-called quench action method [54]. In particular, it is proved that the long time behaviour of local observables is given by the expectation value with respect to a single representative state of the (post-quenched) Hamiltonian. This representative state can be obtained as the result of a variational method, exact in the thermodynamic limit, which requires the knowledge of the overlap of any eigenstate of the post-quenched Hamiltonian with the initial state. As a consequence, this method is potentially very useful for interacting integrable models [49,[55][56][57][58][59] because it does not require having the expression of the conserved charges, differently from the GGE approach. Nonetheless, the expression for the overlaps is often highly non trivial [60,61] and always initial state dependent.
In the present paper, by making use of the quench action method, we solve two quench problems that have recently attracted the attention of both experimentalists and theoreticians [5,[62][63][64][65][66][67][68][69][70][71]. We obtain results for the steady-state reached by an initially free bosonic system after switching on at t = 0 repulsive contact interaction (Lieb-Liniger model [72][73][74][75][76]). As initial states, we consider two experimentally relevant configurations, that are both eigenstates of the Lieb-Liniger model at zero interaction strength. The first is a rotating BEC and the second consists of many couples of counter propagating BECs with opposite momentum, moving on a ring. In this latter case, our calculations represent a close theoretical description of the cornerstone Quantum Newton's Cradle [5]. The only differences are that we do not take into account the presence of the trap [77] and that we actually perform a quench of the interaction strength c, from c = 0 at t = 0 to c > 0 at t > 0, while in the actual experiment the two counter-propagating wavepackets are prepared and let evolved at a fixed value of c.
We note that these also represent physical examples of a quench where the initial state is an excited state of the pre-quenched Hamiltonian [39,40], depending on the value of the momentum. This is different from the majority of quenches studied in literature, where the initial state is the ground state of the pre-quenched Hamiltonian.
In the case of the initially rotating BEC, we have found that the global shape of the rapidity root distribution in the asymptotic state is not affected by the rotation of the BEC and it is just shifted uniformly with respect to the static BEC case by the momentum of each boson. For the initially counter propagating BECs, instead, the entire structure of the root distribution is modified and it is very different in the large and small quench limit.
The paper is organized as follows: in Sects. 2 and 3 we review the main features of the Lieb-Liniger model and the quench action method while in the remaining sections we present our original results. In Sect. 4 we present the results for the rapidity distribution of the asymptotic state after a quench from the initial rotating BEC. In Sect. 5 we obtain the stationary root distribution after a quench from the colliding BECs and we discuss the approximate solutions in the case of a large and small quench. In Sect. 6 we draw conclusions.

The Model
In this section we review the basics of the Lieb-Liniger model [78,79] and set our conventions, following [55]. We consider a single component bosonic quantum gas of N particles in a 1D box of length L with periodic boundary conditions. From t = 0 the system evolves unitarily according to the Lieb-Liniger Hamiltonian [78] (setting = 2m = 1) The parameter c describes the strength of interaction and we will consider c > 0 in the rest of the paper, so as to depict repulsion between particles. This is a realistic model for interacting bosons in quantum degenerate gases with s-wave scattering potential. The model can be solved by means of the Bethe Ansatz [80]. The solution of (1) can be found taking into account the symmetry of the bosonic wavefunction under two particles' exchange and the fact that the N particles are free unless two of them occupy the same position. Accordingly, we can consider all the possible N ! domains defined by the ordered non-coinciding particles' positions (P) : 0 ≤ x P 1 < x P 2 < · · · < x P N ≤ L, where P is the permutation of the number set {1, 2, . . . , N }. The total wavefunction can be written as the sum of the wavefunctions in each domain where ψ(x) is the bosonic wavefunction written in first quantization. Due to the Bose statistics the wavefunction is the same in all domains: For the model (1), ψ 1 can be written as a superposition of N ! plane waves with , (5) The norm of the wavefunction is given by the Gaudin determinant [81].
Imposing periodic boundary conditions, the rapidities get quantized and have to satisfy the Bethe equations [80], a set of N coupled algebraic equations where I ≡ {I 1 , I 2 , . . . , I N } are the quantum numbers of the rapidities, which label an eigenstate uniquely and are c-independent, differently from λ. The full set of allowed quantum numbers is the union of the occupied ones I and the unoccupied ones, defined as I h . If N is odd, I are the integers and I h the semi-integers; if N is even the vice-versa. The solution of the Bethe equations (6) provides an exhaustive knowledge of the spectrum of the Lieb-Liniger model; for a given set of rapidities λ, the total momenta and energy of the system are obtained as Higher conserved charges Q n , n ∈ N, are of the form [82,83] Q n = N j=1 λ n j .
As already mentioned in Sect. 1, in our quench protocol we are interested in the thermodynamic limit of this model (N → ∞, L → ∞, with N /L fixed). It is useful to introduce the (particle) root density ρ(λ), i.e. the density distribution function of the rapidities defined by the particle numbers I in the interval (λ, λ + λ) Accordingly we can define a hole root density ρ h (λ), i.e. a density distribution function of the rapidities defined by I h , which are the unoccupied quantum numbers of the allowed set, in the interval (λ, λ + λ). Hence the Bethe equations (6) can be rewritten in the thermodynamic limit as where The momentum and energy of the system per unit length (7) thus become

Quench Action Method: Overview
In this section, following [54], we give an overview of the quench action method that we will use to obtain the results for the steady states. This method provides a way of computing the expectation value of a local observable at any time during the unitary evolution following a quench, in a thermodynamically large system. We will show that the long-time behaviour is given by the expectation value on a single eigenstate of the post-quenched Hamiltonian. This representative state can be simply computed extremizing a particular functional of the root density ρ(λ), called quench action. The aim is to compute the expectation value of a generic local observable O at time t with |ψ(t) being the time evolved state with the post-quenched Hamiltonian H from the initial state |ψ(0) |ψ(t) = e −i Ht |ψ(0) .
Let us first focus on the numerator of (13); expanding it on the post-quenched eigenbasis |I we get where S I ≡ − ln I |ψ(0) and H |I = I |I . The practical difficulty in evaluating such a double sum can be overcome by going to the thermodynamic limit. Any sum over a discrete set of quantum numbers can be transformed into a functional integral over all the microscopic configurations sharing the same root density ρ(λ), each of which weighted with the Yang- The where the hole density ρ h (λ) is related to the particle density ρ(λ) by the Bethe equations (10). To extract the leading behaviour of (15) in the thermodynamic limit, it is sufficient to transform only one of the two sums in its continuum form with eq. (16). In fact, since O is local, I |O|I is non vanishing only provided that the states |I , |I have the same (macroscopic) root density in the thermodynamic limit. After enforcing (16), (15) becomes where the integration and functional dependence on ρ h is kept implicit from now on. The states |I and |ρ can only differ for a small number of microscopic particle-hole excitations. This yields I I | = e ρ + e|, I ρ + δω e and S * I + S ρ 2 Re(S ρ ) + δs * e . Notice that δω e and δs e are the microscopic differences in energy and overlap between |ρ + e and |ρ , while S ρ and ρ are respectively the extensive part of the overlap and energy, in the thermodynamic limit. Eq. (18) can thus be rewritten as This integral can be computed with a saddle point evaluation in the thermodynamic limit, since S Y Y [ρ] and Re(S ρ ) are extensive in the system's size. The most significant contribution is given by the saddle point distribution ρ s that makes the exponent in (19) stationary. Defining the quench action as Let us now focus on the denominator of (13); its leading contribution is given by the same root density ρ s of Eq. (21) Putting everything together, the expectation value (13) becomes Eq. (23) means that in the thermodynamic limit O(t) is fully determined by states whose root density distribution differs only by few microscopic particle-hole excitations from the saddle point one.
When t → ∞ a stationary phase approximation yields the additional condition so that (23) becomes lim The conclusion is that only one state matters for the stationary behaviour: it is the one fixed by the stationarity of the quench action (20).

Quench from a Rotating BEC to the Lieb-Liniger
In this section we present the first of our original contributions: the solution of the stationary state reached after a quench from H L L (c = 0) to H L L (c > 0) starting from a rotating BEC.
We represent the rotating BEC as a system of N bosons on a ring of circumference L, each of which with a given value of the momentum k, depicted in Fig. 1.
The state is given by . At t = 0 we suddenly turn on the interaction to an arbitrary positive value and then let the system evolve unitarily until a steady state is reached. The first step for obtaining the saddle point root distribution in the stationary state is to construct the quench action (20). It entails the computation of the overlap between a generic Bethe state, eigenfunction of the post-quenched Lieb-Liniger Hamiltonian, and the initial state (see Appendix 1 for the detailed computation). Due to momentum conservation, a non-vanishing overlap arises only for the Bethe states with rapidities that are symmetrically distributed around k. The saddle point evaluation (19), on which the quench action is based, implies that only the extensive part of the overlap S ρ is relevant for identifying ρ s in the thermodynamic limit. As explicitly computed in Appendix 1, it is given by Comparing it to the overlap between a Bethe state and a non-rotating BEC [55], i.e. (26) with k = 0, S ρ has the same expression, with the only difference that the root distribution ρ(λ) is shifted of k.

Let us now impose δS Q A [ρ]
δρ ρ s = 0. Given (17), (20) and (27) the ρ-dependent part of the quench action per unit length reads We have to impose the normalization condition +∞ −∞ dλ ρ(λ) = n; this can be done modifying our functional measure by adding a Lagrange multiplier μ [81] in the following way Imposing that the variation of the quench action under an infinitesimal transformation of the Bethe roots be vanishing we get This is an equation for the root density ρ ≡ ρ(λ, c, k, μ). In particular, * denotes the convolution product defined as and the Bethe equations have already been imposed at the level of the relation between δρ h (λ) and δρ(λ). Switching to dimensionless variables x ≡ λ c and defining τ ≡ e μ/2 and a(x, τ, α) as the ratio of the particle to hole density the saddle point condition (29) can be rewritten as where K(x) ≡ 1 1+x 2 . The connection between a(x, τ, α) and the root density is obtained by applying on the previous equation the operator and comparing it to the Bethe equation written in dimensionless variables Eventually we get that the root density expressed in terms of the rescaled variables, ρ(x, τ, α) has the formρ In conclusion, finding the saddle point distribution ρ s amounts to solving the non-linear integral equation (31) for the function a(x, τ, α) and then plugging it into (34). A perturbative study of the solution ρ s can be used as a guideline towards its full analytic solution. Let us focus on Eq. (31) in the τ → 0 limit. As τ → 0, for any fixed x > τ, ln is a negative, divergent quantity, so the convolution integral gives a subdominant contribution. Hence the zero-th order for a is As it could be argued from the analogy between (27) and the expression of the overlap for the initial BEC state, (35) coincides with the zero-th order term of the function a for the initial BEC state [55], upon adding a constant shift of α to the dimensionless variables. We have checked that this argument holds also for higher order terms in τ , hence the full analytic solution for a(x) is Eq. (36) is the main result of this section. After applying (34), the saddle point distribution is derived, yieldingρ Being its exact expression mathematically involved, we prefer showing the plots in Fig. 2 for the saddle point distribution as a function of x − α, for several values of τ . As we can see from (36), α only acts as a reshift of the centre of the function a; the same happens for the distributionρ s . The parameter τ can be related to the dimensionless interaction constant of the Lieb-Liniger, namely γ [85][86][87][88]. Defined as the ratio of the interaction strength c and the particle density n = N /L, γ is the key-parameter of the Lieb-Liniger model, in terms of which all the physical properties of the system can be described (at zero temperature). Thanks to (9), it is related to the root distribution by Integrating numerically (38) for several fixed values of τ , we find that for the specific case of the rotating BEC γ B EC is independent of α and is related to τ as follows We would like to highlight that, in this quench protocol from c = 0 to c finite, γ is also directly a measure of the quench amplitude. We have a small quench for γ 1 [89][90][91], an intermediate quench for γ 1 and a large quench for γ 1 [86,[92][93][94]. In conclusion we have found that the shape of the root distribution in the stationary state does not change if the initial state, evolved with the Lieb-Liniger, is a non-rotating or rotating BEC. The rotating BEC can be thought of as a state with a steady current, given by the bosons all moving in the same direction with the same value of the momentum. The effect of this

Quench from a State with Oppositely Moving BECs to the Lieb-Liniger
In this section we will present the solution of the saddle point distribution for the same quench of the interaction strength from H L L (c = 0) to H L L (c > 0) on an initial state given by two colliding BECs.
The initial configuration is a product state in which each of the N bosons is in a superposition of right and left moving plane wave states with a definite value of the momentum (see Fig. 3), ie.
with the k-dependent normalization factor f (k) = k L+sin(k L) k L ensuring that x|0 n is normalized to 1. Using periodic boundary conditions on the initial state, we have that k = 2πn/L, hence the normalization constant f (k) reduces to 1.
To obtain the saddle point root distribution, we follow the same path outlined in the previous section. It can be shown (see Appendix 2) that, in the thermodynamic limit, the leading contribution to the overlap between this initial state and a Bethe state is given by parity invariant Bethe states. The extensive part of the overlap is (see Appendix 2 for the detailed computation) We note that, irrespectively of the boundary conditions used, the normalization factor f (k) of the initial state can be dropped in S ρ because it does not contribute to the saddle point distribution, being ρ-independent. Using (17), (20) and (41) to construct the quench action, extremizing it with respect to ρ(λ) and switching to dimensionless variables x = λ c yields the following non-linear integral equation ln(a(x, τ, α)) = ln a(y, τ, α)), (42) with a(x, τ, α) and τ defined as in Sect. 4. Eq. (42) differs from (32), obtained for the initial rotating BEC, only in the contribution of the overlap. Once it is solved for a(x, τ, α), the root distribution can be derived using the definition of a(x, τ, α) and the dimensionless Bethe root equation (33) written in terms of ρ(x, τ, α) ≡ ρ(λ = cx, c, k, τ ). Since the analytic solution is highly non trivial for generic values of τ and α, we will discuss the solution obtained from the numerical integration of Eq. (42), for various values of α and τ .
In Fig. 4 we show the numerical plots for the saddle point distribution for α = 1 as a function of x, obtained for τ = 0.1, 1, 100. Computing γ from (38) for each of the saddle point distributions, we see that these three plots correspond to the three possible quench regimes: Fig. 4(a) represents a large quench (γ 1), Fig. 4(b) an intermediate quench (γ 1) and Fig. 4(c) a small quench (γ 1). Differently from the rotating BEC, in this case γ depends on both α and τ but there is still the correspondence between small(large) values of τ and large(small) values of γ . As an example, we show in Fig. 5 the plot of γ as a function of τ in the case α = 1.
Computing γ for several values of α, τ and imposing it to be of order one, we see that the separation between small/large quenches sits approximately at τ max{α, 1}. The classification in terms of the parameter γ turns out to be useful, in fact the saddle point distributions obtained for different values of α and τ in a specific γ regime show the same Let us now focus on the analysis of the two opposite regime of large and small quenches, where approximate solutions were found, and comment Fig. 4. For the large quench case, we can use a perturbative expansion in the small parameter τ since the large quench regime γ 1 corresponds to τ 1; up to order τ 3 , the solution is It is represented by the continuous line in Fig. 4(a), which agrees extremely well with the numerical plots and is derived in Appendix 3. The behaviour for large x, for τ small but finite, is given bỹ whose leading order is 1 x 4 , the same as for the case of the initial BEC state. To appreciate the difference between the root distribution of the state with oppositely moving BECs and the rotating BEC in the large quench limit we can compareρ B EC s (x − α, τ ) withρ (1) s (x, τ, α) (see Fig. 6). The difference between the two plots increases with increasing α; it is clear that the saddle point distribution for the oppositely moving BECs can not be obtained trivially from reshifting the root density distribution of the BEC case.
Let us now focus on the small quench limit. In this regime, τ > max{α, 1}, so α is negligible with respect to τ ; this explains why the root distribution in Fig. 4(c), for a given value of τ , is the same for any values of α ranging over three order of magnitudes, from 0.01 to 10. In Fig. 4(c), in addition to the root distribution in the small quench limit, we have also plotted the Wigner-semicircular law (represented by the continuous line) Indeed in [55], for the initial BEC state, it was rigorously proved that the Wigner-semicircular law describes the stationary root distribution in the case of a small quench (more precisely in the limit c → 0 with n fixed). Interestingly enough, (45) also represents the leading term of the ground state root distribution of the Lieb-Liniger model in the low-γ expansion [95], i.e. of the ground state of the post-quenched Hamiltonian. Actually, as can be seen in Fig. 4(c), we note that also starting from two counter propagating BECs the stationary root distribution turns out to be very close to the Wigner law, apart from a tail for |x| > 1 π √ γ , probably due to initial excitations.

Conclusions
In this paper we have studied in detail the asymptotic state reached after a quench to arbitrary values of the interaction strength of the Lieb-Liniger model, from two experimentally relevant initial configurations. In the first case, when considering bosons forming a rotating BEC state, the steady state root distribution was found to be the same as when the bosons are initially in a BEC configuration, simply shifted by the momentum of each boson. In the second case, when the initial state can be represented as two oppositely moving BECs, we showed that the non vanishing momentum of the bosons radically modifies the stationary root distribution. We were able to identify three regimes for the saddle point root distribution, characterized by the parameter γ (large, intermediate and small quench regime) and to find approximate solutions in both the small and large quench limit that fit very well with the numerical data.
It would be very interesting to analyze the dependence on the quench amplitude γ of the local two and three point function of the stationary state. Such analysis could allow us to compare with the well known results for the correlation functions in highly excited thermal states of the Lieb-Liniger [87,88,96].
Acknowledgments Open access funding provided by Max Planck Society (Institute for the Physics of Complex Systems). I acknowledge support and hospitality from SISSA and the Max Planck Institute for the Physics of Complex Systems, Dresden, where this work was completed. I would like to thank Pasquale Calabrese for useful comments and for reading the final version of the manuscript.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Overlap Between a Bethe Eigenstate and the Rotating BEC
In this appendix we give the detailed derivation of formula (27), which is the extensive part of the overlap between the rotating BEC and a generic eigenstate of H L L (c > 0), in the thermodynamic limit.

Overlap for Generic N, L
Due to momentum conservation, the overlap between a Bethe state and the rotating BEC is non-vanishing if the sum of the rapidities of the Bethe state is equal to the total momentum of the initial state, N k. This class of states includes as a special case the Bethe states with rapidities symmetrically distributed around k. Nevertheless we can safely consider only this subset of states because the Bethe equations are not mutually consistent for Bethe states that are not parity invariant (with respect to k).
To obtain the overlap, let us first explicitly write the generic normalized (subscript n ) Bethe eigenstate. It is given by (4) divided by its norm, the Gaudin determinant [81] where G is a N × N matrix with elements A Y Y is the Yang-Yang action for the Lieb-Liniger model [84], defined as with (λ) ≡ λ 0 dλ 2arctan( λ c ). Since the Hessian of the Yang-Yang action depends only on the difference between rapidities, which are symmetric with respect to k, the Gaudin determinant is independent of k and has the same expression as for the initial (non-rotating) BEC state, i.e.
The overlap can thus be written as where we have defined and we recall that 1 = {0 ≤ x 1 < x 2 < · · · < x N ≤ L}. To compute I → (L) let us explicitly write the expression for the not-normalized eigenfunction hence we get where λ l ≡ λ l − k. This is exactly of the same form as the one computed for the initial BEC case [55], provided we substitute the rapidities λ with the distance of the rapidities from the centre k of the distribution λ . By exploiting the properties of the Laplace transform of the overlap [97,98], (53) can thus be rewritten generalizing to the present case the expression for the BEC of [55] where Res s indicates to take the residue with respect to the s-variable.

Extensive Part of the Overlap in the Thermodynamic Limit
Deriving the thermodynamic limit of the overlap is a non trivial task (see for example [99]) but, for the computation of the quench action (20), we just need the extensive part of it. In [55] the full expression of the overlap in the thermodynamic limit was analytically computed for the initial BEC state and its extensive part turned out to coincide with the leading contribution of its zero-density limit. In the present case, we first compute the zero density limit of the overlap and then we check that it coincides with the leading order in L of the overlap computed in the finite density case N = 2, L finite. Based on the analogy with the exact calculation in the BEC case, this gives us the sought expression for the extensive part of the overlap in the thermodynamic limit.
To compute the zero-density limit of (54), we have to isolate the leading term in L when the number of particles N is kept finite. In this case the computation is simplified with respect to the thermodynamic limit because the rapidities are not quantized, hence we do not have to enforce any Bethe equations. Since, for any N , L, the overlap with the rotating BEC has the same expression (in terms of λ ) as the overlap of the BEC, this holds also in the zero-density limit; thus we can straightforwardly consider the expression obtained in [55] and adapt it to our rotating case Since the Gaudin determinant, to the leading order in L, is of the form the zero-density limit of (50) becomes The expansion in 1/L can be rewritten as an expansion in n, n being the density of particles n = N /L, which reads (58) The only difference from the non rotating case is that here ρ is shifted of k.
Let us compute the overlap in a finite density case, i.e. N = 2 with L finite. Working out calculations from (54) with the Bethe state |λ 1 , λ 2 = |k + β, k − β we get Enforcing the Bethe equation the overlap reads The k-dependent term in (61) is subleading in L. Given the non trivial form of the overlap already for the finite size case with N = 2, we can extrapolate that any k-dependent term will give subleading contributions for any finite size case. The extensive part of the overlap in the thermodynamic limit is thus given by (58)

Appendix 2: Overlap Between a Bethe Eigenstate and a State with Oppositely Moving BECs
In this appendix we derive formula (41). We first compute the expression for the overlap between a state with oppositely moving BECs and a generic Bethe state for generic N , L.
Then we derive its zero-density limit which we identify with the extensive part of the overlap in the thermodynamic limit similarly to what done in Appendix 1.

Overlap for Generic N, L
The overlap is 1 Defining the expression for the overlap is Let us compute I (L). By inserting the expression for the Bethe wavefunction (52) in (65) we get which can be written in a more compact form as In the previous equation we have defined λ P = {λ P 1 , λ P 2 , . . . , λ P N } as the set of rapidities associated to the particles at positions x = {x 1 , x 2 , . . . , x N } within the permutation P and k = {k 1 , k 2 , . . . , k N } as the set of momenta associated to the particles at positions x, with k i = ±k, i = 1 . . . N ; lastly we have defined Differently from the rotating case, a state with two oppositely moving BECs is not an eigenstate of the momentum operator hence, in principle, it could have non vanishing overlap with non parity symmetric Bethe states. Nevertheless it is possible to show that in the thermodynamic limit the leading contribution to the overlap is obtained just by considering the parity symmetric Bethe states. First of all, let us stress again that, since the extensive part of the overlap in the thermodynamic limit coincides with the extensive part of the overlap in the zero density limit, it is enough to prove the previous statement in the zero density limit. We have explicitly computed the overlap for a generic N = 2 Bethe state, which is an already complicated expression. Considering opposite rapidities, which is only possible realization of parity symmetric state for N = 2, we get . (68) In the L → ∞ limit, considering that λ 1 + λ 2 = 0 and enforcing the Bethe equations, ( 68) is extensive in the system's size We explicitly checked that the non parity invariant states give a subleading contribution to the overlap, which can be neglected in the thermodynamic limit. Hence, we can safely compute (67) only for parity invariant Bethe states.
Let us now go back to (67). If we consider a single term in the summation over k, (67) represents the overlap between a Bethe state with rapidities λ P and a state where the bosons have momenta k. For N = 2, given a Bethe state {β, −β}, the only relevant combinations are {β + k, −β − k}, {β − k, −β + k}. We will label Q as an element of the permutation group acting on a parity symmetric k state and we will denote λ Q P ≡ λ P + k Q . Exploiting the properties of the Laplace transform [55,97,98], (67) can be rewritten as where Res s indicates to take the residue with respect to the s-variable.

Extensive Part of the Overlap in the Thermodynamic Limit
We will follow the same path outlined in Appendix 1, which consists in deriving first the zero density limit of the overlap and then computing the leading contribution in L in the finite density case. To compute the zero density limit of the overlap, we have to isolate the leading term in L keeping N finite. Let us first consider how this is done for the k = 0 case and then discuss the general k = 0 case.
Setting k = 0 in (70), the highest order in L is carried by the pole in s = 0; it is of order Each of the possible target configurations can be identified with a set of N 2 numbers that are collected in a vector σ whose j-th element indicates if the j-th pair (λ j , −λ j ) in the target configuration is reversed (−λ j , λ j ) (σ j = 1) or not (σ j = 0). Each σ corresponds to N 2 ! configurations of rapidities; instead of summing on P we can sum on σ and insert a factor of N 2 !. Let us now consider the case with k = 0 given by (70). In order for the pole in s = 0 to be of maximal order, i.e. N 2 + 1, k should not only be parity invariant, but also parity invariant within each pair (λ j , −λ j ). Each of these configurations is thus indicated by a vector with N 2 elements, ζ , whose j-th element indicates if associated to the j-th pair of λ (λ j , −λ j ) one should sum (−k, k) (ζ j = 1) or (k, −k) (ζ j = 0). Then we get Taking into account that and putting everything together we find that (70), in the zero density limit, is given by where we recall that α = k c . According to (65) As L → ∞, it reduces to which is (74) for N = 2.
For N = 4, with a state of the form |λ 1 , −λ 1 , λ 2 , −λ 2 = |β, −β, γ , −γ , the leading in L term of the overlap turns out to be of the form which is (74) for N = 4. Since the leading term in L of the overlap, in the finite density cases we have examined, coincides with its zero density limit, we can safely state that the extensive part of the overlap in the thermodynamic limit is given by (74).

Appendix 3: Saddle Point Distribution for a State with Oppositely Moving BECs in the Large Quench Limit
In this section we derive expression (43) for the saddle point solutionρ (1) s (x) in the case of oppositely moving BECs in the large quench limit.