# Universality and two-body losses: lessons from the effective non-Hermitian
dynamics of two particles

Alice Marché  Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Hironobu Yoshida  Department of Physics, Graduate School of Science, The
University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan  Alberto Nardin
Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France  Hosho Katsura
Department of Physics, Graduate School of Science, The University of Tokyo,
7-3-1 Hongo, Tokyo 113-0033, Japan  Institute for Physics of Intelligence, The
University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan  Trans-scale Quantum
Science Institute, The University of Tokyo, 7-3-1, Hongo, Tokyo 113-0033,
Japan  Leonardo Mazza  Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay,
France

(August 28, 2024)

###### Abstract

We study the late-time dynamics of two particles confined in one spatial
dimension and subject to two-body losses. The dynamics is exactly described by
a non-Hermitian Hamiltonian that can be analytically studied both in the
continuum and on a lattice. The asymptotic decay rate and the universal power-
law form of the decay of the number of particles are exactly computed in the
whole parameter space of the problem. When in the initial state the two
particles are far apart, the average number of particles in the setup decays
with time $t$ as $t^{-1/2}$; a different power law, $t^{-3/2}$, is found when
the two particles overlap in the initial state. These results are valid both
in the continuum and on a lattice, but in the latter case a logarithmic
correction appears.

##  I Introduction

Universality is a key notion in physics related to the possibility of
performing the transition from the microscopic to the macroscopic scale by
systematically discarding irrelevant information [1, 2, 3]. In recent years
significant interest has been sparked by the discovery of emergent universal
behaviors in open quantum many-body systems [4, 5]. Among the many setups
where they have been identified, we list here polariton condensates [6, 7],
lossy gases [8], quantum reactive systems [9, 10, 11], Rydberg atoms [12, 13]
or quantum spin chains [14, 15]. Our focus here is on two-body losses, one of
the most investigated mechanisms that can induce a many-body correlated open
quantum dynamics in a bosonic or fermionic degenerate gas [16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41], and that is also of particular experimental relevance in the
framework of ultra-cold gases [42, 43, 44, 45, 46, 47, 48, 49, 50, 51].
Additionally, some articles have highlighted potential metrological and
quantum information applications [19, 33].

So far, theoretical studies have focused on two universal features: (i) the
scaling of the asymptotic exponential decay rate (ADR) with the linear system
size, $\lambda_{\rm ADR}\sim L^{-\alpha}$, and (ii) the late-time decay of the
gas density, $n(t)\sim t^{-\beta}$. The calculation of the exponent $\alpha$
benefits from several simplifications that are inherent to lossy gases, and
that have allowed for studies in wide ranges of parameters [31]. On the other
hand, the exponent $\beta$ has proven more difficult to compute, and it has
been characterized only under certain restrictive assumptions. In the bosonic
context, for instance, analytical expressions have been derived only for
weakly-interacting gases with adiabatic losses, or for hard-core gases [25,
35, 52, 36, 53]; remarkably, it has been possible to perform the study also in
the presence of a harmonic confinement [54, 55, 53, 41]. However, identifying
the value of $\beta$ for generic values of interparticle interactions or for
arbitrary loss rates remains an open and challenging task.

In this article, we address the problem of computing $\alpha$ and $\beta$ for
generic model parameters by considering a simplified initial condition, namely
that of a quantum system consisting of just two particles. In this situation,
only one two-body loss process can take place and this allows for important
technical simplifications. Indeed, for open quantum systems, it is always
possible to consider the no-click dynamics by performing a postselection, so
that the Lindblad master equation maps onto a simpler non-Hermitian dynamics
[56]. In the situation that we are considering, the number of particles can be
deduced exactly from this effective non-Hermitian dynamics [35].

Since we consider local losses, in order to have a nontrivial loss dynamics we
need to consider an orbital wavefunction that is symmetric under exchange, so
that two particles can occupy the very same position in space. Notice that
such a state can be realized both with bosons (be they spinless or spinful, in
which case the spin part of the wavefunction needs to be symmetric under spin
exchange) or spin-$1/2$ fermions (provided the spin part of the wavefunction
is a singlet). We will discuss two spatial geometries in which the two atoms
are confined: a translationally invariant one-dimensional (1D) tube and a
uniform tight-binding chain with a single energy level per spatial site. In
both cases, we revisit the known results for the ADR exponent $\alpha$ and
compute analytically the value of the density power-law exponent $\beta$ for
different initial conditions.

Our approach is based on an effective “linearization” of the eigenvalue
equations of the non-Hermitian Hamiltonian around the most long-lived modes in
the limit $L\to\infty$. This allows us to describe correctly the late-time
dynamics of the gas in the whole space of parameters, also for intermediate
interaction strengths or loss rates. We present numerical exact results that
are well described by our analytical approximation, that exhibits a universal
scaling collapse. The approach has the additional merit of highlighting the
link existing between the exponents $\alpha$ and $\beta$, which has been
overlooked so far.

In the continuum limit, we compute the value $\alpha=3$; on a lattice, not
only we find a set of modes with $\alpha=3$, but also additional branches
characterized by $\alpha=5$. The value of the exponent $\beta$ depends non-
trivially on the initial conditions. We consider the case of two overlapping
particles, for which $\beta=3/2$, and the case of two well-separated
particles, for which $\beta=1/2$; we will also demonstrate that by changing
initial conditions we can create situations where the dynamics of the mean
number of particles crosses over between these two algebraic decays; no other
exponents $\beta$ appear in the problem. These universal results are valid for
the entire parameter space and provide clear scaling functions that allow one
to collapse the data onto universal curves. The branches characterized by
$\alpha=5$ are responsible for the appearance of a logarithmic correction in
the lattice. It is important to observe that it is difficult to identify the
logarithmic correction by pure numerical means, as in the situations that we
studied it could be interpreted as a different and non-universal scaling
exponent $\beta^{\prime}$. The possibility of studying the problem in a limit
where analytical calculations could be performed has the merit of clarifying
the difficulties associated with numerical studies of $\beta$ in lattice
setups.

These analytical insights result from the simplified initial condition with
two particles. We believe, however, that the linearization technique that we
introduce here, based on the approximation of the dynamics around the most
long-lived modes, can be employed in an appropriate form also for the more
demanding situations in which a many-particle gas at finite density is
initially prepared. One such technique would allow to investigate whether the
results presented in this article capture the dynamics also in finite-density
situations when the time is sufficiently late for the system to be close to
the stationary state.

The article is organized as follows. In Sec. II we introduce the models for a
continuum setup and on a lattice, and we give an overview of the mapping to a
non-Hermitian Hamiltonian that we employ in the study of the dynamics. In Sec.
III, we present a brief summary of our results and describe the general
behavior of the lossy dynamics. In Sec. IV, we discuss the spectrum of the
non-Hermitian Hamiltonian and the linearization procedure. In Sec. V, we
derive universal expressions for the mean number of particles at late times;
we comment also on the appearance of the quantum Zeno effect in this system.
The conclusions and perspectives of our results are presented in Sec. VI. The
article is concluded by nine appendices.

##  II Models and effective non-Hermitian dynamics

We begin by presenting the two models (in the continuum and on a lattice) that
will be studied in the article; we will only consider the case of local
losses, meaning that a loss event can occur only if the two particles have a
non-zero probability of being at the same position. Our focus is on initial
states with two identical spinless bosons and in the final part of this
section, we describe the exact mapping of its Lindblad master-equation time-
evolution onto a non-Hermitian dynamics.

###  II.1 Two spinless bosons in the continuum

We consider a one-dimensional continuous system of length ${\ell}$ (we assume
periodic boundary conditions) populated by bosons with contact interactions
described by the Lieb-Liniger model. We denote by $M$ the mass of the bosons
and by $\sigma$ the typical spatial extension of their individual initial
wavefunctions; the typical initial energy of each boson is
$E_{0}=\hbar^{2}/(2M\sigma^{2})$. In the following, lengths are measured in
units of $\sigma$, energies in units of $E_{0}$, and times in units of
$\hbar/E_{0}$.

Let $\Psi^{\dagger}(x)$ and $\Psi(x)$ be the creation and annihilation
operators of a boson at position $x$. In this case, the Lindblad master
equation [57] for the density matrix of the system $\rho$ reads [16, 24]

$\frac{d\rho}{dt}=-i[H,\rho]+\int_{0}^{\ell}\Bigl{[}L(x)\rho
L^{\dagger}(x)-\frac{1}{2}\\{L^{\dagger}(x)L(x),\rho\\}\Bigr{]}dx,$ (1)

where $L(x)=\sqrt{\frac{\gamma}{2}}\Psi(x)^{2}$ are the jump operators
modeling the local two-body losses and $\gamma>0$ is the loss rate. The Lieb-
Liniger Hamiltonian $H$ is defined by [58, 59]

$H=\int_{0}^{\ell}\Psi^{\dagger}(x)\left(-\frac{\partial^{2}}{\partial
x^{2}}+\frac{g}{2}\Psi^{\dagger}(x)\Psi(x)\right)\Psi(x)dx,$ (2)

where $g>0$ is the repulsive interaction strength; this Hamiltonian can be
diagonalized by the Bethe ansatz [60, 61, 62, 63]. The number operator in this
context reads $N=\int_{0}^{\ell}\Psi^{\dagger}(x)\Psi(x)dx$.

###  II.2 Two spinless bosons on a lattice

We will also consider a one-dimensional lattice of $L$ sites with periodic
boundary conditions populated by bosons that can hop between neighboring
lattice sites and feature contact interaction; this system is described by the
Bose-Hubbard model. We set the lattice spacing $a$ to $1$. We denote
$b_{j}^{{\dagger}}$ and $b_{j}$ the creation and annihilation operators of a
boson at site $j\in\\{1,\ldots,L\\}$. In the following, energies are measured
in units of $J$, the hopping amplitude between neighboring sites, and times
are measured in units of $\hbar/J$. In this case, the master equation reads
[17, 25]

$\frac{d\rho}{dt}=-i\left[H,\rho\right]+\sum_{j=1}^{L}\left(L_{j}\rho
L_{j}^{\dagger}-\frac{1}{2}\\{L_{j}^{\dagger}L_{j},\rho\\}\right).$ (3)

Here $H$ is the single-band Bose-Hubbard Hamiltonian,
$H=H_{\mbox{\scriptsize{hop}}}+H_{\mbox{\scriptsize{int}}}$, which comprises a
hopping term and a local repulsive interaction with strength $U>0$

$H_{\mbox{\scriptsize{hop}}}=-\sum_{j=1}^{L}\left(b_{j+1}^{{\dagger}}b_{j}+h.c.\right),\quad
H_{\mbox{\scriptsize{int}}}=\frac{U}{2}\sum_{j=1}^{L}b_{j}^{{\dagger}2}b_{j}^{2}.$
(4)

On the other hand, $L_{j}$ is a jump operator modeling a two-body loss event
at site $j$, $L_{j}=\sqrt{\frac{\gamma}{2}}b_{j}^{2}$. The operator counting
the total number of particles is $N=\sum_{j=1}^{L}n_{j}$ where
$n_{j}=b_{j}^{{\dagger}}b_{j}$ counts the particles at site $j$. We note in
passing that, unlike the Lieb-Liniger model, the Bose-Hubbard model is no
longer solvable via the Bethe ansatz for $N>2$ [64].

###  II.3 Remarks on the fermionic problem

Dealing with two bosonic particles forces the wavefunction of the system
written in first quantization to be symmetric under exchange of particles.
Two-fermion states are characterized by a symmetric orbital wavefunction under
the condition that the spin wavefunction is antisymmetric under exchange. In
the case of spin-$1/2$ fermions, this is what happens when we consider the two
particles in a singlet state. As we further elaborate in Appendix A, the
results that we are going to derive in the rest of the paper also apply to
fermionic problems with two particles in a singlet state and local two-body
losses. We note in passing that some of our results also apply to the two-
point correlation function in the tight-binding chain of spinless fermions
with dephasing [65, 66].

###  II.4 Effective non-Hermitian dynamics

In general terms, the Lindblad master equations that we introduced in Eqs. (1)
and (3) can be written using $\mathcal{L}$, the Liouvillian superoperator in
Lindblad form acting on the density matrix of the system $\rho$:

$\frac{d\rho}{dt}=\mathcal{L}[\rho],\qquad\mathcal{L}=\mathcal{K}+\mathcal{J}.$
(5)

We divide $\mathcal{L}$ into two superoperators $\mathcal{K}$ and
$\mathcal{J}$: the former conserves the number of particles whereas the latter
decreases it by enforcing the two-body loss. The action of $\mathcal{K}$ can
be written by introducing an effective non-hermitian Hamiltonian
$H_{\mbox{\scriptsize{eff}}}$

$\mathcal{K}[\rho]=-i\left(H_{\mbox{\scriptsize{eff}}}\rho-\rho
H_{\mbox{\scriptsize{eff}}}^{\dagger}\right),$ (6)

which commutes with $N$, the total number operator,
$[H_{\mbox{\scriptsize{eff}}},N]=0$.

In the continuum case, by comparing Eqs. (1) and (2) with Eqs. (5) and (6), we
identify

$H_{\mbox{\scriptsize{eff}}}=\int_{0}^{\ell}\Psi^{\dagger}(x)\left(-\frac{\partial^{2}}{\partial
x^{2}}+c\Psi^{\dagger}(x)\Psi(x)\right)\Psi(x)dx$ (7)

and $\mathcal{J[\rho]}=\int_{0}^{\ell}L(x)\rho L^{\dagger}(x)dx$. The
effective Hamiltonian $H_{\mbox{\scriptsize{eff}}}$ is a Lieb-Liniger
Hamiltonian with a complex two-body interaction strength $c=g/2-i\gamma/4$.
Analogously, in the lattice case we identify

$H_{\mbox{\scriptsize{eff}}}=H-\frac{i}{2}\sum_{j=1}^{L}L_{j}^{\dagger}L_{j}=H_{\mbox{\scriptsize{hop}}}+2u\sum_{j=1}^{L}b_{j}^{{\dagger}2}b_{j}^{2}$
(8)

and $\mathcal{J}[\rho]=\sum_{j=1}^{L}L_{j}\rho L_{j}^{\dagger}$. The effective
Hamiltonian is a single-band Bose-Hubbard Hamiltonian with an effective
complex on-site interaction parameter $2u=U/2-i\gamma/4$. In both the lattice
and continuum cases, we observe that $\mathcal{J}$ satisfies:

$N\mathcal{J}[\rho]=\mathcal{J}[N\rho]-2\mathcal{J}[\rho].$ (9)

When $\rho$ has a well-defined number of particles, this equality can be
interpreted as a mathematical statement of the fact that $\mathcal{J}$
decreases $N$ by two.

Let us now focus on the dynamics of the mean number of particles in the system

$\mathcal{N}(t)=\Tr[\rho(t)N].$ (10)

We assume that initially the system is prepared in a pure two-particle state
$\ket{\Psi_{0}}$: only one loss event can take place and once it occurs the
system is described by the vacuum, a state that contributes trivially to the
observable of our interest. Therefore, $\mathcal{N}(t)$ can be computed via
the non-Hermitian and number-conserving dynamics of
$H_{\mbox{\scriptsize{eff}}}$ (see Ref. [35] and Appendix B for a proof):

$\mathcal{N}(t)=2\innerproduct{\Psi(t)}{\Psi(t)},\qquad\ket{\Psi(t)}=e^{-iH_{\mbox{\scriptsize{eff}}}t}\ket{\Psi_{0}}.$
(11)

The effective non-Hermitian Hamiltonian $H_{\mbox{\scriptsize{eff}}}$ is
responsible for reducing the squared norm of $\ket{\Psi(t)}$ in Eq. (11), and
thus gives a nontrivial time dependence to the average particle number, even
though it commutes with the total particle number, i.e., $[H_{\rm eff},N]=0$.
On the other hand, as we anticipated, the operator $\mathcal{J}$ enforces a
two-body loss event since the jump operators satisfy $[N,L(x)]=-2L(x)$ in the
continuum and $[N,L_{j}]=-2L_{j}$ in the lattice setup.

The dynamics of $\mathcal{N}(t)$ depends on the spectral properties of
$H_{\mbox{\scriptsize{eff}}}$. Since we are focusing on models with local
contact interactions or local losses, we can label the eigenstates with two
momenta $k_{1}$ and $k_{2}$; however, since $H_{\rm eff}$ is non-Hermitian,
these momenta are complex. We denote $\ket{\Psi^{R}_{k_{1}k_{2}}}$ and
$\ket{\Psi^{L}_{k_{1}k_{2}}}$ the right and the left eigenvectors of
$H_{\mbox{\scriptsize{eff}}}$, respectively (see Ref. [56] for a broad
overview of non-Hermitian Hamiltonians):

$\displaystyle H_{\mbox{\scriptsize{eff}}}\ket{\Psi_{k_{1}k_{2}}^{R}}=$
$\displaystyle E_{k_{1}k_{2}}\ket{\Psi_{k_{1}k_{2}}^{R}},$ (12a)
$\displaystyle
H_{\mbox{\scriptsize{eff}}}^{\dagger}\ket{\Psi_{k_{1}k_{2}}^{L}}=$
$\displaystyle E_{k_{1}k_{2}}^{*}\ket{\Psi_{k_{1}k_{2}}^{L}}.$ (12b)

The eigenenergies are complex and will be parametrized as:

$E_{k_{1}k_{2}}=\epsilon_{k_{1}k_{2}}-\frac{i}{2}\Gamma_{k_{1}k_{2}}$ (13)

with $\epsilon_{k_{1}k_{2}}$ and $\Gamma_{k_{1}k_{2}}$ real-valued. In our
model, which includes only lossy events and no gain, $\Gamma_{k_{1},k_{2}}\geq
0$.

We remark here that in the rest of the article we will also make use of the
alternative parametrization in terms of the total momentum $K=k_{1}+k_{2}$ and
of the relative momentum $\frac{\delta}{2}=\frac{k_{1}-k_{2}}{2}$ [67, 68, 69,
35]. Here $K$ will be real because we are considering Hamiltonians
$H_{\mbox{\scriptsize{eff}}}$ that are translationally invariant; it is
expected to take the form $2\pi\mathbb{Z}/\ell$ in the continuum and $2\pi
p/L$ on the lattice, with $p\in\\{0,\cdots,L-1\\}$.

With a few algebraic passages that make use of the completeness relation
$\sum_{(k_{1},k_{2})}\ket{\Psi^{L}_{k_{1}k_{2}}}\bra{\Psi^{R}_{k_{1}k_{2}}}=\mathbb{I}$,
we obtain

$\mathcal{N}(t)=\sum_{(k_{1},k_{2})}\sum_{({q_{1}},{q_{2}})}f\left(k_{1},k_{2},{q_{1}},{q_{2}}\right)e^{-i\left(E_{{q_{1}}{q_{2}}}-E_{k_{1}k_{2}}^{*}\right)t}$
(14)

with

$f\left(k_{1},k_{2},{q_{1}},{q_{2}}\right)=2\bra{\Psi_{0}}\ket{\Psi^{L}_{k_{1}k_{2}}}\bra{\Psi^{R}_{k_{1}k_{2}}}\ket{\Psi^{R}_{{q_{1}}{q_{2}}}}\bra{\Psi^{L}_{{q_{1}}{q_{2}}}}\ket{\Psi_{0}}.$
(15)

This equation cannot be further simplified because, unlike in Hermitian
systems, the right eigenvectors are not orthogonal.

In the long-time limit, the expression (14) can be further simplified: complex
exponentials dephase in time and one can simply retain a single sum:

$\mathcal{N}(t)\simeq\sum_{(k_{1},k_{2})}\tilde{f}\left(k_{1},k_{2}\right)e^{-\Gamma_{k_{1}k_{2}}t},$
(16)

where we have introduced the notations
$\tilde{f}\left(k_{1},k_{2}\right)=f\left(k_{1},k_{2},k_{1},k_{2}\right)$. The
decay of the particles is the sum of many exponential decays, and the most
long-lived states correspond to the energy eigenvalues of
$H_{\mbox{\scriptsize{eff}}}$ with the smallest absolute value of the
imaginary part, i.e., closest to zero; this is the definition of the ADR:

$\lambda_{\rm ADR}=\min_{\begin{subarray}{c}k_{1},k_{2}\\\ \text{such that
}\Gamma_{k_{1}k_{2}}\neq 0\end{subarray}}\frac{\Gamma_{k_{1}k_{2}}}{2}.$ (17)

The ADR is often called the Liouvillian gap in the context of master
equations. Our goal is to show that there is an intermediate time regime where
the decaying exponentials conjure an effective universal power-law decay
susceptible to a scaling collapse for all parameters.

##  III Summary of results

Before entering into the details of our calculations, we present here an
overview of our results.

First of all, in this article we revisit the study of the spectral properties
of the non-Hermitian Hamiltonian $H_{\mbox{\scriptsize{eff}}}$, which has
complex eigenvalues with negative imaginary part. Eigenvalues can be
characterized by two quantum numbers, and for this discussion the best choice
is to consider the total and relative momenta of the two particles.

| continuum, $\forall K$ | lattice, $K\sim 0\;$ | lattice, $K\sim\pi\;$
---|---|---|---
$\lambda_{\rm ADR}$ | $\ell^{-3}$ | $L^{-3}$ | $L^{-5}$
$\alpha$ | 3 | 3 | 5
Table 1:  Universal scaling of the asymptotic decay rate $\lambda_{\rm ADR}$
with the system size in the continuum $\sim\ell^{-\alpha}$ or on a lattice
$L^{-\alpha}$ in the symmetry sector with total momentum $K$; on a lattice,
the scaling depends on whether $K$ is close to $0$ or $\pi$.

In the continuum limit, for any value of $K$ we find states with $\delta\sim
0$ which have an imaginary part of the energy, $-\Gamma_{K,\delta}/2$, scaling
as $\ell^{-3}$. In the lattice, depending on the value of $K$, two scalings
appear: for $K\sim 0$ the scaling is $L^{-3}$, for $K\sim\pi$ the scaling is
$L^{-5}$; this is summarized in Table 1. To the best of our knowledge, the
existence of a $\lambda_{\rm ADR}$ in the lattice scaling as $L^{-5}$ has
never been observed before. In Appendix C we show that this scaling persists
also for larger numbers of particles, $N>2$, in the context of the non-
Hermitian fermionic Hubbard model, that can be studied using the techniques of
integrability [31].

Figure 1: Sketch of the spectrum of $H_{\mbox{\scriptsize{eff}}}$ plotted in
the complex plane, $E_{K,\delta}=\varepsilon_{K,\delta}-i\Gamma_{K,\delta}/2$.
We consider the continuum Hamiltonian (7) and fix $K=0$. The long-lived modes
correspond to the smallest $\Gamma_{K=0,\delta}$ while the fast-decaying modes
correspond to the largest $\Gamma_{K=0,\delta}$. The linearized spectrum drawn
with a dark blue line captures well the long-lived modes.

The general form of the spectrum is sketched in Fig. 1, where we plot a branch
of complex eigenvalues $E_{K,\delta}$ for the continuum non-Hermitian
Hamiltonian in Eq. (7). The sketch highlights the different role played by the
eigenvalues in the dynamics, as it can be intuitively deduced from Eq. (16).
Those characterized by a large $\Gamma_{K=0,\delta}$ are labeled as fast-
decaying modes and they influence the initial transient dynamics. The
universal long-time dynamics is instead determined by the long-lived modes
with small $\Gamma_{K=0,\delta}$. The exact description of the universal
dynamics therefore does not necessitate the use of the correct expression for
all eigenvalues. The technique that we propose in this paper is to approximate
the spectrum with a linearized branch that coincides with the exact one only
in the region with $\Gamma_{K,\delta}\sim 0$. The poor approximation of the
fast decaying modes is solely responsible for the fact that we produce a poor
description of the initial transient dynamics, which is non-universal and is
not of interest here.

Figure 2: Sketch of the four initial states studied in this article. (a)
Continuum case. In the top sketch both particles are localized around
$\frac{{\ell}}{2}$; in the bottom one, one particle is localized around
$\frac{{\ell}}{4}$ and the other one around $\frac{3{\ell}}{4}$. We set
$\sigma=1$ as unit of length. (b) Lattice case. In the top sketch both
particles are at the site $x_{0}$; in the bottom one, one particle is at the
site $x_{0}$ and the other one is at the site $x_{0}+\frac{L-1}{2}$ (in this
article $L$ is odd).

Another set of results concerns the dynamics of the system. Figure 2 shows a
sketch of the initial conditions that we consider, both in the continuum
(panel (a)) and on the lattice (panel (b)). In each case, we consider two
scenarios: two particles spatially localized at the same position (top sketch)
and two particles separated by half the system size (bottom sketch).

Figure 3: Sketch of the dynamics of the mean number of particles as a function
of time $\mathcal{N}(t)$ in log-log scale. The dynamics in the green area is
predicted analytically for any value of interaction strength, loss rate and
system size.

For these initial conditions, we find that the behavior of the average number
of particles $\mathcal{N}(t)$ in Eq. (11) evolves through three different
regimes, as we sketch in Fig. 3 and explain below. First, a non-universal
transient dynamics takes place; this transient depends both on the initial
condition and on the specific parameters of the master equation. Since the
goal of this article is to look for universal behaviors of the lossy gas, in
this study we will not consider it.

After this initial transient, the system enters a regime in which the
universal features of $\mathcal{N}(t)$ do emerge. The main result of this
paper is the characterization of such an intermediate universal regime as
summarized in Table 2. In the continuum, we find a power-law decay
$\mathcal{N}(t)\sim t^{-\beta}$ where $\beta=3/2$ if the two particles are
initially at the same position and $\beta=1/2$ if the two particles are
initially separated by half of the system size. On the lattice, a logarithmic
correction appears: we find $\mathcal{N}(t)\sim
t^{-\frac{3}{2}}\left(\log(\kappa t^{-\frac{1}{2}})+{C}\right)$ if the two
particles are initially at the same position and $\mathcal{N}(t)\sim
t^{-\frac{1}{2}}\left(\log(\kappa t^{-\frac{1}{2}})+{C^{\prime}}\right)$ if
the two particles are initially separated by half of the system size. Here
$\kappa$ is a known function of $U$, $\gamma$ and $L$ while $C$ and
${C^{\prime}}$ are constants; their explicit expressions are presented in Sec.
V.2.

Finally, after this universal long-time dynamics, an exponential decay
$\mathcal{N}(t)\sim e^{-2\lambda_{\rm ADR}t}$ appears that is related to the
ADR. This takes place at a time scaling as $\lambda_{\rm ADR}^{-1}$ that in
our problem increases polynomially with the system size.

|  Initially at the same position |  Initially far apart
---|---|---
Continuum | $t^{-\frac{3}{2}}$ | $t^{-\frac{1}{2}}$
Lattice | $t^{-\frac{3}{2}}\left(\log(\kappa t^{-\frac{1}{2}})+C\right)$ | $t^{-\frac{1}{2}}\left(\log(\kappa t^{-\frac{1}{2}})+{C^{\prime}}\right)$
Table 2: Summary of the universal long-time dynamics of $\mathcal{N}(t)$. Here
$\kappa$ is a known function of $U$, $\gamma$ and $L$ while $C$ and
$C^{\prime}$ are known constants. In the continuum, the number of particles
decays as a power-law. In the lattice, a log correction appears.

##  IV Structure of the spectrum, linearization procedure, and asymptotic
decay rate

In this Section, we diagonalize the effective non-Hermitian Hamiltonian in the
$N=2$ sector. We characterize its eigenvectors, analyze the structure of its
eigenvalues, and linearize the Bethe equations around the most long-lived
modes to compute the ADR.

###  IV.1 Right and left eigenvectors of $H_{\mbox{\scriptsize{eff}}}$ and
associated eigenvalues

We begin by discussing the eigenvectors of $H_{\mbox{\scriptsize{eff}}}$. In
the continuum, we write [60, 61, 62, 63]

$\ket{\Psi^{R}_{k_{1}k_{2}}}=\int_{0}^{\ell}\int_{0}^{\ell}dx_{1}dx_{2}\Psi^{R}_{k_{1}k_{2}}(x_{1},x_{2})\Psi^{\dagger}(x_{1})\Psi^{\dagger}(x_{2})\ket{v}$
(18a) where $\ket{v}$ denotes the vacuum state. On the lattice, we analogously
have [70, 71]
$\ket{\Psi_{k_{1}k_{2}}^{R}}=\sum_{x_{1}=1}^{L}\sum_{x_{2}=1}^{L}\Psi^{R}_{k_{1}k_{2}}(x_{1},x_{2})b_{x_{1}}^{\dagger}b_{x_{2}}^{\dagger}\ket{v}.$
(18b)

In both cases, we define

$\displaystyle\Psi^{R}_{k_{1}k_{2}}$ $\displaystyle(x_{1},x_{2})\propto$
$\displaystyle\begin{cases}e^{i\left(k_{1}x_{1}+k_{2}x_{2}\right)}+s_{k_{1}k_{2}}e^{i\left(k_{2}x_{1}+k_{1}x_{2}\right)},&x_{1}\leq
x_{2}\\\
e^{i\left(k_{1}x_{2}+k_{2}x_{1}\right)}+s_{k_{1}k_{2}}e^{i\left(k_{1}x_{1}+k_{2}x_{2}\right)},&x_{1}>x_{2}\end{cases}$
(19)

where the two-particle S-matrix $s_{k_{1}k_{2}}$ in the continuum case takes
the value:

$s_{k_{1}k_{2}}=\dfrac{k_{1}-k_{2}-ic}{k_{1}-k_{2}+ic},$ (20a) while on a
lattice it reads: $s_{k_{1}k_{2}}=\dfrac{\sin k_{1}-\sin k_{2}-2iu}{\sin
k_{1}-\sin k_{2}+2iu}.$ (20b)

Note that $\Psi^{R}_{k_{1}k_{2}}(x_{1},x_{2})$ in Eq. (19) can be expressed as
a separable function of the center-of-mass position $\frac{x_{1}+x_{2}}{2}$
and of the relative distance $x_{1}-x_{2}$. One then obtains a natural writing
in terms of the center-of-mass momentum $K=k_{1}+k_{2}$ and the relative
momentum $\frac{\delta}{2}=\frac{k_{1}-k_{2}}{2}$, as we anticipated in Sec.
II.4:

$\begin{split}\Psi^{R}_{k_{1}k_{2}}&(x_{1},x_{2})\propto
e^{i\left(k_{1}+k_{2}\right)\frac{x_{1}+x_{2}}{2}}\times\\\
&\Bigl{(}e^{-i\frac{k_{1}-k_{2}}{2}|x_{1}-x_{2}|}+s_{k_{1}k_{2}}e^{i\frac{k_{1}-k_{2}}{2}|x_{1}-x_{2}|}\Bigr{)}.\end{split}$
(21)

Because of the periodic boundary conditions, any pair of quasi-momenta
$(k_{1},k_{2})$ should satisfy the Bethe equations, which in the continuum
read:

$e^{ik_{1}{\ell}}=s_{k_{2}k_{1}}=-\frac{c-i\left(k_{1}-k_{2}\right)}{c+i\left(k_{1}-k_{2}\right)},\qquad
k_{1}+k_{2}=\frac{2\pi}{{\ell}}\mathbb{Z}.$ (22a) On the lattice, they read:
$e^{ik_{1}L}=s_{k_{2}k_{1}}=-\frac{2u-i\left(\sin{k_{1}}-\sin{k_{2}}\right)}{2u+i\left(\sin{k_{1}}-\sin{k_{2}}\right)},\quad
k_{1}+k_{2}=\frac{2\pi}{L}p;$ (22b)

with $p\in\\{0,\cdots,L-1\\}$. Once the momenta $(k_{1},k_{2})$ are determined
with these equations, one eventually obtains the eigenenergies of the
Hamiltonian, whose right eigenvalues in the continuum are:

$E_{k_{1}k_{2}}=k_{1}^{2}+k_{2}^{2},$ (23a) and on a lattice read:
$E_{k_{1}k_{2}}=-2\left(\cos k_{1}+\cos k_{2}\right).$ (23b)

For both setups, we can show that the left eigenvectors of $H_{\rm eff}$ can
be related to the functional form of the right ones (see Appendix D and Ref.
[72] for more details):

$\ket{\Psi^{L}_{k_{1}k_{2}}}=\frac{1}{\bra{\Psi^{R}_{k_{1}k_{2}}}\ket{\Psi^{R}_{k_{1}^{*}k_{2}^{*}}}}\ket{\Psi^{R}_{k_{1}^{*}k_{2}^{*}}}.$
(24)

In the discrete case, when $L$ is even, an eigenstate of total momentum
$K=\pi$ appears that does not take the Bethe-ansatz form in Eq. (19). It is of
the form $\frac{1}{\sqrt{2L}}\sum_{j=1}^{L}e^{i\pi
j}(b_{j}^{\dagger})^{2}\ket{v}$. The state can be interpreted as a bosonic
analogue of the $\eta$-pairing states introduced in Ref. [73] in the context
of the fermionic Hubbard model. This state has energy $U-i\gamma/2$ and
corresponds to a state that loses particles at a rate that does not depend on
$L$ and contributes to the initial non-universal transient dynamics.

When $L$ is even and $K=\pi$, the Bethe equation (22b) has solutions
$k_{1}=(2m-1)\pi/L$ and $k_{2}=\pi-k_{1}$ $(m\in\mathbb{Z})$. The
corresponding energy and wavefunction are $E=0$ and
$\Psi^{R}_{k_{1}k_{2}}(x_{1},x_{2})\propto
e^{i\pi\frac{x_{1}+x_{2}}{2}}\sin\left(\frac{k_{1}-k_{2}}{2}|x_{1}-x_{2}|\right)$,
respectively. Note that the wavefunction is nonvanishing only when
${4m-2\neq(2n-1)L}$ $(n\in\mathbb{Z})$. Such a state is an eigenstate of the
Bose-Hubbard model with zero interaction energy; as a consequence, it will not
experience two-body losses and it will play the role of a nontrivial steady
state of the dynamics (the trivial steady state is the vacuum). Since we are
interested in the universal features of the decaying density, we will mostly
work with odd $L$, so that the system does not have any nontrivial steady
state and we are certified that, asymptotically, the system will be empty
because the vacuum is the only stationary state.

###  IV.2 Linearization of the spectrum in the continuum and asymptotic decay
rate

If we use the center-of-mass and relative momenta in the continuum, Eq. (22a)
becomes

$e^{\frac{i}{2}\delta{\ell}}=(-1)^{p+1}\frac{c-i\delta}{c+i\delta},\qquad
K=\frac{2\pi}{{\ell}}p,\qquad p\in\mathbb{Z},\delta\in\mathbb{C};$ (25)

and accordingly the energy in Eq. (23a) reads

$E_{p\delta}=\frac{1}{2}\delta^{2}+\frac{1}{2}\left(\frac{2\pi}{{\ell}}p\right)^{2}.$
(26)

Figure 4: Eigenenergies $E_{p\delta}$ of $H_{\mbox{\scriptsize{eff}}}$ for two
bosons in the continuum for $p=0$, $1$ and $c=1-i$ and ${\ell}=100$. The other
values of $p$ can be obtained with a rigid translation along the
$\real[E]$-axis. The spectrum continues towards infinite $\real[E]$, with the
$\imaginary[E]$ tending asymptotically towards the value
$4\imaginary[c]/\ell$.

The sets of relative momenta $\delta$ which are solutions of Eq. (25) depend
on $p$ only via its parity. Therefore, the full spectrum can be deduced from
the two branches corresponding to $p=0$ and $p=1$ only; the other branches are
obtained by translation along the $\real[E]$-axis. We solve Eq. (25) with an
iterative numerical method described in Ref. [74] and we use Eq. (26) to
obtain the spectrum of $H_{\mbox{\scriptsize{eff}}}$. In Fig 4, we show the
energy branches for $p=0,1$.

Figure 5: Eigenenergies $E_{p\delta}$ of $H_{\mbox{\scriptsize{eff}}}$ for two
bosons in the continuum for $p\in\\{0,\ldots,7\\}$ and $c=1-i$ and
${\ell}=100$. The colored points are obtained by a numerical solution of Eq.
(25) and the black crosses are the linearized spectrum $\mathcal{E}_{p,q}$.

By inspecting the previous equations in certain regimes, we are also able to
extract some analytic asymptotic behavior. From Eq. (26), and recalling the
definition in Eq. (13), one can easily show that
$\Gamma_{p\delta}=-2\real[\delta]\imaginary[\delta]$ and hence in order to
characterize the properties of the spectrum for long-lived states we perform a
Taylor expansion of Eq. (25) for small $\delta$, namely $|\delta/c|\ll 1$. We
obtain the following linearized momenta (see Appendix E)

$\delta=\begin{cases}\dfrac{4\pi q+4\pi}{{\ell}+4/c}\quad&\mbox{if }p\mbox{ is
odd},\\\ \\\ \dfrac{4\pi q+2\pi}{{\ell}+4/c}\quad&\mbox{if }p\mbox{ is
even};\end{cases}\qquad q\in\mathbb{N}.$ (27)

We substitute this analytical approximate expression into Eq. (26) to obtain
an approximate expression of the spectrum that we call _linearized spectrum_
$\mathcal{E}_{p,q}$. We observe that in the complex energy plane, the
linearized spectrum occupies a line because (see Appendix E)

$\frac{\imaginary[\mathcal{E}_{p=0,q}]}{\real[\mathcal{E}_{p=0,q}]}=\frac{8\imaginary[c]\left(|c|^{2}\ell+4\real[c]\right)}{\left(|c|^{2}\ell+4\real[c]\right)^{2}+16\imaginary[c]^{2}}$
(28)

is a constant that does not depend on $q$. The linearized spectrum is compared
to the exact spectrum in Fig. 5, where it is plotted with black crosses.

Only the eigenvalues corresponding to the smallest $\Gamma_{p\delta}$, and
thus to the largest decay times, are well captured by the linearized spectrum.
The linearization is accurate in the limit ${\ell}|c|\gg 1$ because the most
long-lived modes become accumulation points for the eigenvalues of
$H_{\mbox{\scriptsize{eff}}}$. By substituting Eq. (27) for $p$ even and $q=0$
into Eq. (26) we find an analytical expression for the asymptotic decay rate
$\lambda_{\rm ADR}$ that we defined in Eq. (17):

$\lambda_{\rm
ADR}=\frac{4\pi^{2}\gamma}{\left(\frac{g^{2}}{4}+\frac{\gamma^{2}}{16}\right){\ell}^{3}}+o({\ell}^{-3}).$
(29)

A second regime that we can characterize analytically is the opposite one of
$|\delta/c|\gg 1$, corresponding to the short-lived states which will
ultimately be responsible for the initial transient dynamics. By doing an
analysis of Eq. (25) for $|\delta/c|\gg 1$, we show that asymptotically the
imaginary part of the energy tends to the value $4|\imaginary[c]|/{\ell}$ (see
Appendix F).

###  IV.3 Linearization of the spectrum on the lattice and asymptotic decay
rate

Figure 6: Eigenenergies $E_{p\delta}$ of $H_{\mbox{\scriptsize{eff}}}$ for two
bosons on a lattice of $L=21$ sites for $U=1$ and $\gamma=1$. The branches
characterized by a total momentum $K=2\pi p/L$ are represented with the same
color. The plot highlights the symmetry $p\leftrightarrow L-p$, for which the
eigenvalues coincide exactly.

Using the total and relative momenta for the lattice problem, Eq. (22b)
becomes

$\displaystyle e^{\frac{i}{2}\delta L}=$
$\displaystyle(-1)^{p+1}\frac{2u-2i\sin(\frac{\delta}{2})\cos(\frac{\pi}{L}p)}{2u+2i\sin(\frac{\delta}{2})\cos(\frac{\pi}{L}p)},\qquad\delta\in\mathbb{C};$
$\displaystyle\quad K=$ $\displaystyle\frac{2\pi}{L}p,\qquad
p\in\\{0,\ldots,L-1\\}.$ (30)

The energy reads

$E_{p\delta}=-4\cos(\frac{\delta}{2})\cos(\frac{\pi}{L}p).$ (31)

We first solve numerically the Eq. (30) for each $p$ and we substitute these
solutions into Eq. (31); the resulting spectrum is plotted in Fig. 6. We
observe the appearance of a short-lived branch that does not have an analogue
in the continuum setup and that is composed of states with different total
momenta. This branch is well-separated from the other eigenstates, which have
a smaller imaginary part, and that can be grouped into branches with a fixed
$p$. Figure 7 illustrates that the eigenenergies having the smallest $\Gamma$
are associated to a $\delta$ close to $0$ on one side of the branch and
$\delta$ close to $2\pi$ on the other side.

Figure 7: Eigenenergies $E_{p\delta}$ of $H_{\mbox{\scriptsize{eff}}}$ for two
bosons on a lattice of $L=21$ sites for $U=1$ and $\gamma=1$; the modes
corresponding to the same value of $|\delta|$ are represented with the same
color. We do not display the branch of short-lived states and focus only on
the branches with $0\leq p\leq 10$ (for the other values of $p$ the value of
$|\delta|$ is reverted).

Since our goal is to produce a simple theory that captures well the long-lived
modes, we proceed in a way that is analogous to the previous subsection and
perform a Taylor expansion of Eq. (30) for $|\delta|\ll 1$ or
$|\delta-2\pi|\ll 1$. From a simple observation of Fig. 7, we should expect a
significant dependence on $p$ of the linearized momenta $\delta$, and indeed
we find the following expressions for $\delta\sim 0$:

$\delta=(q+D_{p})\frac{4\pi}{L+\frac{2}{u}\cos(\frac{\pi}{L}p)},\qquad
q\in\mathbb{N};$ (32a) and for $\delta\sim 2\pi$
$\delta=\left(q+D_{p}-\frac{1}{u}\cos(\frac{\pi}{L}p)\right)\frac{4\pi}{L-\frac{2}{u}\cos(\frac{\pi}{L}p)},\quad
q\in\mathbb{N}.$ (32b)

Here, $D_{p}=\frac{1}{2}$ if $p$ is even and $D_{p}=1$ if $p$ is odd. In the
case of $\delta\sim 0$, the most long-lived modes are described by $q=0$. In
the case of $\delta\sim 2\pi$, the most long-lived modes correspond to $q\sim
L/2$, and more precisely to
$q=\left(L-(L\bmod{2})\right)/2-1-\left(p\bmod{2}\right)\left(1-(L\bmod{2})\right)$;
this latter expression is valid both for even and odd $L$. In Fig. 8 we plot
the linearized spectrum $\mathcal{E}_{p,\delta}$ obtained by substituting the
expressions in Eqs. (32) into Eq. (31) for two values of $p$.

Finally, we notice that the asymptotic decay rate $\lambda_{\rm ADR}$ takes
the following simple form:

$\lambda_{\rm
ADR}=\frac{4\pi^{4}\Delta_{L}^{2}\gamma}{\left(\frac{U^{2}}{4}+\frac{\gamma^{2}}{16}\right)L^{5}}+o(L^{-5})$
(33)

with

$\Delta_{L}=\begin{cases}\frac{1}{2}&\text{ for }L\bmod{2}=1;\\\ 2&\text{ for
}L\bmod{4}=0;\\\ 1&\text{ for }L\bmod 4=2.\end{cases}$ (34)

The ADR corresponds to an eigenvalue with $p=\pm(\frac{L}{2}-\Delta_{L})$,
hence to $K\sim\pi$; Eq. (33) is obtained by substituting this value of $p$
and Eq. (32a), with $D_{p}=\frac{1}{2}$ and $q=0$, into Eq. (31). In contrast
to the continuum case, here the scaling is $\lambda_{\rm ADR}\sim L^{-5}$.

Figure 8: Eigenenergies $E_{p\delta}$ of $H_{\mbox{\scriptsize{eff}}}$ for two
bosons on a lattice of $L=21$ sites for $U=1$ and $\gamma=1$. The colored
points are obtained with a numerical solution of Eq. (30) whereas the black
crosses correspond the linearized spectrum represented only for $p=9,10$.

However, if we consider the branches characterized by $p$ close to $0$, and
hence by $K\sim 0$, we obtain:

$\lambda_{\rm
ADR}^{(p)}=\min\limits_{\delta}|\imaginary[E_{p\delta}]|=\frac{4\pi^{2}\gamma}{\left(\frac{U^{2}}{4}+\frac{\gamma^{2}}{16}\right)L^{3}}+o(L^{-3}).$
(35)

This latter equation is the exact analogue of the expression for $\lambda_{\rm
ADR}$ found in the continuum, reported in Eq. (29). In both cases, we have a
scaling with the inverse system size to the third power. The ADR scaling as
$L^{-5}$ is associated to a total momentum at the edge of the Brillouin zone,
and it is thus a lattice effect. The formulas that we just discussed on the
ADR are summarized in Table 1. For completeness, we mention that in Ref. [31]
the authors compute the ADR for a fermionic lattice system at finite density
and show that $\lambda_{\rm ADR}\sim L^{-2}$.

##  V Universal long-time scalings

In this section, we derive the universal expressions for the mean number of
particles at late times for the two initial conditions illustrated in Fig. 2,
both for the continuum system and the lattice.

###  V.1 Universal long-time scaling in the continuum

By using the total and relative momenta, Eq. (16) can be rewritten as

$\mathcal{N}(t)=\sum_{p=-\infty}^{\infty}\sum_{\delta\in\mathcal{S}_{p}}\tilde{f}\left(p,\delta\right)e^{-\frac{i}{2}({\delta}^{2}-{\delta^{*}}^{2})t},$
(36)

where $\mathcal{S}_{p}$ is the set of solutions of the Eqs. (25). If $\delta$
is a solution, then $-\delta$ is also so; thus, when performing actual
calculations, in order to avoid double counting, we always take
$\real[\delta]>0$. If we consider a time much larger than the typical fast-
decay time, $\tau_{\rm fd}=\left(4\frac{|\imaginary[c]|}{\ell}\right)^{-1}$,
then the dynamics is well captured by the linearized spectrum. Indeed,
$\tau_{\rm fd}$ is the typical decay-time of the majority of eigenmodes, which
are characterized by $|\delta/c|\gg 1$ and that are solely responsible for the
short-time transient dynamics. We can thus simply sum over the $\delta$ in the
set $\mathcal{S}^{(0)}_{p}$, corresponding to the values given in Eq. (27);
for consistency, $\tilde{f}(p,\delta)$ has to be expanded to the leading order
in $\delta$.

We detail this procedure for the case where the initial state is
$\Psi_{0}(x_{1},x_{2})\propto
e^{-\left(x_{1}-\frac{{\ell}}{2}\right)^{2}}e^{-\left(x_{2}-\frac{{\ell}}{2}\right)^{2}}$,
illustrated in the top sketch of Fig. 2(a), with both particles localized
around the position ${\ell}/2$ and corresponding to two particles initially at
the same position. The explicit expression for $\tilde{f}(p,\delta)$ is
obtained using the definitions given in Sec. IV and is reported in Appendix G;
an expansion for $\delta$ close to zero, corresponding to the most long-lived
modes, gives:

$\tilde{f}_{0}(p,\delta)\approx
g_{0}({\ell},c)\frac{|\delta|^{4}}{{\delta}^{2}-{\delta^{*}}^{2}}e^{-\frac{1}{4}(\frac{2\pi
p}{{\ell}})^{2}}.$ (37)

The explicit expression of $g_{0}(\ell,c)$ is in Eq. (68) (see Appendix H for
details). We obtain:

$\mathcal{N}(t)\approx
g_{0}({\ell},c)\sum_{p=-\infty}^{\infty}\sum_{\delta\in\mathcal{S}_{p}^{(0)}}\frac{|\delta|^{4}\;e^{-\frac{1}{4}(\frac{2\pi
p}{{\ell}})^{2}}}{{\delta}^{2}-{\delta^{*}}^{2}}e^{-\frac{i}{2}({\delta}^{2}-{\delta^{*}}^{2})t}.$
(38)

By transforming the sums into integrals, it can be shown that long after the
time $\tau_{\rm fd}$ has passed, corresponding to the typical decay time of
the modes characterized by $|\delta/c|\gg 1$, but before the asymptotic
exponential decay is reached (the Liouvillian gap being small but finite) the
average number of particles decays as a power law according to

$\mathcal{N}(t)\approx
h_{0}({\ell},c)\left(\Gamma_{c}t\right)^{-3/2},\qquad\Gamma_{c}=\frac{|\imaginary[c]|}{|c|^{2}}.$
(39)

The explicit expression of the universal rescaling function $h_{0}(\ell,c)$ is
given in Eq. (75) (see Appendix H for details). This is one of the main
results of this article and it shows the existence of a universal behavior in
the time-dependence of the density of the system. Remarkably, the power law
$t^{-3/2}$ and the exponent $\beta=3/2$ do not depend on the complex
interaction constant $c$. This result was not easily anticipated as in
previous works looking for universal behaviors in gases at finite density, the
gas was shown to have two different universal decays for weak and strong
interactions (for instance, for the fermionic case, see Refs. [30, 52]).

In Fig. 9, we compare the approximate result we just derived with the exact
expression of $\mathcal{N}(t)$ computed using the exact formula in Eq. (14).
The comparison is excellent and shows that the linearization of the spectrum
is a powerful technique for extracting the universal decay of the density as
$t^{-3/2}$. Moreover, thanks to the rescaling function $h_{0}(\ell,c)$, we
observe a perfect collapse of the exact numerical data.

Figure 9: Dynamics of the mean number of particles $\mathcal{N}(t)$ rescaled
by the scaling function $h_{0}({\ell},c)$ for the initial state
$\Psi_{0}(x_{1},x_{2})\propto
e^{-\left(x_{1}-\frac{{\ell}}{2}\right)^{2}}e^{-\left(x_{2}-\frac{{\ell}}{2}\right)^{2}}$
and for ${\ell}=500$. The solid lines correspond to the exact computations
from Eq. (14) for three different values of $c$. The dashed black line
corresponds to the universal expression given by Eq. (39); the collapse of the
data is excellent.

We now consider a different initial state $\Psi_{1}(x_{1},x_{2})\propto
e^{-(x_{1}-\frac{{\ell}}{4})^{2}}e^{-(x_{2}-\frac{3{\ell}}{4})^{2}}+e^{-(x_{2}-\frac{{\ell}}{4})^{2}}e^{-(x_{1}-\frac{3{\ell}}{4})^{2}}$,
illustrated in the bottom sketch of Fig. 2(a). One particle is localized
around $\frac{{\ell}}{4}$ and the other one is localized around
$\frac{3{\ell}}{4}$, so that the particles are separated by half of the system
size. For this initial condition, we obtain the following expansion of
$\tilde{f}(p,\delta)$:

$\tilde{f}_{1}\left(p,\delta\right)=g_{1}({\ell},c)\frac{|\delta|^{2}}{{\delta}^{2}-{\delta^{*}}^{2}}e^{-\frac{1}{4}(\frac{2\pi
p}{{\ell}})^{2}}.$ (40)

The different behavior of $\tilde{f}_{1}$ for small $\delta$ with respect to
$\tilde{f}_{0}$ is responsible for a different universal long-time decay.
Following the same steps outlined previously, we find, for the same range of
$t$ corresponding to the intermediate algebraic regime,

$\mathcal{N}(t)\approx
h_{1}({\ell})\left(\Gamma_{c}t\right)^{-1/2},\qquad\Gamma_{c}=\frac{|\imaginary[c]|}{|c|^{2}}.$
(41)

The explicit expressions of $g_{1}(\ell,c)$ and $h_{1}(\ell)$ are respectively
presented in Eqs. (76) and (77) in Appendix H. In Fig. 10 we plot the function
$\mathcal{N}(t)$ given by the exact expression in Eq. (14) and we compare it
with the approximate analytical expression given by Eq. (41) to highlight the
universal power-law scaling $t^{-1/2}$. The collapse of the data onto the
universal curve is apparent.

Figure 10: Dynamics of the mean number of particles $\mathcal{N}(t)$ rescaled
by the scaling function $h_{1}({\ell},c)$ for the initial state
$\Psi_{1}(x_{1},x_{2})\propto
e^{-(x_{1}-\frac{{\ell}}{4})^{2}}e^{-(x_{2}-\frac{3{\ell}}{4})^{2}}+x_{1}\leftrightarrow
x_{2}$ and for ${\ell}=500$. The solid lines correspond to the exact
computations from Eq. (14) for three different values of $c$. The dashed black
line corresponds to the universal expression given by Eq. (41); the collapse
of the data is excellent.

We can briefly estimate for which range of times the universal scalings that
we just presented are valid. In general, the decay of particles is
characterized by a first non-universal dynamics that is determined by the
fast-decaying modes, which decay on a typical time scale $\tau_{\rm fd}$, and
that constitute a lower bound for the appearance of the universal behavior. At
very late times, the universal power-law scaling is then followed by an
exponential decay $\mathcal{N}(t)\sim e^{-2\lambda_{\rm ADR}t}$. Concluding,
we can thus estimate that the expressions in Eqs. (39) and (41) are valid for
$\tau_{\rm fd}\ll t\ll\frac{1}{\lambda_{\rm ADR}}$.

####  V.1.1 Interpolating the two initial conditions

Figure 11: Dynamics of the mean number of particles $\mathcal{N}(t)$ for the
initial state $\Psi_{1-s}(x_{1},x_{2})\propto
e^{-(x_{1}-\frac{{\ell}}{4}-s\frac{{\ell}}{4})^{2}}e^{-(x_{2}-\frac{3{\ell}}{4}+s\frac{{\ell}}{4})^{2}}+x_{1}\leftrightarrow
x_{2}$ for ${\ell}=500$. We consider one interaction strength $c=1-i$ and
various $s\in\\{0.2,0.5,0.8,0.95\\}$. The solid curves correspond to the exact
computation from Eq. (14); the dashed and dashed-dotted lines correspond to
the algebraic decays discussed in the main text. For $s\sim 1$ the density of
bosons $\mathcal{N}(t)$ crosses over between two algebraic decays.

As it can be understood intuitively, the scaling $\mathcal{N}(t)\sim t^{-1/2}$
is not peculiar to the situation where the two particles are initially
separated by exactly half of the system size. The same scaling emerges for all
initial states consisting of two non-overlapping Gaussians, i.e. separated by
a distance greater than their typical width. In Fig. 11, we show the mean
number of particles computed with the exact expression given by Eq. (14) for
the initial states

$\Psi_{1-s}(x_{1},x_{2})\propto
e^{-(x_{1}-\frac{{\ell}}{4}-s\frac{{\ell}}{4})^{2}}e^{-(x_{2}-\frac{3{\ell}}{4}+s\frac{{\ell}}{4})^{2}}+x_{1}\leftrightarrow
x_{2}$ (42)

for various values of the parameter $s$ in the range $[0,1]$. This
wavefunction places the two particles at an initial distance of
$\frac{\ell}{2}\left(1-s\right)$ and thus interpolates between the two
particles at the same position ($s=1$) and the two particles at their largest
distance ($s=0$).

At small $s$, the particles are initially far apart, their overlap is
negligible and the initial condition to any purpose completely
indistinguishable from the $s=0$ one. We observe numerically that
$\mathcal{N}(t)$ decays as $t^{-1/2}$. We can thus conclude that this power
law characterizes the algebraic decay of $\mathcal{N}(t)$ as long as the two
particles are initially far apart.

The situation is different when $s$ is close to $1$; initially the two
particles are not far apart (as compared to their typical widths); due to the
unitary dynamics, the two wave packets spread and can develop a large overlap,
leading to a situation which will not be distinguishable from the one
consisting of two perfectly overlapping Gaussian wave packets at the initial
time. We numerically observe this to be the case in Fig. 11: after an initial
$\mathcal{N}(t)\propto t^{-1/2}$ decay (which gets shorter as $s$ increases),
a second algebraic region $\mathcal{N}(t)\propto t^{-3/2}$ emerges.

Importantly, even though as we discussed there may be a competition between
the two asymptotic algebraic decays, no new power law emerges in this case,
showcasing the robustness of the power laws we characterized to a perturbation
of the system’s initial conditions.

####  V.1.2 Experimental relevance

So far we studied very large system sizes in order to let the universal decay
emerge in full clarity; however, since we are considering only two particles,
a smaller ${\ell}$ could be more experimentally relevant. For the
linearization to be valid, the inequality ${\ell}|c|\gg 1$ should be
satisfied; we can therefore reduce the system size by introducing larger
dissipation rates. In Fig. 12, we show that the scaling $\mathcal{N}(t)\sim
t^{-1/2}$ still emerges for a smaller system of length ${\ell}=30$; the
universal behavior is more visible for large values of $|c|$, as expected.
Since in this situation the universal decay appears for significantly shorter
times and for higher values of the mean number of particles with respect to
the examples studied so far, we believe that the universal decay that we are
proposing is not completely outside the possibilities of an experimental
observation.

Figure 12: Dynamics of the mean number of particles $\mathcal{N}(t)$ for
$\Psi_{0}(x_{1},x_{2})\propto
e^{-(x_{1}-\frac{{\ell}}{4})^{2}}e^{-(x_{2}-\frac{3{\ell}}{4})^{2}}+x_{1}\leftrightarrow
x_{2}$ for a small setup of length ${\ell}=30$. The solid lines correspond to
the exact computation from Eq. (14). The dashed black line corresponds to the
universal expression given by Eq. (41). Even if the size of the system is
rather small, the universal behavior can be clearly distinguished, especially
for large values of $|c|$.

###  V.2 Universal long-time scaling on the lattice

We now move to the lattice geometry and compute $\mathcal{N}(t)$ for two
initial states paralleling those discussed in the continuum: in the first
initial condition, the two particles are at the same site; in the second
initial condition, the two particles are at distant sites. As we already said,
we take $L$ odd in order to avoid the presence of nontrivial steady states.

In order to see the universal long-time scaling, we assume that $L|u|\gg 1$
and we consider times such that $L^{3}\ll\Gamma_{u}t\ll L^{5}$ with

$\Gamma_{u}=\frac{256\pi^{2}|\imaginary(2u)|}{|2u|^{2}}=\frac{1024\pi^{2}\gamma}{4U^{2}+\gamma^{2}}.$
(43)

Under this assumption, the dynamics is well captured by the long-lived modes
and the relative momenta can be linearized as in Eqs. (32a) and (32b).

We first consider the state in which two particles are initially on the same
site. In this case, we compute

$\mathcal{N}(t)\approx\sqrt{\frac{2L^{3}}{\pi\left(\Gamma_{u}t\right)^{3}}}\frac{2\Gamma_{u}}{\gamma}\times\left[\log{\left(2\sqrt{\frac{L^{5}}{\Gamma_{u}t}}\right)}+C\right].$
(44)

We find that $\mathcal{N}(t)$ scales as $t^{-3/2}$, exactly as in the
continuum setup that we presented in Eq. (39). However, unlike the continuum
case, we find a logarithmic correction to the scaling. $C$ is a universal
constant that evaluates approximately to $1.75$, whose exact expression is in
Eq. (94) in Appendix I.

On the other hand, if we consider an initial state where one particle is at
site $x_{0}$ and the other is at $x_{0}+(L-1)/2$, we obtain

$\mathcal{N}(t)\approx
16\sqrt{\frac{L}{2\pi\Gamma_{u}t}}\times\left[\log{\left(2\sqrt{\frac{L^{5}}{\Gamma_{u}t}}\right)}+C^{\prime}\right].$
(45)

In this case $\mathcal{N}(t)$ scales as $t^{-1/2}$, analogously to the
continuum problem presented in Eq. (41). Here $C^{\prime}$ is a universal
constant that evaluates approximately to $1.52$, we refer the interested
reader to Eq. (93). in Appendix I for the detailed expression. We sketch the
derivation of Eq. (45) in Appendix I; Eq. (44) can be derived in a completely
analogous way, but the calculations are slightly more involved.

The presence of the logarithmic correction is robust and is a direct
consequence of the fact that the spectrum features two different ADR for the
symmetry sectors with total momentum $K\sim 0$ and for $K\sim\pi$, as we
discussed in Sec. IV.3 and in Table 1. It can be misleading if one tries to
fit the numerical data without prior knowledge of the universal form of
$\mathcal{N}(t)$. Indeed, a fit using a power-law decay of the form
$\mathcal{N}(t)\sim t^{-\beta}$ might seem to properly describe the data;
however, by doing so, one recovers a $\beta$ that depends on $U$, $\gamma$ and
$L$. This blind approach misses the universal behavior of the decay of the
density of particles and could be relevant also in other contexts, such as
those discussed in Refs. [52, 15].

Figure 13: Dynamics of the mean number of particles $\mathcal{N}(t)$ for the
initial state
$\ket{\Psi_{2}}=b_{x_{0}}^{\dagger}b_{x_{0}}^{\dagger}\ket{v}/\sqrt{2}$ and
for $L=121$ sites. The solid curves correspond to the exact computation. The
dashed black line corresponds to the universal expression given by Eq. (44);
the collapse of the data is excellent.  Figure 14: Dynamics of the mean number
of particles $\mathcal{N}(t)$ for the initial state
$\ket{\Psi_{3}}=b_{x_{0}}^{\dagger}b_{x_{0}+(L-1)/2}^{\dagger}\ket{v}$ and for
$L=121$ sites. The solid curves correspond to the exact computation. The
dashed black line corresponds to the universal expression given by Eq. (45);
the collapse of the data is excellent.

In Figs. 13 and 14, we compare the analytical universal scaling functions in
Eqs. (44) and (45) with numerical results. For both initial conditions, we
have a non-universal dynamics at first. The universal power-law scaling with
logarithmic correction can be seen up to times scaling as $\sim
L^{5}/\Gamma_{u}$, which is the scaling of $\lambda_{\rm ADR}^{-1}$ for total
momentum $K\sim\pi$. For longer times, the density decays exponentially with a
time scale determined by this ADR. The collapse of the exact numerical data
onto the universal scaling functions that we just derived is excellent and
further proves the usefulness of the linearization of the spectrum for the
long-lived modes. The results that we just discussed are summarized in Table
2.

####  V.2.1 Crossover from the transient dynamics to the universal long-time
behavior in the weakly-interacting and weakly-lossy regime

The results that we found on the universal late-time properties of
$\mathcal{N}(t)$ are particularly surprising when we consider the weakly-
interacting and weakly-lossy region, for which a previous article has found an
exponential decay in time for a similar setup characterized by two-body losses
[35]:

$\mathcal{N}(t)=2e^{-2\gamma t/L}.$ (46)

The applicability of this result to our problem is corroborated by the study
of the spectrum of $H_{\mbox{\scriptsize{eff}}}$. As $\gamma$ and $U$
decrease, the eigenmodes of the effective non-Hermitian Hamiltonian accumulate
on the horizontal line $\imaginary[E]=-\gamma/L$ as shown in Fig. 15 for the
lattice case. The four panels refer to different values of $U$ and $\gamma$
and the tendency that we described appears rather clearly. The presence of
other eigenvalues with an imaginary part of the energy that is closer to zero
suggests that at late times the dynamics of $\mathcal{N}(t)$ will not be
described by Eq. (46).

Figure 15: Eigenergies of $H_{\mbox{\scriptsize{eff}}}$ for $L=21$ and various
values of $U$ and $\gamma$. The dashed black line on which the modes
accumulate for $U,\gamma\sim 0$ corresponds to $\imaginary[E]=-\gamma/L$. The
collapse of the eigenvalues onto the black dashed line is evident, even if
outlying eigenvalues describing modes decaying more slowly are always present.

In Fig. 16, we compare the exact dynamics to the analytical form given by Eq.
(46) for the initial state consisting of two particles separated by half of
the system size. As $\gamma$ and $U$ decrease, not only the analytical curve
of Eq. (46) approaches the exact computation, but also the time interval over
which this agreement is good increases and the crossover to the universal
behavior appears later. These plots explain the compatibility of our universal
result with other approaches that could be developed to discuss the transient
dynamics of the gas. The techniques presented in this article focus on the
asymptotically late-time dynamics that precede the exponential decay in time
dictated by the ADR: this does not exclude the possibility of finding other
analytical behaviors as a function of time for the initial stages of the
dynamics.

Figure 16: Crossover in time from the early-time transient dynamics in Eq.
(46) to the universal late-time behaviour for both $U$ small and $\gamma$
small. Simulations are carried out on the lattice for $L=121$ sites and the
initial state is
$\ket{\Psi_{3}}=b_{x_{0}}^{\dagger}b_{x_{0}+(L-1)/2}^{\dagger}\ket{v}$; the
solid curves correspond to the exact computation. The black dashed lines
correspond to the universal expression given by Eq. (45); the red dashed-
dotted line corresponds to the analytical form given by Eq. (46), that only
describes the early-time behaviour of $\mathcal{N}(t)$.

###  V.3 Quantum Zeno effect

For the lattice case, the emergence of the time scale $\Gamma_{u}^{-1}$ in
Eqs. (44) and (45) can be seen as a manifestation of the continuous quantum
Zeno (QZ) effect [25, 44, 52], a phenomenon that is here invoked because
strong losses suppress coherent processes (in our case, the hoppings that
change the number of double occupancies). The suppression of motion within the
lattice often slows down particle losses, and in our discussion we will speak
of QZ behavior whenever this takes place. For instance, if $\gamma\ll U$ then
$\Gamma_{u}\sim\gamma$, but if $\gamma\gg U$, then $\Gamma_{u}\sim 1/\gamma$.

In Fig. 17, we illustrate this effect by studying $\mathcal{N}(t)$ at a fixed
time $t$ and varying the two-body loss rate; we speak of QZ behaviour when the
curve is increasing for increasing loss rate. We begin our analysis by
considering the lattice situation with the two particles that are initially
well-separated; as illustrated in Fig. 17 (d), the number of particles has a
non-monotonic behavior as a function of $\gamma$. This behavior is completely
inherited from the properties of $\Gamma_{u}$, that is the only time scale
appearing in Eq. (45).

Figure 17: Top panels: mean number of particles $\mathcal{N}(t)$ for the
continuum setup as a function of $|\imaginary[c]|$ for ${\ell}=500$ and at
time $t={\ell}^{2}$; IC stands for “initial condition”. (a) The particles are
initially well-separated and the blue curve is obtained from Eq. (41);
$\mathcal{N}(t)$ is a non-monotonic function of $|\imaginary[c]|$, which is a
manifestation of the quantum Zeno (QZ) effect. (b, c) The particles are
initially at the same position; the orange and green curves are obtained from
Eq. (39). On the panel (b) we have $|c|\ll 1$ and no QZ effect emerges; on the
panel (c) we have $|c|\gg 1$ and the QZ effect can be seen for
$|\imaginary[c]|\gg\real[c]$. Bottom panels: mean number of particles
$\mathcal{N}(t)$ for the lattice setup as a function of $\gamma$ for $L=121$
and at time $t=L^{4}$. (d) The particles are initially well-separated and the
blue curve is obtained from Eq. (45). As in panel (a), $\mathcal{N}(t)$ is a
non-monotonic function of $\gamma$, which is a manifestation of the quantum
Zeno effect. (e) The particles are initially at the same site and the orange
curve is obtained from Eq. (44); $\mathcal{N}(t)$ is a decreasing function of
$\gamma$ which means that no quantum Zeno effect emerges.

However, if the two particles are initially on the same site, an initial
exponential decay with a rate $\gamma$ is expected and the number of particles
decreases rapidly at first. The rapid loss of particles cannot be prevented by
considering large values of $\gamma$. Indeed, the initial state with a double
occupancy is sensitive to loss and even if in the regime of large $\gamma$ the
hopping is suppressed, for this initial condition no QZ behaviour is expected.
Mathematically, the absence of QZ behaviour is enforced by the prefactor
$\Gamma_{u}/\gamma$ in Eq. (44), that decreases rapidly for increasing
$U^{2}+\gamma^{2}/4$. With a few algebraic passages, one finds that this
enforces the strictly decreasing behavior of $\mathcal{N}(t)$ as a function of
$\gamma$ at a given time, as shown in Fig. 17 (e). In this situation, thus,
increasing $\gamma$ makes the gas more dissipative and no QZ behaviour
appears.

For the continuum case, if the particles are initially well-separated,
$\Gamma_{c}^{-1}$ is clearly the unique time-scale of the dynamics since the
scaling function $h_{1}({\ell})$ is independent of $c$, as we can see in Eq.
(41). Thus, the situation is analogous to the lattice case; the QZ behavior
occurs when $|\imaginary[c]|\gg\real[c]$, as we show in Fig. 17 (a).

However, for the initial state consisting of two perfectly overlapping
Gaussian wave packets, the scaling function $h_{0}({\ell},c)$, appearing in
Eq. (39) and explicitly written in Eq. (75), depends on $c$; this makes the
discussion slightly more complicated. If $|c|\ll 2\pi$, the number of
particles always decreases faster when $|\imaginary[c]|$ increases, similarly
to the lattice case. On the contrary, if $|c|\gg 2\pi$ then the QZ behavior
emerges in the usual limit $|\imaginary[c]|\gg\real[c]$. We illustrate these
various cases for two particles in the continuum at the same position
initially in Fig. 17 (b) and (c).

##  VI Conclusions and perspectives

In this article, we characterized the universal properties of the dynamics of
a quantum system composed of two bosons that undergo a two-body loss. We have
focused on two quantities, the scaling of the asymptotic decay rate with the
linear system size, $\lambda_{\rm ADR}\sim L^{-\alpha}$, and the late-time
dynamics of the mean number of particles, $\mathcal{N}(t)\sim t^{-\beta}$. Our
results have been carried out with analytical techniques and have shown a
nontrivial dependence both on the geometry of the system (a continuum one-
dimensional setup or a lattice) and, for what concerns $\beta$, on the initial
state (with overlapping or separated particles).

Our main results are summarized in Tables 1 and 2 and we briefly recapitulate
them here. Concerning the $\lambda_{\rm ADR}$, in the continuum $\alpha=3$;
this result is obtained also for lattice setups if we restrict the study to
the symmetry sector with total momentum equal to zero. On the lattice,
however, eigenmodes with total momentum close to the edge of the Brillouin
zone display an asymptotic decay rate that decreases faster with the linear
system size, characterized by $\alpha=5$.

The universal properties of the decay of the number of particles in the system
depend on the initial conditions. When particles initially occupy the same
position, we compute $\beta=3/2$, whereas when the particles are well
separated we obtain $\beta=1/2$. No other exponent has been found by changing
the relative position of the two particles, but rather a crossover between the
two: when the particles are separated but not too distant, the decay time
follows $\beta=1/2$ at short times and then turns into $\beta=3/2$.

The study of $\beta$ on the lattice has revealed the presence of a logarithmic
correction which is the consequence of the aforementioned two-fold spectral
structure, characterized by an asymptotic decay rate whose scaling depends on
the momentum symmetry sector. The logarithmic correction cannot be easily
fitted from the numerical data, which are compatible with a power-law decay
characterized by a non-universal exponent. Similar non-universal behaviors
have been numerically found in other setups [52] and this work suggests
investigating whether these results could be interpreted as logarithmic
corrections.

The main outlook of this work is the possibility of performing these
calculations for an initial number of particles that has finite density. In
this case, the dynamics cannot be described using the tools of a non-Hermitian
dynamics but one has to employ the full Lindblad dynamics. Yet, our work
suggests that the study of the spectral properties of the Lindblad
superoperator could be linked to the universal exponent $\beta$. The
investigation of such a connection using analytical or numerical tools is an
exciting perspective that we leave for future work.

Finally, we conclude by commenting on the fact that this work is an example of
the rich properties that can be associated to a many-particle non-Hermitian
dynamics. It would be interesting to extend this study to the many-particle
case using the techniques of Bethe ansatz and to find the universal decay of
population in the system due to two-body losses using a non-Hermitian
dynamics. Even if it is known that the study cannot describe the real dynamics
of many particles under two-body losses, it remains of the utmost interest to
assess which properties can be correctly captured with this approach, that is
technically less demanding.

###### Acknowledgements.

We acknowledge insightful discussions with R. Fazio, B. Laburthe-Tolra, B.
Pasquiou and M. Robert-de-Saint-Vincent. L.M. thanks A. Biella, J. De Nardis
and L. Rosso for enlightening discussions on previous works on related topics.
This work was carried out in the framework of the joint Ph. D. program between
the CNRS and the University of Tokyo. This work is supported by the ANR
project LOQUST ANR-23-CE47-0006-02. H.Y. was supported by JSPS KAKENHI Grant-
in-Aid for JSPS fellows Grant No. JP22J20888, the Forefront Physics and
Mathematics Program to Drive Transformation, and JSR Fellowship, the
University of Tokyo. H.K. was supported by JSPS KAKENHI Grants No. JP23H01086,
No. JP23H01093, and MEXT KAKENHI Grant-in-Aid for Transformative Research
Areas A “Extreme Universe” (KAKENHI Grant No. JP21H05191). A.N. acknowledges
financial support from Fondazione Angelo dalla Riccia, LoCoMacro 805252 from
the European Research Council, and LPTMS for warm hospitality.

##  Appendix A Details about two fermions in a spin singlet on a lattice

In this Appendix, we examine in more detail the case of two spin-$1/2$
fermions in a singlet trapped on a lattice. We denote
$c_{j,\sigma}^{{\dagger}}$ and $c_{j,\sigma}$ the creation and annihilation
operators of a fermion with spin $\sigma\in\\{\uparrow,\downarrow\\}$ at site
$j\in\\{1,\ldots,L\\}$; $n_{j,\sigma}=c_{j,\sigma}^{{\dagger}}c_{j,\sigma}$
counts the number of particles at site $j$ with spin $\sigma$ and the operator
counting the total number of particles is
$N=\sum_{j=1}^{L}\sum_{\sigma=\uparrow,\downarrow}n_{j,\sigma}$.

The master equation of this fermionic problem is given by Eq. (3) with $H$
being the ${\rm SU}(2)$ Fermi-Hubbard Hamiltonian [67]
$H=H_{\mbox{\scriptsize{hop}}}+H_{\mbox{\scriptsize{int}}}$ with

$\displaystyle H_{\mbox{\scriptsize{hop}}}$
$\displaystyle=-J\sum_{j=1}^{L}\sum_{\sigma=\uparrow,\downarrow}\left(c_{j+1,\sigma}^{{\dagger}}c_{j,\sigma}+h.c.\right),$
(47) $\displaystyle H_{\mbox{\scriptsize{int}}}$
$\displaystyle=U\sum_{j=1}^{L}n_{j,\uparrow}n_{j,\downarrow}.$ (48)

The jump operators are

$L_{j}=\sqrt{\gamma}c_{j,\downarrow}c_{j,\uparrow}.$ (49)

Thus, the effective non-Hermitian Hamiltonian is

$H_{\mbox{\scriptsize{eff}}}=H-\frac{i}{2}\sum_{j=1}^{L}L_{j}^{\dagger}L_{j}=H_{\mbox{\scriptsize{hop}}}+H_{\mbox{\scriptsize{int}}}^{\prime}$
(50)

with

$H_{\mbox{\scriptsize{int}}}^{\prime}=\left(U-i\frac{\gamma}{2}\right)\sum_{j=1}^{L}n_{j,\uparrow}n_{j,\downarrow}.$
(51)

We also note that $H_{\mbox{\scriptsize{eff}}}$ is SU(2) symmetric, and in
particular, commutes with the total spin operator $\hat{S}^{2}$.

For this fermionic system, the dynamics of the mean number of particle
$\mathcal{N}(t)$ is the same as for the bosonic lattice system described in
Sec. II.2. Indeed, we can always decompose a fermionic two-body initial state
into a spin singlet part and a spin triplet part $\ket{\Psi_{0}}$ =
$\ket{\Psi_{0}^{T}}$ \+ $\ket{\Psi_{0}^{S}}$ with
$\bra{\Psi_{0}^{T}}\ket{\Psi_{0}^{S}}=0$ because the two states are
eigenvectors of the Hermitian operator $S^{2}$ with different eigenvalues. The
spin triplet part is associated with an antisymmetric orbital part, which
cannot contain any double occupancy and thus does not yield any two-fermion
decay, as exemplified by
$H_{\mbox{\scriptsize{int}}}^{\prime}\ket{\Psi_{0}^{T}}=0$. As a consequence,
using the spin symmetry of the Hamiltonian, we can single out two
contributions to the time evolution of the average number of particles in Eq.
(11):

$\innerproduct{\Psi(t)}{\Psi(t)}=\innerproduct{\Psi_{0}^{T}}{\Psi_{0}^{T}}+\innerproduct{\Psi^{S}(t)}{\Psi^{S}(t)}$
(52)

where we defined

$\ket{\Psi^{S}(t)}=e^{-iH_{\mbox{\scriptsize{eff}}}t}\ket{\Psi_{0}^{S}}.$ (53)

This shows that the spin triplet part of $\ket{\Psi_{0}}$ has a trivial, time-
independent, contribution to $\innerproduct{\Psi(t)}{\Psi(t)}$. For this
reason, it is more interesting to consider $\ket{\Psi_{0}}$ to be a spin
singlet state, meaning that $\ket{\Psi_{0}^{T}}=0$ and
$\ket{\Psi_{0}}=\ket{\Psi_{0}^{\mbox{\scriptsize{orb}}}}\otimes(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})/\sqrt{2}$,
with the orbital part $\ket{\Psi_{0}^{\mbox{\scriptsize{orb}}}}$ being
symmetric under particle exchange. This makes the fermionic problem completely
equivalent to the study of two spinless bosons.

##  Appendix B Derivation of Eq. (11)

In this Appendix, we derive Eq. (11). Using the decomposition (5),
$\mathcal{N}(t)$ can be calculated as

$\mathcal{N}(t)=\Tr[Ne^{\mathcal{L}t}[\rho(0)]]=\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\Tr\left[N(\mathcal{K}+\mathcal{J})^{n}[\rho(0)]\right].$
(54)

Then, $\mathcal{K}$ conserves the number of particles and $\mathcal{J}$
decreases the number of particles by two. Since $N\rho(0)=2\rho(0)$ for the
initial state, only the terms in the form
$N\mathcal{K}^{n}[\rho(0)]=2\mathcal{K}^{n}[\rho(0)]$ are nonvanishing in the
right-hand side of Eq. (54). Therefore,

$\displaystyle\mathcal{N}(t)$
$\displaystyle=2\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\Tr\left[\mathcal{K}^{n}[\rho(0)]\right]=2\Tr\left[e^{\mathcal{K}t}[\rho(0)]\right]$
$\displaystyle=2\Tr\left[e^{-iH_{\mathrm{eff}}t}\ket{\Psi_{0}}\bra{\Psi_{0}}e^{iH^{\dagger}_{\mathrm{eff}}t}\right]$
$\displaystyle=2\innerproduct{\Psi(t)}{\Psi(t)},$ (55)

where $\ket{\Psi(t)}=e^{-iH_{\mathrm{eff}}t}\ket{\Psi_{0}}$.

##  Appendix C The $L^{-5}$ scaling in the case of $N>2$ fermions on a lattice

In this Appendix, we show the existence of eigenvalues whose imaginary part
scales with $L^{-5}$ in the case of $N>2$ spin-$1/2$ fermions with $M$ down
spins on $L$ lattice sites. We consider the dilute case: we fix $N$, $M$ and
take $L$ large. We also assume that $|u|=|U/4-i\gamma/8|$ is large. When $|u|$
is large enough, the spin-charge separation occurs [75], and the eigenvalues
of $H_{\mathrm{eff}}$ are given by [67]

$\displaystyle E$
$\displaystyle=-2\sum_{j=1}^{N}\cos\left(k_{j}\right)+\frac{2e_{s}}{u}\left[\frac{N}{L}\sum_{j=1}^{N}\sin^{2}\left(k_{j}\right)-\frac{1}{L}\left(\sum_{j=1}^{N}\sin\left(k_{j}\right)\right)^{2}\right]+O(1/u^{2}),$
(56)

where

$e_{s}=-\frac{2}{N}\sum_{l=1}^{M}\frac{1}{\Lambda_{l}^{2}+1}$ (57)

and $k_{j}$ and $\Lambda_{l}$ are charge momenta and spin rapidities,
respectively, which satisfy

$\displaystyle k_{j}=\frac{2\pi
n_{j}}{L}+\frac{P_{s}}{L},\qquad\prod_{l=1}^{M}\frac{\Lambda_{l}-i}{\Lambda_{l}+i}=e^{iP_{s}}$
(58)

with $n_{j}\in\\{0,\ldots,L-1\\}$ and $P_{s}=\frac{2\pi m}{N}$
$(m\in\\{0,\ldots,L-1\\})$. Then, the imaginary part of $E$ reads

$\displaystyle\imaginary E$
$\displaystyle=\imaginary\left(\frac{2e_{s}}{uL}\right)\sum_{1\leq j<l\leq
N}(\sin k_{j}-\sin k_{l})^{2}+O(1/u^{2})$ (59)
$\displaystyle=\imaginary\left(\frac{8e_{s}}{uL}\right)\sum_{1\leq j<l\leq
N}\cos^{2}\left(\frac{k_{j}+k_{l}}{2}\right)\sin^{2}\left(\frac{k_{j}-k_{l}}{2}\right)+O(1/u^{2}).$
(60)

If we choose $k_{j}$ as $k_{j}=\frac{\pi}{2}-p_{j}\frac{\pi}{L}$, where all
$p_{j}$ are distinct constants, we have

$\displaystyle\imaginary E$
$\displaystyle=\imaginary\left(\frac{8e_{s}}{uL}\right)\sum_{1\leq j<l\leq
N}\sin^{2}\left[\frac{\pi(p_{j}+p_{l})}{2L}\right]\sin^{2}\left[\frac{\pi(p_{j}-p_{l})}{2L}\right]+O(1/u^{2})$
(61)
$\displaystyle=\frac{1}{L^{5}}\imaginary\left(\frac{\pi^{4}e_{s}}{2u}\right)\sum_{1\leq
j<l\leq N}\left(p_{j}^{2}-p_{l}^{2}\right)^{2}+O(1/u^{2})+O(1/L^{7}).$ (62)

Since $p_{j}$, $N$ and $e_{s}$ are constants with respect to $L$, $\imaginary
E$ scales with $L^{-5}$. Note that this scaling does not hold when the
fermionic density $N/L$ is fixed, in which case $\imaginary E\propto L^{-2}$
[31].

##  Appendix D Derivation of Eq. (24)

In this Appendix, we derive Eq. (24); we make also reference to Ref. [72]
where this link between left and right eigenvectors was already present.
First, we note that $(k_{1},k_{2})$ is a solution of the Bethe equations
associated to the diagonalization of $H_{\mbox{\scriptsize{eff}}}$, if and
only if $(k_{1}^{*},k_{2}^{*})$ is a solution of the Bethe equations
associated to the diagonalization of $H_{\mbox{\scriptsize{eff}}}^{\dagger}$.
In the continuum case, for instance, starting from a solution of the Bethe
equation

$e^{ik_{1}\ell}=-\frac{c-i\left(k_{1}-k_{2}\right)}{c+i\left(k_{1}-k_{2}\right)},\qquad
k_{1}+k_{2}=\frac{2\pi}{\ell}\mathbb{Z},$ (63)

we obtain a solution for the Bethe equation for
$H_{\mbox{\scriptsize{eff}}}^{\dagger}$

$\begin{split}e^{ik_{1}^{*}\ell}&=-\frac{c^{*}-i\left(k_{1}^{*}-k_{2}^{*}\right)}{c^{*}+i\left(k_{1}^{*}-k_{2}^{*}\right)},\qquad
k_{1}^{*}+k_{2}^{*}=\frac{2\pi}{\ell}\mathbb{Z}.\end{split}$

The solution $(k_{1}^{*},k_{2}^{*})$ corresponds to an energy
$E_{k_{1}k_{2}}^{*}$ and a total momentum $k_{1}^{*}+k_{2}^{*}=k_{1}+k_{2}$.
Therefore, we have the proportionality relation
$\ket{\Psi_{k_{1}k_{2}}^{L}}=A\ket{\Psi_{k_{1}^{*}k_{2}^{*}}^{R}}$. The
normalization constant $A$ is set by the constraint
$\bra{\Psi_{k_{1}k_{2}}^{R}}\ket{\Psi_{k_{1}k_{2}}^{L}}=1$.

##  Appendix E Linearized momenta in the continuum

In this Appendix, we first derive Eq. (27). If $|\delta/c|\ll 1$ and $p$ is
even, Eq. (25) becomes
$e^{\frac{i}{2}\delta{\ell}}e^{-i\pi}=\frac{1-i\frac{\delta}{c}}{1+i\frac{\delta}{c}}\approx
1-2i\frac{\delta}{c}$. Thus, we have
$i\frac{\delta{\ell}}{2}-i\pi\approx\ln(1-2i\frac{\delta}{c})+2\pi
iq\approx-2i\frac{\delta}{c}+2\pi iq$ with $q\in\mathbb{Z}$. Finally, we
obtain $\delta\approx\frac{4\pi q+2\pi}{{\ell}+\frac{4}{c}}$. In Eq. (27), we
only consider $q\geq 0$ to avoid double counting of the eigenstates and
eigenenergies. The case $p$ odd is treated in a similar way.

Let us now discuss Eq. (28). If we consider $p=0$ as an example, we have
$\frac{\imaginary[\mathcal{E}_{p=0,q}]}{\real[\mathcal{E}_{p=0,q}]}=\frac{\imaginary[\delta^{2}]}{\real[\delta^{2}]}=\frac{\imaginary[\tilde{\delta}^{2}]}{\real[\tilde{\delta}^{2}]}$
where $\tilde{\delta}=\frac{\delta}{4\pi q+2\pi}=({\ell}+4/c)^{-1}$. Thus,
$\frac{\imaginary[\mathcal{E}_{p=0,q}]}{\real[\mathcal{E}_{p=0,q}]}$ is
independent of $q$ and we obtain Eq. (28).

##  Appendix F Asymptotic value of $|\imaginary[E]|$ for large $\delta$ in the
continuum

In this Appendix, we derive the expression of the asymptotic value of
$|\imaginary[E]|$ in the limit $|c/\delta|\ll 1$ in the continuum case. If
$|c/\delta|\ll 1$, Eq. (25) becomes
$e^{\frac{i}{2}\delta\ell}=(-1)^{p}\frac{1+i\frac{c}{\delta}}{1-i\frac{c}{\delta}}\approx(-1)^{p}\left(1+2i\frac{c}{\delta}\right)$.
Thus, we have $e^{i\delta\ell}\approx\left(1+4i\frac{c}{\delta}\right)$, and
$i\delta\ell\approx\ln(1+4i\frac{c}{\delta})+2\pi iq\approx
4i\frac{c}{\delta}+2\pi iq$ where $q\in\mathbb{Z}$. We obtain
$\delta^{2}\ell-2\pi q\delta-4c=0$ which is solved by $\delta_{\pm}=\frac{\pi
q}{\ell}\left(1\pm(1+\frac{2c\ell}{(\pi q)^{2}})\right)$ in the limit of large
integers $q\gg|c|\ell$; in order to be consistent with the numerical solution
we should select the solution $\delta_{+}$. Finally, we find
$|\imaginary[E]|\approx|\real[\delta_{+}]\imaginary[\delta_{+}]|\approx\frac{4|\imaginary[c]|}{\ell}$.

##  Appendix G Exact expressions of $\tilde{f}(p,\delta)$

In this Appendix, we give the exact expressions of $\tilde{f}(p,\delta)$ for
the four initial states studied in the main text. In the continuum, assuming
${\ell}\gg 1$,

  1. (i)

if the initial wave-function is
$\Psi_{0}(x_{1},x_{2})=\frac{1}{\sqrt{\pi}}\left(e^{-\left(x_{1}-\frac{{\ell}}{2}\right)^{2}}e^{-\left(x_{2}-\frac{{\ell}}{2}\right)^{2}}\right)$,
then

$\tilde{f}(p,\delta)=\frac{64\pi(c-c^{*})e^{-\frac{1}{8}(2(\frac{2\pi
p}{{\ell}})^{2}+{\delta}^{2}+{\delta^{*}}^{2})}{|\delta|^{2}}|{\delta}+c\operatorname{erfi}(\frac{{\delta}}{2\sqrt{2}})|^{2}}{{\ell}|4c+{\ell}\left(c^{2}+{\delta}^{2}\right)|^{2}\left({\delta}^{2}-{\delta^{*}}^{2}\right)};$
(64)

  2. (ii)

if the initial wave-function is
$\Psi_{1}(x_{1},x_{2})=\frac{1}{\sqrt{2\pi}}\left(e^{-(x_{1}-\frac{{\ell}}{4})^{2}}e^{-(x_{2}-3\frac{{\ell}}{4})^{2}}+e^{-(x_{2}-\frac{{\ell}}{4})^{2}}e^{-(x_{1}-3\frac{{\ell}}{4})^{2}}\right)$
, then

$\tilde{f}(p,\delta)=\frac{32\pi(c-c^{*})e^{-\frac{1}{8}(2(\frac{2\pi
p}{{\ell}})^{2}+\delta^{2}+{\delta^{*}}^{2})}|\delta|^{2}}{{\ell}|\delta^{2}+c^{2}|\left(\delta^{2}-{\delta^{*}}^{2}\right)|4c+{\ell}(c^{2}+\delta^{2})|^{2}}\left|(-1)^{p}(c+i\delta)\left(-i\delta+c\erf{\frac{{\ell}-i\delta}{2\sqrt{2}}}\right)+(c-i\delta)\left(i\delta+c\erf{\frac{{\ell}+i\delta}{2\sqrt{2}}}\right)\right|^{2},$
(65)

where by definition, $\erf(z):=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt$
and $\operatorname{erfi}(z):=-i\erf(iz)$.

On the lattice,

  1. (i)

if the initial state is
$\ket{\Psi_{2}}=b_{x_{0}}^{\dagger}b_{x_{0}}^{\dagger}\ket{v}/\sqrt{2}$, then

$\tilde{f}(p,\delta)=\frac{8|\sin^{4}(\frac{\delta}{2})|\cos^{3}(\frac{p\pi}{L})\left(u-u^{*}\right)}{L|Lu^{2}+2u\cos(\frac{p\pi}{L})\cos(\frac{\delta}{2})+L\cos^{2}(\frac{p\pi}{L})\sin^{2}(\frac{\delta}{2})|^{2}\left(\cos(\frac{\delta^{*}}{2})-\cos(\frac{\delta}{2})\right)};$
(66)

  2. (ii)

if the initial state is
$\ket{\Psi_{3}}=b_{x_{0}}^{\dagger}b_{x_{0}+(L-1)/2}^{\dagger}\ket{v}$, then

$\tilde{f}(p,\delta)=\frac{8\cos(\frac{p\pi}{L})\left(u-u^{*}\right)|\sin(\frac{\delta}{2})|^{2}|u^{2}+\cos^{2}(\frac{p\pi}{L})\sin^{2}(\frac{\delta}{2})|\left(\cos(\frac{{\delta^{*}}-\delta}{4})+(-1)^{p}\cos(\frac{{\delta^{*}}+\delta}{4})\right)}{L|Lu^{2}+2u\cos(\frac{p\pi}{L})\cos(\frac{\delta}{2})+L\cos^{2}(\frac{p\pi}{L})\sin^{2}(\frac{\delta}{2})|^{2}\left(\cos(\frac{\delta^{*}}{2})-\cos(\frac{\delta}{2})\right)},$
(67)

where by definition, $\cos(z):=\frac{e^{iz}+e^{-iz}}{2}$ and
$\sin(z):=\frac{e^{iz}-e^{-iz}}{2i}$ for $z\in\mathbb{C}$.

##  Appendix H Derivation of Eq. (39)

In this Appendix, we derive Eq. (39). Assuming ${\ell}|c|\gg 1$ and ${\ell}\gg
1$, we find

$g_{0}({\ell},c)=\frac{32(c-c^{*})(2\pi+|c|^{2}+\sqrt{2\pi}(c+c^{*}))}{|c|^{4}{\ell}^{3}}.$
(68)

We first rewrite the sum over $p$ in Eq. (38) as a sum over even $p$ and a sum
over odd $p$ because the expression for the linearized momenta $\delta$ in Eq.
(27) depends only on the parity of $p$ and never on its actual value. This
means that instead of defining $\mathcal{S}_{p}^{(0)}$, the set of allowed
values of $\delta$ for a given $p$, as we did in the main text, we can simply
define $\mathcal{S}_{e}^{(0)}$ and $\mathcal{S}_{o}^{(0)}$, where $e$ stands
for “even” and $o$ for “odd”. The double sum now takes a factorized
expression:

$\displaystyle\mathcal{N}(t)$ $\displaystyle\approx
g_{0}({\ell},c)\left(\sum_{p\text{ even}}e^{-\frac{1}{4}\left(\frac{2\pi
p}{{\ell}}\right)^{2}}\right)\left(\sum_{\delta\in\mathcal{S}^{(0)}_{e}}\frac{|\delta|^{4}e^{-\frac{i}{2}(\delta^{2}-\delta^{*2})t}}{\delta^{2}-\delta^{*2}}\right)$
$\displaystyle+g_{0}({\ell},c)\left(\sum_{p\text{
odd}}e^{-\frac{1}{4}\left(\frac{2\pi
p}{{\ell}}\right)^{2}}\right)\left(\sum_{\delta\in\mathcal{S}^{(0)}_{o}}\frac{|\delta|^{4}e^{-\frac{i}{2}(\delta^{2}-\delta^{*2})t}}{\delta^{2}-\delta^{*2}}\right).$
(69)

In the large $\ell$ limit, we transform the sums over $p$ into integrals

$\begin{split}\sum_{p\mbox{ \scriptsize{even}}}e^{-\frac{1}{4}(\frac{2\pi
p}{{\ell}})^{2}}&\approx\sum_{p\mbox{
\scriptsize{odd}}}e^{-\frac{1}{4}(\frac{2\pi p}{{\ell}})^{2}}\\\
&\approx\int_{-\infty}^{\infty}e^{-\frac{1}{4}(\frac{4\pi
p}{{\ell}})^{2}}\,dp=\frac{{\ell}}{2\sqrt{\pi}}.\end{split}$ (70)

Thus, we have

$\mathcal{N}(t)\approx
g_{0}({\ell},c)\frac{{\ell}}{2\sqrt{\pi}}\sum_{a=e,o}\underbrace{\sum_{\delta\in\mathcal{S}^{(0)}_{a}}\frac{|\delta|^{4}}{{\delta}^{2}-{\delta^{*}}^{2}}e^{-\frac{i}{2}({\delta}^{2}-{\delta^{*}}^{2})t}}_{\mathcal{J}_{a}}.$
(71)

We use the definitions of $\mathcal{S}^{(0)}_{e}$ and $\mathcal{S}^{(0)}_{o}$
and the fact that ${\ell}|c|\gg 1$ to find the following expression, valid for
$a=e,o$:

$\mathcal{J}_{a}=\frac{i128\pi^{4}}{\tilde{\Gamma}_{c}{\ell}}\sum_{q=0}^{\infty}(q+D_{a})^{2}e^{-(q+D_{a})^{2}\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t},$
(72)

with $D_{e}=1/2$ and $D_{o}=1$; we denote
$\tilde{\Gamma}_{c}=128\pi^{2}\frac{|\imaginary[c]|}{|c|^{2}}=128\pi^{2}\Gamma_{c}$.

When $|\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t|\ll 1$, we can approximate the
sum by an integral, because
$(q+D_{a})^{2}e^{-(q+D_{a})^{2}\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t}$ is
slowly varying when $q$ varies by $1$:

$\mathcal{J}_{a}=\frac{i128\pi^{4}}{\tilde{\Gamma}_{c}{\ell}}\int_{0}^{\infty}dq(q+D_{a})^{2}e^{-(q+D_{a})^{2}\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t}.$
(73)

Changing variable and using
$y=(q+D_{a})^{2}\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t$, we obtain

$\mathcal{J}_{a}=\frac{i64\pi^{4}}{\tilde{\Gamma}_{c}{\ell}}\left(\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t\right)^{-3/2}\int_{D_{a}^{2}\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t}^{\infty}\,dy\sqrt{y}e^{-y}$
(74)

The lower bound of the integral can be approximated by zero when
$|\frac{\tilde{\Gamma}_{c}}{{\ell}^{3}}t|\ll 1$; in this limit the integral
evaluates to $\sqrt{\pi}/2$. This leads us to Eq. (39) with

$h_{0}({\ell},c)=\frac{{\ell}^{3/2}\left(2\pi+|c|^{2}+2\sqrt{2\pi}\real[c]\right)}{64\sqrt{2}\pi|c|^{2}}.$
(75)

Eq. (41) can be derived along similar lines; here, we just report the explicit
expressions of $g_{1}({\ell},c)$ and $h_{1}({\ell})$ that we obtained assuming
${\ell}|c|\gg 1$ and ${\ell}\gg 1$:

$\displaystyle g_{1}({\ell},c)=$
$\displaystyle\begin{cases}\dfrac{128\pi(c-c^{*})}{{\ell}^{3}|c|^{2}},\quad&\mbox{for
}p\mbox{ is even};\\\ \\\ 0,\quad&\mbox{for }p\mbox{ is odd};\end{cases}$ (76)
$\displaystyle h_{1}({\ell})=$ $\displaystyle\frac{\sqrt{{\ell}}}{2\sqrt{2}}.$
(77)

##  Appendix I Derivation of Eq. (45)

In this Appendix, we derive Eq. (45). By using the total and relative momenta,
Eq. (16) can be rewritten as

$\mathcal{N}(t)=\sum_{p=0}^{L-1}\sum_{\delta\in\mathcal{S}_{p}^{\prime}}\tilde{f}\left(p,\delta\right)e^{-\frac{i}{2}({\delta}^{2}-{\delta^{*}}^{2})\cos(\frac{p\pi}{L})t},$
(78)

where $\mathcal{S}^{\prime}_{p}$ is the set of solutions of the equations
(30); in order to avoid double counting of solutions, we always take
$0<\real[\delta]<2\pi$. The explicit expression for $\tilde{f}(p,\delta)$,
obtained using the definitions given in Sec. IV and an initial state of the
form

$\ket{\Psi_{3}}=b_{x_{0}}^{\dagger}b_{x_{0}+(L-1)/2}^{\dagger}\ket{v},$ (79)

is reported in Eq. (67), see Appendix G. Similarly to the continuum case, if
we consider a sufficiently large time, then the dynamics is well captured by
the linearized spectrum, and we have to sum over the $\delta$ in the set
$\mathcal{S}^{\prime(0)}_{p}$ and $\mathcal{S}^{\prime(2\pi)}_{p}$,
corresponding to the values given in Eq. (32a) and Eq. (32b) respectively. We
assume the following symmetry between the sums over $\delta$ close to $0$ and
over those close to $2\pi$:

$\begin{split}\sum_{\delta\in\mathcal{S}_{p}^{\prime(0)}}&\tilde{f}\left(p,\delta\right)e^{-\frac{i}{2}({\delta}^{2}-{\delta^{*}}^{2})\cos(\frac{p\pi}{L})t}\\\
&\approx\sum_{\delta\in\mathcal{S}_{p}^{\prime(2\pi)}}\tilde{f}\left(p,\delta\right)e^{-\frac{i}{2}({\delta}^{2}-{\delta^{*}}^{2})\cos(\frac{p\pi}{L})t},\end{split}$
(80)

which we verified numerically. This observation significantly simplifies the
calculation, because we can sum only over $\mathcal{S}_{p}^{\prime(0)}$ and
approximate $\tilde{f}\left(p,\delta\right)$ by its expansion for $\delta\sim
0$, that we call $\tilde{f}_{3}\left(p,\delta\right)$. Assuming $L|u|\gg 1$,
we have

$\tilde{f}_{3}\left(p,\delta\right)=\frac{16\left(1+(-1)^{p}\right)\left(u-u^{*}\right)\cos(\frac{p\pi}{L})}{|u|^{2}L^{3}}\frac{|\delta|^{2}}{{\delta}^{2}-{\delta^{*}}^{2}}.$
(81)

Note that $\tilde{f}_{3}\left(p,\delta\right)=0$ if $p$ is odd. This technical
simplification does not appear for the initial state consisting of two bosons
at the same site and it makes the derivation of Eq. (44) more lengthy. By
using the analytical expression of the linearized momenta Eq. (32a), we find

$\mathcal{N}(t)\approx\frac{16}{L^{2}}S$ (82)

with

$S=\sum_{p=0}^{(L-1)/2}\sum_{q=0}^{(L-1)/4-1}e^{-\frac{1}{2}\left(q+\frac{1}{2}\right)^{2}\cos^{2}(\frac{2\pi
p}{L})\xi},$ (83)

and with $\xi=\frac{\Gamma t}{L^{3}}$. Without loss of generality (provided
that $L$ is odd) we consider $L\equiv 1\pmod{4}$. Since we are interested in
the intermediate time regime that is located after the initial transient
dynamics but before the exponential asymptotic decay, we will consider
$1\ll\xi\ll L^{2}$ to further simplify the expression for $S$. Under this
assumption, we can perform a series expansion of the cosine at the exponent
for $p$ close to $L/4$ and extend the ranges of the two sums towards infinity
to arrive at

$S\approx\sum_{p=-\infty}^{\infty}\sum_{q=0}^{\infty}e^{-\frac{2\pi^{2}\xi}{L^{2}}\left(q+\frac{1}{2}\right)^{2}\left(p-\frac{L}{4}\right)^{2}}.$
(84)

We now change variable $p\rightarrow-p+(L-1)/4$, so that $p-L/4$ becomes
$p+1/4$:

$S\approx\sum_{p=-\infty}^{\infty}\sum_{q=0}^{\infty}e^{-\frac{2\pi^{2}\xi}{L^{2}}\left(q+\frac{1}{2}\right)^{2}\left(p+\frac{1}{4}\right)^{2}}.$
(85)

Using the Fourier representation of the Gaussian
$e^{-\sigma_{q}(p+\frac{1}{4})^{2}}=\sqrt{\frac{\pi}{\sigma_{q}}}\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{4\sigma_{q}}+i\frac{x}{4}}e^{ipx}\frac{dx}{2\pi}$
with $\sigma_{q}=\frac{2\pi^{2}\xi}{L^{2}}\left(q+\frac{1}{2}\right)^{2}$, we
introduce a Dirichlet kernel
$\sum_{p=-\infty}^{\infty}e^{ipx}=2\pi\sum_{n=-\infty}^{\infty}\delta(x-2\pi
n)$:

$\begin{split}S\approx&\sum_{q=0}^{\infty}\frac{1}{2\sqrt{\sigma_{q}\pi}}\int_{-\infty}^{\infty}e^{-\frac{x^{2}}{4\sigma_{q}}+i\frac{x}{4}}\sum_{p=-\infty}^{\infty}e^{ipx}dx\\\
=&\sum_{q=0}^{\infty}\sqrt{\frac{\pi}{\sigma_{q}}}\sum_{n=-\infty}^{\infty}e^{-\frac{\pi^{2}n^{2}}{\sigma_{q}}}i^{n}.\end{split}$
(86)

The terms corresponding to $n$ odd are canceled out by those corresponding to
$-n$. We obtain

$\begin{split}S\approx\sum_{q=0}^{\infty}\underbrace{\sqrt{\frac{1}{2\pi\xi(\frac{q}{L}+\frac{1}{2L})^{2}}}\theta_{4}\left(e^{-\frac{1}{\frac{\xi}{2}\left(\frac{q}{L}+\frac{1}{2L}\right)^{2}}}\right)}_{F(q)}\end{split}$
(87)

where $\theta_{4}(q)=\sum_{n=-\infty}^{\infty}(-1)^{n}q^{n^{2}}$ is one of the
Jacobi theta functions [76]. In the large $L$ limit, we substitute the
summation with an integral and include the first and second Euler-Maclaurin
correction [76]:

$S\approx\int_{0}^{\infty}F(q)dq+\frac{F(\infty)+F(0)}{2}+\frac{F^{\prime}(\infty)-F^{\prime}(0)}{12}$
(88)

In the $\xi\ll L^{2}$ limit, the correction
$\frac{F(\infty)+F(0)}{2}+\frac{F^{\prime}(\infty)-F^{\prime}(0)}{12}$ reads
$\frac{4}{3}\sqrt{\frac{L^{2}}{2\pi\xi}}$. Changing the variable and using
$x=\frac{L}{\sqrt{\xi}\left(q+1/2\right)}$, we obtain

$S\approx\sqrt{\frac{L^{2}}{2\pi\xi}}\left(\frac{4}{3}+\underbrace{\int_{0}^{2L/\sqrt{\xi}}\frac{dx}{x}\theta_{4}\left(e^{-2x^{2}}\right)}_{I}\right).$
(89)

Notice that $2L/\sqrt{\xi}\gg 1$; at large $x$,
$\theta_{4}\left(e^{-2x^{2}}\right)\approx 1$. The integral is therefore
logarithmically divergent. We can rewrite it by introducing a cutoff scale
$1\ll\Lambda\ll 2L/\sqrt{\xi}$ as

$\begin{split}I\approx\int_{0}^{\Lambda}\frac{dx}{x}\theta_{4}\left(e^{-2x^{2}}\right)+\int_{\Lambda}^{2L/\sqrt{\xi}}\frac{dx}{x}.\end{split}$
(90)

In the second integral we approximated
$\theta_{4}\left(e^{-2x^{2}}\right)\approx 1$, since $\Lambda\gg 1$. We thus
get

$\begin{split}\int_{0}^{2L/\sqrt{\xi}}\frac{dx}{x}\theta_{4}\left(e^{-2x^{2}}\right)\approx\zeta+\log\left(\frac{2L}{\sqrt{\xi}}\right)\end{split}$
(91)

where

$\zeta=\lim_{\Lambda\to\infty}\int_{0}^{\Lambda}\frac{dx}{x}\theta_{4}\left(e^{-2x^{2}}\right)-\log(\Lambda)$
(92)

which is just a real constant. We evaluated it numerically and obtained
$\zeta\approx 0.183599$. We finally arrive at

$\begin{split}S=\sqrt{\frac{L^{2}}{2\pi\xi}}\left(\underbrace{\frac{4}{3}+\zeta}_{C^{\prime}}+\log\left(\frac{2L}{\sqrt{\xi}}\right)\right)\end{split}$
(93)

which, together with Eq. (82), gives Eq. (45) we reported in the main text.

Equation (44) can be derived along similar lines; we here just report the
constant $C$ we obtained

$:\begin{split}C=&\lim_{\Lambda\to\infty}\left[\int_{0}^{\Lambda}\frac{dx}{x}\theta_{4}(0,e^{-2x^{2}})-\log\Lambda\right]-\frac{1}{2}\log(2)+\frac{23}{12}\\\
\approx&1.75\end{split}$ (94)

## References

  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49, 435 (1977).
  * Goldenfeld [1992] N. Goldenfeld,  _Lectures On Phase Transitions And The Renormalization Group_ (CRC Press, 1992).
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum Field Theory and Critical Phenomena_ (Oxford University Press, 2002).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Reports on Progress in Physics 79, 096001 (2016).
  * Sieberer  _et al._ [2023] L. M. Sieberer, M. Buchhold, J. Marino, and S. Diehl, Universality in driven open quantum matter (2023), arXiv:2312.03073 [cond-mat.stat-mech] .
  * Altman  _et al._ [2015] E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J. Toner, Two-dimensional superfluidity of exciton polaritons requires strong anisotropy, Phys. Rev. X 5, 011017 (2015).
  * Fontaine  _et al._ [2022] Q. Fontaine, D. Squizzato, F. Baboux, I. Amelio, A. Lemaître, M. Morassi, I. Sagnes, L. Le Gratiet, A. Harouri, M. Wouters, I. Carusotto, A. Amo, M. Richard, A. Minguzzi, L. Canet, S. Ravets, and J. Bloch, Kardar–Parisi–Zhang universality in a one-dimensional polariton condensate, Nature 608, 687–691 (2022).
  * Barontini  _et al._ [2013] G. Barontini, R. Labouvie, F. Stubenrauch, A. Vogler, V. Guarrera, and H. Ott, Controlling the dynamics of an open many-body quantum system with localized dissipation, Phys. Rev. Lett. 110, 035302 (2013).
  * Perfetto  _et al._ [2023a] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, Reaction-limited quantum reaction-diffusion dynamics, Phys. Rev. Lett. 130, 210402 (2023a).
  * Perfetto  _et al._ [2023b] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, Quantum reaction-limited reaction-diffusion dynamics of annihilation processes, Phys. Rev. E 108, 064104 (2023b).
  * Rowlands  _et al._ [2024] S. Rowlands, I. Lesanovsky, and G. Perfetto, Quantum reaction-limited reaction–diffusion dynamics of noninteracting bose gases, New J. Phys. 26, 043010 (2024).
  * Buchhold  _et al._ [2017] M. Buchhold, B. Everest, M. Marcuzzi, I. Lesanovsky, and S. Diehl, Nonequilibrium effective field theory for absorbing state phase transitions in driven open quantum spin systems, Phys. Rev. B 95, 014308 (2017).
  * Wintermantel  _et al._ [2021] T. M. Wintermantel, M. Buchhold, S. Shevate, M. Morgado, Y. Wang, G. Lochead, S. Diehl, and S. Whitlock, Epidemic growth and Griffiths effects on an emergent network of excited atoms, Nature Communications 12, 103 (2021).
  * Cai and Barthel [2013] Z. Cai and T. Barthel, Algebraic versus exponential decoherence in dissipative many-particle systems, Phys. Rev. Lett. 111, 150403 (2013).
  * Begg and Hanai [2024] S. E. Begg and R. Hanai, Quantum criticality in open quantum spin chains with nonreciprocity, Phys. Rev. Lett. 132, 120401 (2024).
  * Dürr  _et al._ [2009] S. Dürr, J. J. García-Ripoll, N. Syassen, D. M. Bauer, M. Lettner, J. I. Cirac, and G. Rempe, Lieb-Liniger model of a dissipation-induced Tonks-Girardeau gas, Phys. Rev. A 79, 023614 (2009).
  * García-Ripoll  _et al._ [2009] 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 J. Phys. 11, 013053 (2009).
  * Baur and Mueller [2010] S. K. Baur and E. J. Mueller, Two-body recombination in a quantum-mechanical lattice gas: Entropy generation and probing of short-range magnetic correlations, Phys. Rev. A 82, 023626 (2010).
  * Foss-Feig  _et al._ [2012] M. Foss-Feig, A. J. Daley, J. K. Thompson, and A. M. Rey, Steady-state many-body entanglement of hot reactive fermions, Phys. Rev. Lett. 109, 230501 (2012).
  * Grišins  _et al._ [2016] P. Grišins, B. Rauer, T. Langen, J. Schmiedmayer, and I. E. Mazets, Degenerate bose gases with uniform loss, Phys. Rev. A 93, 033634 (2016).
  * Johnson  _et al._ [2017] A. Johnson, S. S. Szigeti, M. Schemmer, and I. Bouchoule, Long-lived nonthermal states realized by atom losses in one-dimensional quasicondensates, Phys. Rev. A 96, 013623 (2017).
  * Yamamoto  _et al._ [2019] K. Yamamoto, M. Nakagawa, K. Adachi, K. Takasan, M. Ueda, and N. Kawakami, Theory of Non-Hermitian Fermionic Superfluidity with a Complex-Valued Interaction, Phys. Rev. Lett. 123, 123601 (2019).
  * Goto and Danshita [2020] S. Goto and I. Danshita, Measurement-induced transitions of the entanglement scaling law in ultracold gases with controllable dissipation, Phys. Rev. A 102, 033316 (2020).
  * Bouchoule  _et al._ [2020] 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, 44 (2020).
  * Rossini  _et al._ [2021] 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, Phys. Rev. A 103, L060201 (2021).
  * Nakagawa  _et al._ [2020] M. Nakagawa, N. Tsuji, N. Kawakami, and M. Ueda, Dynamical Sign Reversal of Magnetic Correlations in Dissipative Hubbard Models, Phys. Rev. Lett. 124, 147203 (2020).
  * Booker  _et al._ [2020] C. Booker, B. Buča, and D. Jaksch, Non-stationarity and dissipative time crystals: spectral properties and finite-size effects, New J. Phys. 22, 085007 (2020).
  * Bouchoule and Dubail [2021] I. Bouchoule and J. Dubail, Breakdown of Tan’s Relation in Lossy One-Dimensional Bose Gases, Phys. Rev. Lett. 126, 160603 (2021).
  * Bouchoule  _et al._ [2021] I. Bouchoule, L. Dubois, and L.-P. Barbier, Losses in interacting quantum gases: Ultraviolet divergence and its regularization, Phys. Rev. A 104, L031304 (2021).
  * Rosso  _et al._ [2021] 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 104, 053305 (2021).
  * Nakagawa  _et al._ [2021] M. Nakagawa, N. Kawakami, and M. Ueda, Exact Liouvillian Spectrum of a One-Dimensional Dissipative Hubbard Model, Phys. Rev. Lett. 126, 110404 (2021).
  * Yamamoto  _et al._ [2021] K. Yamamoto, M. Nakagawa, N. Tsuji, M. Ueda, and N. Kawakami, Collective excitations and nonequilibrium phase transition in dissipative fermionic superfluids, Phys. Rev. Lett. 127, 055301 (2021).
  * Rosso  _et al._ [2022a] L. Rosso, L. Mazza, and A. Biella, Eightfold way to dark states in SU(3) cold gases with two-body losses, Phys. Rev. A 105, L051302 (2022a).
  * Yamamoto  _et al._ [2022] K. Yamamoto, M. Nakagawa, M. Tezuka, M. Ueda, and N. Kawakami, Universal properties of dissipative Tomonaga-Luttinger liquids: Case study of a non-Hermitian XXZ spin chain, Phys. Rev. B 105, 205125 (2022).
  * Yoshida and Katsura [2023] H. Yoshida and H. Katsura, Liouvillian gap and single spin-flip dynamics in the dissipative Fermi-Hubbard model, Phys. Rev. A 107, 033332 (2023).
  * Huang  _et al._ [2023] C.-H. Huang, T. Giamarchi, and M. A. Cazalilla, Modeling particle loss in open systems using Keldysh path integral and second order cumulant expansion, Phys. Rev. Res. 5, 043192 (2023).
  * Mazza and Schirò [2023] G. Mazza and M. Schirò, Dissipative dynamics of a fermionic superfluid with two-body losses, Phys. Rev. A 107, L051301 (2023).
  * Hamanaka  _et al._ [2023] S. Hamanaka, K. Yamamoto, and T. Yoshida, Interaction-induced Liouvillian skin effect in a fermionic chain with a two-body loss, Phys. Rev. B 108, 155114 (2023).
  * Nagao  _et al._ [2023] K. Nagao, I. Danshita, and S. Yunoki, Semiclassical descriptions of dissipative dynamics of strongly interacting Bose gases in optical lattices (2023), arXiv:2307.16170 [cond-mat.quant-gas] .
  * Xiao and Watanabe [2023] T. Xiao and G. Watanabe, Local non-hermitian hamiltonian formalism for dissipative fermionic systems and loss-induced population increase in fermi superfluids (2023), arXiv:2306.16235 [cond-mat.quant-gas] .
  * Maki  _et al._ [2024] J. Maki, L. Rosso, L. Mazza, and A. Biella, Loss induced collective mode in one-dimensional Bose gases (2024), arXiv:2402.05824 [cond-mat.quant-gas] .
  * Syassen  _et al._ [2008] 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 320, 1329 (2008).
  * Yan  _et al._ [2013] 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 501, 521–525 (2013).
  * Zhu  _et al._ [2014] 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. 112, 070404 (2014).
  * Rauer  _et al._ [2016] 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. 116, 030402 (2016).
  * Tomita  _et al._ [2017] 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, Sci. Adv. 3, e1701513 (2017).
  * Schemmer and Bouchoule [2018] M. Schemmer and I. Bouchoule, Cooling a bose gas by three-body losses, Phys. Rev. Lett. 121, 200401 (2018).
  * Sponselee  _et al._ [2019] 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 Sci. Technol. 4, 014002 (2019).
  * Tomita  _et al._ [2019] T. Tomita, S. Nakajima, Y. Takasu, and Y. Takahashi, Dissipative Bose-Hubbard system with intrinsic two-body loss, Phys. Rev. A 99, 031601 (2019).
  * Bouchoule and Schemmer [2020] I. Bouchoule and M. Schemmer, Asymptotic temperature of a lossy condensate, SciPost Phys. 8, 060 (2020).
  * Honda  _et al._ [2023] 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. 130, 063001 (2023).
  * Rosso  _et al._ [2023] 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 107, 013303 (2023).
  * Gerbino  _et al._ [2023] F. Gerbino, I. Lesanovsky, and G. Perfetto, Large-scale universality in Quantum Reaction-Diffusion from Keldysh field theory (2023), arXiv:2307.14945 [cond-mat.stat-mech] .
  * Rosso  _et al._ [2022b] L. Rosso, A. Biella, and L. Mazza, The one-dimensional Bose gas with strong two-body losses: the effect of the harmonic confinement, SciPost Phys. 12, 44 (2022b).
  * Riggio  _et al._ [2024] F. Riggio, L. Rosso, D. Karevski, and J. Dubail, Effects of atom losses on a one-dimensional lattice gas of hard-core bosons, Phys. Rev. A 109, 023311 (2024).
  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Non-hermitian physics, Advances in Physics 69, 249–435 (2020).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Lieb and Liniger [1963] E. H. Lieb and W. Liniger, Exact Analysis of an Interacting Bose Gas. I. The General Solution and the Ground State, Phys. Rev. 130, 1605 (1963).
  * Franchini [2017] F. Franchini,  _An Introduction to Integrable Techniques for One-Dimensional Quantum Systems_ (Springer International Publishing, 2017).
  * Gaudin [1983] M. Gaudin,  _La fonction d’onde de Bethe_ (Masson, Paris, 1983) The Bethe Wavefunction (translation by J.-S. Caux), Cambridge University Press, 2014.
  * Korepin  _et al._ [1997] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin,  _Quantum inverse scattering method and correlation functions_ , Vol. 3 (Cambridge University Press, 1997).
  * Zvonarev [2010] M. Zvonarev,  _Notes on Bethe Ansatz_ (2010).
  * Šamaj and Bajnok [2013] L. Šamaj and Z. Bajnok,  _Introduction to the statistical physics of integrable many-body systems_ (Cambridge University Press, 2013).
  * Choy and Haldane [1982] T. Choy and F. Haldane, Failure of Bethe-Ansatz solutions of generalisations of the Hubbard chain to arbitrary permutation symmetry, Phys. Lett. A 90, 83 (1982).
  * Medvedyeva  _et al._ [2016] M. V. Medvedyeva, F. H. Essler, and T. Prosen, Exact Bethe ansatz spectrum of a tight-binding chain with dephasing noise, Phys. Rev. Lett. 117, 137202 (2016).
  * Alba [2023] V. Alba, Free fermions with dephasing and boundary driving: Bethe ansatz results (2023), arXiv:2309.12978 [cond-mat.stat-mech] .
  * Essler  _et al._ [2005] F. H. Essler, H. Frahm, F. Göhmann, A. Klümper, and V. E. Korepin,  _The one-dimensional Hubbard model_ (Cambridge University Press, 2005).
  * Zhi and Hong [1991] S. Zhi and S. Hong, Pseudospin Wave State in One-Dimensional Half-Filled Hubbard Model, Commun. Theor. Phys. 16, 485 (1991).
  * Zhang and Song [2021] X. Zhang and Z. Song, $\eta$-pairing ground states in the non-Hermitian Hubbard model, Phys. Rev. B 103, 235153 (2021).
  * Oelkers and Links [2007] N. Oelkers and J. Links, Ground-state properties of the attractive one-dimensional Bose-Hubbard model, Phys. Rev. B 75, 115119 (2007).
  * Li  _et al._ [2022] Y. Li, D. Schneble, and T.-C. Wei, Two-particle states in one-dimensional coupled Bose-Hubbard models, Phys. Rev. A 105, 053310 (2022).
  * Tanaka  _et al._ [2013] A. Tanaka, N. Yonezawa, and T. Cheon, Exotic quantum holonomy and non-Hermitian degeneracies in the two-body Lieb–Liniger model, J. Phys. A: Math. Theor. 46, 315302 (2013).
  * Yang [1989] C. N. Yang, $\eta$ pairing and off-diagonal long-range order in a Hubbard model, Phys. Rev. Lett. 63, 2144 (1989).
  * Giamarchi [2003] T. Giamarchi,  _Quantum Physics in One Dimension_ (Oxford University Press, 2003).
  * Ogata and Shiba [1990] M. Ogata and H. Shiba, Bethe-ansatz wave function, momentum distribution, and spin correlation in the one-dimensional strongly correlated Hubbard model, Phys. Rev. B 41, 2326 (1990).
  * Abramowitz [1974] M. Abramowitz,  _Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables,_ (Dover Publications, Inc., USA, 1974).