The effect of atom losses on the distribution of rapidities in the one-dimensional Bose gas

We theoretically investigate the effects of atom losses in the one-dimensional (1D) Bose gas with repulsive contact interactions, a famous quantum integrable system also known as the Lieb-Liniger gas. The generic case of K-body losses (K = 1,2,3,...) is considered. We assume that the loss rate is much smaller than the rate of intrinsic relaxation of the system, so that at any time the state of the system is captured by its rapidity distribution (or, equivalently, by a Generalized Gibbs Ensemble). We give the equation governing the time evolution of the rapidity distribution and we propose a general numerical procedure to solve it. In the asymptotic regimes of vanishing repulsion -- where the gas behaves like an ideal Bose gas -- and hard-core repulsion -- where the gas is mapped to a non-interacting Fermi gas -- we derive analytic formulas. In the latter case, our analytic result shows that losses affect the rapidity distribution in a non-trivial way, the time derivative of the rapidity distribution being both non-linear and non-local in rapidity space.

We theoretically investigate the effects of atom losses in the one-dimensional (1D) Bose gas with repulsive contact interactions, a famous quantum integrable system also known as the Lieb-Liniger gas. The generic case of K-body losses (K = 1, 2, 3, . . . ) is considered. We assume that the loss rate is much smaller than the rate of intrinsic relaxation of the system, so that at any time the state of the system is captured by its rapidity distribution (or, equivalently, by a Generalized Gibbs Ensemble). We give the equation governing the time evolution of the rapidity distribution and we propose a general numerical procedure to solve it. In the asymptotic regimes of vanishing repulsion -where the gas behaves like an ideal Bose gas -and hard-core repulsion -where the gas is mapped to a non-interacting Fermi gas -, we derive analytic formulas. In the latter case, our analytic result shows that losses affect the rapidity distribution in a non-trivial way, the time derivative of the rapidity distribution being both non-linear and non-local in rapidity space.
Introduction. Nowadays, trapped ultracold atoms offer versatile platforms for the study of isolated quantum many-body dynamics. The system is never perfectly isolated though, and it is always weakly coupled to its environment. The main effect breaking unitary evolution is the atom losses. Although of primary interest to understand the limitations of the simulation of quantum manybody physics, a complete theoretical description of losses is still lacking. Different loss processes occur in cold atom experiments: one-body losses might be non negligible; two-body losses due to inelastic two-body collisions are sometimes present [1] or engineered [2]; three-body losses, where a deeply bound diatomic molecule is formed, are always present and are usually dominant [3,4]. All those K-body loss processes -where K is the number of atoms involved, and lost, in each loss event -are local and their effect on the mean atomic density n reads dn/dt = −Gn K Kg K . Here g K is the normalized zerodistance K-body correlation function and G, which has units of length d(K−1) .time −1 for a gas of dimension d, is the constant quantifying the loss process. However, the many-body state at time t is not characterized solely by its atom density n. In particular, the above equation is not a closed equation for n, because the determination of the correlation function g K itself requires knowledge of the many-body state. For chaotic systems, an important simplification occurs if the rate Gn K−1 is much smaller than the relaxation time of the system. Then, as far as local observables are concerned, the system is described at any time by a thermal state which is entirely determined by the mean atomic density n and mean energy density e. The time derivative of the mean atom density n, a local observable, can be computed once g K is calculated for the thermal state [5]. Computing the time evolution of e is more difficult. This can be done for a system where interactions are weak, since the mean energy of a lost atom is simply e/n. The time-evolution of the system can then be computed [6,7].
However, such an analysis is not valid for onedimensional (1D) bosons with point-like repulsive interactions, also known as the Lieb-Liniger gas [8]. This quantum gas is a famous integrable system [9,10], and in the past fifteen years it has been established both experimentally and theoretically that it does not thermalize [11,12], see the reviews in the special issue [13]. Relaxation is still meaningful but, owing to the infinite number of conserved quantities, the system after relaxation is described not only by two quantities n and e, but by a whole function, known as the rapidity distribution [14][15][16][17]. Several early works on atom losses -predating the ones on the absence of thermalization -have focused on the calculation of g 3 in the ground state of the Lieb-Liniger gas [18,19] and g 2 in thermal states [20][21][22][23]. These results were soon extended to excited states of the gas [24,25], culminating in general expressions for g K valid for arbitrary rapidity distributions [26,27]. However, these results are not sufficient to fully describe the evolution of the system under atom losses. Recent attempts in that direction have concentrated on the quasicondensate regime. This regime is characterized by weak correlations between atoms, with g K 1, and is well modeled using a Bogoliubov approximation where the system is described by a collection of independent collective modes. The evolution of the energy in each collective mode under the loss process was computed in Refs. [28,29] and, in particular, it was shown that the system evolves towards a non-thermal state [30]. Predictions concerning phonons were verified experimentally [31][32][33].
In this paper, we revisit the problem of atom losses in the Lieb-Liniger gas. We assume that the loss process is slow compared to the intrinsic dynamics, in such a way that the system at any time is described, locally, by a rapidity distribution ρ(k), or equivalently by a Generalized Gibbs Ensemble (GGE). We give a complete description of the effects of losses by computing the evolution of the rapidity distribution. We devise a numerical pro-cedure valid for any initial state. In the two asymptotic cases where the gas lies in the ideal Bose gas and Tonks-Girardeau [34] (hard core) regimes, we obtain analytical expressions for the evolution of the rapidity distribution.
The rapidity distribution is a key notion for this paper. Its basic definition does not require integrability, and is based on the notion of asymptotic states in scattering theory. Imagine that a homogeneous Lieb-Liniger gas is confined in a flat box potential of length L and that the box potential is suddenly released so that the gas expands freely in 1D. The rapidity distribution ρ(k) of the original state, is simply the density, per unit asymptotic velocity k and per unit length of the original box, of asymptotic particles obtained after an infinite time of this 1D expansion of the gas. What is special about integrable models is that, by elastic and factorised scattering [35,36], the rapidity distribution ρ(k) is conserved by the dynamics. Thus the rapidity distribution is a good parametrisation of any state after relaxation. Crucially, by this definition, ρ(k) is also a measurable quantity. The above thought experiment typical of scattering theory is nothing else but a 1D expansion [37][38][39][40][41]. After a sufficiently long expansion time t 1D the atoms are propagating freely and their velocities are nothing but the rapidities. Thus, measuring the velocity distribution at t 1D amounts to measuring the rapidity distribution, as has been very recently done for a Lieb-Liniger gas in the hard core regime [42]. Note that the atoms' velocity distribution is not conserved by the dynamics and the initial velocity distribution is different from the rapidity distribution. The velocity distribution at t 1D can be measured by time-of-flight, for instance letting the gas expand further in 1D and measuring the density profile, homothetic to that of velocities. The rapidity distribution is an object of central experimental relevance in the investigation of out-of-equilibrium 1D gases, and understanding how it is affected by atom losses is therefore of paramount importance. This is what we do in this paper.
Pinpointing the problem. In principle, one would like to describe the gas by a density matrixρ -not to be confused with the distribution of rapidities ρ(k) -evolving according to the Markovian Lindblad equation (1) Here Ψ(x) is the bosonic atom annihilation operator, and H =´L 0 Ψ + −∂ 2 x /2 + (g/2)Ψ + Ψ Ψ dx is the Hamiltonian of the Lieb-Liniger gas. We set = m = 1. The dimensional parameter G is the same as in the introduction; it quantifies the loss process.
As the size of the density matrixρ is exponential in the number of atoms N , the Lindblad equation is not tractable in physically relevant setups where N ∼ 10 2 − 10 5 . Fortunately, in the asymptotic limit of small G, where the dynamics generated by losses is slow, the complexity of the problem is dramatically reduced.
The key notion that permits simplification of the prob-lem is the notion of relaxation alluded to above. Recall that an isolated Lieb-Liniger gas relaxes at long times, in the sense that averages of local observables approach asymptotic values. For lengths of systems typically found in experiments, τ is much smaller than the Poincaré recurrence time, thus relaxation is a good approximation. As argued above, the asymptotic values of local observables are all determined by a single intensive function, the rapidity distribution ρ(k). How do we evaluate the averages of local observables in terms of ρ(k)? For this purpose, we use the techniques of integrability. The eigenstates of Bethe-Ansatz form naturally encode the asymptotic velocities of the scattering states. An eigenstate is parametrised by a set of rapidities |{λ i } , and its associated rapidity distribution is sim- At large L, the rapidities form a continuum. To compute mean values of local quantities, we consider a subsystem of length where is much larger than the correlation lengths of the system but much smaller than L. The rest of the system acts as a reservoir of rapidities, and the subsystem relaxes to a GGE. Up to corrections of order 1/ , the reduced density matrix of this subsystem is diagonal in the basis of its Bethe ansatz states and takes the form , is determined by entropy maximisation with the constraint that, for all k, . This is a simple generalization [17] of the thermodynamics calculations done in [14].
We assume that the loss process occurs on times much longer than the intrinsic relaxation time of the system, τ . More precisely, the typical rate Gn K−1 is assumed to be much smaller than 1/τ . Then one can assume that, at any time t, the system has relaxed with respect to its Hamiltonian H. Thus, according to the above considerations, the system at any time t is completely determined by its rapidity distribution ρ(k). The reduction of the complexity of the problem stems from having replaced the full density matrixρ by the one-dimensional function ρ(k). To lowest order in Gn K−1 τ , the Lindblad equation leads to Since´ρ(k)dk = n and dn/dt The key problem is to determine the functional F , which is the main goal of this paper.
The functional F as an expectation value of a local observable. In order to determine the functional F , the Lindblad equation (1) is translated into an evolution equation for averages of local quantities q(x), obtained by inserting Eq. (1) into d q(x) /dt = Tr(q(x)dρ/dt).
Eq. (1) is translational invariant so we can assume thatρ also is, and we omit the variable x in q(x) . In Eq. (1), the contribution of the Hamiltonian term can be written, using cyclic invariance of the trace, as the mean value of the commutator [q, H]. The latter is a local quantity since q is local. [Although H is extensive, its commutator with the local operator q is local, as can be seen expanding H and q with respect to the field operators Ψ, Ψ + and their derivatives.] One can therefore use the GGE density matrixρ GGE to represent its average, and we then find that the contribution of this term vanishes since [H,ρ GGE ] = 0. Thus only the non-hermitian term contributes to d q /dt. Using translational invariance of ρ, the integral over x can be recast into a form which involes the total charge Q =´L 0 q(x)dx and we obtain The notation [ρ] means that the mean values are computed for the GGE corresponding to ρ(k), which is justified since the operators inside the brackets are local operators, Q appearing only inside a commmutator with a local operator. For pedagogical purposes, we rederive (3) using a toy model where lost atoms are absorbed by an environment of oscillators, in Appendix B.
The evolution of the distribution of rapidities ρ(k) is obtained by choosing q k such that the total charges Q k are the conserved quantities of eigenvalues Here δ σ (λ) is any approximation of the delta function with a rapidity spread of order σ (for instance a Gaussian of width σ), where we choose σ 1/L; as a consequence, the density q k , of extent σ −1 in position space [43], is local. Choosing σ ∆k, where ∆k is the scale over which ρ(k) varie, one has q k ρ(k) and (3) gives where we used the fact thatρ GGE commutes with the conserved quantity Q k . The formulation (2) of the problem, and the definition (4) of the functional F , is the first main result of our paper. In the following, we use Eq. (4) to compute F . General case: numerical summation over Bethe states. To evaluate Eq. (4), one must be able to calculate expectation values . . . [ρ] . Analytically, this is a very hard problem, and at present there exists no general method to solve it -at least, not for arbitrary repulsion strength g -. Therefore, in this paragraph, we turn to numerics and design a general numerical method to evaluate F .
Eq. (4) can be evaluated numerically in finite size by computing a double sum over Bethe states. The first sum comes from the expectation value . [ρ] , taken with respect to the GGE parameterized by the rapidity distribution ρ(k). This is a diagonal density matrix in the basis of Bethe states |{λ i } , with entries p({λ i }), Z is a normalization factor such that {λi} p({λ i }) = 1, and the weight W [ρ](λ) is related to ρ by the (generalized) thermodynamics Bethe-Ansatz equation of Yang and Yang [14] -with the differential scattering phase The second sum comes from inserting a set of intermediate states, 1 = {µj } |{µ j } {µ j }|, between the observables Ψ +K (0) and [Q k , Ψ K (0)] in Eq. (4). The commutator action is evaluated by using the eigenvalue equation for Q k , and we obtain This expression of F in terms of the form factors {µ i }| Ψ(0) K |{λ j } of the Bethe states is the second main result of our paper. The physical meaning of this equation is clear. If the initial state of the system is |{λ i } , the probability to have a loss event during the time interval dt and that the system at t + dt is found in the In such a case, the final rapidity distribution The probability for the system to stay in the initial state, and thus that q k stays equal to (1/ ) . Computing q k summing over the different cases detailed above, we arrive at Eq. (7).
as well as the K-body correlation in a given Bethe state g K ({λ j }) = {λ j }|Ψ +K Ψ K |{λ j } /n K , we can rewrite Eq. (7) as This enables us to evaluate F numerically, by measuring the expectation value of g K ({λ i }) ( i δ σ (k − λ i ) − j δ σ (k − µ j )) with respect to the probability distribution p({λ i }) p({µ j }|{λ i }). To do this, we sample pairs of Bethe states |{λ i } , |{µ j } , using two Markov chains which have equilibrium distributions p({λ i }) and p({µ j }|{λ i }) respectively. Our two Markov chains are constructed using a Metropolis-Hastings algorithm, by implementing moves of Bethe integers and solving the Bethe equations [8][9][10] with a Newton-Raphson method to find the associated configurations of rapidities {λ i } and {µ j } (see Refs. [44][45][46] for similar numerical summations over Bethe states, and especially Refs. [47,48] for similar samplings of GGEs). Crucially for our method, exact analytical formulas are available for the form factors of Ψ(0) K thanks to recent work by Piroli and Calabrese [49], and for g K ({λ j }) thanks to work by Pozsgai [26]. Our numerical procedure heavily relies on these exact formulas. The procedure is computationally costly, however at present we do not know of any realistic alternative to evaluate the function F numerically. In Fig. (1), we show the results obtained with this method, for the rapidity distribution ρ(k) of a thermal state at T = 0.2n 2 and g/n = 1, for K = 1, 2, 3. The function δ σ (k) approximating the delta function in rapidity space is a Gaussian of width σ = 0.06. To obtain these results we work with N 30 particles in average (the number of particles is let to fluctuate around some fixed mean value in our code), and sum over 10 5 independent pairs of states |{λ i } , |{µ j } . We have checked that the two Markov chains are long enough so that the pairs are truly independent, and that the results do not significantly change as we increase N (see Fig. 1). In total, the computation shown in Fig. 1 takes about 10 hours on a laptop. Notice that, since it is essentially a Monte Carlo integration method, our procedure can be trivially parallelized. This could be important for practical purposes.
Finally, we would like to stress that, in principle, for a sufficiently large system size , a single Bethe state would be sufficient to evaluate the expectation value of a local observable, according to ideas of typicality of eigenstates (see e.g. Ref. [50]). Thus, instead of a double sum, a single sum needs to be evaluated, in principle (see for instance Ref. [51] where this property is exploited to numerically evaluate the dynamical density-density correlation). However, because our observable is the rapidity distribution itself, we find that for our purposes this idea does not work in practice. The discrete nature of the rapidities induces a rugosity of their distribution at finite . Huge system sizes , hardly tractable numerically, would be necessary to mitigate this effect. Numerically, we find that the above method works, while using a single Bethe state doesn't.
Ideal Bose gas limit. When the typical energy per atom E is much larger than both the scattering energy g 2 (= mg 2 / 2 ) and the interaction energy gn, then interactions play a negligible role and the gas is well described by an ideal Bose gas. E g 2 ensures that a collision event between two atoms leads to negligible reflexion [52], while E gn ensures that the bosons are far from the quasicondensate regime. The rapidity distribution is then simply the momentum distribution of the ideal Bose gas: because the gas is free, this momentum distribution is what would be measured by a 1D time of flight. More precisely, defining the canonical bosonic operators Ψ k =´ 0 e −ikx Ψ(x)dx/ √ , which annihilate an atom of momentum k (with values in 2πZ/ ), the Hamiltonian reduces to H = k (k 2 /2)Ψ † k Ψ k and the rapidity distribution is ρ(k) = Ψ + k Ψ k /(2π). In order to compute the evolution of the rapidity distribution due to losses, we use again (4). In this expression, the delta-function approximation for the operator Q k may be chosen simply as δ σ (k) = δ k,0 /(2π). Therefore, Q k = Ψ + k Ψ k /(2π). The density matrix isρ GGE = kρ k whereρ k is gaussian, such that one can use Wick theorem. The commutator in (4) is immediately obtained as and the expression on the right-hand side of (4) is evaluated using Wick's theorem, with contractions Ψ + (0)Ψ(0) [ρ] = n and √ 2π Ψ + (0)Ψ k [ρ] = ρ(k). The result is therefore We see that the rapidity distribution keeps the same form, being simply rescaled in amplitude as time goes  7), for K = 1, 2 respectively. The black dashed line corresponds to KK!ρ(k). The results for K = 1 (resp. K = 2) are obtained by summing over 10 4 (resp. 2.10 3 ) independent pairs of Bethe states for N 30 atoms. We use a Gaussian of width σ = 0.15 to smoothen the rapidity distribution.
on. For a thermal state, ρ(k) 1 2π /(e (k 2 /2−µ)/T − 1) is close to the Bose-Einstein distribution. However, a rescaled Bose-Einstein distribution is no longer a Bose-Einstein distribution, unless the gaz is in the classical regime where T |µ| (i.e. T n 2 ). Thus losses bring the system to a non-thermal state.
In Fig. 2, we display our numerical results for the thermal state at temperature T = 5n 2 and repulsion strength g = 0.1n. This is close to the ideal Bose gas regime, although deviations of ρ(k) from the Bose-Einstein distribution are clearly visible at small k. We compute F numerically using Eq. (7), and we find that the results are in good agreement with formula (10). On the technical side, the numerical evaluation of the double-sum (7) is more difficult than in the regime shown in Fig. 1, because the auto-correlation time of our Markov chain is much longer. We thus work with smaller samples of pairs of Bethe states (10 4 pairs for K = 1 and 2.10 3 for K = 2), which explains the fluctuations visible in Fig. 2 (especially for K = 2).
Tonks-Giradeau limit. The hard-core regime is obtained when the typical energy per atom fulfills E g 2 . The probability to find more than one atom at a given position is vanishing in this regime, so that only the case K = 1 is relevant. We restrict to this case here. It is well-known that, in this regime, the Lieb-Liniger gas is mapped to non-interacting fermions by the Jordan- With the vacuum |0 , both sets of states {c + λ |0 } and {c + µ |0 }, for λ ∈ Z p and µ ∈ Z ap , form a basis of the one-particle Hilbert space L 2 ( ). They are related to each other as and vice-versa. [The fact that this is an involutive change of basis follows from the identity ] The Bethe states, which diagonalise the Hamiltonian, are of the form c † λ1 c † λ2 . . . c † λ N |0 , where λ j ∈ Z p if N is odd and λ j ∈ Z ap if N is even. Thus, while the boundary conditions are always periodic for bosons, the Hamiltonian is diagonalised on fermionic fields with either periodic or anti-periodic boundary conditions depending on the parity of atom number. The rapidities of the Bethe states are the momenta of the fermions and, on each parity sector, the Hamiltonian reduces to that of a non-interacting Fermi gas, H = k (k 2 /2)c + k c k . The GGE density matrix likewise separates into the parity sectors. With the projectors P ± = (1 ± (−1) N )/2 it is of the formρ GGE = P + k∈Zapρ k + P − k∈Zpρ k , whereρ k is diagonal in the Fock basis, so that Wick theorem applies. Using the mixed projector P k = P + δ k∈Zap + P − δ k∈Zp , the conserved charges may be taken as Q k = P k c + k c k /(2π).
One has ρ(k) = P k c + k c k /(2π). Expanding the commutator in Eq. (4), we obtain two terms: The first term is easily computed using Ψ(0) = q c q / √ and Wick's theorem: The second term amounts to computing P k c + k c k on a state obtained from the initial state by the removal of one atom. Crucially, after one loss the parity of the atom number has changed, which induces a sudden change of boundary conditions for the fermions. Thus, in order to use the basis that diagonalises the density matrix, in P k c + k c k one must use the change of basis equation (11). Hence, we have where λ and k are in different sectors. Using Wick's theorem and c(0) = 1 √ λ c λ , one gets Recall that the subsystem size is assumed to be large enough so that ρ(k) varies slowly on the scale 1/ . Then the sums can be replaced by integrals. In order to avoid a divergence in λ ρ(λ) (k−λ) 2 , we rewrite this term as λ ρ(λ)−ρ(k) (k−λ) 2 + 2 4 ρ(k). This leads to ffl ' is the Cauchy principal value of the integral. Combining (13) with (12), we arrive at our final result for the functional F . The evolution of the rapidity distribution under one-body losses in the Tonks-Girardeau gas is determined by Eq. (2) with This formula is the third main result of this paper. It shows that, in the hard-core regime, losses affect the rapidity distribution in a very non-trivial way. The functional F is both non-local in rapidity space -i.e. F [ρ](k) depends on ρ(λ) for any λ, not just on ρ(k) -and nonlinear in ρ(λ). This is in stark contrast with the ideal Bose gas regime, see formula (10). We stress that it is also remarkable because, even though the Tonks-Girardeau gas is mapped to a non-interacting Fermi gas, the effect of losses completely differs from the one of fermionic losses in such a gas. This is of course coming from the non-locality of the Jordan-Wigner transformation. In Fig. 3 we compare formula (14) with numerical evaluation using the above procedure, for the rapidity distribution of a thermal state at T = 1.02n 2 and g = 10 5 n. The agreement is excellent, which further validates our numerical method.
Remarkably, we find that the evolution equation (2) with the loss term (14) can be solved analytically. For an initial distribution ρ 0 (k) at time t = 0, the distribution at time t is given by the exact formula where n 0 =´ρ 0 (k)dk is the initial particle density and i = √ −1. We defer the derivation of that formula to Appendix A.
Numerical time integration The time evolution of the rapidity distribution is obtained by numerical timeintegration of Eq. (2). Fig. 4 shows the resulting evolution of ρ(k) for a gas whose initial rapidity distribution is the thermal state of Fig. 1. The time step is chosen such that the decrease of atom number is 5% of the initial atom number at each step. We perform the same calculation for a gas that lies deep into the hard-core regime, and compare to the exact result (15). The agreement is excellent, which shows that the time step is sufficiently small to provide accurate predictions.
Inhomogeneous profiles and Generalized Hydrodynamics. Our equation (2) is readily generalized to account for the evolution of an inhomogeneous Lieb-Liniger gas, for instance in the presence of an external potential V (x). At large scales the gas is described by a positiondependent distribution of rapidities ρ(x, k), which evolves according to  Fig. 3. To benchmark our method we compare our numerical results with the analytical formula (15) (black dashed line): the agreement is excellent. Right: three-body losses, away from the hard-core regime. The initial state is the same as in Fig. 1. where the "effective velocity" v eff is a functional of ρ defined in Refs. [53][54][55]. The nonzero term on the righthand side extends the Euler-scale Generalized Hydrodynamics equations [54][55][56][57][58][59] to include atom losses; other types of integrability breaking situations have been studied recently [41,[60][61][62][63][64]. Solving equation (16) is a challenging numerical problem, because it requires the calculation of F for many different rapidity distributions. It may be doable with the methods we presented, using a large amount of parallelization. Analytical progress on the evaluation of the sum (7)-(8) would also be desirable and could possibly lead to drastic reduction of the computational time needed to evaluate F , thus facilitating a numerical solution of Eq. (16). This would lead to important improvements in the theoretical modeling of out-of-equilibrium cold atom experiments by GHD [65].
Conclusion. This work on the effect of losses in the one-dimensional Bose gas, an integrable system, calls for further studies.
First, in the hard-core regime, it is possible to show from Eq. (14) that F [ρ](k) generically behaves as 1/k 4 at large k. Thus, in sharp contrast with its thermal equilibrium distribution, the gas typically develops 1/k 4 tails in its rapidity distribution because of atom losses. This stems from the short-range correlations between atoms in the Lieb-Liniger model: right after a loss event, the many-body wavefunction presents a cusp at the position of the lost atom. We expect this feature to exist also beyond the hard-core regime, and we hope to come back to this in future works. This observation is of prime importance for experimental simulations of the Lieb-Liniger model, where proper characterisation of the initial state is required. It could explain the experimental evidence for non-Gibbs ensembles in cold atoms experiment [30].
Second, it is clear that quantitative comparisons with experimental data, where the gas is usually inhomogeneous (see Eq. (16)), requires further analytical and numerical developments, in order to facilitate the evaluation of the functional F . Analytical progress on the thermodynamic limit of the form factors of ψ K (0) [49] would be needed (see e.g. Refs. [66][67][68][69][70] for such studies of form factors of other operators), as well as methods for resumming those form factors (see e.g. Refs. [71][72][73][74][75]). These are very challenging tasks. A simpler problem, which might serve as a good starting point for further analytical developments, would be to compute the functional F at low temperature using an effective Luttinger liquid approach, or more generally in states close to zero-entropy states where this approach can be generalized [76][77][78][79].
Finally, the link between the results of this paper and those previously obtained in the quasi-BEC regime, both theoretically [28,30] and experimentally [32,33], remains to be made. toy model with unitary Hamiltonian evolution representing the loss processes. This is done by connecting the Lieb-Liniger model to an external environment, or reservoir, able to absorb the atoms. Such toy models are in fact a standard way of deriving the full Lindblad equation, and the derivation we present follows the textbook discussions on this subject, see for instance [80,81]. This will also make the connection of our work with recent works on perturbation of integrable models, especially [41,60,61,63], more explicit. From the prespective of the toy model, Eq. (3) is nothing else but the secondorder perturbation theory evolution equation discussed in these works.
Let H be the Lieb-Liniger Hamiltonian, and consider H = H + γV + H E where V represents the interaction with the environment, and H E is the environment's Hamiltonian. The environment interacting with the atoms at any given position, may be thought of as being composed of a family of harmonic oscillators of all frequencies. As we want to describe loss processes occurring at every point in space, the total environment is composed of one such family for every point x. Thus, we have canonical oscillators c(x, ω), where x represents the position and ω the frequency, with The environment's Hamiltonian acts as e iH E t c + (x, ω)e −iH E t = e iωt c + (x, ω).
The interaction describing losses is simply that where atoms are exchanged between the Lieb-Liniger gas at position x, and the family of oscillators at position x. In order to describe exchanges which are instantaneous in time (that is, assuming that the time taken for the loss to occur is much smaller than the typical evolution timescales of the gas), the interaction -or equivalently the distribution of oscillators in the environment -is taken to be flat in frequency space, with an infinite band of frequencies. Further, in order to ensure sufficient decoherence, the frequency ω = 0 is not coupled. This gives V =ˆω =0 dxdω Ψ +K (x)c(x, ω) + c + (x, ω)Ψ K (x) .

(B3)
Similarly to what is done in the main text, we make the assumptions of homogeneity of the full system, of a small interaction strength γ, and of relaxation between interaction events (which are, here, loss events). We are interested in the evolution of the Lieb-Liniger gas and the environment with respect to the full Hamiltonian H under these assumptions. In these assumptions, the relaxation is towards the GGEs with respect to the subsystem H 0 = H + H E ; like H, this subsystem is also integrable. This GGE is described by a rapidity distribution ρ(λ) for the Lieb-Liniger gas, and a distribution of environment's oscillators f (ω) defined as c + (x, ω)c(x , ω ) [ρ,f ] = δ(x − x )δ(ω − ω )f (ω). The evolution of any conserved density q(x) (conserved with respect to H 0 ) under these assumptions can be obtained from a standard second-order perturbation theory, and takes the form where V (s) = e iH0t V e −iH0t and v = dω Ψ + (0)c(0, ω) + c + (0, ω)Ψ(0) . This general equation, written, with Q = Q k , as an evolution equation for the rapidity density, is at the basis of the Boltzmann kinetic formulation of perturbed integrable models [41,60,61,63].
As we wish to describe loss events, and not events where particles are re-absorbed by the Lieb-Liniger gas, we take the initial state of the environment to be the vacuum, f (ω) = 0. By looking at the evolution of the conserved densities q ω (x) = c + (x, ω)c(x, ω), one can show that the environment's state stays the vacuum throughout time under the evolution equation (B4). Using the vacuum property and the canonical commutation relations (B1), one then finds, for all conserved densities q k of the Lieb-Liniger gas, where the environment's contribution has been factored out. This is indeed Eq. (3) where Q = Q k is taken to be a conserved quantity, and where we identify γ 2 = G.