# Effects of atom losses on a one-dimensional lattice gas of hardcore bosons

François Riggio, Lorenzo Rosso†, Dragi Karevski and Jérôme Dubail  †Université
Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France
Université de Lorraine, CNRS, LPCT, F-54000 Nancy, France
<EMAIL_ADDRESS>
###### Abstract

Atom losses occur naturally during cold atoms experiments. Since this
phenomenon is unavoidable, it is important to understand its effect on the
remaining atoms. Here we black investigate a toy model: a lattice gas of hard-
core bosons subject to $K$-body losses (where $K=1,2,3,\dots$ is the number of
atoms lost in each loss event) on $K$ neighboring sites. In particular we
investigate the effect of losses on the rapidity distribution $\rho(k)$ of the
atoms. Under the assumption that losses are weak enough so that the system
relaxes between two loss events, we are able to determine the loss functional
$F[\rho](k)$ encoding the loss process for $K$-body losses. We derive closed
expressions for the cases of one- and two-body losses, and show their effects
on the evolution of the total number of particles. Then we add a harmonic
trapping potential and study the evolution of the position-dependent rapidity
distribution of this system by solving numerically the evolution equation for
one-, two- and three-body losses.

##  I Introduction

During the past two decades, cold atom experiments have become a key
simulation platform to investigate the physics of low-dimensional quantum
many-body systems [1]. While these experiments are typically well isolated and
their dynamics at short times is well approximated by unitary dynamics, cold
atom setups are always subject to atom losses. The effects of the latter on
the dynamics of the gas can become important at longer times, and they are
typically hard to describe theoretially. Cold atom experiments involve various
loss processes, for instance one-body losses resulting from scattering with
background thermal atoms [2], which can be significant. Inelastic two-body
collisions can occur naturally or be intentionally engineered, leading to two-
body losses [3, 4, 5, 6, 7, 8, 9, 10]. Three-body losses, where a highly bound
diatomic molecule is formed, are invariably present and typically dominate the
overall loss process [11, 12, 13]. In principle, loss processes involving more
than three atoms also exist. In particular, losses involving four atoms have
been reported in Refs. [14, 15]. Loss processes can sometimes be controlled
and engineered [16] to bring out physical phenomena such as cooling [17, 18,
19, 20, 21, 22, 23, 24] and the quantum Zeno effect [25, 26, 27, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39].

Despite the relevance of loss processes, a consistent general theory
describing the effects of losses on the dynamics does not exist yet. One
standard way to model atom losses is to use the Lindblad equation [26, 40, 25,
41, 42]. Several studies have inspected the interplay between the unitary
dynamics and the lossy one in both bosonic [26, 25, 40, 43, 44] and fermionic
gases [45, 46, 47, 48, 49, 50, 51, 44].

black In one spatial dimension, several ultracold gases are integrable or
nearly integrable, including continuous gases of bosons [52, 53, 54, 55, 56]
or fermions [57, 58], sometimes with multiple components [59, 60], but also
lattice gases like hard-core bosons [61, 62] or the Fermi-Hubbard model [63].
The question of losses in those integrable gases is particularly interesting.
An integrable model admits a macroscopic number of conserved charges in
comparison with non-integrable models where, typically, only the energy and
the number of particle are conserved. Due to this large set of conserved
quantities, an isolated integrable system has a singular property: the
stationary states of the system are modeled by Generalized Gibbs Ensembles
(GGE). A GGE is constructed by maximizing the entropy under the constraints
imposed by all conservations laws [64]. However, under dissipative processes
like atom losses, the conserved quantities of an integrable system are not
conserved anymore, as the coupling between the system and its environnement
typically causes integrability breaking [65]. In principle, one then expects
that integrability breaking leads to thermalization at very long times,
possibly with a pre-thermalization phenomenon at intermediate time scales
[66]. Some studies seem to confirm the thermalization but other works suggest
otherwise [17].

black Motivated by the outstanding challenge of developping a general theory
of the effects of losses on the dynamics of integrable gases, here we
introduce a toy model where we can theoretically investigate the effects of
$K$-body losses ($K\geq 1$) on the dynamics of the simplest such gas. We
propose to investigate hard-core bosons on a lattice —a famous model of
integrable quantum gas dynamics that has been extensively studied
theoretically [67, 68, 69, 70, 71, 72] and realized experimentally [61, 62]—
subject to weak $K$-body atom losses: whenever $K$ atoms occupy $K$
neighboring sites, they can escape the system with some rate (Fig. 1).

We stress that experimentally realistic loss processes in optical lattices
usually involve atoms on the same lattice site, so, strictly speaking, only
the $K=1$ version of our toy model is of experimental relevance. For $K\geq
2$, we are not aware of existing experimental setups that could engineer our
$K$-site loss processes. However, from a theoretical perspective this model is
very natural. Here our goal is to understand how integrability breaking caused
by such loss events affects the dynamics of the gas. Of course, under atom
losses, the stationary state at very long times is always the vacuum. However,
what we want to investigate is how the vacuum is approached, and in particular
whether or not to the gas is in a quasi-stationary thermal state. We
investigate for instance the mean density $n(t)=\braket{N(t)}/L$ (where $N$ is
the number of atoms and $L$ is the number of lattice sites) for different
initial states and observe different non-trivial behaviors depending on the
number $K$ of atoms lost in each loss event. In particular, for two-body
losses ($K=2$), we analytically show that $n(t)\propto 1/t$ or $\propto
1/t^{1/2}$, depending on a specific property of the initial state (see below).
More generally, for other values of $K\geq 2$ our numerical results are in
agreement with a power-law decay $n(t)\propto t^{\alpha}$ where the exponent
$\alpha$ typically depends on the initial state, and generically differs from
the mean-field result $1/(K-1)$. This observation holds also if we add a
harmonic potential; however in that case the exponent $\alpha$ changes and is
different from the one found in the homogeneous case.

The paper is organized as follows. In Section II we define the model and
discuss our assumptions. Our main hypothesis of slow losses, and fundamental
concepts such as the rapidity distribution and the loss functional are
introduced. In Section III we present our calculation of the loss functional
for $K$-body losses, which crucially depends on the parity of $K$. In Section
IV we investigate the effect of atom losses on the rapidity distribution and
on the mean density in the homogeneous gas, by combining analytical and
numerical calculations. In Section V, we add a harmonic trapping potential. We
use a hydrodynamic-like approach similar to Generalized Hydrodynamics [73, 74]
where we incorporate the loss functional (following similar proposals in Refs.
[40, 43]), and design a numerical method for solving the resulting evolution
equation for the rapidity distribution in the gas. We conclude in Section VI.

Figure 1: We consider a lattice gas of hardcore bosons (pink dots) that can
jump from site to site. Two bosons cannot be on the same site. Losses are
represented by the large black arrow which connects the bosons chain to its
environment (blue background). Here we show the situation for one-body losses
($K=1$), two-body losses ($K=2$) and three-body losses ($K=3$). For each
cases, $K$ consecutive atoms are removed and are lost in the ‘environment’. In
this paper losses are modeled by the Lindblad equation (2).

##  II The model

###  II.1 Definition: lattice Tonks-Girardeau gas with $K$-body losses

We consider a lattice Tonks-Girardeau gas subject to atom losses. Each site
$j\in\mathbb{Z}$ is occupied by either zero or one boson. We write
$\sigma^{+}_{j}$/$\sigma^{-}_{j}$ for the operator that creates/annihilates a
boson on site $j$. Because of the hard-core constraint, these operators do not
satisfy the usual bosonic canonical commutation relations. Instead, they
satisfy the algebra of Pauli matrices,
$(\sigma^{+}_{j})^{2}=(\sigma^{-}_{j})^{2}=0$, and
$[\sigma^{+}_{i},\sigma^{-}_{j}]=\delta_{i,j}\sigma^{z}_{j}$.

We consider the hard-core boson (HCB) Hamiltonian with nearest-neighbor
hopping,

$\displaystyle H_{\rm
HCB}=-\frac{1}{2}\sum_{j\in\mathbb{Z}}(\sigma^{+}_{j}\sigma^{-}_{j+1}+\sigma^{+}_{j+1}\sigma^{-}_{j}).$
(1)

This Hamiltonian generates the unitary part of the evolution of the gas. In
addition, we assume that the gas is subject to incoherent $K$-body loss
processes, with $K$ a positive integer. To describe the losses, we assume that
the dynamics is Markovian, and we consider the following Lindblad equation for
the density matrix $\hat{\rho}$,

$\displaystyle\dot{\hat{\rho}}(t)$ $\displaystyle=$ $\displaystyle-i[H_{\rm
HCB},\hat{\rho}(t)]$ (2)
$\displaystyle+\Gamma\sum_{j\in\mathbb{Z}}\left(L_{j}\,\hat{\rho}(t)\,L^{\dagger}_{j}-\dfrac{1}{2}\\{L^{\dagger}_{j}L_{j},\hat{\rho}(t)\\}\right),$

Here $\Gamma$ is a constant that sets the loss rate, while the Lindblad
operators

$\displaystyle L_{j}=\prod^{K-1}_{l=0}\sigma^{-}_{j+l}$ (3)

remove $K$ bosons from the $K$ consecutive sites $j,j+1,\dots,j+K-1$.

We stress that, with these loss terms, the model is not exactly solvable, so
it is necessary to develop some effective approaches to tackle it.

###  II.2 Adiabatic losses, effective description by slow motion of the
charges

To simplify the description of the system, we follow the approach of Ref.
[40]. There the approach was developed for losses in a continuous gas, and it
is based on the assumption that the losses are slow, so that the gas remains
in a Generalized Gibbs Ensemble (GGE) with parameters that slowly drift in
time (see e.g. [75] for a review). A similar description of slowly-evolving
nearly integrable systems with weak integrability breaking had previously
appeared in Refs. [76, 77, 78]. Following these ideas, here we also assume
that the loss processes occur on very long times compared to the relaxation
time scale due to the unitary evolution of the gas. In that limit, the gas has
time to reach a local stationary state after each loss event. To efficiently
exploit that assumption, we look at the slow dynamics of the conserved
charges.

For convenience, from now on we focus on a finite system of $L\gg 1$ sites,
with periodic boundary conditions. The Hamiltonian $H_{\rm HCB}$ commutes with
an infinite set of hermitian operators $Q_{a}$, $a=0,1,2\dots$ that can be
constructed using the Jordan-Wigner mapping to non-interacting fermions (see
Subsec. II.3 for details). These operators also commute among themselves,
$[Q_{a},Q_{b}]=[H_{\rm HCB},Q_{a}]=0$. Moreover, they are local in the sense
that $Q_{a}=\sum_{j=1}^{L}q_{a,j}$ where $q_{a,j}$ is a charge density
operator that has compact support (i.e. it acts on a finite number of sites
around $j$).

The time evolution of the expectation value $\left<Q_{a}\right>(t)={\rm
tr}[\hat{\rho}(t)Q_{a}]$, is obtained from Eq. (2),

$\displaystyle\dot{\expectationvalue{Q_{a}}}(t)=\frac{\Gamma}{2}\sum_{j=1}^{L}\expectationvalue{L^{\dagger}_{j}\,[Q_{a},\,L_{j}]+[L^{\dagger}_{j},\,Q_{a}]L_{j}}(t).$
(4)

Moreover, the hermiticity of $Q_{a}$ implies
$\langle[L_{j}^{\dagger},Q_{a}]L_{j}\rangle=\langle
L^{\dagger}_{j}\,[Q_{a},\,L_{j}]\rangle^{*}$. Thus,

$\displaystyle\dot{\expectationvalue{Q_{a}}}(t)=\Gamma\sum_{j=1}^{L}{\rm
Re}\,\expectationvalue{L^{\dagger}_{j}[Q_{a},L_{j}]}(t).$ (5)

This equation is exact and it is a direct consequence of Eq. (2). Notice that,
in this form, it is not particularly useful, because to evaluate the r.h.s one
needs to know the exact density matrix $\hat{\rho}(t)$.

Now comes the crucial step. Importantly, the operator
$L^{\dagger}_{j}[Q_{a},L_{j}]$ is local, because both the operator $L_{j}$ and
the charge density $q_{a,j}$ have compact support. This, together with the
assumption of slow losses, allows us to use the idea of local relaxation in
the system. Namely, we expect that, under unitary evolution, the density
matrix of a small subsystem quickly relaxes to a Generalized Gibbs Ensemble
(GGE). Expectation values of local observables can then be evaluated with
respect to the GGE density matrix,

$\displaystyle\hat{\rho}_{{\rm GGE},\\{\langle
Q_{a}\rangle\\}}\,\propto\,e^{-\sum_{a}\beta_{a}Q_{a}},$ (6)

where the Lagrange multipliers $\beta_{a}$ are fixed by the expectation values
of the charges $\langle Q_{a}\rangle$, which must be equal to ${\rm
tr}[\hat{\rho}_{{\rm GGE},\\{\langle Q_{b}\rangle\\}}\,Q_{a}]$. Evaluating the
r.h.s of Eq. (4) in the GGE density matrix leads to a closed evolution
equation for the slow motion of the charges induced by the losses,

$\displaystyle\frac{d}{dt}\langle Q_{a}\rangle=\Gamma\sum_{j=1}^{L}{\rm
Re}\,\expectationvalue{L^{\dagger}_{j}[Q_{a},L_{j}]}_{{\rm GGE},\\{\langle
Q_{b}\rangle\\}}.$ (7)

It is this evolution equation that we study in great detail in this paper. For
lattice HCB, the description is further simplified by specifying the form of
the conserved charges $Q_{a}$. This is what we do next, by introducing the
distribution of rapididites.

###  II.3 Slow evolution of the rapidity distribution

Hard-core bosons can be mapped to free fermions by a Jordan-Wigner
transformation,

$\displaystyle\sigma^{+}_{j}=\prod_{i=1}^{j-1}(-1)^{c^{\dagger}_{i}c_{i}}\,c^{\dagger}_{j},\qquad\sigma^{-}_{j}=\prod_{i=1}^{j-1}(-1)^{c^{\dagger}_{i}c_{i}}\,c_{j}.$
(8)

Here the operators $c^{\dagger}_{j}/c_{j}$ create/annihilate a fermion on site
$j$. They satisfy the canonical anticommutation relations
$\\{c_{i},c^{\dagger}_{j}\\}=\delta_{ij}$. Under the Jordan-Wigner mapping,
the Hamiltonian (1) becomes

$H_{\rm
HCB}\,=\,-\frac{1}{2}\sum_{j=1}^{L}(c^{\dagger}_{j}c_{j+1}+c^{\dagger}_{j+1}c_{j}).$
(9)

Moreover, the fermions satisfy antiperiodic (resp. periodic) boundary
conditions if the number of particles $N$ in the system is even (resp. odd):

$c_{L+1}^{\dagger}=(-1)^{N-1}c_{1}^{\dagger}.$ (10)

The Fourier modes are

$c^{\dagger}(k)=\frac{1}{\sqrt{L}}\sum_{j=1}^{L}e^{ikj}c_{j}^{\dagger}$ (11)

with $k\in\frac{2\pi}{L}(\mathbb{Z}+\frac{1}{2})$ if $N$ is even, and
$k\in\frac{2\pi}{L}\mathbb{Z}$ if $N$ is odd. Either way the Hamiltonian reads

$H_{\rm HCB}\,=\,\sum_{k}\varepsilon(k)c^{\dagger}(k)c(k),$ (12)

with $\varepsilon(k)=-\cos k$.

It is clear from the form (12) that any operator of the form

$Q[f]=\sum_{k}f(k)c^{\dagger}(k)c(k),$ (13)

for any function $f(k)$, commutes with the Hamiltonian $H_{\rm HCB}$. Moreover
these conserved charges also commute among themselves. Convenient choices for
$f(k)$ are $\cos(nk)$ or $\sin(nk)$ for $n\in\mathbb{N}$, which leads to a
hermitian basis set of charges, where each charge has a charge density that is
compactly supported.

However, for the purposes of this paper, rather than to work with a specific
choice of basis for the space of conserved charges $Q_{a}$ (or $Q[f]$), it is
more convenient to work directly with the occupation number, or ‘rapidity
distribution’,

$\rho(k)\,\underset{L\rightarrow\infty}{=}\,\left<c^{\dagger}(k)c(k)\right>,\quad\rho(k)\,\in\,[0,1].$
(14)

It is clear that if we know the rapidity distribution $\rho(k)$, then we also
know the expectation values of any charge $Q[f]$, because
$\left<Q[f]\right>=\sum_{k}f(k)\rho(k)$.

Following Ref. [40], we can turn the evolution equation for the slow motion of
the charges (4) into an equation for the slow evolution of the rapidity
distribution itself,

$\displaystyle\dot{\rho}(k)=-\Gamma F[\rho](k),$ (15)

where the loss functional

$\displaystyle F[\rho](k)=\sum_{j=1}^{L}{\rm
Re}\,\expectationvalue{L^{\dagger}_{j}[L_{j},c^{\dagger}(k)c(k)]}_{{\rm
GGE},\rho},$ (16)

and the GGE density matrix itself is parameterized by the rapidity
distribution $\rho(k)$. More precisely, the GGE density matrix is Gaussian for
the fermions $c_{j}^{\dagger},c_{j}$, and it is characterized by its two-point
function $\expectationvalue{c^{\dagger}(k)c(k^{\prime})}_{{\rm
GGE},\rho}=\rho(k)\delta_{k,k^{\prime}}$. All higher-order correlations can be
computed using Wick’s theorem for fermionic operators.

The functional (16) is the central object of this paper. In the next section
we compute it explicitly for K-body loss processes. For one- and two-boson
loss processes we get simple closed expressions. For loss events involving
larger numbers $K$ of bosons, we will see that we can express the loss
functional as a small determinant, which follows from applying Wick’s theorem
to Eq. (16).

##  III Deriving the loss functional

In this section we compute the functional $F[\rho](k)$ explicitly.
Importantly, in our calculation we uncover a different structure depending on
the parity of the number $K$ of bosons lost in each loss event, which can be
traced back to the Jordan-Wigner string appearing in the mapping (8) to non-
interacting fermions.

###  III.1 One-body losses

For $K=1$ the Lindblad dissipators are $L_{j}=\sigma^{-}_{j}$. Using
translational invariance, the loss functional (16) that we need to compute is

$\displaystyle
F[\rho](k)\,=\,L\expectationvalue{\sigma^{+}_{1}[\sigma^{-}_{1},c^{\dagger}(k)c(k)]}_{{\rm
GGE},\rho}$
$\displaystyle=\,L\expectationvalue{\sigma^{+}_{1}\sigma^{-}_{1}c^{\dagger}(k)c(k)}_{{\rm
GGE},\rho}-L\expectationvalue{\sigma^{+}_{1}c^{\dagger}(k)c(k)\sigma^{-}_{1}}_{{\rm
GGE},\rho},$

where $L$ is the length of the system and the loss operator acts on the site
$j=1$. Both terms in the second line of (III.1) can be computed using the fact
that the GGE is a gaussian state for the fermions, which allows us to use
Wick’s theorem. For the first term we have (using
$c_{1}=\frac{1}{\sqrt{L}}\sum_{q}e^{iq}\,c(q)$):

$\displaystyle
L\expectationvalue{\sigma^{+}_{1}\sigma^{-}_{1}c^{\dagger}(k)c(k)}_{{\rm
GGE},\rho}\,=\,L\expectationvalue{c^{\dagger}_{1}c_{1}c^{\dagger}(k)c(k)}_{{\rm
GGE},\rho}$
$\displaystyle=\,\sum_{qq^{\prime}}e^{i(q-q^{\prime})}\expectationvalue{c^{\dagger}(q^{\prime})c(q)c^{\dagger}(k)c(k)}_{{\rm
GGE},\rho}$
$\displaystyle=\,\sum_{qq^{\prime}}e^{i(q-q^{\prime})}\left(\langle
c^{\dagger}(q^{\prime})c(q)\rangle\langle c^{\dagger}(k)c(k)\rangle+\right.$
$\displaystyle\qquad\qquad\left.\langle
c^{\dagger}(q^{\prime})c(k)\rangle\langle c(q)c^{\dagger}(k)\rangle\right)$
$\displaystyle=\,\langle N\rangle\rho(k)+\rho(k)(1-\rho(k)).$ (18)

The second term requires more care, because the operator $c^{\dagger}(k)c(k)$
is inserted between $\sigma_{1}^{+}$ and $\sigma_{1}^{-}$, and the latter
change the parity of the number of particles in the system. The boundary
conditions for the fermions are modified according to Eq. (10). Thus we need
to relate the Fourier modes of the fermions with periodic boundary conditions
to the ones with anti-periodic boundary conditions. For conciseness, let us
introduce the two corresponding sets of momenta,

$\displaystyle Q^{\rm p}=\frac{2\pi}{L}\times\\{1,2,\dots,L\\},$ (19)
$\displaystyle Q^{\rm
ap}=\frac{2\pi}{L}\times\\{\frac{1}{2},\frac{3}{2},\dots,L-\frac{1}{2}\\}.$
(20)

Then we have the following identities,

$\displaystyle(k\in Q^{\rm p})\qquad c(k)=\dfrac{i}{L}\sum_{q\in Q^{\rm
ap}}\dfrac{e^{i(q-k)/2}}{\sin((q-k)/2)}\,c(q),$ (21) $\displaystyle\qquad(k\in
Q^{\rm ap})\qquad c(k)=\dfrac{i}{L}\sum_{q\in Q^{\rm
p}}\dfrac{e^{i(q-k)/2}}{\sin((q-k)/2)}\,c(q).$ (22)

We can insert them into the second term of Eq. (III.1), which leads to

$\displaystyle\expectationvalue{\sigma^{+}_{1}c^{\dagger}(k)c(k)\sigma^{-}_{1}}=\dfrac{1}{L^{2}}\sum_{q,q^{\prime}}\dfrac{e^{i(q-q^{\prime})/2}\,\expectationvalue{c^{\dagger}_{1}c^{\dagger}(q^{\prime})\,c(q)c_{1}}}{\sin((q-k)/2)\sin((q^{\prime}-k)/2)}.$
(23)

This correctly implements the change of boundaries of the fermions. Next we
can apply Wick’s theorem to evaluate the four-fermion correlator
$\expectationvalue{c^{\dagger}_{1}c^{\dagger}(q^{\prime})\,c(q)c_{1}}$. This
leads to

$\displaystyle
L\expectationvalue{\sigma^{+}_{1}c^{\dagger}(k)c(k)\sigma^{-}_{1}}=$
$\displaystyle\frac{\expectationvalue{N}}{L^{2}}\sum_{q}\dfrac{\rho(q)}{\sin^{2}(\frac{q-k}{2})}-\left(\frac{1}{L}\sum_{q}\cot(\frac{q-k}{2})\,\rho(q)\right)^{2}-\dfrac{\expectationvalue{N}^{2}}{L^{2}}.$
(24)

The first term in the above equation has a pole of order 2 at
$q=k+2\pi\mathbb{Z}$, and it is convenient to reduce its degree using the
identity $\sum_{q}1/\sin^{2}(\frac{q-k}{2})=L^{2}$. This leads to the
equivalent expression

$\displaystyle
L\expectationvalue{\sigma^{+}_{1}c^{\dagger}(k)c(k)\sigma^{-}_{1}}$
$\displaystyle=\frac{\expectationvalue{N}}{L^{2}}\sum_{q}\dfrac{\rho(q)-\rho(k)}{\sin^{2}(\frac{q-k}{2})}+\expectationvalue{N}\rho(k)$
$\displaystyle-\left(\frac{1}{L}\sum_{p}\cot(\frac{p-k}{2})\,\rho(p)\right)^{2}-\dfrac{\expectationvalue{N}^{2}}{L^{2}}.$
(25)

Putting the two terms (III.1)-(III.1) together and taking the thermodynamic
limit $L\rightarrow\infty$, we arrive at the following form of the one-body
loss functional

$\displaystyle F[\rho](k)$
$\displaystyle=\rho(k)-\rho^{2}(k)+\left(\fint^{\pi}_{-\pi}\frac{dp}{2\pi}\,\cot(\frac{k-p}{2})\,\rho(p)\right)^{2}$
$\displaystyle+n\left(n+\fint^{\pi}_{-\pi}\frac{dq}{2\pi}\,\dfrac{\rho(k)-\rho(q)}{\sin^{2}(\frac{k-q}{2})}\right),$
(26)

where $n=\expectationvalue{N}/L$ is the density of particle and $\fint$ means
the Cauchy principal value of the integral. This is our main result for one-
body losses. It is similar to (but different from) the formula given in Ref.
[40] for the Tonks-Girardeau gas in the continuum. Notice that the functional
is non-linear in $\rho(k)$, and also non-local in rapidity space.

###  III.2 Two-body losses

For $K=2$, the dissipators are $L_{j}=\sigma^{-}_{j}\sigma^{-}_{j+1}$. Under
the Jordan-Wigner mapping they become

$\displaystyle
L_{j}=\sigma^{-}_{j}\sigma^{-}_{j+1}=c_{j}(-1)^{c^{\dagger}_{j}c_{j}}c_{j+1}=-c_{j}c_{j+1}.$
(27)

Then, to compute the functional $F$, we simply need to insert the dissipator
(27) in the definition (16),

$\displaystyle
F[\rho](k)=\sum_{j}\expectationvalue{\sigma^{+}_{j+1}\sigma^{+}_{j}[\sigma^{-}_{j}\sigma^{-}_{j+1},c^{\dagger}(k)c(k)]}$
$\displaystyle=\sum_{q,q^{\prime},p,p^{\prime}}\frac{e^{i(2p^{\prime}+p-2q-q^{\prime})}}{M}\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})[c(p)c(p^{\prime}),c^{\dagger}(k)c(k)]}.$
(28)

Expanding the commutator in the braket leads to two terms

$\displaystyle\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})[c(p)c(p^{\prime}),c^{\dagger}(k)c(k)]}$
$\displaystyle=\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})c(p)c(p^{\prime})c^{\dagger}(k)c(k)}$
$\displaystyle-\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})c^{\dagger}(k)c(k)c(p)c(p^{\prime})}.$
(29)

The first term can be expressed as

$\displaystyle\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})c(p)c(p^{\prime})c^{\dagger}(k)c(k)}$
$\displaystyle=\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})c^{\dagger}(k)c(k)c(p)c(p^{\prime})}$
$\displaystyle+\delta_{p^{\prime}k}\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})c(p)c(k)}-\delta_{pk}\expectationvalue{c^{\dagger}(q)c^{\dagger}(q^{\prime})c(p^{\prime})c(k)}.$
(30)

The second term is the expectation value of $c^{\dagger}(k)c(k)$ in a state
where two atoms have been removed, the parity of the initial number of
particle is then unchanged. Hence, in contrast to the $K=1$ case treated in
the previous subsection, here the parity of the number of atoms does not
change.

One can then apply Wick’s theorem, and take the thermodynamic limit to obtain
the loss functional,

$\displaystyle F[\rho](k)$
$\displaystyle=\dfrac{2}{\pi}\int^{\pi}_{-\pi}dq\,\sin^{2}\left(\frac{k-q}{2}\right)\;\rho(q)\,\rho(k).$
(31)

This is our main result for two-body losses. That functional presents some
similarities with the loss functional found in the relation (5) of Ref. [25].

We now generalise this calculation to the case of $K$-body losses for $K$ an
arbitrary even integer.

###  III.3 $K$-body losses with K even

In this subsection we show that is possible to find a closed formula for the
loss functional defined in (16) where the Lindblad operator is given by
$L_{j}=\sigma_{j}\sigma_{j+1}\dots\sigma_{j+K-1}$. Taking the Fourier
transform of $L_{j}$ and $L^{\dagger}_{j}$ in (16), the loss functional reads

$\displaystyle F^{\rm even}[\rho](k)=$
$\displaystyle\frac{1}{L^{K-1}}\sum_{\begin{subarray}{c}q_{1},\dots,q_{K}\\\
q^{\prime}_{1},\dots,q^{\prime}_{K}\end{subarray}}\exp{i\sum^{K}_{l=1}(q_{l}-q^{\prime}_{l})l}\times$
$\displaystyle\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})[c(q_{1})\dots c(q_{K}),c^{\dagger}(k)c(k)]}.$ (32)

As mentionned in III.2, in the case of even $K$-body losses the commutator in
(III.3) reduces to $K$ terms

$\displaystyle\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})[c(q_{1})\dots c(q_{K}),c^{\dagger}(k)c(k)]}$
$\displaystyle=\delta_{kq_{K}}\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c(q_{1})\dots c(q_{K-1})c(k)}$
$\displaystyle-\delta_{kq_{(K-1)}}\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c(q_{1})\dots c(q_{K-2})c(q_{K})c(k)}$
$\displaystyle+\dots$ (33)

and applying Wick’s theorem on each terms leads to a product of $K$ terms of
the form $\expectationvalue{c^{\dagger}(q^{\prime})c(q)}$. Using the property
$\expectationvalue{c^{\dagger}(q^{\prime})c(q)}=\delta_{qq^{\prime}}\,\rho(q)$
and taking the thermodynamic limit, the loss functional in (III.3) can be
expressed as a sum of $K$ terms each consisting of a $K$ by $K$ matrix
determinant. Let us introduce the $K\times K$ matrix $A^{(j)}_{[\rho]}$ with
matrix elements

$\displaystyle[A^{(j)}_{[\rho]}]_{ab}=\left\\{\begin{array}[]{ll}\frac{1}{2\pi}\int^{\pi}_{-\pi}dq\,e^{i(b-a)q}\,\rho(q)\quad{\rm
if}\;b\neq j\\\ e^{i(b-a)k}\,\rho(k)\quad{\rm if}\;b=j,\end{array}\right.$
(36)

for indices $a,b=1,\dots,K$. The superscript $j$ indicates which column
depends on the rapidity $k$. Apart from the $j^{\rm th}$ column, the matrix
essentially contains Fourier transforms of the rapidity distribution
$\rho(k)$. The loss functional then reduces to

$\displaystyle F^{\rm even}[\rho](k)=\sum^{K}_{j=1}\det(A^{(j)}_{[\rho]}).$
(37)

The relation (37) is another fundamental result of this paper.

Figure 2:  black Effect of one-body losses on different rapidity
distributions. The initial distribution is a thermal (Fermi-Dirac)
distribution $\rho(k)=(1+\exp(-\cos(k)/T))^{-1}$, and we show the results for
different initial temperatures $T=0.1,1.0,10.0,100.0$. The time evolution is
obtained by solving Eq. (15) with the loss functional (III.1) with the Runge-
Kutta method. Notice that the solution of Eq. (15) depends on the loss rate
$\Gamma$ only through the dimensionless time $\Gamma t$. The black dashed
lines correspond to formula (54) and show a perfect agreement with the
numerical solution of the evolution equation.

###  III.4 $K$-body losses with K odd

In this last subsection we investigate the case of losses process with an odd
number $K$ of lost atoms. The reasoning is similar to the one we developped in
the previous subsection, however like in the $K=1$ case, we need to be careful
about the change of boundary conditions for the fermions. We start by taking
the Fourier transform of the Lindblad operators in the definition (16), which
leads to the relation (III.3). The commutator in the loss functional gives two
terms

$\displaystyle\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})[c(q_{1})\dots c(q_{K}),c^{\dagger}(k)c(k)]}$
$\displaystyle=\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c(q_{1})\dots c(q_{K})c^{\dagger}(k)c(k)}$
$\displaystyle-\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c^{\dagger}(k)c(k)c(q_{1})\dots c(q_{K})}.$ (38)

As we have already discussed in the subsection III.1, the first term is the
expectation value of $c^{\dagger}(k)c(k)$ in a state where the initial number
of particles is preserved. However the second term corresponds to the
expectation value of $c^{\dagger}(k)c(k)$ in a state where $K$ atoms have been
removed. Since $K$ is an odd number, the parity of the number of particle is
changed and one needs to use the relations (21) to express
$c^{\dagger}(k)c(k)$ in the appropriate parity sector. Inserting the relations
(21) in the second term, one has

$\displaystyle\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c^{\dagger}(k)c(k)c(q_{1})\dots c(q_{K})}=$
$\displaystyle\sum_{q,q^{\prime}}\dfrac{e^{i(q-q^{\prime})/2}\,\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c^{\dagger}(q^{\prime})c(q)c(q_{1})\dots
c(q_{K})}}{L^{2}\sin((q-k)/2)\sin((q^{\prime}-k)/2)}.$ (39)

Before using Wick’s theorem on the above formula, one can notice that the
first term in the right-hand side of (III.4) can be written as

$\displaystyle\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c(q_{1})\dots c(q_{K})c^{\dagger}(k)c(k)}$
$\displaystyle=\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c^{\dagger}(k)c(k)c(q_{1})\dots c(q_{K})}$
$\displaystyle+\delta_{kq_{K}}\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c(q_{1})\dots c(q_{K-1})c(k)}$
$\displaystyle-\delta_{kq_{(K-1)}}\expectationvalue{c^{\dagger}(q^{\prime}_{K})\dots
c^{\dagger}(q^{\prime}_{1})c(q_{1})\dots c(q_{K-2})c(q_{K})c(k)}$
$\displaystyle+\dots,$ (40)

where we used the anti-commutation relation for the fermionic operators. As we
proceed in the previous subsection, the Wick’s contractions of (III.4) and of
the first term in the right-hand side of (III.4) can be written as two
determinants of two matrices $B$ and $C$. The matrices $B$ and $C$ are
$(K+1)\times(K+1)$ hermitian matrices and their matrix elements depend on the
Fourier and Hilbert transforms [79] of $\rho(k)$

$\displaystyle[B_{[\rho]}]_{ab}=\left\\{\begin{array}[]{ll}\frac{1}{2\pi}\int^{\pi}_{-\pi}dq\,e^{i(b-a)q}\,\rho(q)\;\;{\rm
if}\;a,\,b<K+1\\\ e^{-iak}\,\rho(k)\;\;{\rm if}\;b=K+1\\\ 0\;{\rm
if}\;a=b=K+1,\end{array}\right.$ (44)
$\displaystyle[C_{[\rho]}]_{ab}=\left\\{\begin{array}[]{ll}\frac{1}{2\pi}\int^{\pi}_{-\pi}dq\,e^{i(b-a)q}\,\rho(q)\;{\rm
if}\;a,\,b<K+1\\\
\frac{1}{2\pi}\fint^{\pi}_{-\pi}dq\,e^{-i(a-1)q}\rho(q)(\cot(\dfrac{k-q}{2})+i)\\\
{\rm if}\;b=K+1\\\
\frac{1}{2\pi}\fint^{\pi}_{-\pi}dq\,\dfrac{\rho(q)-\rho(k)}{\sin^{2}(\dfrac{k-q}{2})}\;{\rm
if}\;a=b=K+1.\end{array}\right.$ (49)

The loss functional for odd $K$ takes the final form

$\displaystyle F^{\rm odd}_{K}[\rho](k)$
$\displaystyle=\left(\sum^{K}_{j=1}\det(A^{(j)}_{[\rho]})\right)$
$\displaystyle+\left[\det(B_{[\rho]})-\det(C_{[\rho]})\right].$ (50)

It is possible to write a general expression valid both for even and odd $K$
by introducing the factor $\dfrac{1-(-1)^{K}}{2}$ which vanishes for $K$ even,
so that our final result, valid in all cases, reads

$\displaystyle F_{K}[\rho](k)$
$\displaystyle=\left(\sum^{K}_{j=1}\det(A^{(j)}_{[\rho]})\right)$
$\displaystyle+\dfrac{1-(-1)^{K}}{2}\left[\det(B_{[\rho]})-\det(C_{[\rho]})\right].$
(51)

##  IV Evolution of the rapidity distribution in a homogeneous gas

Having established the general form of the loss functional $F[\rho]$ for
$K$-body losses, Eqs. (III.1)-(31)-(37)-(III.4)-(III.4), we now turn to the
time evolution of the rapidity distribution. We solve the evolution equation

$\displaystyle\dot{\rho}(k)=-\Gamma F[\rho](k)$ (52)

numerically (and analytically for the special case $K=1$), and we focus in
particular on the time evolution of the atom density
$n=\int_{-\pi}^{\pi}\rho(k)\frac{dk}{2\pi}$.

###  IV.1 Results for $K=1$

For one-body losses the atom density always decays exponentially,

$n(t)=n(0)e^{-\Gamma t}.$ (53)

This simply follows from Eq. (4) applied to the total particle number
$Q_{a}=N$ and to $L_{j}=\sigma^{-}_{j}$: it gives
$\langle\dot{N}\rangle=-\Gamma\langle N\rangle$, which implies Eq. (53) for
the atom density $n=N/L$.

It turns out the evolution equation (52) for the rapidity distribution can be
solved exactly for the loss functional for $K=1$ (see Eq. (III.1)). The
solution is derived in Appendix A; it reads

$\displaystyle\rho(t,k)=n_{0}e^{-\Gamma
t}\real\left(\dfrac{\tanh(n_{0}(e^{-\Gamma
t}-1))+\frac{i}{n_{0}}\,I(t,k)}{1+\frac{i}{n_{0}}\tanh(n_{0}(e^{-\Gamma
t}-1))\;I(t,k)}\right),$ (54)

where $n_{0}=n(0)$ is the initial atom density, $\rho_{0}(k)=\rho(t=0,k)$ is
the initial rapidity distribution, and $I(t,k)$ is the integral

$\displaystyle
I(t,k)=\int^{\pi}_{-\pi}\frac{dq}{2\pi}\dfrac{\rho_{0}(q)}{\tan(\frac{k-q}{2}+in_{0}(1-e^{-\Gamma
t}))}.$ (55)

A related expression for the continuous Tonks-Girardeau gas with one-body
losses was obtained in Ref. [40].

In Fig. (2) we show the evolution of the rapidity distribution from thermal
initial states at different temperatures. We display the analytical result
(54), as well as the numerical solution of Eq. (52) obtained from the Runge-
Kutta method; they are in perfect agreement, as they should.

We see that loss processes spread the distribution in rapidity space. In the
limit of large temperature for the initial state, the rapidity distribution is
flat, and remains flat at all times. For smaller initial temperatures, it
evolves into a bell-shape distribution at late times, which is close to a
Boltzmann distribution $\rho(k)\propto\exp[\cos(k)/T]$ with a density going to
zero according to Eq. (53). This is further illustrated in Fig. 3.(a), where
we plot the ratio $\rho(k,t)/n(t)$ in the limit $t\rightarrow\infty$, and fit
the result with a Boltzmann distribution. The agreement is very good, even
though it is clear from the exact formulas (54)-(55) that the distribution
$\rho(k,t)$ never becomes exactly thermal, even at infinite time.

To find a more striking signature of the fact that the system never goes to a
low-density thermal distribution, we consider the case of an oscillating
initial rapidity distribution $\rho_{0}(k)=(1-\cos(sk))/2$, where $s$ is an
integer. In that case the long-time limit of the integral (55) can be
evaluated analytically,
$\lim\limits_{t\rightarrow\infty}I(t,k)=(1-e^{-2sn_{0}}e^{iks})/2$, and when
injected in Eq. (54) it leads to a late-time rapidity distribution of the form

$\displaystyle\frac{\rho(k,t)}{n(t)}\,\underset{t\rightarrow\infty}{=}\,\frac{\alpha+\beta\,\cos(sk)}{\gamma+\delta\,\cos(sk)},$
(56)

where the coefficients $\alpha,\beta,\gamma,\delta$ depend on the initial
density $n_{0}$ and on the integer $s$. Thus, even at long time, the rescaled
rapidity distribution is sensitive to the structure of the inital distribution
(see also Fig.3.(b)). We conclude that in general the rapidity distribution
does not go to a low-density thermal distribution at long times.

Figure 3: black Long time behavior of different rapidity distributions
rescaled by the corresponding density. a) The initial distribution (blue
curve) is a Fermi-Dirac distributions $\rho_{0}(k)=(1+\exp(-\cos(k)/T))^{-1}$
with $T$=0.1. The red curve is the distribution at long time according to Eq.
(54). The black dashed line is a numerical fit with a Boltzmann distribution
$a_{1}\,\exp(\cos(k)/a_{2})$. b) The initial distribution (blue curve) is a
cosine function $\rho_{0}(k)=(1-\cos(sk))/2$. The red curve is the
distribution for $t\rightarrow\infty$. Here it is clear that the distribution
does not become thermal. Figure 4: black Effect of two-body losses on the
rapidity distribution. The initial rapidity distribution
$\rho(k,0)=1/(1+\exp(-\cos(k)/T))$ is chosen with $T=0.1,1,10,100$ from left
to right. The black dashed curve is the analytic solution (see Eq. (58)).
After several loss events the rapidity distribution takes the form of a
gaussian centered at $k=0$. As in the one-body case, if initially the rapidity
distribution is flat then it remains flat under lossy evolution.

###  IV.2 Results for $K=2$

We now consider the time evolution equation for the rapidity distribution for
the $K=2$ case, which is characterized by the functional given by Eq. (31).
For simplicity, here we focus on initial rapidity distributions $\rho_{0}(k)$
that are symmetric under reflection $k\rightarrow-k$. Since the master
equation is also invariant under $k\rightarrow-k$, this property is conserved
throughout the entire evolution. Then Eq. (31) simplifies to the following
expression,

$F[\rho]=2\,\rho(k)\,n(t)-\dfrac{1}{\pi}\cos{(k)}\rho(k)\int^{\pi}_{-\pi}dq\,\cos{(q)}\;\rho(q).$
(57)

Eq. (57) highlights the two distinct contributions to the time evolution of
$\rho(k,t)$: the first term in the right-hand side represents a mean-field
contribution, as it does not introduce any structure in rapidity space, while
the second term is responsible for generating quantum correlations and,
consequently, introducing structure in $k-$space.

After some algebra presented in App. B, one can derive an exact (implicit)
expression for the rapidity distribution at all times:

$\displaystyle\rho(k,t)=\rho_{0}(k)$
$\displaystyle\times\exp{-2\Gamma\int_{0}^{t}\left(1-\sigma_{0}\cos{(k)}\sqrt{1+\frac{\partial_{\tau}n(\tau)}{2\Gamma\,n(\tau)^{2}}}\right)n(\tau)d\tau},$
(58)

where $\sigma_{0}=\text{sgn}(\int_{-\pi}^{\pi}\cos{(k)}\rho_{0}(k)dk)$, with
$\text{sgn}(x)=\pm 1$ the sign function. Eq. (58) shows that $\rho(k,t)$ is
entirely determined by $\rho_{0}(k)$ and $n(t)$. Notably, it reveals that
rapidities are distributed according to a cosine law, with the $k=0$ mode that
has the longest lifetime.

In Fig. 4 we show the evolution of the rapidity distribution for initial
thermal states at different temperatures, and in Fig. 5 we show the
corresponding evolution of the mean atom density $n(t)$. It appears that,
except for an initial infinite temperature state, the mean density always
decays as $n(t)\propto 1/\sqrt{t}$ at very long times, while the density
decreases as $1/t$ for an infinite temperature. These two behaviors follow
from Eqs. (57)-(58), as we now explain.

Eq. (57) reveals that initial rapidity distributions that have a vanishing
first Fourier mode $\int_{-\pi}^{\pi}dq\,\cos(q)\rho(q)$ always follow the
exact same dynamics as the mean density, namely

$n(t)=\dfrac{1}{1+2n(0)\Gamma t},$ (59)

characterized by a long-term decay as $\sim 1/t$. This power law is then
always found for initial rapidity distributions with vanishing first Fourier
mode, including infinite temperature states.

Figure 5: black The mean density under two-body losses for different thermal
initial rapidity distributions $\rho(k,0)=1/(1+\exp(-\cos(k)/T))$. Colored
curves are obtained by solving numerically the time evolution equation (15)
for the loss functional (31) with the Runge-Kutta method. The red dashed line
is the mean density associated to an initial rapidity distribution which is
flat (i.e. infinite temperature), see Eq. (66). The dashed line is the
expected long-time behavior $\sim t^{-1/2}$, see Eq.(64). Figure 6: blackLong
time behavior of different rapidity distributions rescaled by the
corresponding density $n(t)$ under two-body losses. In a), b) and c), the blue
curve represents the initial rapidity distribution while the red curve is the
rapidity distribution at long time. The violet curve is the rescaled rapidity
distribution at an intermediate time. a) The initial distribution is a Fermi-
Dirac distributions $\rho_{0}(k)=(1+\exp(-\cos(k)/T))^{-1}$ with $T$=0.1. The
green dashed curve is a fit with a Boltzmann distribution
$a_{1}\,\exp(a_{2}\cos(k))$ where $a_{2}$ is positive. b) The initial rapidity
distribution is $(1-\cos(k))/2$ which has a non-vanishing first Fourier mode.
The green dashed curve is a fit with a Boltzmann distribution
$a_{1}\,\exp{a_{2}\cos(k)}$ where $a_{2}$ is negative. c) The initial
rapididty distribution is $(1-\cos(2k))/2$ which has no first Fourier mode.
Figure 7: black Evolution of the rapidity distribution under three-body
losses. We numerically solve (15) for the loss functional $F[\rho](k)$
corresponding to $K=3$, using the Runge-Kutta method. The initial rapidity
distribution is the Fermi-Dirac distribution
$\rho(k,0)=1/(1+\exp(-\cos(k)/T))$ at half filling and for different
temperatures $T=0.10,1.0,100.0$ (from left to right).

In contrast, initial rapidity distributions that have a non-vanishing first
Fourier mode decay as $\sim 1/\sqrt{t}$ at long times. This can be understood
by looking at the long time limit of Eq. (58). Let us introduce the two time-
dependent functions $g(t)=\int_{0}^{t}n(\tau)\,d\tau$ and
$f(t)=\int_{0}^{t}\sqrt{1+\frac{\partial_{\tau}n(\tau)}{2\Gamma\,n(\tau)^{2}}}n(\tau)d\tau$.
Numerically we observe that $|\partial_{\tau}n(\tau)|\ll 2\Gamma n(\tau)^{2}$
at long times, as soon as the initial rapidity distributions has a non-
vaishing first Fourier mode. Then, expanding at first order in
$|\partial_{\tau}n(\tau)|/(2\Gamma n(\tau)^{2})$, $f(t)$ becomes

$\displaystyle f(t)\simeq g(t)+\frac{1}{4\Gamma}\ln(\dfrac{n(t)}{n(0)}),$ (60)

implying that, at large $t$, the difference between $f(t)$ and $g(t)$ grows as
$\ln(t)$. Integrating Eq. (58) over $k$ leads to the mean density

$\displaystyle n(t)=\frac{e^{-2\Gamma
g(t)}}{2\pi}\int^{\pi}_{-\pi}dk\,\rho_{0}(k)\,e^{2\Gamma\sigma_{0}\cos{(k)}f(t)}.$
(61)

Since the function $f(t)$ diverges at large $t$, the latter integral can be
evaluated by the saddle-point approximation; we denote $k^{*}_{\sigma_{0}}$
the saddle point: $k^{*}_{\sigma_{0}}=0$ if $\sigma_{0}=+1$ and
$k^{*}_{\sigma_{0}}=\pi$ if $\sigma_{0}=-1$. We are thus left with:

$\displaystyle\int^{\pi}_{-\pi}dk\,\rho_{0}(k)\,$ $\displaystyle
e^{2\Gamma\sigma_{0}\cos{(k)}f(t)}$
$\displaystyle\underset{t\rightarrow\infty}{\simeq}e^{2\Gamma\,f(t)}\int^{\infty}_{-\infty}dk\,\rho_{0}(k)\,e^{-\Gamma\,k^{2}f(t)}$
$\displaystyle=\sqrt{\frac{\pi}{\Gamma
f(t)}}\,\rho_{0}(k^{*}_{\sigma_{0}})\,e^{2\Gamma\,f(t)}.$ (62)

Using Eq. (60), we find

$\displaystyle n(t)\underset{t\rightarrow\infty}{\simeq}\frac{1}{\Gamma
f(t)}\frac{\rho_{0}(k^{*}_{\sigma_{0}})^{2}}{4\pi n(0)}\simeq\frac{1}{\Gamma
g(t)}\frac{\rho_{0}(k^{*}_{\sigma_{0}})^{2}}{4\pi n(0)},$ (63)

where, in the second identity, we have used the fact that the logarithmic term
in Eq. (60) is subleading. Since $n(t)=\partial_{t}g(t)$, we arrive at an
ordinary differential equation of the form $\partial_{t}g(t)\propto 1/g(t)$.
Consequently, $g(t)\propto\sqrt{t}$, and then

$n(t)\propto t^{-1/2},$ (64)

as expected from our numerical results, see Fig. 5. We note that a similar
result was found recently in a lattice gas with similar but different two-body
loss term [25] as well as in its continuous analog [43], although we stress
that the loss functionals and rate equations for these models are different
from the ones of this paper.

We conclude this subsection with an investigation of the long time behavior of
the rapidity distribution $\rho(k,t)$, which is determined by the mean density
$n(t)$ according to Eq. (58). We have just we etablished that the first
Fourier mode of the initial rapidity distribution strongly influences the long
time behavior. In the case of a vanishing first Fourier mode, the rapidity
density at time $t$ is simply given (see Eqs. (59)-(58)) by
$\rho(k,t)=\rho_{0}(k)\,n(t)/n(0)$. The ratio $\rho(k,t)/n(t)$ is then time-
independent, as illustrated in Fig. 6(c). In contrast, when the first Fourier
mode of the initial distribution $\rho_{0}(k)$ is non-zero, the ratio
$\rho(k,t)/n(t)$ loses its dependence on the initial rapidity distribution at
very long times. Indeed, in that case the rapidity distribution goes to a low-
density, low-temperature, Boltzmann distribution of the form
$\rho(k,t)/n(t)\simeq e^{\beta(t)\cos k}$ with the effective inverse
temperature $\beta(t)=2\sigma_{0}\Gamma f(t)$. This is illustrated for the
case of an initial thermal rapidity distribution in Fig. 6(a), where we see
that the ration $\rho(k,t)/n(t)$ gets concentrated around $k=0$ and is very
close to a Boltzmann distribution. Notice also that the effective temperature
$\beta(t)$ is negative when the sign of the first Fourier mode of the initial
rapidity distribution is negative. This is illustrated in Fig. 6(b), where we
display the ratio $\rho(k,t)/n(t)$ at late time for the far-from-thermal
initial rapidity distribution $\rho_{0}(k)=(1-\cos(k))/2$. We observe that, at
late times, the distribution gets concentrated around $k=\pi$ and corresponds
to a Boltzmann distribution at negative temperature.

Remarkably, these observations are in stark contrast with our findings for the
$K=1$ case. While we found that, for $K=1$, the rapidity distribution never
goes to a thermal distribution at late time, here for $K=2$ the distribution
goes to a low-density, low-temperature (possibly negative), thermal
distribution. This is always true, except in the special case where the first
Fourier mode of the rapidity distribution vanishes; in that case the
rapididity distribution is simply rescaled by a factor $n(t)/n(0)$ under lossy
evolution.

Figure 8: black Time evolution of the boson density under three-body losses
for different initial rapidity distributions at half filling. Colored curves
are obtained by solving numerically the time evolution equation of $\rho(k)$
using the Runge-Kutta method with an non-regular time step. From blue to red,
the simulation is performed with an initial distribution which is a Fermi-
Dirac distribution $\rho(k,0)=(1+\exp(-\cos(k)/T))^{-1}$. The dark and light
green curves are respectively obtained from initial rapidity distributions
$\rho(k,0)=(1-\cos(k))/2$ and $\rho(k,0)=(1-\cos(2k))/2$. The three dashed
lines are a guide to the eye showing power-law decay for three different
exponents. Figure 9: blackLong time behavior of different rapidity
distributions rescaled by the corresponding density under three-body losses.
For the subplots a, b and c, the blue curve represents the initial rapidity
distribution while the red curve is the rapidity distribution at long time. a)
The initial distribution is a Fermi-Dirac distributions
$\rho_{0}(k)=1/(1+\exp(-\cos(k)/T))$ with $T$=0.1. b) Non-thermal initial
rapidity distribution is $(1-\cos(k))/2$. c) Non-thermal initial rapididty
distribution is $(1-\cos(2k))/2$.

###  IV.3 Generic observations for arbitrary $K$

We now turn to the case of higher $K$, and draw some general conclusions.

Numerically, we solve the time evolution equation of the rapidity distribution
for three-body losses ($K=3$), see Fig. 7 for the evolution of the rapidity
distribution from an initial thermal state, and Fig. 8 for the atom density
$n(t)$. In Fig. 7 we see that the effect of three-body losses is to spread the
rapidity distribution in rapidity space, as already observed for one-body and
two-body losses. We expect that this is a generic effect caused by $K$-body
losses for any $K$. In Fig. 8, we observe that the mean density decays as
$t^{-1/2}$ for an initial infinite temperature state, while for any non-zero
initial temperature it crosses over to a $t^{-\alpha}$ decay at long times
with an exponent $\alpha\simeq 0.21$. This exponent seems to be independent of
the initial temperature as long as it is non-zero, see Fig. 8. However, for an
initial rapidity distribution that is far from thermal, such as for instance
$\rho(k,t=0)=(1-\cos k)/2$ or $(1-\cos(2k))/2$, we find that the density also
decays as a power-law at late time, although with a different exponent
$\alpha$: the exponent is close to $0.21$ for $\rho(k,t=0)=(1-\cos k)/2$, and
close to $0.38$ for $\rho(k,t=0)=(1-\cos(2k))/2$.

We have not been able to analytically derive the observed generic power-law
decay for $K=3$ or for higher $K$, beyond the special case of the initial
infinite temperature state. The latter case is easily understood because, for
an infinite temperature state the rapidity distribution is constant,
$\rho(k)=n$, and the equation (15) can be solved analytically. Then the
determinant of the matrix $B$ is equal to the determinant of $C$, and the
matrices $A$ reduce to identical and diagonal matrices. Therefore the loss
functional is simply given by

$\displaystyle F_{K}[\rho](k)=Kn^{K},$ (65)

which is the result expected from the mean-field approach. Then the solution
of the evolution equation (15) gives the mean density

$\displaystyle n(t)=\dfrac{n(0)}{(1+n(0)^{K-1}\,K(K-1)\,\Gamma t)^{1/(K-1)}}.$
(66)

Beyond that simple case, we have not been able to express the loss functional
in a simple form so as to derive the long-time decay of the mean density.

Similarly to the $K=1$ and $K=2$ cases, we have investigated the behavior of
the rescaled rapidity distribution $\rho(k,t)/n(t)$ at late times. Recall that
this ratio reveals that the gas generically (i.e. unless the first Fourier
mode of $\rho(k)$ is tuned to zero) goes to a ,low-density, low-temperature
thermal state for $K=2$, while for $K=1$ it never does. In Fig. 9 we display
this ratio at late time for $K=3$. Fig. 9.(a) corresponds to a thermal initial
rapidity distribution, and Figs. 9.(b)-(c) to non-thermal initial rapidity
distributions $\rho(k,t=0)=(1-\cos k)/2$ and $(1-\cos(2k))/2$ respectively. We
observe that the rescaled rapidity distribution concentrates around the maxima
of the initial rapidity distribution at long times. Even for an initial
thermal distribution (Fig. 9 (a)), the long time behavior of the density
profile can not be described by a Boltzmann distribution, as it looks like a
bell-shaped distribution that has a small dip at $k=0$. A similar conclusion
holds for Fig. 9.(b). Finally Fig. 9.(c) shows the emergence of peaks
localised at $k=\pm\pi/2$. We conclude that, in contrast with the $K=2$ case,
the late-time rapidity distribution is generically non-thermal.

##  V Harmonically trapped gas

In many cold atom experiments, the gas lies in a longitudinal trapping
potential. This prompts us to study the influence of the trapping potential on
the dynamics of our lossy lattice hard-core gas. For simplicity we restrict to
a harmonic potential $V(x)=\omega^{2}x^{2}/2$.

We adopt a coarse-grained perspective of the gas: we assume that the gas can
be divided into fluid cells which contain a large number of bosons, and that
the state of the gas within each fluid cell $[x,x+dx]$ is a certain macrostate
represented by the local density of rapidities $\rho(x,k)$. Such coarse-
grained descriptions have been very successful lately in describing the out-
of-equilibrium quantum many-body dynamics of nearly integrable gases [74, 73,
64, 80, 81, 75, 82]. Here we investigate the effect of losses on our lattice
hard-core gas within that coarse-grained description.

The equation satisfied by the position-dependent rapidity distribution is

$\displaystyle\partial_{t}\rho(x,k,t)+\sin(k)\,\partial_{x}\rho(x,k,t)-\omega^{2}x\,\partial_{k}\rho(x,k,t)$
$\displaystyle=-\Gamma F[\rho(x,.,t)](k).$ (67)

In the first line, the term $\partial_{x}\sin(k)\rho(x,k)$ corresponds to the
gradient of the current of quasi-particles with rapidity $k$,
$j(x,k)=\sin(k)\rho(x,k)$. Here $\sin(k)$ is the group velocity of quasi-
particles with lattice dispersion relation $\varepsilon(k)=-\cos(k)$. The term
$-\omega^{2}x\,\partial_{k}\rho(x,k)$ in Eq. (V) corresponds to Newton’s
second law, and encodes the fact that the quasi-particles feel the harmonic
potential and are accelerated according to
$\dot{k}=-\partial_{x}V(x)=-\omega^{2}x$. Finally, the r.h.s of Eq. (V) is the
loss term at position $x$, which follows from the assumption that the gas is
locally homogeneous so that we can apply the formalism developed in previous
sections, this time within each fluid cell $[x,x+dx]$.

Figure 10: black Evolution of the position-dependent distribution $\rho(x,k)$
in a harmonic trap and under $K$-body losses. The initial distribution is
$\rho(x,k)=1/(1+\exp((-\cos(k)+\omega^{2}x^{2}/2)/T))$ with harmonic trap
frequency $\omega$=1/2 and temperature $T$=0.1. The loss rate is $\Gamma$=0.1.
The different snapshots show increasing times in units of $\tau=2\pi/\omega$.
Below each phase portrait we add the corresponding density profiles $n(x,t)$
obtained by integrating $\rho(x,k,t)$ over $k$. The density profiles are
rescaled by the initial density at the trap center $n_{0}=n(x=0,t=0)=0.5$. The
two first rows show the effects of one-body losses ($K$=1) on the distribution
in phase-space and the corresponding density profiles. The two next rows show
the effects of two-body losses ($K$=2) and the last two rows show the effect
of three-body losses ($K$=3). Numerically, the $x$ and $k$-axis are
discretized with 300 points on each axis, and we use a time step of $\Delta
t=0.05$; we provide a check of convergence of the numerical method in Appendix
E.

###  V.1 Numerical method

Our main goal in this section is to solve numerically the evolution equation
(V). For this we use a split-step method. Assuming that we know the rapidity
distribution $\rho_{t}(x,k)$, from time $t$ to $t+\Delta t$ we first compute
the new rapidity distribution $\rho^{\prime}_{t+\Delta t}(x,k)$ generated by
the transport of quasi-particles, and then compute $\rho_{t+\Delta t}(x,k)$
from $\rho^{\prime}_{t+\Delta t}(x,k)$ by implementing localized lossy
evolution during a time step $\Delta t$.

black This gives the following scheme. For the transport part we rely on the
fact that the underlying dynamics is the one of non-interacting quasi-
particles. Each quasi-particle at position $(x,k)$ in phase space evolves as

$\left\\{\begin{array}[]{rcl}\dfrac{dx}{dt}=\sin(k)\\\ \\\
\dfrac{dk}{dt}=-\omega^{2}x.\end{array}\right.$ (68)

Let us call $X(x,k,t)$ and $K(x,k,t)$ the analytical solution of that
equation. Then, the first step of our numerical scheme is

$\displaystyle\rho^{\prime}_{t+\Delta t}(X(x,k,t),K(x,k,t))$ $\displaystyle=$
$\displaystyle\rho_{t}(x,k),$ (69)

and the second step is

$\rho_{t+\Delta t}(x,k)=y(\Delta t,k),$ (70)

where $y(\tau,k)$ is the solution of the differential equation

$\partial_{\tau}y(\tau,k)=-\Gamma F[y(\tau,.)](k),\qquad
y(0,k)=\rho^{\prime}_{t+\Delta t,x}.$ (71)

The rapidity distribution $\rho(x,k)$ is discretized on a regular grid in
phase space, and the two steps in Eq. (69) and Eq. (70) are implemented as
follows.

For the transport step in (69), we numerically integrate Eq. (68) for each
point of the grid. Then, starting from the values $\rho_{t}(x,k)$ on the
regular grid in phase space, we move each node of the grid according to
$(x,k)\rightarrow(X(x,k,\Delta t),K(x,k,\Delta t))$. This gives us the
rapidity distribution after transport over a time $\Delta t$, which however is
defined on a new grid, different from the initial regular grid. To get the new
rapidity distribution $\rho^{\prime}_{t+\Delta t}(x,k)$ on the initial regular
grid, we use linear interpolation. Notice that the whole procedure, including
both the numerical integration of (68) and the linear interpolation, needs to
be performed only once. To benchmark this transport part, we have checked that
this method gives excellent numerical precision for the simulated transport in
the absence of losses.

The second step consists in solving numerically the differential equation (71)
for each column of the grid in phase-space. The solution of each differential
equation is obtained by the Runge-Kutta method with the initial condition
$y(0,k)=\rho^{\prime}_{t+\Delta t,x}$. We do this for each column of the grid,
and we thus get the new space-dependent rapidity distribution $\rho_{t+\Delta
t}(x,k)$ according to Eq. (70).

The combination of the two steps allows us to go from a rapidity distribution
$\rho_{t}(x,k)$ to the new distribution $\rho_{t+\Delta t}(x,k)$, both defined
on the same phase-space grid. We then repeat this procedure many times with a
small time step $\Delta t$ to simulate the lossy evolution of the gas in the
trap.

###  V.2 Results

We have performed numerical simulations of the evolution of the position-
dependent rapidity distribution $\rho(x,k)$ under $K$-body losses using the
algorithm presented in the previous section, see Fig. 10. For the initial
state, we use a thermal (Fermi-Dirac) rapidity distribution
$\rho(x,k)=1/(1+\exp{(-\cos(k)+\omega^{2}x^{2}/2-\mu)/T})$, where $\mu$ is the
chemical potential and $T$ the temperature.

Our numerical study allows us to make the following general observations,
illustrated in Fig. 10. For all $K$ (we have simulated $K=1,2,3$) the
distribution typically spreads in phase-space, similarly to the homogeneous
case. However, while the number of particles decays exponentially for $K=1$,
the loss dynamics is much slower for higher $K$, and this has visible effects
on the distribution of rapidities after a given percentage of lost atoms.

For $K=1$ we observe that the edges of the phase-space distribution get
depopulated very fast and quickly results in a halo around the origin. The
spreading is also clearly visible in real space in the particle density
profile $n(x)=\int\rho(x,k)dk/(2\pi)$, see the second line of Fig. 10.

For $K=2$, the situation is a little bit different, see Fig. 10. Until
$\sim$20$\%$ of the atoms have been lost, the dynamics is similar to the one
for $K=1$, but after that we observe the formation of spirals in the bulk of
the phase-space distribution. The spirals become more visible as the ratio
$\Gamma/\omega$ is increased, see Appendix D. Compared to $K=1$, the bulk of
the distribution also gets depopulated, leading to an approximately uniform
circular droplet in phase-space, with a density that decays with time. This is
visible in Fig. 10 after $\sim$50$\%$ of the atoms have been lost; we note
that, compared to the case of one-body losses, the distribution spreads less
significantly in phase-space.

Figure 11: black Time evolution of the total particle number under two-body
losses. The initial distribution is
$1/(1+\exp((-\cos(k)+\omega^{2}x^{2}/2-\mu)/T)$, where $T=0.1$ and $\mu$, the
chemical potential, is chosen to fix the initial mean particle density to 0.5
and $N(0)$ represents the initial number of particle. The loss rate is
$\Gamma$=0.1. (a) Results for two-body losses. The red dashed line is the long
time behavior of the total particle number for the homogeneous case with an
initial rapidity distribution which has no first Fourier mode (see Fig.(4)).
(b) Three-body losses. The black dashed line represents the asymptotic
behavior of the mean density and the red dashed line is the mean-field
prediction.

For $K=3$, we also observe a small spiral appearing at the center of the
distribution after $\sim$30$\%$ of the atoms have been lost. This spiral
remains localised at the center during the dynamics. Interestingly, as one can
see from the fifth line in Fig. 10, this time it is the center of the phase-
space distribution that decreases faster compared to the edges. After
$\sim$60$\%$ of atoms lost, a hole starts developing at the center of the
phase-space distribution, and the latter looks more and more like a ring. This
is a clear signature of a strongly out-of-equilibrium gas in the trap, with a
population inversion: the higher-energy single-particle orbitals get more
populated than the low-energy ones. This effect is reflected in the
corresponding real-space particle density $n(x)$ then acquires a doubly-peaked
shape, with a local minimum at $x=0$, similarly to what we observed also in
the homogeneous case.

Finally, having simulated the dynamics of the position-dependent rapidity
distribution $\rho(x,k)$, one can easily get the evolution of the total
particle number $N(t)$ by integrating $\rho(x,k,t)$ over $x$ and $k$ at fixed
time $t$. For one-body losses, one always finds an exponential decay
$N(t)=e^{-\Gamma t}N(0)$. For two-body and three-body losses the result is
more interesting. In Fig. 11.(a), we show the evolution of the mean density
under two-body losses for different trap frequencies $\omega$. Starting with a
thermal distribution at temperature $T$ and chemical potential $\mu$, we
observe that the mean density decreases at long times as $\sim 1/t$, a result
that coincides with the mean-field (or infinite temperature) decay for the
homogeneous gas. blackMoreover, we see that the higher the trap frequency the
faster the decay of the total particle number. Indeed, as the trap frequency
is increasing, the distribution concentrates around the center of the
confinement potential. Since the density of particles is higher, the
probability of finding two neighboring bosons increases and the depopulation
is then accelerated.

For three-body losses, we see in Fig. 11.(b) that, like in the two-body case,
a stronger confinement speeds up the decrease of total particle number.
However, this time we do not observe a clear convergence towards the
expectation from the mean-field (or infinite temperature) result for the
homogeneous case, which would be $N(t)\sim t^{-1/2}$. We observe a number of
particles that decreases approximately as a power-law $\sim t^{-\alpha}$ with
an exponent $\alpha\simeq 0.6$. Our numerics does not allow us to draw a clear
conclusion as to whether or not this would go to the mean-field exponent $1/2$
at longer times. blackNevertheless, let us stress that, qualitatively, the
effect of three-body losses is the same as in the two-body case: compared to
the homogeneous case, the trap dramatically speeds up the losses. This can be
explained by the same mechanism as in the two-body case: since the density of
particles is higher in the center, the probability of finding three
neighboring bosons increases and the losses are enhanced.

##  VI Conclusion

We have studied the effects of $K$-body losses on a gas of lattice hardcore
bosons, in particular their effect on the thermalisation of the gas at late
time. For this, we have relied on the hypothesis of adiabatic losses used
previously in Refs. [40, 76, 77, 78, 25]. We derived analytical results for
the loss functional for any integer $K$ in the form of a small finite
determinant, and closed expressions in the cases $K=1,2$.

For $K=1$ and $K=2$, we solved analytically the time evolution equation of the
rapidity distribution of the spatially homogeneous gas. In the case of one-
body losses, our formula (54) shows that the loss functional is in general
non-linear and non-local in rapidity space, as already observed for the
continuous Lieb-Liniger gas in Ref. [40]. After investigating the long time
behavior of the rapidity distribution, we concluded that one-body losses do
not drive the gas to a low-density thermal equilibrium state at long times. In
the case of two-body losses, our formula (58) gives an implicit expression for
the rapidity distribution and using a similar method as in Ref. [25], we were
able to investigate the long time behavior of the rapidity distribution and
the mean particle density. In particular, we found that it decays generically
as $\sim 1/\sqrt{t}$, except when the first Fourier mode of the initial
distribution vanishes; in that case the particle density decays as $\sim 1/t$.
A similar conclusion was drawn for a different loss process in Refs.[25, 43].

Finally, we considered the inhomogeneous system consisting in a lattice
hardcore bosons gas in harmonic potential. We provided a numerical method to
solve the dynamics combining the effects of the losses and of the trapping
potential. We observed that for $K\geq 2$ the trap generically speeds up the
decay of the total particle number. We also found that the gas typically
evolves towards a highly non-thermal state; in particular for $K=3$ we observe
a striking ring-shaped distribution in phase space, which signals a inversion
of population (i.e. higher energy single-particle orbitals get more populated
than the lower energy ones), see Fig. 10. Further investigations are needed to
draw more quantitative conclusions about the dynamics in the trap.

Note: while we were finishing this paper, a preprint by Perfetto, Carollo,
Garrahan and Lesanovsky appeared [83], where a similar model of spinless
fermions with $K=2,3,4$-body losses is studied within the context of quantum
reaction-diffusion dynamics of annihilation processes. Our hard-core model
coincides with their model for $K$ even, but not for $K$ odd because of the
Jordan-Wigner mapping, as we explained in detail in Sec. III. Perfetto et al.
do not focus on the effect of losses on the rapidity distributions of the gas,
but rather on the evolution of the number of particles in the homogeneous
setting. For the case of $K$ even, where our models coincide, their findings
about the evolution of the number of particles in the homogeneous setting are
in agreement with ours.

###### Acknowledgements.

We thank Alberto Biella, Isabelle Bouchoule, Mario Collura and Leonardo Mazza
for very useful discussions and for joint work on closely related topics. The
work of JD, FR and DK is supported by the Agence Nationale de la Recherche
through ANR-20-CE30-0017-01 project ‘QUADY’ and ANR-22-CE30-0004-01 project
‘UNIOPEN’. L.R. acknowledges hospitality from LPCT during the completion of
this work.

##  Appendix A Deriving solution (54)

For one-body losses, the time evolution of the rapidity distribution
$\rho(t,k)$ is given by

$\displaystyle\partial_{t}\rho=-\Gamma\left(\rho-(\rho^{2}-\mathcal{H}(\rho)^{2}-n^{2}(t))+2n(t)\mathcal{H}^{\prime}(\rho)\right),$
(72)

where $\Gamma$ is the loss rate. We introduced the Hilbert transform
$\mathcal{H}(f(x))=\frac{1}{2\pi}\int^{\pi}_{-\pi}dy\frac{f(y)}{\tan(\frac{x-y}{2})}$
with $f(x)$ a periodic function. The mean density $n(t)$ is known:
$n(t)=n_{0}e^{-\Gamma t}$.

Here the rapidity distribution is a $2\pi$-periodic real-valued function. From
the rapidity distribution, we can construct a complex-valued function whose
the imaginary part is the Hilbert transform of the real part:
$Q=\rho(k)+i\mathcal{H}(\rho(k))$. Such a function is called an analytic
signal and can be analytically continued to the upper half-plane: the function
$Q(z)=\frac{i}{2\pi}\int^{\pi}_{-\pi}\frac{dq\;\rho(q)}{\tan(\frac{z-q}{2})}$
is well-defined for $\imaginary(z)>0$ and $\real(z)\in[-\pi,\pi]$ and reduces
to $\rho(k)$ on the real axis.

Taking the Hilbert transform of (72)

$\displaystyle\partial_{t}\mathcal{H}(\rho)=-\Gamma\left(\mathcal{H}(\rho)-\mathcal{H}(\rho^{2}-\mathcal{H}(\rho)^{2})+2n(t)\mathcal{H}(\mathcal{H}^{\prime}(\rho))\right)$
(73)

and adding (72) and i (73), one has

$\displaystyle\partial_{\tau}Q(\tau,z)=-(Q(\tau,z)-i2n\partial_{z}Q(\tau,z)-Q^{2}(\tau,z)+n^{2}(\tau)).$
(74)

We used some properties of the Hilbert transform: i)
$\mathcal{H}(\mathcal{H}(f))=-f$, ii)
$\mathcal{H}^{\prime}(f)=\mathcal{H}(f^{\prime})$. Moreover, since $Q^{2}(z)$
is analytic for $\imaginary(z)>0$, the function
$\rho^{2}-\mathcal{H}(\rho)^{2}+i2\rho\mathcal{H}(\rho)$ is an analytic signal
if and only if
$\mathcal{H}(\rho^{2}-\mathcal{H}(\rho)^{2})=2\rho\mathcal{H}(\rho)$.
Introducing the function $Y(\tau,z)=Q(\tau,z+i2n(\tau))$, one gets

$\displaystyle\partial_{\tau}Y(\tau,z)=Y^{2}(\tau,z)-Y(\tau,z)-n^{2}(\tau)$
(75)

This equation can be solved if one assumes
$Y(\tau,z)=\alpha(\tau,z)\;e^{-\tau}$. Indeed, thanks to this trick the above
equation reduces to

$\displaystyle\partial_{\tau}\alpha(\tau,z)=(\alpha^{2}(\tau,z)-n^{2}_{0})\;e^{-\tau}$
(76)

Putting all terms depending on $\alpha$ in the left-hand side, one has

$\displaystyle\int\dfrac{d\alpha}{\alpha^{2}-n^{2}_{0}}=-e^{-\tau}+C_{1},$
(77)

which leads to

$\displaystyle\alpha(\tau,z)=n_{0}\tanh(n_{0}\,e^{-\tau}+C_{2}).$ (78)

The initial condition $Y(0,z)=Y_{0}$ sets the constant:
$C_{2}=\tanh^{-1}(Y_{0}/n_{0})-n_{0}$. Thus, one can write

$\displaystyle Y(\tau,z)$
$\displaystyle=n(\tau)\tanh(n_{0}(e^{-\tau}-1)+\tanh^{-1}(Y_{0}/n_{0}))$
$\displaystyle=n(\tau)\left(\dfrac{\tanh(n_{0}(e^{-\tau}-1))+Y_{0}/n_{0}}{1+\tanh(n_{0}(e^{-\tau}-1))Y_{0}/n_{0}}\right).$
(79)

Finally, the rapidity distribution reads

$\displaystyle\rho(t,k)=n_{0}e^{-\Gamma
t}\real\left(\dfrac{\tanh(n_{0}(e^{-\Gamma t}-1))+\frac{i}{2\pi
n_{0}}\int^{\pi}_{-\pi}dq\dfrac{\rho_{0}(q)}{\tan(\frac{k-q}{2}+in_{0}(1-e^{-\Gamma
t}))}}{1+\frac{i}{2\pi n_{0}}\tanh(n_{0}(e^{-\Gamma
t}-1))\int^{\pi}_{-\pi}dq\dfrac{\rho_{0}(q)}{\tan(\frac{k-q}{2}+in_{0}(1-e^{-\Gamma
t}))}}\right).$ (80)

##  Appendix B Derivation of Eq. (58) of the main text

In this section we present the main steps to derive the exact expression of
the rapidity distribution in the homogeneous case for $K=2$. We consider the
time evolution of $n(t)$ written as:

$\partial_{t}n(t)=\dfrac{1}{2\pi}\int_{-\pi}^{\pi}\partial_{t}\rho(k,t)dk.$
(81)

We now insert the evolution equation (57) obtaining:

$\displaystyle\partial_{t}n(t)=$ $\displaystyle-\dfrac{2\Gamma
n(t)}{2\pi}\int_{-\pi}^{\pi}\rho(k,t)dk+\dfrac{\Gamma}{2\pi^{2}}\left(\int_{-\pi}^{\pi}\cos{(k)}\rho(k,t)dk\right)^{2}$
$\displaystyle=$ $\displaystyle-2\Gamma
n(t)^{2}+\dfrac{\Gamma}{2\pi^{2}}\left(\int_{-\pi}^{\pi}\cos{(k)}\rho(k,t)dk\right)^{2}$
(82)

By inverting the latter relation, one obtains:

$\left\lvert\int_{-\pi}^{\pi}\cos{(k)}\rho(k,t)dk\right\rvert=\pi\sqrt{\dfrac{2}{\Gamma}\left(\partial_{t}n(t)+2\Gamma
n(t)^{2}\right)},$ (83)

where $\lvert\bullet\rvert$ denotes the absolute value. At this point of the
derivation it is useful to introduce the following variable:
$\sigma(t)=\text{sgn}(\int_{-\pi}^{\pi}\cos{(k)}\rho(k,t)dk)$, with
$\text{sgn}(x)$ being the sign function. We now claim that the function
$\sigma(t)$ is solely determined by its value at initial time, i.e.
$\sigma(t)=\sigma(0)\doteqdot\sigma_{0}$, the argument goes as follows. (i) If
the first Fourier mode vanishes at some time $t$, then it must vanish also at
any later time. This follows from Eq. (57). (ii) This implies that the sign of
the first Fourier mode is continuous in time. Since it can take only discrete
values, it is in fact a constant.

By inserting the latter equation in Eq. (57) the time evolution of the
rapidity distribution can be then recasted into the following form:

$\dot{\rho}(k,t)=-2\Gamma\,\rho(k)\,n(t)+\Gamma\sigma_{0}\cos{(k)}\rho(k)\sqrt{\dfrac{2}{\Gamma}\left(\partial_{t}n(t)+2\Gamma
n(t)^{2}\right)}.$ (84)

We now divide both sides by $\rho(k,t)$ and then integrate:

$\ln{\left(\dfrac{\rho(k,t)}{\rho_{0}(k)}\right)}=-2\Gamma\,\int_{0}^{t}n(t^{\prime})dt^{\prime}+\Gamma\sigma_{0}\cos{(k)}\int_{0}^{t}\sqrt{\dfrac{2}{\Gamma}\left(\partial_{t^{\prime}}n(t^{\prime})+2\Gamma
n(t^{\prime})^{2}\right)}dt^{\prime}.$ (85)

By exponentiating the latter equation we get:

$\rho(k,t)=\rho_{0}(k)\exp{-2\Gamma\int_{0}^{t}n(t^{\prime})dt^{\prime}+\sigma_{0}\cos{(k)}\int_{0}^{t}\sqrt{\dfrac{2\left(\partial_{t^{\prime}}n(t^{\prime})+2\Gamma
n(t^{\prime})^{2}\right)}{\Gamma}}dt^{\prime}},$ (86)

this concludes the derivation of Eq. (58).

##  Appendix C Additional data for the homogeneous $K=2$ case

Figure 12: Left panel: the mean density under two-body losses for two
different non-thermal rapidity distributions. Solid colored curves are
obtained by solving numerically the time evolution equation of $\rho(k)$ for
the loss functional (31) with a time step $dt$=0.05 and a loss rate
$\Gamma=$0.1. blackThe red dashed line is the long time limit of the mean
density given by Eq. (59). The long time behavior for non vanishing first
Fourier mode distributions is presented by the black dashed line (see
Eq.(64)). Right panel: difference between the two function $f(t)$ and $g(t)$
defined in the main text for a thermal and non-thermal distributions.

In this section we present additional data concerning the homogeneous $K=2$
case. In particular, we show in Fig. 12 (left panel) the dynamics of the mean
density for two different initial rapidity distributions which are not
thermal. Firstly, we consider an initial distribution given by
$\rho_{0}(k)=\dfrac{1}{2}(1-\cos{(k)})$, which has first Fourier mode
different from zero. Secondly, we consider
$\rho_{0}(k)=\dfrac{1}{2}(1-\cos{(5k)})$, whose first Fourier mode vanishes.
We see that the dynamics induced by the former distribution has a longtime
behaviour guven by $\sim 1/\sqrt{t}$, whereas the latter, due to its
vanishining first Fourier mode, is described by Eq. (59). This corroborates
our findings for initial thermal distributions presented in the main text.
Moreover, we show in Fig. 12 (right panel) the quantity $f(t)-g(t)$, where
$g(t)=\int_{0}^{t}n(\tau)\,d\tau$ and
$f(t)=\int_{0}^{t}\sqrt{1+\frac{\partial_{\tau}n(\tau)}{2\Gamma\,n(\tau)^{2}}}n(\tau)d\tau$
for two different rapidity distributions. In the main text we took the first
order expansion of $f(t)$ resulting in a logarthmic growth for the quantity
$f(t)-g(t)$, which is thus corroborated by the numerical data here presented.
As such, given an initial rapidity distribution whose first Fourier mode is
non-vanishing, one has a longtime decay of the mean density given by $n(t)\sim
1/\sqrt{t}$.

##  Appendix D Brief discussion on the spiral in figure Fig. 10

During the evolution of the position-dependent rapidity distribution in phase
space (see Fig. 10) , the distribution exhibits a spiral, which is clearly
visible after 40$\%$ of atoms lost for two-body losses. In principle, in a
regime where the trap frequency $\omega$ is highly dominating the loss rate
$\Gamma$, one expects that the distribution remains rotation invariant at all
times. This implies that the spiral vanishes for $\omega\gg\Gamma$. To check
this statement we compare the quantity $\rho(x,k)$ for two distincts values of
$\omega$ (see the figure below). In figure Fig. 13 We can see that the spiral
appearing for $\omega=5\Gamma$ covers entirely the distribution, while for
$\omega=20\Gamma$ the spiral is localised at the distribution’s center.
Moreover, in the case $\omega=20\Gamma$ we observe small oscillations between
the edges and the center of the distribution. The frequency of these
oscillations is high compared to the center of the distribution.

Figure 13: Phase portraits showing the spiral for two different values of
$\omega$. The initial distribution is identical to the one used in figure Fig.
10 and here we only consider the distribution after 40$\%$ of atoms lost.

##  Appendix E Convergence of numerical results

black In this Appendix we provide evidence for numerical convergence of the
results displayed in Fig. 10. We have performed the same simulation as the one
in Fig. 10 with a bigger numerical time step $\Delta=0.25$ (instead of
$\Delta=0.05$). The physical parameters we use are (as in Fig. 10):

$\displaystyle\rho(x,k,t=0)=\frac{1}{1+\exp{(-\cos(k)+\omega^{2}x^{2}/2)T}},$
(87)

with $T=0.1$ and $\omega=1/2$. We show results for one-body losses ($K=1$),
with a loss rate $\Gamma$=0.1.

In Fig. 14, we compare the results obtained for two different time steps
$\Delta t=0.25$ (first row) and $\Delta t=0.05$. We see that the results are
converged.

Figure 14: black Dynamics of the local rapidity distribution $\rho(x,k)$ in
phase space under one-body losses. The two first lines show the evolution of
the distribution in phase space for two different time steps $\Delta t=0.25$
and $\Delta t=0.05$. The trap frequency is fixed to $\omega=1/2$, $\Gamma=0.1$
is the loss rate and $\tau$ is defined as $\tau=2\pi/\omega$. The last line
corresponds to the density profils obtained by integrating over $k$ the local
rapidity distribution. The numerical results are converged.

## References

  * [1] I. Bloch, J. Dalibard, and W. Zwerger, “Many-body physics with ultracold gases,” Rev. Mod. Phys., vol. 80, pp. 885–964, Jul 2008.  10.1103/RevModPhys.80.885.
  * [2] S. Knoop, J. S. Borbely, R. van Rooij, and W. Vassen, “Nonexponential one-body loss in a bose-einstein condensate,” Phys. Rev. A, vol. 85, p. 025602, Feb 2012.  10.1103/PhysRevA.85.025602.
  * [3] Browaeys, A., Poupard, J., Robert, A., Nowak, S., Rooijakkers, W., Arimondo, E., Marcassa, L., Boiron, D., Westbrook, C. I., and Aspect, A., “Two body loss rate in a magneto-optical trap of metastable he,” Eur. Phys. J. D, vol. 8, no. 2, pp. 199–203, (2000).  10.1007/s100530050027.
  * [4] L. Franchi, L. F. Livi, G. Cappellini, G. Binella, M. Inguscio, J. Catani, and L. Fallani, “State-dependent interactions in ultracold 174yb probed by optical clock spectroscopy,” New Journal of Physics, vol. 19, p. 103037, nov 2017.  10.1088/1367-2630/aa8fb4.
  * [5] T. Tomita, S. Nakajima, I. Danshita, Y. Takasu, and Y. Takahashi, “Observation of the mott insulator to superfluid crossover of a driven-dissipative bose-hubbard system,” Science Advances, vol. 3, no. 12, p. e1701513, 2017\.  10.1126/sciadv.1701513.
  * [6] N. Syassen, D. M. Bauer, M. Lettner, T. Volz, D. Dietze, J. J. García-Ripoll, J. I. Cirac, G. Rempe, and S. Dürr, “Strong dissipation inhibits losses and induces correlations in cold molecular gases,” Science, vol. 320, no. 5881, pp. 1329–1331, 2008.  10.1126/science.1155309.
  * [7] T. Tomita, S. Nakajima, Y. Takasu, and Y. Takahashi, “Dissipative bose-hubbard system with intrinsic two-body loss,” Phys. Rev. A, vol. 99, p. 031601, Mar 2019.  10.1103/PhysRevA.99.031601.
  * [8] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, “Observation of dipolar spin-exchange interactions with lattice-confined polar molecules,” Nature, vol. 501, pp. 521–525, Sep 2013\.  10.1038/nature12483.
  * [9] K. Sponselee, L. Freystatzky, B. Abeln, M. Diem, B. Hundt, A. Kochanke, T. Ponath, B. Santra, L. Mathey, K. Sengstock, and C. Becker, “Dynamics of ultracold quantum gases in the dissipative fermi–hubbard model,” Quantum Science and Technology, vol. 4, p. 014002, sep 2018.  10.1088/2058-9565/aadccd.
  * [10] K. Honda, S. Taie, Y. Takasu, N. Nishizawa, M. Nakagawa, and Y. Takahashi, “Observation of the sign reversal of the magnetic correlation in a driven-dissipative fermi gas in double wells,” Phys. Rev. Lett., vol. 130, p. 063001, Feb 2023.  10.1103/PhysRevLett.130.063001.
  * [11] J. Söding, D. Guéry-Odelin, P. Desbiolles, F. Chevy, H. Inamori, and J. Dalibard, “Three-body decay of a rubidium bose–einstein condensate,” Applied Physics B, vol. 69, pp. 257–261, Oct 1999.  10.1007/s003400050805.
  * [12] T. Weber, J. Herbig, M. Mark, H.-C. Nägerl, and R. Grimm, “Three-body recombination at large scattering lengths in an ultracold atomic gas,” Phys. Rev. Lett., vol. 91, p. 123201, Sep 2003.  10.1103/PhysRevLett.91.123201.
  * [13] B. L. Tolra, K. M. O’Hara, J. H. Huckans, W. D. Phillips, S. L. Rolston, and J. V. Porto, “Observation of reduced three-body recombination in a correlated 1d degenerate bose gas,” Phys. Rev. Lett., vol. 92, p. 190401, May 2004.  10.1103/PhysRevLett.92.190401.
  * [14] F. Ferlaino, S. Knoop, M. Berninger, W. Harm, J. P. D’Incao, H.-C. Nägerl, and R. Grimm, “Evidence for universal four-body states tied to an efimov trimer,” Phys. Rev. Lett., vol. 102, p. 140401, Apr 2009.  10.1103/PhysRevLett.102.140401.
  * [15] J. H. Gurian, P. Cheinet, P. Huillery, A. Fioretti, J. Zhao, P. L. Gould, D. Comparat, and P. Pillet, “Observation of a resonant four-body interaction in cold cesium rydberg atoms,” Phys. Rev. Lett., vol. 108, p. 023005, Jan 2012.  10.1103/PhysRevLett.108.023005.
  * [16] D. M. Harber, J. M. McGuirk, J. M. Obrecht, and E. A. Cornell, “Thermally Induced Losses in Ultra-Cold Atoms Magnetically Trapped Near Room-Temperature Surfaces,” Journal of Low Temperature Physics, vol. 133, pp. 229–238, Nov. (2003).  10.1023/A:1026084606385.
  * [17] A. Johnson, S. S. Szigeti, M. Schemmer, and I. Bouchoule, “Long-lived nonthermal states realized by atom losses in one-dimensional quasicondensates,” Physical Review A, vol. 96, p. 013623, July (2017).  10.1103/PhysRevA.96.013623.
  * [18] I. Bouchoule, M. Schemmer, and C. Henkel, “Cooling phonon modes of a Bose condensate with uniform few body losses,” SciPost Phys., vol. 5, p. 043, 2018.  10.21468/SciPostPhys.5.5.043.
  * [19] L. H. Dogra, J. A. P. Glidden, T. A. Hilker, C. Eigen, E. A. Cornell, R. P. Smith, and Z. Hadzibabic, “Can three-body recombination purify a quantum gas?,” Phys. Rev. Lett., vol. 123, p. 020405, Jul 2019.  10.1103/PhysRevLett.123.020405.
  * [20] P. Grišins, B. Rauer, T. Langen, J. Schmiedmayer, and I. E. Mazets, “Degenerate bose gases with uniform loss,” Phys. Rev. A, vol. 93, p. 033634, Mar 2016.  10.1103/PhysRevA.93.033634.
  * [21] B. Rauer, P. Grišins, I. E. Mazets, T. Schweigler, W. Rohringer, R. Geiger, T. Langen, and J. Schmiedmayer, “Cooling of a one-dimensional bose gas,” Phys. Rev. Lett., vol. 116, p. 030402, Jan 2016.  10.1103/PhysRevLett.116.030402.
  * [22] M. Schemmer, A. Johnson, R. Photopoulos, and I. Bouchoule, “Monte carlo wave-function description of losses in a one-dimensional bose gas and cooling to the ground state by quantum feedback,” Phys. Rev. A, vol. 95, p. 043641, Apr 2017.  10.1103/PhysRevA.95.043641.
  * [23] M. Schemmer and I. Bouchoule, “Cooling a bose gas by three-body losses,” Phys. Rev. Lett., vol. 121, p. 200401, Nov 2018.  10.1103/PhysRevLett.121.200401.
  * [24] I. Bouchoule and M. Schemmer, “Asymptotic temperature of a lossy condensate,” SciPost Phys., vol. 8, p. 060, 2020.  10.21468/SciPostPhys.8.4.060.
  * [25] D. Rossini, A. Ghermaoui, M. B. Aguilera, R. Vatré, R. Bouganne, J. Beugnon, F. Gerbier, and L. Mazza, “Strong correlations in lossy one-dimensional quantum gases: from the quantum Zeno effect to the generalized Gibbs ensemble,” Physical Review A, vol. 103, p. L060201, June (2021).  10.1103/PhysRevA.103.L060201.
  * [26] J. J. García-Ripoll, S. Dürr, N. Syassen, D. M. Bauer, M. Lettner, G. Rempe, and J. I. Cirac, “Dissipation-induced hard-core boson gas in an optical lattice,” New Journal of Physics, vol. 11, p. 013053, Jan. (2009).  10.1088/1367-2630/11/1/013053.
  * [27] E. Misra and C. Sudarshan, “The zeno’s paradox in quantum theory,” J. Math. Phys., vol. 18, p. 756, 1977.
  * [28] W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J. Wineland, “Quantum zeno effect,” Phys. Rev. A, vol. 41, pp. 2295–2300, 1990.  10.1103/PhysRevA.41.2295.
  * [29] A. Beige, D. Braun, and P. L. Knight, “Driving atoms into decoherence-free states,” New J. Phys., vol. 2, pp. 22–22, sep 2000.  10.1088/1367-2630/2/1/322.
  * [30] A. Beige, D. Braun, B. Tregenna, and P. L. Knight, “Quantum computing using dissipation to remain in a decoherence-free subspace,” Phys. Rev. Lett., vol. 85, pp. 1762–1765, Aug 2000.  10.1103/PhysRevLett.85.1762.
  * [31] J. Kempe, D. Bacon, D. A. Lidar, and K. B. Whaley, “Theory of decoherence-free fault-tolerant universal quantum computation,” Phys. Rev. A, vol. 63, p. 042307, Mar 2001.  10.1103/PhysRevA.63.042307.
  * [32] P. Facchi, H. Nakazato, and S. Pascazio, “From the quantum zeno to the inverse quantum zeno effect,” Phys. Rev. Lett., vol. 86, pp. 2699–2703, Mar 2001\.  10.1103/PhysRevLett.86.2699.
  * [33] P. Facchi and S. Pascazio, “Quantum zeno subspaces,” Phys. Rev. Lett., vol. 89, p. 080401, Aug 2002.  10.1103/PhysRevLett.89.080401.
  * [34] R. Schützhold and G. Gnanapragasam, “Quantum zeno suppression of three-body losses in bose-einstein condensates,” Phys. Rev. A, vol. 82, p. 022120, Aug 2010.  10.1103/PhysRevA.82.022120.
  * [35] K. Stannigel, P. Hauke, D. Marcos, M. Hafezi, S. Diehl, M. Dalmonte, and P. Zoller, “Constrained dynamics via the zeno effect in quantum simulation: Implementing non-abelian lattice gauge theories with cold atoms,” Phys. Rev. Lett., vol. 112, p. 120406, Mar 2014.  10.1103/PhysRevLett.112.120406.
  * [36] Z. Gong, S. Higashikawa, and M. Ueda, “Zeno hall effect,” Phys. Rev. Lett., vol. 118, p. 200401, May 2017.  10.1103/PhysRevLett.118.200401.
  * [37] H. Fröml, C. Muckel, C. Kollath, A. Chiocchetta, and S. Diehl, “Ultracold quantum wires with localized losses: Many-body quantum zeno effect,” Phys. Rev. B, vol. 101, p. 144301, Apr 2020.  10.1103/PhysRevB.101.144301.
  * [38] K. Snizhko, P. Kumar, and A. Romito, “Quantum zeno effect appears in stages,” Phys. Rev. Research, vol. 2, p. 033512, Sep 2020.  10.1103/PhysRevResearch.2.033512.
  * [39] A. Biella and M. Schiró, “Many-Body Quantum Zeno Effect and Measurement-Induced Subradiance Transition,” Quantum, vol. 5, p. 528, Aug. 2021.  10.22331/q-2021-08-19-528.
  * [40] I. Bouchoule, B. Doyon, and J. Dubail, “The effect of atom losses on the distribution of rapidities in the one-dimensional Bose gas,” SciPost Phys. 9, 044, June (2020).  10.21468/SciPostPhys.9.4.044.
  * [41] I. Bouchoule, L. Dubois, and L.-P. Barbier, “Losses in interacting quantum gases: Ultraviolet divergence and its regularization,” Phys. Rev. A, vol. 104, p. L031304, Sep 2021.  10.1103/PhysRevA.104.L031304.
  * [42] F. Minganti, V. Savona, and A. Biella, “Dissipative phase transitions in $n$-photon driven quantum nonlinear resonators,” 2023.  10.48550/arXiv.2303.03355.
  * [43] L. Rosso, A. Biella, and L. Mazza, “The one-dimensional Bose gas with strong two-body losses: the effect of the harmonic confinement,” SciPost Physics, vol. 12, p. 044, Jan. (2022).  10.21468/SciPostPhys.12.1.044.
  * [44] C.-H. Huang, T. Giamarchi, and M. A. Cazalilla, “Modeling particle loss in open systems using keldysh path integral and second order cumulant expansion,” 2023.  10.48550/arXiv.2305.13090.
  * [45] B. Zhu, B. Gadway, M. Foss-Feig, J. Schachenmayer, M. L. Wall, K. R. A. Hazzard, B. Yan, S. A. Moses, J. P. Covey, D. S. Jin, J. Ye, M. Holland, and A. M. Rey, “Suppressing the loss of ultracold molecules via the continuous quantum zeno effect,” Phys. Rev. Lett., vol. 112, p. 070404, Feb 2014.  10.1103/PhysRevLett.112.070404.
  * [46] L. Rosso, D. Rossini, A. Biella, and L. Mazza, “One-dimensional spin-1/2 fermionic gases with two-body losses: Weak dissipation and spin conservation,” Phys. Rev. A, vol. 104, p. 053305, Nov 2021.  10.1103/PhysRevA.104.053305.
  * [47] L. Rosso, L. Mazza, and A. Biella, “Eightfold way to dark states in su(3) cold gases with two-body losses,” Phys. Rev. A, vol. 105, p. L051302, May 2022\.  10.1103/PhysRevA.105.L051302.
  * [48] L. Rosso, A. Biella, J. De Nardis, and L. Mazza, “Dynamical theory for one-dimensional fermions with strong two-body losses: Universal non-hermitian zeno physics and spin-charge separation,” Phys. Rev. A, vol. 107, p. 013303, Jan 2023.  10.1103/PhysRevA.107.013303.
  * [49] M. Nakagawa, N. Tsuji, N. Kawakami, and M. Ueda, “Dynamical sign reversal of magnetic correlations in dissipative hubbard models,” Phys. Rev. Lett., vol. 124, p. 147203, Apr 2020.  10.1103/PhysRevLett.124.147203.
  * [50] M. Nakagawa, N. Kawakami, and M. Ueda, “Exact liouvillian spectrum of a one-dimensional dissipative hubbard model,” Phys. Rev. Lett., vol. 126, p. 110404, Mar 2021.  10.1103/PhysRevLett.126.110404.
  * [51] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, “Reaction-limited quantum reaction-diffusion dynamics,” Phys. Rev. Lett., vol. 130, p. 210402, May 2023.  10.1103/PhysRevLett.130.210402.
  * [52] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, “Tonks–girardeau gas of ultracold atoms in an optical lattice,” Nature, vol. 429, no. 6989, pp. 277–281, 2004.
  * [53] T. Kinoshita, T. Wenger, and D. S. Weiss, “Observation of a one-dimensional tonks-girardeau gas,” Science, vol. 305, no. 5687, pp. 1125–1128, 2004\.
  * [54] E. Haller, M. Gustavsson, M. J. Mark, J. G. Danzl, R. Hart, G. Pupillo, and H.-C. Nägerl, “Realization of an excited, strongly correlated quantum gas phase,” Science, vol. 325, no. 5945, pp. 1224–1227, 2009.
  * [55] A. Van Amerongen, J. Van Es, P. Wicke, K. Kheruntsyan, and N. Van Druten, “Yang-yang thermodynamics on an atom chip,” Physical review letters, vol. 100, no. 9, p. 090402, 2008.
  * [56] T. Jacqmin, J. Armijo, T. Berrada, K. V. Kheruntsyan, and I. Bouchoule, “Sub-poissonian fluctuations in a 1d bose gas: from the quantum quasicondensate to the strongly interacting regime,” Physical review letters, vol. 106, no. 23, p. 230405, 2011.
  * [57] H. Moritz, T. Stöferle, K. Günter, M. Köhl, and T. Esslinger, “Confinement induced molecules in a 1d fermi gas,” Physical review letters, vol. 94, no. 21, p. 210401, 2005.
  * [58] Y.-a. Liao, A. S. C. Rittner, T. Paprotta, W. Li, G. B. Partridge, R. G. Hulet, S. K. Baur, and E. J. Mueller, “Spin-imbalance in a one-dimensional fermi gas,” Nature, vol. 467, no. 7315, pp. 567–569, 2010.
  * [59] G. Pagano, M. Mancini, G. Cappellini, P. Lombardi, F. Schäfer, H. Hu, X.-J. Liu, J. Catani, C. Sias, M. Inguscio, et al., “A one-dimensional liquid of fermions with tunable spin,” Nature Physics, vol. 10, no. 3, pp. 198–201, 2014.
  * [60] P. Wicke, S. Whitlock, and N. Van Druten, “Controlling spin motion and interactions in a one-dimensional bose gas,” arXiv preprint arXiv:1010.4545, 2010.
  * [61] J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, I. P. McCulloch, F. Heidrich-Meisner, I. Bloch, and U. Schneider, “Expansion dynamics of interacting bosons in homogeneous lattices in one and two dimensions,” Physical review letters, vol. 110, no. 20, p. 205301, 2013\.
  * [62] L. Vidmar, J. P. Ronzheimer, M. Schreiber, S. Braun, S. S. Hodgman, S. Langer, F. Heidrich-Meisner, I. Bloch, and U. Schneider, “Dynamical quasicondensation of hard-core bosons at finite momenta,” Physical review letters, vol. 115, no. 17, p. 175301, 2015.
  * [63] F. H. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin, The one-dimensional Hubbard model.  Cambridge University Press, 2005.
  * [64] B. Doyon, “Lecture notes on Generalised Hydrodynamics,” SciPost Physics Lecture Notes, p. 18, Aug. (2020).  10.21468/SciPostPhysLectNotes.18.
  * [65] A. Hutsalyuk and B. Pozsgay, “Integrability breaking in the one dimensional Bose gas: Atomic losses and energy loss,” Physical Review E, vol. 103, p. 042121, Apr. (2021).  10.1103/PhysRevE.103.042121.
  * [66] B. Bertini, F. H. Essler, S. Groha, and N. J. Robinson, “Prethermalization and Thermalization in Models with Weak Integrability Breaking,” Physical Review Letters, vol. 115, p. 180601, Oct. (2015).  10.1103/PhysRevLett.115.180601.
  * [67] M. Rigol and A. Muramatsu, “Fermionization in an expanding 1d gas of hard-core bosons,” Physical review letters, vol. 94, no. 24, p. 240403, 2005.
  * [68] M. Rigol, A. Muramatsu, and M. Olshanii, “Hard-core bosons on optical superlattices: Dynamics and relaxation in the superfluid and insulating regimes,” Physical Review A, vol. 74, no. 5, p. 053616, 2006.
  * [69] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, “Relaxation in a completely integrable many-body quantum system: an ab initio study of the dynamics of the highly excited states of 1d lattice hard-core bosons,” Physical review letters, vol. 98, no. 5, p. 050405, 2007.
  * [70] L. Vidmar, S. Langer, I. McCulloch, U. Schneider, U. Schollwöck, and F. Heidrich-Meisner, “Sudden expansion of mott insulators in one dimension,” Physical Review B, vol. 88, no. 23, p. 235117, 2013.
  * [71] L. Vidmar, W. Xu, and M. Rigol, “Emergent eigenstate solution and emergent gibbs ensemble for expansion dynamics in optical lattices,” Physical Review A, vol. 96, no. 1, p. 013608, 2017.
  * [72] W. Xu and M. Rigol, “Expansion of one-dimensional lattice hard-core bosons at finite temperature,” Physical Review A, vol. 95, no. 3, p. 033617, 2017\.
  * [73] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, “Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium,” Physical Review X, vol. 6, p. 041065, Dec. (2016).  10.1103/PhysRevX.6.041065.
  * [74] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, “Transport in Out-of-Equilibrium X X Z Chains: Exact Profiles of Charges and Currents,” Physical Review Letters, vol. 117, p. 207201, Nov. (2016).  10.1103/PhysRevLett.117.207201.
  * [75] I. Bouchoule and J. Dubail, “Generalized hydrodynamics in the one-dimensional bose gas: theory and experiments,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2022, p. 014003, jan 2022.  10.1088/1742-5468/ac3659.
  * [76] F. Lange, Z. Lenarčič, and A. Rosch, “Time-dependent generalized Gibbs ensembles in open quantum systems,” Physical Review B, vol. 97, p. 165138, Apr. (2018).  10.1103/PhysRevB.97.165138.
  * [77] F. Lange, Z. Lenarčič, and A. Rosch, “Pumping approximately integrable systems,” Nature Communications, vol. 8, p. 15767, Jun 2017\.  10.1038/ncomms15767.
  * [78] Z. Lenarčič, F. Lange, and A. Rosch, “Perturbative approach to weakly driven many-particle systems in the presence of approximate conservation laws,” Phys. Rev. B, vol. 97, p. 024302, Jan 2018.  10.1103/PhysRevB.97.024302.
  * [79] P. L. Butzer and R. J. Nessel, Fourier analysis and approximation.  Lehrbücher und Monographien aus dem Gebiete der exakten Wissenschaften. Mathematische Reihe, Bd. 40, Basel, Stuttgart: Birkhäuser, 1971.
  * [80] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, “Generalized hydrodynamics on an atom chip,” Phys. Rev. Lett., vol. 122, p. 090601, Mar 2019.  10.1103/PhysRevLett.122.090601.
  * [81] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, “Generalized hydrodynamics in strongly interacting 1d bose gases,” Science, vol. 373, no. 6559, pp. 1129–1133, 2021.  10.1126/science.abf0147.
  * [82] M. Coppola and D. Karevski, “Some speculations about local thermalization of nonequilibrium extended quantum systems,” Condensed Matter Physics, vol. 26, no. 1, p. 13502, (2023).  10.5488/CMP.26.13502.
  * [83] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, “Quantum reaction-limited reaction-diffusion dynamics of annihilation processes,” arXiv preprint arXiv:2305.06944, 2023.  10.48550/arXiv.2305.06944.