The one-dimensional Bose gas with strong two-body losses: the effect of the harmonic confinement

We study the dynamics of a one-dimensional Bose gas in presence of strong two-body losses. In this dissipative quantum Zeno regime, the gas fermionises and its dynamics can be described with a simple set of rate equations. Employing the local density approximation and a Boltzmann-like dynamical equation, the description is easily extended to take into account an external potential. We show that in the absence of confinement the population is depleted in an anomalous way and that the gas behaves as a low-temperature classical gas. The harmonic confinement accelerates the depopulation of the gas and introduces a novel decay regime, which we thoroughly characterise.


Introduction
The study of quantum many-body physics with ultra-cold gases has recently attracted considerable attention [1]. The unprecedented possibility of studying a quantum many-body system in isolated conditions and with negligible effects from the environment has made possible the experimental characterisation of its closed-system dynamics in well-controlled settings, see for instance the experiments reported in Refs. [2][3][4]. Restricting our focus to one-dimensional Bose gases that at low temperature are described by the Lieb-Liniger model [5,6], several developments, among which the generalised hydrodynamics [7,8], have produced a theoretical framework that successfully describes the experiments where bosonic gases are put out of equilibrium via a quantum quench of the external potential [9][10][11].
Yet, even these quantum systems cannot be considered as perfectly isolated, and the most relevant form of environment to which they are coupled is represented by the surrounding vacuum into which atoms or molecules can leak [12][13][14]. Even if current experiments do not suffer from significant loss effects, the community has recently witnessed a renewed interest in the dynamics that they induce in correlated quantum gases. On the one hand, they need to be taken into account in order to develop a quantitative description of experiments [9,[12][13][14]. On the other hand, it was early recognised that open-system dynamics (such as losses) can be extremely interesting and can lead to the appearance of phenomena where quantum mechanics and quantum coherence play an important role [15][16][17][18][19][20]. One well-explored phenomenon is that of cooling by atom losses [21][22][23][24]. Moreover, the recent possibility of inducing and engineering losses in well-controlled conditions has finally opened the path to the experimental study and quantum simulation of many-body dissipative dynamics [25].
In this article we focus on the one-dimensional Bose gas with two-body losses. From an experimental viewpoint, our work has several motivations. First, the recent advances in the cooling to quantum degeneracy of novel atomic species have produced ultra-cold gases characterised by strong inelastic interactions. For instance, ytterbium atoms are characterised by an excited and long-lived metastable state that suffers from two-body recombination processes, so that a quantum-degenerate gas in this metastable state is a natural setup to study the interplay between two-body losses and quantum dynamics [26][27][28]. Second, photoassociation processes can induce two-body losses with controllable strength in standard alkaline atoms, and have already been used in several contexts, such as in the study of dissipative phase transition [25]. Last, the use of molecules gives the possibility of studying a quantum degenerate gas subject to two-body losses and it is one of the setups where these experimental studies were first performed [29,30]. For completeness we mention that experiments exist also in the fermionic case [31][32][33].
From the theoretical viewpoint, a first important step in the description of the stronglycorrelated and lossy Bose gas has been presented in Refs. [34,35], focusing primarily on the limit of weak losses (the case of strong interactions and fermionisation is discussed only for one-body losses). Here, instead, we are interested in the case of strong two-body losses, and our goal is to characterise the simplest experimental observable, namely the dynamics of the total number of particles composing the gas. This problem has already been discussed for a lattice gas described by the Bose-Hubbard model without any external confinement in Ref. [36]. There, it is shown that the gas effectively fermionises, namely that it becomes a gas of hard-core bosons, thus developing strong particle-particle correlations and a peculiar momentum distribution function that does not resemble any equilibrium situation. Concerning the continuum situation, in Ref. [37] the authors present a solution of the dissipative Lieb-Liniger model, where the interaction constant has also an imaginary part. The study, that has the merit of highlighting in a transparent way the fermionisation induced by strong two-body losses, does not produce any prediction for the dynamics of the population of the gas. Finally, we mention the existence of a conspicuous literature focusing on the problem of two-body losses in one-dimensional fermionic gases [38][39][40][41].
In this article we present a theoretical description of the effect of a harmonic confinement on the dynamics of the one-dimensional Bose gas with strong two-body losses. As a starting point, and for pedagogical reasons, we show that in the absence of harmonic confinement the longtime dynamics deviates from naive mean-field expectations; in particular we show that the gas is turned into a low-temperature classical gas described by a Maxwell-Boltzmann momentum distribution function. The harmonic confinement has a qualitative impact on such behaviour. Two regimes are identified: that of weak confinement, where the period of an oscillation is long with respect to the loss rate, and of strong confinement, with opposite properties. Whereas in the former case the dynamics has all qualitative features of the homogeneous system, in the latter case we find that the harmonic trap accelerates the depletion of the gas. The methodology that we employ relies on the fact that strong two-body losses fermionise the gas; the dynamical properties are indeed obtained by characterising the population of the emergent fermionic momenta (or rapidities) using a Boltzmann-like equation.
The article is organised as follows. In Sec. 2 we introduce the theoretical model, we define the limit of strong two-body losses and present the set of rate equations that describe its dynamics in the absence of any external potential. In Sec. 3 we solve the aforementioned equations assuming an equilibrium initial configuration. In Sec. 4 we consider the presence of the harmonic confinement: we generalise the equations describing the homogeneous gas and discuss the dynamics of an initial equilibrium state. Our conclusions are presented in Sec. 5. Five appendices conclude the article.
2 Homogeneous Bose gas with strong two-body losses: the Quantum Zeno effect We consider a bosonic gas trapped in one dimension and subject to strong two-body losses; for simplicity, we consider a homogeneous sample of length L with periodic boundary conditions. We introduce a field Ψ(x) satisfying canonical commutation relations [Ψ(x), Ψ(x )] = 0 and [Ψ(x), Ψ † (x )] = δ(x − x ), and describe the dynamics using the Lindblad master equation, which make sense only in 1D (see Ref. [20] with: In Eq. (1a) we recognize the Lieb-Liniger model; the term in Eq. (1b) accounts for twobody losses. A study of this model in the weakly-dissipative limit has been presented in Ref. [34,35]; an ab-initio derivation of the master equation and critical discussion of the underlying approximations can be found in Ref. [37].

Strong dissipation and fermionisation
We are interested in studying this model in the strongly-dissipative regime, which is often called quantum Zeno regime. In order to better clarify this limit, it is useful to consider the non-Hermitian Hamiltonian associated to the master equation; we thus introduce the non-Hermitian version of the Lieb-Liniger model that has been thoroughly discussed in Ref. [37]: The strength of the non-linear term which accounts for the two-body elastic and inelastic interaction, is quantified by the dimensionless parameter: where n is the one-dimensional density of the gas (H NHLL is number-conserving and thus the density n = N/L is well-defined). According to Ref. [37], the model can be mapped to a Tonks-Girardeau gas whenever |ξ| → ∞. It is well-known that the Bose gas becomes a gas of hard-core particles when g 2 n/m; this is true also when g = 0 and the following inequality is satisfied: This inequality defines the regime of strong dissipation. If we consider a 87 Rb gas (atomic mass: m = 1.43×10 −25 Kg) with typical density n 1×10 7 m −1 as in Ref. [9], we obtain that in order to enter the strongly-dissipative regime, the inequality reads: Although the dissipation does not conserve the number of particles, it conserves the parity of the number of particles. Thus, in order to study the problem in the sector with an even number of particles, we introduce the field operators in momentum space via Fourier transform: c k = L − 1 2 c(x)e −ikx dx where k = π(2n + 1)/L and n ∈ Z. The Hamiltonian in momentum space reads: Summarizing, in the limit of infinitely-strong dissipation |ξ| → ∞, Hamiltonian (2) is welldescribed by a non-interacting fermionic model, whose eigenstates and eigenenergies are defined by the occupation numbers of the modes k. This information is of great help also for discussing the stable modes of the master equation (1). Let us for instance consider a density matrix that is diagonal in the basis of the modes k and that is solely parametrised by a set of {λ k } coefficients related to the occupation of the modes: In order to prove that the latter density matrix is a stationary state of the master equation (1) we recast Eqs. (1) as: Within the fermionised approximation, H NHLL is given by Eq. (6), which commutes with the density matrix (7); we consequently conclude that the first part of the r.h.s. of Eq. (8) vanishes. The part proportional to γ vanishes as well because Ψ(x) 2 is equal to zero for a fermionised state. We thus conclude that the fermionic modes identified above with the non-Hermitian Lieb-Liniger model are stable modes of the master equation (1). The operators {c † k c k } k constitute an infinite set of operators that commute with the Hamiltonian; the state (7) is the associated generalised Gibbs ensemble (GGE) [42]. The fermionic momenta, also called rapidities in the Bethe-ansatz context, can be observed in experiments by letting the gas expand in the one-dimensional tube [10,[43][44][45][46]. This is different from performing a standard time-of-flight experiment on a one-dimensional gas: in this case the gas expands freely in three-dimensional space, and one instead probes the bosonic momentum distribution function [1].

Rate equations for the dissipative dynamics
In reality, in any experimental setting, the parameter ξ is never going to have infinite modulus. The main consequence of this fact is that the fermionic excitations described by the c k acquire a finite decay time and become quasi-stable. From an experimental viewpoint, this means that particles leak from the sample and the setup is depleted. It is interesting to observe that if the initial density satisfies the inequality (4), the gas will remain in the stronglydissipative regime at all times because n(t) is a non-increasing function (the system can only lose particles). The curious thing about the system that we are studying is that even if we are working in the strongly-dissipative regime, we have identified a set of modes (or quasiparticles) which are almost stable. The dynamics of these modes is thus weakly dissipative: as we will see, their decay rate scales as 1/γ and this counterintuitive effect is known in the literature as environment-induced quantum Zeno effect [15,47].
It has already been observed that in the presence of weak losses, the dynamics can be described with a time-dependent state that is a simple generalisation of the GGE proposed in (7) (see Refs. [36,48,49]): Ref. [48] presents the equations that determine the dynamics of the λ k (t) in general. However, there is a one-to-one correspondence between λ k (t) and in this work, instead of working with λ k (t) we prefer to work with n k (t).
The first result of this article is the identification of the differential equations obeyed by the n k (t) for this specific problem: As anticipated, Γ eff , which dictates the typical decay rate of the fermionic modes, scales as γ −1 . Considering again the abovementioned 87 Rb gas discussed in Ref. [9], we estimate g ∼ 8.9 × 10 −36 kg·m 3 ·s −2 and take a value for γ that is in the strongly-dissipative regime, γ = 2 × 10 −2 m·s −1 . We obtain: Γ eff = 4.58 × 10 −19 m 3 ·s −1 . The rate equations (11) generalise those already obtained for a one-dimensional bosonic gas trapped in a lattice and described by the Bose-Hubbard model [36]. The derivation of the equations is given here below and it is obtained by regularising the problem onto a lattice. The uninterested reader can jump directly to Sec. 3 without compromising the understanding of the rest of the article.

Derivation of the rate equations
In order to derive the rate equations (11), we regularise the original problem introducing a short-length cutoff, a, that is the shortest length-scale of the problem. In this way, the original master equation is turned into a lattice problem. After introducing the lattice bosonic operators b † m , the master equation reads: The relation between the original parameters and the lattice ones is given by: Ja 2 ↔ 2 /(2m); U a ↔ g, and γ L a ↔ γ. Concerning the relation between the physical length of the original system L and the number of lattice points M it holds that M a ↔ L.
Following the methods employed in Ref. [36], we fermionise the problem, introduce the momenta k = 2π L (n + 1), with n = 1, 2, . . . M and obtain the following rate equations for the occupation number of each fermionic mode: We now take the limit a → 0 + and obtain the rate equations that describe our system in the continuum limit. First, we observe that: Since L is fixed and a is going to zero, the number of lattice points M is increasing; when M 1 we can change the sum into an integral via k → L 2π dk and using L/M = a we obtain: Let us now take the limit a k −1 max , where k max is the maximal wavevector that is going to be occupied during the dynamics. Since each function n k (t) is monotously decreasing (its derivative in Eq. (13) is never positive), it is enough to consider the largest wavevector that is occupied at time t = 0; this quantity is well defined. In this limit: sin(ka) ∼ ka for all relevant k; moreover the integration limits can be safely expanded to ±∞. We obtain: Note that Γ eff has a well-defined continuum limit and the good dimensions of L 3 × T −1 , since n k is an adimensional quantity. This concludes our derivation.

Depletion of an initial equilibrium state
We discuss the depletion dynamics of a Bose gas prepared in the ground state of the Tonks-Girardeau gas (|ξ| → ∞) with density n in . We consider the rate equations (11) with the following initial conditions: n k (0) = 1 for |k| < πn in and n k (0) = 0 otherwise. The initial fermionic momentum profile, n k (0), is symmetric with respect to k → −k transformations, and thus (k − q) 2 n q dq = (k 2 + q 2 )n q dq. Since the master equation is invariant under k → −k exchanges, this property is preserved during the whole dynamics and thus Eq. (11) is substituted by: In order to get a better understanding of the problem we introduce the rescaled momentum k = k/n in and the rescaled timet = n 3 in Γ eff t so that the equations read: The adimensional normalised density of the gas readsñ(t) = 1 2π nk(t)dk and it satisfies n(0) = 1. It is linked to the physical density via the simple relation n = n in ×ñ.

"Heating" to a low-temperature classical gas
With the aim of studying the dynamics of the gas, we employ Eq. (18) in order to write the two following equations (see App. A for details on the derivation): If we divide the latter equation by nk(t) we obtain the following relation: which is solved by (see again App. A for the details): This latter relation shows that nk(t) is fully determined byñ(t), and in particular that momenta distribute according to a Gaussian function centred atk = 0 and truncated at k = ±π. The variance of the Gaussian is 4π t 0ñ (t )dt −1 , which decays to zero ast − 1 2 in the long-time limit (see below Sec. 3.2). Thus, at sufficiently long times, we can disregard the truncation atk = ±π and interpret the distribution as a Maxwell-Boltzmann classical 2mk B T with time-dependent temperature and chemical potential. Reintroducing for clarity the dimensional units, they read: The initial quantum gas at zero-temperature is turned, at long times, into a low-temperature and low-density classical gas. Note that at long times µ k B T , which is consistent with the proposed Maxwell-Boltzmann interpretation.

Long-time behaviour
We define the function ν(t) = t 0ñ (t )dt and taking the integral onk ∈ [−π, π] of Eq. (21) we obtain: in terms of the function ν(t) we have: A discussion of the properties of the system at short times is in Appendix B. At long time, we expect that ν(t) → ∞ and consequently we insert in Eq. (24) the value lim x→∞ Erf[x] = 1. We obtain: The statement used in the previous section to discuss the variance of nk(t) is thus proved. By differentiating with respect tot, we find the asymptotic behaviour of the density, which for clarity we present here also in dimensional units: We find a long-time behaviour characterised by n(t) ∼ t −1/2 , similarly to what has already been found in the analogous lattice problem [36]. In order to appreciate the interest of this result, it is useful to compare it with the mean-field prediction n(t) ∼ t −1 obtained from the equation ∂ t n = −κn 2 , which follows from the approximation of uncorrelated bosons. Indeed, the correct equation that describes the population of the homogeneous gas under local twobody losses reads [34,37] In weakly-correlated bosonic gases, it is often the case that g (2) (0) = 1 is a good approximation (uncorrelated bosons), and the mean-field behaviour is recovered. The typical decay time scales as (n 3 in Γ eff ) −1 , so that a dense gas has a shorter lifetime than a dilute one. To get a more concrete idea about this time-scale, we consider once more the example of the 87 Rb gas discussed above: for an initial density of 1 × 10 7 m −1 we obtain a typical time of 2 × 10 −3 s. The Zeno effect can be fully appreciated if one compares this value with the original typical decay time of the gas: (γn in ) −1 ∼ 5 × 10 −6 s.

Numerical solution
In order to test the previous predictions, we have solved the Eqs. (18) with a 4 th -order Runge-Kutta numerical algorithm; the integration step is dt = 10 −3 and N step = 2 · 10 4 , while we have discretised thek space in 10 3 points in the interval [−π, π]. Comparing with the results obtained employing dt = 10 −4 , we estimate a relative error of 10 −6 , completely negligible for our purposes. The results are summarized in Fig. 1. We first compare the numericallycomputed densityñ(t) with the long-time behaviour given by Eq. (26); the latter faithfully describes the behaviour ofñ(t) even for values oft that are of order 10 −1 . Concerning the nk, we observe that the initial Fermi-sea att = 0 first evolves into a Gaussian that is truncated atk = ±π; then, for sufficiently long times, it is possible to disregard this truncation and interpret the distribution as a Maxwell-Boltzmann classical distribution.

An analytical formula forñ(t)
We conclude this section proposing the following formula to describe the dynamics ofñ(t): where A and B are two coefficients to be determined. Notice that the same formula has been proposed in Ref. [36]. This functional form is motivated by the fact that it can capture both the short-time and the long-time behaviours. The short-time behaviour is discussed in Appendix B, and in particular in Eq. (66) we show thatñ(t) ∼ 1 − αt. The long-time behaviour is discussed in Eq. (26), whereñ(t) ∼ (4π) −1t−1/2 . By comparing our ansatz with the two exact limits, we obtain two conditions for the coefficients A and B that are solved by the following values (details are in Appendix C): In Fig. 1 left panel we compare the numerically-computedñ(t) with Eq. (28); the plot shows that our formula describes very well the numerical data.

The harmonic confinement
We now include in our discussion the presence of an external potential, focusing on the experimentally-relevant case of a harmonic confinement: In order to generalise the previous approach to an inhomogeneous situation, we employ the local-density approximation (LDA) and promote each fermionic momentum occupation number to a space-dependent quantity: n k (t) → f (x, k, t), where f (x, k, t) should be intended as an average over a small macroscopic region centered around x. The spatial density is obtained integrating over all possible momenta: The time evolution of f (x, k, t) is dictated by a Boltzmann-like equation that includes the effect of losses: The force field depends on the potential; in the case of the harmonic potential: The initial condition can be determined using the LDA for the equilibrium density profile n in (x) of a gas in a harmonic potential: Here, R represents the LDA radius of the gas, N in is the initial number of particles and This approach is an example of the treatments based on the so-called generalised hydrodynamics that has been recently introduced to discuss the dynamics of one-dimensional integrable models with Bethe-ansatz techniques.

Dimensionless units
In order to get a better understanding of the problem, we introduce the following dimensionless variables:x Notice thatk has a different definition with respect to the one presented in Sec. 3. This choice is motivated by the fact that in these units the initial condition (34) takes the particularly simple form of a circle with unitary radius: f (x,k, 0) = 1 forx 2 +k 2 ≤ 1; 0 otherwise.
Also the Boltzmann-like equation (32) gets an intuitive form when thex andk coordinates are employed. We rescale time introducing: and we obtain: As it is standard in the classical motion of the harmonic oscillator in phase space, in the absence of losses the initial unitary circle is remapped onto itself and rotates with a period T ω = 2π/ω. Note that the losses do not scatter momenta outside the unitary circle, so that the dynamics remains confined within it. The ratio ω/Γ emerges as the relevant parameter that measures the competition between the harmonic confinement and two-body losses.
Similarly to what we have done in the homogeneous case, our focus will be the time evolution of the rescaled number of atomsÑ , defined asÑ = N/N (0); since the area of the initial circle is π, we have:Ñ

Numerical solution
We simulate Eq. (38) with a numerical algorithm; details about the actual implementation are reported in Appendix D.
Our results for f (x,k,t) are presented in Fig. 2 for four different values of ω/Γ. For some values of ω/Γ a link to the animated version of the dynamics is provided in footnote 1 . The time dependence of the integrated and normalised populationÑ (t) is reported in Fig. 3; the observation ofÑ (t) shows the existence of two limiting behaviours, for small and large values of ω/Γ, dubbed weak-confinement and strong-confinement regimes, whose description will be the object of the next sections.
To make the discussion more concrete, we identify the value of the trapping frequency that separates the two regimes. The relation ω/Γ = 1 yields the following equation for the trap frequency: Taking the numerical values obtained with the previous numerical estimates and an initial number of atoms N in = 10 3 we find ω = 0.22 s −1 ; this is an extremely small trapping frequency, given that those implemented in standard labs are of order 10 2 s −1 . The observation of the weak confinement regime thus requires the use of a more dilute and less populated gas. If we consider the same parameters used in the previous estimate but take a lower density value, n in = 5 × 10 6 m −1 , and an initial population of N in = 100 bosons, we obtain ω = 226.42 s −1 .

The limit of weak confinement: ω/Γ 1
We begin with the analysis of the data in Fig. 2 for ω/Γ ≤ 1, plotted in the first two lines. We observe that the gas depletion is stronger for higher values ofk, a phenomenon that was characterising also the dynamics of the homogeneous gas, see Fig. 1. Moreover, in this inhomogeneous scenario, the losses are more effective in the centre of the trap, where the density is higher; a faster dynamics at higher densities was also observed in the homogeneous system, see Eq. (26). At the beginning of the dynamics the most long-lived population is thus located at the edges of the trap, where the density is lower and momenta are smaller. This is the first consequence of the presence of the trap: it changes the density profile of the gas with a faster depletion of the population at the centre. In order to present a quantitative theory of the short-time dynamics in the weak-confinement limit, we present an analytical solution for the case ω = 0. Note that we are setting ω = 0 in Eq. (38) but we are keeping the initial condition (36), which is shaped by the presence of the trap. In order to better visualize the depletion phenomenon in real space, in Fig. 4 we show the behavior of the spatial densitỹ  pative evolution leads to a spatial profile much denser at the boundaries with respect to the center of the trap.
Thanks to the LDA, we can model the gas as composed of several independent and homogeneous subparts located at different points of the trap and with different initial densities. According to our study of the homogeneous lossy gas in Sec. 3, each of these subparts features a long-time decay proportional tot − 1 2 . The long-time behaviour of the trapped gas is obtained by integrating all these contributions; after some algebra reported in Appendix E we obtain: and Γ(x) is the Euler Γ function. This result is plotted in Fig. 3, where we observe that it reproduces quantitatively the long-time behaviour of the gas.

From weak to strong confinement
In the weak-confinement regime, once the inhomogeneous density profile shown in Fig. 4 has been created, the presence of the harmonic trap has an additional effect and induces a rotation in phase space (see for instance Fig. 2 for ω/Γ = 0.1, second line). At time t ∼ T ω /4 the long-lived bosons located initially at the edges of the trap have moved to the centre, where they are quickly lost because at the centre of the trap the loss mechanism is more effective (see for instance Fig. 2 for ω/Γ = 1.0, third line). The time t ∼ T ω /4 thus marks the onset of a novel behaviour that is clearly displayed in the plots ofÑ (t) reported in Fig. 3. Since the time T ω /4 corresponds to the rescaled timet = πΓ/(2ω), we observe that when the initial confinement is weak, ω/Γ 1,Ñ (t) departs from the ω = 0 curve and after some oscillations collapses onto a new curve that describes the depopulation in the strong-confinement regime.
We can make sense of this transition by comparing the frequency of the rotation in phase space with the instantaneous decay-time of the gas. We introduce the following adimensional parameter: For r 1 we have a weakly-confined gas, whereas for r 1 we have a strongly-confined gas. For a gas that is initially weakly-confined, using the formula in Eq. (42) we obtain that r(t) = 1/(2ωt), a decreasing function of time. Thus, there must be a time t tr at which the typical decay-time becomes longer than the trapping frequency. For t t tr the gas becomes effectively strongly-confined and the previous approximation cannot hold anymore: the gas cannot be described by Eq. (42). If we set this characteristic time as a quarter of the harmonic period, i.e. ωt tr = π/2, we find which gives r(t tr ) = π, i.e. a value of order unity for the parameter r at t = t tr . Thus, with completely different arguments, we have found again the importance of the harmonic period in signaling the transition from a weakly-confined behaviour to a strongly-confined one.

The limit of strong confinement: ω/Γ 1
We observe in Fig. 3 that for ω/Γ ≥ 1 all curves collapse onto a universal function. In this limit, the dynamics is strongly determined by the trap, and it results in a behaviour that is well described by the formula:Ñ as it is shown in the plot. We now present a simple model that can explain the formula in Eq. (45). For ω/Γ ≥ 1, the inhomogeneities in the momentum distribution function that are created by the losses are rapidly washed out by the action of the trap, that rotates the f (x,k,t) in phase space (see for instance Fig. 2 for ω/Γ = 10, fourth line). Since this rotation takes place on the shortest time-scale of the problem, we introduce a novel distribution function g(ẽ,t) that does not depend onx andk in a separate way, but only on the dimensionless energyẽ =x 2 +k 2 that is conserved by the harmonic-oscillator dynamics; normalisation is ensured by the following relation:Ñ (t) = 1 π ∞ 0 g(ẽ,t)dẽ. This "radial" symmetry in phase space is well highlighted by the plots in Fig. 2 for ω/Γ = 10. In the next section, we show how to derive from Eqs. (38) the rate equation obeyed by g(ẽ,t), which we anticipate here: We did not find an analytical solution of Eq. (46), but we could easily produce a numerical solution for the initial condition: g(e, 0) = 2πe, which satisfies 1 0 g(e, 0)de = π, the area of the circle. We display in Fig. 5 the numerically-computedÑ (t); the result is well described by Eq. (45). The appearance of this novel behaviour leads to an interesting qualitative observation: the presence of a harmonic confinement accelerates the depletion of the gas. Indeed, in the presence of a trap the depletion is eventually scaling as 1/t, whereas for a homogeneous system it scales as 1/t 1/2 . Clearly, the previous is a statement concerning the long-time asymptotics: the possibility of observing experimentally the decay as 1/t 1/2 is not ruled out if one works with a shallow trap and in the appropriate time regime.
Finally, a last consistency check. If we compute the parameter r(t) for the gas in the strong-confinement limit using the formula in Eq. (45), we obtain r(t) = 1/(ωt). Again, we encounter a decreasing function of time. Thus, if we start with r 1, the system is expected to remain in the strongly-confined regime, and the description is self-consistent at all considered times.
4.6 Rate equations in the strong-confinement limit ω/Γ 1 As discussed above, in the strong-confinement limit we consider the phase-space density function g(ẽ,t), and the normalised population readsÑ (t) = 1 π 1 0 g(ẽ,t)dẽ. We now identify the evolution equation that is satisfied by g(ẽ,t).
First, we observe that the part of the Boltzmann-like equation that is responsible for the classical motion in phase space leaves g(ẽ,t) unchanged. Indeed: We then move to the loss term, which is local in space. Since the particles with energy˜ are assumed to be uniformly spread in the region of the phase space with energy˜ , we need to introduce the probability density p(x,˜ ) that describes the fraction of particles with energỹ that are located at positionx. With simple calculations we obtain: It is a good probability density, and indeed is satisfies √˜ − √˜ p(x,˜ )dx = 1, where we have used that the extrema of the classically-allowed region are ± √˜ . With this notation, the number of particles with energy in [ẽ,ẽ + dẽ] located in the infinitesimal space region [x,x + dx] is: Let us now focus on the loss of particles with energyẽ due to the presence of particles with energy˜ . Recall that losses depend on the squared momentum difference (k −q) 2 at the same positionx. The particles with energyẽ at positionx can have two equally-likely momenta: ±kẽ ,x , withkẽ ,x = + √ẽ −x 2 . Note that π ×kẽ ,x × p(ẽ,x) = 1. The same is true for particles with energy˜ at positionx. Thus, concerning the loss-rate momentum dependence (k −q) 2 , four different combination of momenta are then possible: each with probability density Let us first consider˜ <ẽ. Here, we have to restrict the integration overx to ± √˜ , which is the maximal spatial extension of the particles with smallest energy; for |x| > √˜ no interaction is possible. The loss of particles with energyẽ due to the presence of particles with energy˜ reads: This expression can be simplified, because (kẽ ,x −k˜ ,x ) 2 + (kẽ ,x +k˜ ,x ) 2 = 2k 2 e,x + 2k 2 ,x . The generalisation to the situation˜ >ẽ is simple.
By integrating over˜ , we write: With some algebra, we obtain This equation is amenable to the compact writing which was presented in Eq. (46).

Conclusions and perspectives
In this article we have studied the effect of a harmonic confinement on the dynamics of a one-dimensional Bose gas subject to strong two-body losses. Thanks to the fermionisation that takes place in the regime of strong dissipation, the so-called quantum Zeno regime, we have developed a set of rate equations that describe the depopulation of each fermionic mode.
In the homogeneous case we have predicted both the dynamics of the density of the gas, displaying an anomalous t − 1 2 decay, and the time-dependence of the fermionic momentum distribution function. This latter quantity, which is often called distribution of rapidities, is a measurable quantity, via an expansion in a one-dimensional setting. The main result of the article is that the presence of a harmonic confinement leads to a drastic modification of the dynamics of the gas; in particular, when the trap is strong, the density of the gas decreases as t −1 . We have shown the existence of a crossover time that signals the passage from the behaviour of a homogeneous-gas at short times to that of a strongly-confined gas at long times. Several of our predictions can be tested in state-of-the-art experiments.
Two directions of further study can be devised. The first one concerns the fact that during the dynamics the gas might leave the purely one-dimensional situation, and start to populate excited levels of the transverse confinement. The fact that we are working in the stronglydissipative Zeno regime makes this not implausible because experiments are typically longer than standard ones. Usually, including higher energy bands in a numerical or analytical description can be extremely challenging, as it requires the inclusion of a novel degree of freedom, which enlarges exponentially the Hilbert space to be considered. In our case, instead, the rate-equation model scales linearly with the number of fermionic modes. As such, the inclusion of a higher-energy band within this framework is feasible and can be an important step to quantitatively describe the experimental data.
A second investigation direction concerns the study of strong and finite elastic interactions, which could be responsible for the scattering of the fermionic momenta in phase space. These processes could motivate the introduction of novel terms in the dynamics, so far neglected because they are less important than those discussed in the article. Among the processes that could be anticipated, we mention diffusive terms, originally introduced in Ref. [50], that have been shown to be eventually responsible for the thermalisation of the one-dimensional Bose gas in the presence of an external confinement [51]. The rate equations completely disregard interactions between the fermionic momenta, which could be responsible for some forms of redistribution, and eventually change the decay of the gas. The physical importance of these terms motivates further studies in this direction.
We now insert the rate equations (18) obtaining: For what concerns the time evolution of nk(t), the rate equations (18) can be recasted into the following form: where in the last passage we used the previous result for ∂tñ (t). This concludes the derivation of Eqs. (19).
We want now to show that Eq. (20) is solved by Eq. (21). By integrating Eq. (20) one obtains: By means of logarithm properties and subsequent exponentiation we get: this consludes the derivation of Eq. (21).

B Homogeneous Bose gas: short-time behaviour
We now characterise the behaviour of the gas at short time using Eq. (24), which is exact. At short time ν(t) 1 and we can use the series expansion for the Erf The solution of this differential equation at short times is: By differentiating with respect tot we obtain the short-time behaviour of the population: C Homogeneous Bose gas: Analytical expressions for n(t) and ν(t) We now discuss the derivation of the values of the coefficients A and B in Eq. (29). We consider both the short-and long-time behaviour of Eq. (28): By imposing that f (t) has the same short-and long-time behaviour of n(t), we obtain the following two conditions for A and B: The values in Eq. (29) solve this system of equations; the other solution of the system produce a function f (t) that does not reproduce the numerical data.
We can also present an analytical formula for ν(t): This expression is interesting because we can use it to see whether it satisfies the differential equation (24), and thus to check whether the analytical formula that we found is exact or not. We find that this is not the case.

D The numerical implementation of the Boltzmann equation
In this appendix we present some details related to the numerical simulation of the dynamics induced by Eq. (38). We integrate it exploiting the second-order Suzuki-Trotter expansion of the differential operator governing the dynamics, namely the r.h.s. of Eq. (38): f (x,k,t + dt) = D HO dt/2 • D loss dt • D HO dt/2 f (x,k,t) + O(dt 2 ).
Here, D HO dt/2 describes the dynamics induced by the trap, i.e. first term of the r.h.s. of Eq. (38), and D loss dt/2 is responsible for the two-body losses, i.e. second term of the r.h.s. of Eq. (38). The action of the two operators D HO dt/2 on the density distribution can be implemented exactly since it is a rigid rotation of the density function in the phase spacex−k of and angle (ω/γ)dt = ωdt. The action of the the operator D loss dt is implemented with a fourth-order Runge-Kutta method with variable time-step.
We smoothened the initial condition in order to avoid numerical problems connected to the sharpness of the initial distribution. In particular we used: f (x,k, 0) = 1x 2 +k 2 ≤ 1 − ξ; e −[x 2 +k 2 −(1−ξ)] 2 /ξ otherwise, with ξ = 5d. This smoothening is visible in the plots of Fig. 2. For all the data shown in Fig. 2 we discretized the phase space with a grid of 250 × 250 pixels for the phase spacex,k ∈ [−1.2, 1.2]. This discrtetization implies a resolution d = dx = dk = 9.6 · 10 −3 . We check that the results obtained for the time-evolution of the normalized densityÑ (t) do not depend on the size of the pixels. In Fig.6 we showÑ (t) for two different values of d. The two curves collapse for the time window under consideration.

E Analytical solution for ω/Γ = 0
We solve the Eq. (38) for ω/Γ = 0. Whereas the presence of the trap is taken into account by the initial conditions, it is completely disregarded in the dynamics, so that it describes the situation where losses are extremely more important. The Boltzmann equation factorizes the dynamics at differentx: (k −q) 2 f (x,k,t)f (x,q,t)dq. where we have explicitly written that f (x, q, t) is different from zero only for x 2 + q 2 ≤ 1. By rescaling momenta and time as follows: we can recast the Boltzmann-like equation in a form that we have already encountered, see Eq. (18) : ∂ ∂T f (x, K, T ) = − π π (K 2 + Q 2 )f (x, K, T )f (x, Q, T )dQ.
Using the results in Sec. 3, we obtain the long-time limit: 1 2π π −π f (x, K, T )dK ∼ 1 4π We finally obtain the result reported in the main text: