# Large-scale universality in Quantum Reaction-Diffusion from Keldysh field
theory

Federico Gerbino<EMAIL_ADDRESS>Institut für Theoretische
Physik, Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany
Physics Department, University of Trento, via Sommarive, 14 I-38123 Trento,
Italy  Igor Lesanovsky  Institut für Theoretische Physik, Universität
Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany  School of Physics
and Astronomy and Centre for the Mathematics and Theoretical Physics of
Quantum Non-Equilibrium Systems, The University of Nottingham, Nottingham, NG7
2RD, United Kingdom  Gabriele Perfetto  Institut für Theoretische Physik,
Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany

###### Abstract

We consider the quantum reaction-diffusion dynamics in $d$ spatial dimensions
of a Fermi gas subject to binary annihilation reactions $A+A\to\emptyset$.
These systems display collective nonequilibrium long-time behavior, which is
signalled by an algebraic decay of the particle density. Building on the
Keldysh formalism, we devise a field theoretical approach for the reaction-
limited regime, where annihilation reactions are scarce. By means of a
perturbative expansion of the dissipative interaction, we derive a description
in terms of a large-scale universal kinetic equation. Our approach shows how
the time-dependent generalized Gibbs ensemble assumption, which is often
employed for treating low-dimensional nonequilibrium systems, emerges from
systematic diagrammatics. It also allows to exactly compute — for arbitrary
spatial dimension — the decay exponent of the particle density. The latter is
based on the large-scale description of the quantum dynamics and it differs
from the mean-field prediction even in dimension larger than one. We moreover
consider spatially inhomogeneous setups involving an external potential. In
confined systems the density decay is accelerated towards the mean-field
algebraic behavior, while for deconfined scenarios the power-law decay is
replaced by a slower non-algebraic decay.

Introduction.— The identification of large-scale universal behavior in many-
body nonequilibrium systems is a demanding task. Reaction-diffusion (RD)
systems [1, 2, 3], where particles diffuse and react upon meeting, are ideal
systems for the investigation of dynamical universal behavior. For example,
for binary annihilation reactions $A+A\to\emptyset$, the late time decay of
the particle density takes a universal power-law form. In the “diffusion-
limited” regime [4, 5, 6, 7, 8, 9, 10, 11, 12], where diffusion is weak, the
origin of this dynamical behavior are spatial density fluctuations. Here,
mean-field approaches cannot be applied and field-theorical and
renormalization group analyses [13, 14, 15, 16, 17, 18, 19, 20, 21, 22]
correctly predict the observed power-law. The mean-field approximation is,
however, valid in more than one dimension and/or in the “reaction-limited”
regime of fast hopping mixing [2, 3, 6, 23].

Figure 1: Quantum RD dynamics via Keldysh field theory. (a) Comparison of
classical and quantum RD dynamics: classical incoherent diffusion (top-blue
solid lines) is replaced by quantum coherent motion (bottom-blue wiggly
lines), while in both cases annihilation, $A+A\to\emptyset$, is irreversible.
(b) Time evolution of the density matrix $\rho(t)$ along the closed-time
contour associated to the Lindblad map $\mathcal{L}$. Two time branches,
forward ($+$) and backward ($-$), are required. The particle number decreases
during time evolution, as can be seen by comparing the initial $n(x,t_{0})$
and final $n(x,t_{f})$ density profiles. (c) Annihilation interaction vertices
define the self-energy Keldysh matrix $\hat{\Sigma}$, which dresses the
Green’s functions $\hat{G}$ with respect to their bare values $\hat{G}_{0}$.
In the reaction-limited regime $\hbar n(x,t)\Gamma/D\ll 1$, the quasi-particle
dispertion relation $\epsilon_{k}$ is not modified, while quasi-particles
acquire a large-finite lifetime $\sim\Gamma^{-1}$, given by the energy width
$\epsilon$ of $G^{K}(\vec{x},t,\vec{k},\epsilon)$ (sketched in red). This
lifetime is given by tadpole Feynman diagrams (sketched in the light blue
inset). At leading order in the space-time variations, these diagrams
reproduce the TGGE predictions for $n(\vec{x},\vec{k},t)$ (spectral integral
in $\epsilon$ of $G^{K}(\vec{x},t,\vec{k},\epsilon)$ depicted in light blue).

In quantum many-body systems large-scale properties are even harder to uncover
than in the classical realm, already in the one-dimensional case, since they
entail the simulations of large sizes and long times [24, 25, 26, 27, 28, 29,
30]. In this regard, quantum RD systems have moved into the focus of
attention. They follow simple dynamical rules [31, 32, 33, 34, 35, 36, 37, 38,
39, 40], which connect to cold-atomic experiments [41, 42, 43, 44, 45, 46,
47], and they allow for novel forms of particle-density decay beyond mean-
field. This has been shown in Refs. [31, 32, 33, 36, 48, 49, 50, 51, 52, 53]
for one-dimensional systems in the reaction-limited regime. In this limit,
analytical predictions can be obtained under the assumption that the systems
relaxes to a time-dependent generalized Gibbs ensemble (TGGE) [54, 55, 56,
57]. The systematic classification of these universal properties, and the
study of their dependence on the spatial dimensions $d$, requires, however,
the development of a field-theory.

In this manuscript, we accomplish this by exploiting the Keldysh path integral
representation of the quantum master equation [58, 59, 60, 61, 62, 63, 64]. We
investigate as a paradigmatic example the Fermi gas in $d$ spatial dimensions
subject to binary annihilation reactions $A+A\to\emptyset$. From the Keldysh-
field theory, we perform a diagrammatic expansion of the dissipative
interaction vertices. At leading order in space-time derivatives, this
expansion acquires the universal form of a kinetic Boltzmann equation. In
$d=1$, this systematic derivation naturally yields the TGGE ansatz for the
reaction-limited regime. Our results therefore bridge the gap between the
much-employed TGGE relaxation assumption, diagrammatic expansions and large-
scale dynamics in dissipative quantum systems. Crucially, from the field-
theory description, we exactly compute the density decay exponent in arbitrary
dimensions, which is found to deviate from the mean-field prediction even in
$d>1$. This result is in contrast with the classical case and it is rooted
into the large-scale universal description of the underlying quantum dynamics.
We also consider the case of inhomogeneous systems where we study quenches of
a trapping potential confining the fermions. For a quench from a double to
single well potential, we find an acceleration of the particle decay, which
diverts the decay exponent towards the mean-field one. For a trap-release
quench, we, instead, find a qualitatively different scenario: the algebraic
decay first slows down on an intermediate time window, and then it gets
replaced at long times by a slower non-algebraic decay.

Quantum RD Keldysh action.— The dynamics of the considered Fermi gas in $d$
spatial dimensions is governed by the quantum master equation [65, 66, 67, 68]
with a Lindblad map $\mathcal{L}$

$\dot{\rho}(t)=\mathcal{L}[\rho(t)]=-\frac{i}{\hbar}[H,\rho(t)]+\mathcal{D}[\rho(t)].$
(1)

The Hamiltonian $H$ describes coherent free motion, see Fig. 1(a), which is
the quantum mechanical counterpart of classical diffusion. We further allow
for the presence of an external trapping potential $V(\vec{x})$:

$H=\int
d^{d}\vec{x}\,\psi^{\dagger}(\vec{x})[-D\nabla^{2}+V(\vec{x})]\psi(\vec{x}).$
(2)

Here $D=\hbar^{2}/(2m)$, while $\psi$, $\psi^{\dagger}$ are fermionic field
operators satisfying canonical anticommutation relations
$\\{\psi(\vec{x}),\psi^{\dagger}(\vec{x}^{\prime})\\}\,=\,\delta(\vec{x}-\vec{x}^{\prime})$.
The dissipator $\mathcal{D}[\rho]$ embodies irreversible reaction processes

$\displaystyle\mathcal{D}[\rho(t)]=\sum_{\alpha}\int d^{d}\vec{x}\,\Big{[}\,$
$\displaystyle L_{\alpha}(\vec{x})\rho(t)L^{\dagger}_{\alpha}(\vec{x})$
$\displaystyle-\frac{1}{2}\\{\,L^{\dagger}_{\alpha}(\vec{x})L_{\alpha}(\vec{x}),\rho(t)\\}\Big{]},$
(3)

with $\alpha=1,2\dots d$. We focus on binary annihilation reactions
$A+A\,\rightarrow\,\emptyset$, see Fig. 1(a), modelled by the jump operators

$L_{\alpha}(\vec{x})=\sqrt{\Gamma}\,\psi(\vec{x})\,\partial_{x_{\alpha}}\psi(\vec{x}).$
(4)

The constant $\Gamma$ (units: $\mbox{length}^{d+2}/\mbox{time}$) characterizes
the annihilation reactions. In the Supplemental Material [69], the jump
operator (4) is obtained by taking the continuum limit of nearest neighbours
annihilation. The latter is the natural annihilation decay to consider for
fermionic particles, where onsite reactions are forbidden.

Nonequilibrium universal behavior manifests in the power-law decay of the
density $n(x,t)=\braket{\psi^{\dagger}(\vec{x})\psi(\vec{x})}_{t}$ in time.
Here, we characterise this decay in the reaction-limited regime, where $\hbar
n(x,t)\Gamma/D\ll 1$ throughout the dynamics. In this limit, reactions are
weak and the density slowly changes in time. We further allow for weak space
inhomogeneities due the presence of the trapping potential $V(\vec{x})$, which
we assume to vary on macroscopic length scales. In this limit, the
quasiparticle dispersion relation
$\epsilon_{k}(\vec{x})=D\vec{k}^{2}+V(\vec{x})$ of $H$ in Eq. (2), with
$\vec{k}$ the momentum, is locally modified by the external potential. In this
regard, the quantum reaction-limited regime is particularly suitable for the
study of large space-time universality, i.e., taking place over long times and
large distances, in the density decay.

In order to do this, we exploit the Keldysh quantum field theory description
[60, 61, 62, 64, 63] of the open quantum RD dynamics (1)-(4) [69]. A general
feature of Keldysh field theory is the doubling of the $\psi,\,\bar{\psi}$
fields into four fields — $\psi_{+},\,\bar{\psi}_{+}$ and
$\psi_{-},\,\bar{\psi}_{-}$ — evolving along a “forward” ($+$) and a
“backward” ($-$) contour of the considered time interval, respectively, as
sketched in Fig. 1(b). “Plus” and “minus” fields are usually rewritten in
terms of Keldysh-rotated fields [59, 69, 70]
$\phi_{1},\,\bar{\phi}_{1},\,\phi_{2},\,\bar{\phi}_{2}$. The Keldysh partition
function
$Z(t)\\!\\!=\\!\\!\mathrm{tr}[\rho(t)]\\!\\!=\\!\\!\int\mathbf{D}[\phi_{1},\bar{\phi}_{1},\phi_{2},\bar{\phi}_{2}]\,\mathrm{exp}\\{iS[\phi_{1},\bar{\phi}_{1},\phi_{2},\bar{\phi}_{2}]\\}$
includes full information on the system’s microscopic dynamics. The action
$S=S_{\mathrm{diff}}+S_{\mathcal{D}}$ is composed of two sectors. A quadratic
sector $S_{\mathrm{diff}}$ describes coherent motion (2):

$S_{\mathrm{diff}}=\int
d^{d}x\,dt^{\prime}\,[\bar{\phi}_{1}(G_{0}^{R})^{-1}\phi_{1}+\bar{\phi}_{2}(G_{0}^{A})^{-1}\phi_{2}\,],$
(5)

with
$(G_{0}^{R/A})^{-1}\,=\,i\partial_{t^{\prime}}+(D\nabla^{2}-V(\vec{x}))/\hbar\pm
i\delta$ the inverse retarded/advanced bare propagator, respectively. The bare
Keldysh Green’s function $G_{0}^{K}\equiv-i\braket{\phi_{1}\bar{\phi}_{2}}$
associated to $S_{\mathrm{diff}}$ is a regularization factor. The Green’s
function $\braket{\phi_{2}\bar{\phi}_{1}}=0$ is generically zero as a
consequence of trace preservation. The retarded/advanced propagator
$G_{0}^{R/A}$ (with $V(\vec{x})=0$) in $S_{\mathrm{diff}}$ has a structure
similar to the classical RD quadratic counterpart [14, 15, 16, 17, 21, 20,
22], the difference between the two being in the factor $i$ in front of the
time derivative in (5). The second part $S_{\mathcal{D}}$ of the action
contains the interaction vertices of the theory due to annihilation reactions
(4):

$\displaystyle S_{\mathcal{D}}$ $\displaystyle=\frac{i\Gamma}{4}\int
d^{d}x\,dt^{\prime}\Big{[}2(\vec{\nabla}\bar{\phi}_{1}\bar{\phi}_{2}+\vec{\nabla}\bar{\phi}_{2}\bar{\phi}_{1})\cdot(\phi_{1}\vec{\nabla}\phi_{2}+$
$\displaystyle+\phi_{2}\vec{\nabla}\phi_{1})+(\vec{\nabla}\bar{\phi}_{1}\bar{\phi}_{2}+\vec{\nabla}\bar{\phi}_{2}\bar{\phi}_{1})\cdot(\phi_{1}\vec{\nabla}\phi_{1}+\phi_{2}\vec{\nabla}\phi_{2})-$
$\displaystyle-(\vec{\nabla}\bar{\phi}_{1}\bar{\phi}_{1}+\vec{\nabla}\bar{\phi}_{2}\bar{\phi}_{2})\cdot(\phi_{1}\vec{\nabla}\phi_{2}+\phi_{2}\vec{\nabla}\phi_{1})\Big{]}\,.$
(6)

All the interaction vertices are quartic in the fields, differently from the
classical RD field theory, where both cubic and quartic interaction vertices
are present. Furthermore, spatial gradients of the fields appear as a
consequence of fermionic statistics.

Kinetic equation.— In the reaction-limited regime $\hbar n\Gamma/D\ll 1$, the
dressed Green’s functions $\hat{G}$ can be rewritten in terms of the bare
Green’s functions $\hat{G}_{0}$ as the interaction vertices (6) are expanded
perturbatively around the quadratic sector (5). The sum of all internal one-
particle-irreducible contributions to the Feynman diagrams results in the
entries of the “self-energy” Keldysh-space matrix $\hat{\Sigma}$ [58, 59, 60,
61, 62, 63, 64, 71, 72]. The ensuing Dyson equation
$(\hat{G}_{0}^{-1}-\hat{\Sigma})\,\circ\,\hat{G}\,=\,\mathds{1}$ is
pictorially shown in Fig. 1(c).

The Keldysh component of the Dyson equation for
$G^{K}(\vec{x}_{1},t_{1},\vec{x}_{2},t_{1})$ determines the kinetic equation
[69]. The large-scale universal limit of this equation is best extracted by
performing a Fourier transform in the relative space [time] variable
$\vec{x}^{\prime}=\vec{x}_{1}-\vec{x}_{2}$ [$t^{\prime}=t_{1}-t_{2}$],
$G^{K}(\vec{x},t,\vec{k},\epsilon)$ 111In the whole manuscript, for lightness
of notation, we use the same symbol for a function and its Wigner transform.
The two functions can be distinguished from the corresponding arguments.,
i.e., the so-called “Wigner transform”[74, 75, 76]. The set of Wigner center-
of-mass coordinates $\vec{x}\,=\,(\vec{x_{1}}+\vec{x_{2}})/2$
[$t=(t_{1}+t_{2})/2$], as well as momentum $\vec{k}$ and energy $\epsilon$,
characterise the effective macroscopic evolution of Green’s functions.

Figure 2: Binary annihilation decay in $d$ dimensions. Solution of the
homogeneous Boltzmann equation (8) from the Fermi-sea initial state at density
$n_{0}$. The rescaled density $\tilde{n}=n/n_{0}$ decays algebraically as a
function of the dimensionless time $\tilde{t}=n_{0}^{1+2/d}\Gamma t$. From top
to bottom, algebraic decay $\tilde{n}\sim\tilde{t}^{-\frac{d}{d+1}}$ in
$d=1,2,3$ (solid lines). The dashed line represents the mean-field decay
$\tilde{n}(\tilde{t})\,\sim\,\tilde{t}^{-1}$ valid in infinite $d$.  Figure 3:
Binary annihilation inhomogeneous decay in $d=1$. (a) Double-well to harmonic
confinement potential quench: plot of the rescaled particle number
$\tilde{N}(\tilde{t})$ as a function of rescaled time $\tilde{t}$ for
increasing values of the parameter $\Omega$ (from top to bottom). The
algebraic time decay gets accelerated as $\Omega$ is increased towards the
mean-field prediction $\tilde{N}(\tilde{t})\sim t^{-1}$ (bottom dashed line).
(b) The corresponding rescaled spatial density
$\tilde{n}(\tilde{x},\tilde{t})$ profiles at increasing times (from top to
bottom) are plotted as a function of space $\tilde{x}$, with $\Omega=0.1$. For
both plots $C=A\omega\mu/(m^{2}\omega^{4})=0.8$, with the chemical potential
$\mu$ fixing the initial particle number $N_{0}$. (c) Trap release quench:
particle number decay $\tilde{N}(\tilde{t})$ versus time $\tilde{t}$ for
increasing parameter $\Omega$ (from bottom to top). A decay exponent
$\tilde{N}(\tilde{t})\sim\tilde{t}^{-\xi}$ is approximately observed at
intermediate times, with $\xi$ decreasing with $\Omega$, and $\xi=1/2$ for
$\Omega=0$ (black dashed line). At longer times a non-algebraic slow decay
sets in. (d) The associated rescaled spatial density
$\tilde{n}(\tilde{x},\tilde{t})$ profiles are reported at increasing times
(from top to bottom), with $\Omega=0.1$.

In the regime of weak space-time modulations,
$G^{K}(\vec{x},t,\vec{k},\epsilon)$ is a slow function of its arguments and
hence it is possible to perform a derivative expansion of the kinetic
equation, akin to the Moyal expansion carried on in hydrodynamics [77, 78, 79,
80, 81, 82]. Variations in time are also slow in the reaction-limited regime
$\hbar n\Gamma/D\ll 1$. To leading order in the interaction $\Gamma$, indeed,
the quasi-particle lifetime $\sim(n\Gamma)^{-1}$, given by the width of
$\mbox{Im}(G^{R})$, is much larger than the coherent time scale set by
$(D/\hbar)^{-1}$. This implies that $G^{K}(\vec{x},t,\vec{k},\epsilon)$ is a
sharply peaked function of $\epsilon$ around $\epsilon_{k}(\vec{x})$ and
therefore quasiparticles are well-defined even throughout the time evolution
[see Fig. 1(c)]. This allows to identify the equal-time Keldysh Green’s
function $G^{K}(\vec{x},t,\vec{k},t)$ as the emergent large-scale degree of
freedom:

$iG^{K}(\vec{x},\vec{k},t,t)\,=\,1\,-\,2\,n(\vec{x},\vec{k},t).$ (7)

Here, $n(\vec{x},\vec{k},t)$ is the one-body Wigner function, i.e., the
semiclassical phase-space $(\vec{x},\vec{k})$ occupation function [83, 84,
85]. At leading order in the space time-derivatives, it obeys the quantum
Boltzmann-like equation

$\begin{split}\Big{[}\,\partial_{t}&+\vec{v}_{g}(\vec{k})\cdot\vec{\nabla}_{x}-\vec{\nabla}_{x}V/\hbar\cdot\vec{\nabla}_{k}\,\Big{]}n(\vec{x},\vec{k},t)=\\\
&=-\Gamma\int\frac{d^{d}q}{(2\pi)^{d}}(\vec{k}-\vec{q})^{2}\,n(\vec{x},\vec{k},t)\,n(\vec{x},\vec{q},t).\end{split}$
(8)

Crucially, at first order in $\Gamma$, we find that $\hat{\Sigma}$ contributes
via purely imaginary terms, which determine the r.h.s., named collision
integral. The latter is computed in terms of tadpole Feynman diagrams
[depicted in Fig. 1(c)], which we consider at lowest order in the derivatives
$\vec{\nabla}$. The appearing factor $(\vec{k}-\vec{q})^{2}$ stems from the
fermionic statistics. Moreover, the dispersion relation
$\epsilon_{k}(\vec{x})$ ($\vec{v}_{g}(\vec{k})\,=\,2D\vec{k}/\hbar$ is the
group velocity) is not renormalized. In $d=1$, the r.h.s. of Eq. (8) has the
same form derived in Refs. [31, 32, 48, 49, 50] assuming the systems relaxes
to a TGGE state in-between consecutive reactions. Our Keldysh field-theory
analysis thus shows how the TGGE ansatz naturally emerges as the large-scale
limit of tadpole Feynman diagrams generated by the interaction vertices (6).
At the same time, it allows us to consider higher dimensional systems.

Homogeneous decay in $d$ dimensions.— For homogeneous initial states and no
trapping potential $V(\vec{x})=0$, the Wigner function $n(\vec{x},\vec{k},t)$
reduces to the momentum-occupation function $n(\vec{k},t)$. We consider a
Fermi-sea initial state, with equally populated modes up to some Fermi
momentum $k_{d}^{F}\,\sim\,n_{0}^{1/d}$ and a total initial density $n_{0}$.
The asymptotics of the particle density decay can then be worked out for
generic $d$ by reducing the $d$-dimensional problem to a radial one, i.e.,
$n(\vec{k},t)=n(k,t)$. Moreover, it is convenient to introduce the
adimensional re-scaled density $\tilde{n}\,=\,n/n_{0}$ and time
$\tilde{t}\,=\,n_{0}^{1+2/d}\Gamma t$. The long-time asymptotic for
$\tilde{n}(\tilde{t})$ is given by the power law 222$\alpha_{d}$ is a
geometric factor arising form angular integration. $\alpha_{d}\,=\,1$ for $d$
odd, $\alpha_{d}\,=\,\sqrt{\pi/2}$ for $d$ even.

$\tilde{n}(\tilde{t})\,\sim\,\Bigg{[}\,\frac{(\alpha_{d}\,\Theta_{d})^{2}}{2^{d}(d+1)^{d}(2\pi)^{2}}\,\Bigg{]}^{\frac{1}{d+1}}\,\tilde{t}^{-\frac{d}{d+1}},$
(9)

with $\Theta_{d}$ the $d$-dimensional solid angle. The solution of the
homogeneous Boltzmann equation for $d=1,2,3$ from the Fermi-sea initial state
is shown in Fig. 2. In $d=1$, the decay exponent is $1/2$ in agreement with
the TGGE prediction [31, 32, 48, 49, 50]. We find that the algebraic decay is
different from the mean-field $\tilde{t}^{-1}$ even for $d>1$ and it
approaches the latter only for $d\to\infty$. This result is surprising and
fundamentally different from the classical $A+A\to\emptyset$ RD dynamics.
Therein, non mean-field algebraic decay is possible only in $d=1$ in the
diffusion limited regime $\hbar\Gamma n/D\sim 1$ [14, 15, 16, 17, 21, 20, 22],
as a consequence of spatial fluctuations. Conversely, the different
dimensional dependence of the exponent in Eq. (9) does not emerge due to
spatial fluctuations but from the universal large-scale limit (8) of the
quantum dynamics.

Inhomogeneous decay in one dimension.— To show that our approach allows to
treat inhomogeneous settings, we consider a one-dimensional quantum quench of
the trapping potential from the prequench
$V_{0}(x)=Ax^{4}/4-m\omega^{2}x^{2}/2$ double-well to the postquench harmonic
$V(x)=m\omega^{2}x^{2}/2$ form. A similar setting has been consider in Ref.
[32], where also the prequench potential $V_{0}(x)$ is harmonic. The initial
condition is the local-density approximation of the ground state of the
Hamiltonian $H$ in Eq. (2) with potential $V_{0}(x)$. We set the initial
particle number to $N_{0}$ and we rescale variables through the parameters of
the harmonic potential $V(x)$: for the particle number $\tilde{N}=N/N_{0}$,
time $\tilde{t}=t\Gamma\big{(}2N_{0}m\omega/\hbar\big{)}^{3/2}/(2\pi)$, space
$\tilde{x}=x\sqrt{m\omega/(2N_{0}\hbar)}$, momentum
$\tilde{k}=k\sqrt{\hbar/(2N_{0}m\omega)}$ and density
$\tilde{n}=n\sqrt{\hbar/(2N_{0}m\omega)}$. The adimensional parameter
$\Omega=2\pi(D/\hbar N_{0})^{3/2}/(\Gamma\sqrt{\omega})$ quantifies the
interplay between reactions ($\Gamma)$, diffusion ($\hbar/m$) and confinement
($\omega$).

In Fig. 3(a), the decay of $\tilde{N}$ as a function of $\tilde{t}$ is shown.
The decay accelerates periodically as $\Omega$ is increased since particles
bounce off the potential walls and gather up at the centre of the well, as
shown in Fig. 3(b) for the density $\tilde{n}(\tilde{x},\tilde{t})$. As a
consequence of such breathing motion all decay profiles converge from the
short-time $\tilde{N}(\tilde{t})\sim\tilde{t}^{-1/2}$ algebraic behavior
towards the mean-field asymptotic decay
$\tilde{N}(\tilde{t})\sim\tilde{t}^{-1}$.

Next, we consider the long-time decay for the deconfinement dynamics of an
harmonic trap-release quench of the Fermi gas from
$V_{0}(x)=m\omega^{2}x^{2}/2$ to $V(x)=0$. We rescale also in this case
variables with respect to the harmonic potential $V_{0}(x)$ parameters. In
Fig. 3(c) the decay of $\tilde{N}$ in time $\tilde{t}$ is reported. One first
observes [69] an approximate algebraic decay
$\tilde{N}(\tilde{t})\sim\tilde{t}^{-\xi}$, with an exponent $\xi$
continuously decreasing as $\Omega$ is increased (from the value $\xi=1/2$ at
$\Omega=0$). At longer times the algebraic decay is replaced by a much slower
non-algebraic decay. As the particle density $\tilde{n}(\tilde{x},\tilde{t})$
spreads in free space, shown in Fig. 3(d), the most relevant reactions,
indeed, take place between particles with same momenta. The fermionic
statistics, manifest in the factor $(\tilde{k}-\tilde{q})^{2}$ in Eq. (8),
thereby suppress these reactions and determine the unexpectedly slow large-
scale behavior of Figs. 3(c)-(d).

Summary.— We provided a Keldysh field-theory description of quantum RD
dynamics of binary annihilation $A+A\to\emptyset$. We analytically derived the
universal large-scale Boltzmann equation in arbitrary dimension $d$ describing
the reaction-limited regime of slow reactions $\hbar n\Gamma/D\ll 1$. In
$d=1$, our results are equivalent to the TGGE ansatz, connecting the latter to
the field-theoretical diagrammatic expansions. For homogeneous systems, we
analytically showed that the density algebraic decay exponent features an
unexpected dependency on $d$ and it deviates from mean-field value even in
$d>1$, in contrast with classical RD dynamics. In one-dimensional
inhomogeneous setups involving a trapping potential, we found that the decay
is either accelerated towards the mean-field value (confined systems), or
severely slowed down (deconfined systems). Our results find a natural
application in cold-atomic experiments of fermionic gases involving two-body
losses [51, 52, 53]. From the formulation here proposed, several relevant
questions can be addressed. As an example, one can assess the impact of
elastic-Hamiltonian collisions on the decay exponent [52, 53, 87, 88, 89].
Away from the reaction-limited, it is also crucial to study the quantum
diffusion-limited regime $\hbar n\Gamma/D\sim 1$ via renormalization group
schemes, as done for the classical RD [17, 21, 20, 22].

Aknowledgments.— We aknowledge fruitful discussion with M. Buchhold, S. Diehl
and J.P. Garrahan. F.G. acknowledges financial support from the University of
Trento and Collegio Clesio. G.P. acknowledges support from the Alexander von
Humboldt Foundation through a Humboldt research fellowship for postdoctoral
researchers. We acknowledge financial support in part from EPSRC Grant no.
EP/R04421X/1 and EPSRC Grant no. EP/V031201/1. We are also grateful for
funding from the Deutsche Forschungsgemeinsschaft (DFG, German Research
Foundation) under Project No. 435696605, the Research Unit FOR 5413/1, Grant
No. 465199066 and the Research Unit FOR 5522/1, Grant No. 499180199.

## References

  * Henkel  _et al._ [2008] M. Henkel, H. Hinrichsen, S. Lübeck, and M. Pleimling,  _Non-equilibrium phase transitions_, Vol. 1 (Springer, 2008).
  * Hinrichsen [2000] H. Hinrichsen, Adv. Phys. 49, 815 (2000).
  * Vladimir [1997] P. Vladimir,  _Nonequilibrium statistical mechanics in one dimension_ (Cambridge University Press, 1997).
  * Ovchinnikov and Zeldovich [1978] A. Ovchinnikov and Y. Zeldovich, Chem. Phys. 28, 215 (1978).
  * Kang and Redner [1984a] K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 (1984a).
  * Kang and Redner [1984b] K. Kang and S. Redner, Phys. Rev. A 30, 2833 (1984b).
  * Kang and Redner [1985] K. Kang and S. Redner, Phys. Rev. A 32, 435 (1985).
  * Spouge [1988] J. L. Spouge, Phys. Rev. Lett. 60, 871 (1988).
  * Torney and McConnell [1983] D. C. Torney and H. M. McConnell, J. Phys. Chem. 87, 1941 (1983).
  * Privman [1994] V. Privman, Phys. Rev. E 50, 50 (1994).
  * Toussaint and Wilczek [1983] D. Toussaint and F. Wilczek, J. Chem. Phys. 78, 2642 (1983).
  * Rácz [1985] Z. Rácz, Phys. Rev. Lett. 55, 1707 (1985).
  * Doi [1976] M. Doi, J. Phys. A: Math. Gen. 9, 1465 (1976).
  * Peliti [1985] L. Peliti, J. Phys. France 46, 1469 (1985).
  * Peliti [1986] L. Peliti, J. Phys. A: Math. Gen. 19, L365 (1986).
  * Lee [1994] B. P. Lee, J. Phys. A: Math. Gen. 27, 2633 (1994).
  * Cardy [1996] J. Cardy, arXiv:cond-mat/9607163 (1996).
  * Mattis and Glasser [1998] D. C. Mattis and M. L. Glasser, Rev. Mod. Phys. 70, 979 (1998).
  * Grassberger and Scheunert [1980] P. Grassberger and M. Scheunert, Fortschritte der Phys. 28, 547 (1980).
  * Tauber [2002] U. C. Tauber, arXiv:cond-mat/0205327 (2002).
  * Täuber  _et al._ [2005] U. C. Täuber, M. Howard, and B. P. Vollmayr-Lee, J. Phys. A: Math. Gen. 38, R79 (2005).
  * Täuber [2014] U. C. Täuber,  _Critical dynamics: a field theory approach to equilibrium and non-equilibrium scaling behavior_ (Cambridge University Press, 2014).
  * Privman and Grynberg [1992] V. Privman and M. D. Grynberg, J. Phys. A. Math. Gen. 25, 6567 (1992).
  * van Horssen and Garrahan [2015] M. van Horssen and J. P. Garrahan, Phys. Rev. E 91, 032132 (2015).
  * Gillman  _et al._ [2019] E. Gillman, F. Carollo, and I. Lesanovsky, New J. Phys. 21, 093064 (2019).
  * Gillman  _et al._ [2020] E. Gillman, F. Carollo, and I. Lesanovsky, Phys. Rev. Lett. 125, 100403 (2020).
  * Jo  _et al._ [2021] M. Jo, J. Lee, K. Choi, and B. Kahng, Phys. Rev. Research 3, 013238 (2021).
  * Carollo  _et al._ [2019] F. Carollo, E. Gillman, H. Weimer, and I. Lesanovsky, Phys. Rev. Lett. 123, 100604 (2019).
  * Carollo and Lesanovsky [2022] F. Carollo and I. Lesanovsky, Phys. Rev. Lett. 128, 040603 (2022).
  * Carollo  _et al._ [2022] F. Carollo, M. Gnann, G. Perfetto, and I. Lesanovsky, Phys. Rev. B 106, 094315 (2022).
  * Rossini  _et al._ [2021] D. Rossini, A. Ghermaoui, M. B. Aguilera, R. Vatré, R. Bouganne, J. Beugnon, F. Gerbier, and L. Mazza, Phys. Rev. A 103, L060201 (2021).
  * Rosso  _et al._ [2022] L. Rosso, A. Biella, and L. Mazza, SciPost Phys. 12, 44 (2022).
  * Bouchoule  _et al._ [2020] I. Bouchoule, B. Doyon, and J. Dubail, SciPost Phys. 9, 44 (2020).
  * Bouchoule and Dubail [2021] I. Bouchoule and J. Dubail, Phys. Rev. Lett. 126, 160603 (2021).
  * Bouchoule and Dubail [2022] I. Bouchoule and J. Dubail, J. Stat. Mech.: Theory Exp. 2022, 014003 (2022).
  * Rosso  _et al._ [2023] L. Rosso, A. Biella, J. De Nardis, and L. Mazza, Phys. Rev. A 107, 013303 (2023).
  * 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, New J. Phys. 11, 013053 (2009).
  * Everest  _et al._ [2014] B. Everest, M. R. Hush, and I. Lesanovsky, Phys. Rev. B 90, 134306 (2014).
  * Huang  _et al._ [2023] C.-H. Huang, T. Giamarchi, and M. A. Cazalilla, arXiv:2305.13090 (2023).
  * Hamanaka  _et al._ [2023] S. Hamanaka, K. Yamamoto, and T. Yoshida, arXiv:2305.19697 (2023).
  * Burrows  _et al._ [2017] K. A. Burrows, H. Perrin, and B. M. Garraway, Phys. Rev. A 96, 023429 (2017).
  * Bouchoule and Schemmer [2020] I. Bouchoule and M. Schemmer, SciPost Phys. 8, 60 (2020).
  * Traverso  _et al._ [2009] A. Traverso, R. Chakraborty, Y. N. Martinez de Escobar, P. G. Mickelson, S. B. Nagel, M. Yan, and T. C. Killian, Phys. Rev. A 79, 060702 (2009).
  * Yamaguchi  _et al._ [2008] A. Yamaguchi, S. Uetake, D. Hashimoto, J. M. Doyle, and Y. Takahashi, Phys. Rev. Lett. 101, 233002 (2008).
  * Kinoshita  _et al._ [2005] T. Kinoshita, T. Wenger, and D. S. Weiss, Phys. Rev. Lett. 95, 190406 (2005).
  * Söding  _et al._ [1999] J. Söding, D. Guéry-Odelin, P. Desbiolles, F. Chevy, H. Inamori, and J. Dalibard, Appl. Phys. B 69, 257 (1999).
  * Tolra  _et al._ [2004] B. L. Tolra, K. M. O’Hara, J. H. Huckans, W. D. Phillips, S. L. Rolston, and J. V. Porto, Phys. Rev. Lett. 92, 190401 (2004).
  * Perfetto  _et al._ [2023a] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. Lett. 130, 210402 (2023a).
  * Perfetto  _et al._ [2023b] G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, arXiv:2305.06944 (2023b).
  * Riggio  _et al._ [2023] F. Riggio, L. Rosso, D. Karevski, and J. Dubail, arXiv:2307.02298 (2023).
  * Yan  _et al._ [2013] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Nature 501, 521 (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, Phys. Rev. Lett. 112, 070404 (2014).
  * Sponselee  _et al._ [2018] K. Sponselee, L. Freystatzky, B. Abeln, M. Diem, B. Hundt, A. Kochanke, T. Ponath, B. Santra, L. Mathey, K. Sengstock, _et al._ , Quantum Sci. Technol. 4, 014002 (2018).
  * Lange  _et al._ [2018] F. Lange, Z. Lenarčič, and A. Rosch, Phys. Rev. B 97, 165138 (2018).
  * Mallayya  _et al._ [2019] K. Mallayya, M. Rigol, and W. De Roeck, Phys. Rev. X 9, 021027 (2019).
  * Lange  _et al._ [2017] F. Lange, Z. Lenarčič, and A. Rosch, Nat. Commun. 8, 1 (2017).
  * Lenarčič  _et al._ [2018] Z. Lenarčič, F. Lange, and A. Rosch, Phys. Rev. B 97, 024302 (2018).
  * Schwinger [1961] J. Schwinger, J. Math. Phys. 2, 407 (1961).
  * Keldysh [1965] L. V. Keldysh, Sov. Phys. JETP 20, 1018 (1965).
  * Kamenev [2011] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (Cambridge University Press, 2011).
  * Altland and Simons [2010] A. Altland and B. D. Simons,  _Condensed Matter Field Theory_, 2nd ed. (Cambridge University Press, 2010).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Rep. Prog. Phys. 79, 096001 (2016).
  * Tonielli [2016] F. Tonielli,  _Keldysh Field Theory for dissipation-induced states of Fermions_ (University of Pisa, 2016).
  * Thompson and Kamenev [2023] F. Thompson and A. Kamenev, Ann. Phys. , 169385 (2023).
  * Lindblad [1976] G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer Science & Business Media, 2004).
  * [69] See Supplemental Material at [URL will be inserted by publisher] for the details of the calculations.
  * Larkin and Ovchinnikov [1975] A. Larkin and Y. Ovchinnikov, Sov. Phys. JETP 41, 960 (1975).
  * Buchhold and Diehl [2015] M. Buchhold and S. Diehl, Eur. Phys. J. D 69, 1 (2015).
  * Duval and Cherroret [2023] C. Duval and N. Cherroret, Phys. Rev. A 107, 043305 (2023).
  * Note [1] In the whole manuscript, for lightness of notation, we use the same symbol for a function and its Wigner transform. The two functions can be distinguished from the corresponding arguments.
  * Weyl [1931] H. Weyl,  _The theory of groups and quantum mechanics_, Vol. 1951 (Dover publications, 1931).
  * Wigner [1932] E. Wigner, Phys. Rev. 40, 749 (1932).
  * Case [2008] W. B. Case, Am. J. Phys. 76, 937 (2008).
  * Moyal [1949] J. E. Moyal, Proc. Cambridge Philos. Soc. 45, 99–124 (1949).
  * Fagotti [2017] M. Fagotti, Phys. Rev. B 96, 220302 (2017).
  * Fagotti [2020] M. Fagotti, SciPost Phys. 8, 048 (2020).
  * Essler [2022] F. H. Essler, Phys. A: Stat. , 127572 (2022).
  * Dean  _et al._ [2019] D. S. Dean, P. L. Doussal, S. N. Majumdar, and G. Schehr, EPL 126, 20006 (2019).
  * De Nardis and Doyon [2023] J. De Nardis and B. Doyon, J. Phys. A: Math. 56, 245001 (2023).
  * Bettelheim  _et al._ [2006] E. Bettelheim, A. G. Abanov, and P. Wiegmann, Phys. Rev. Lett. 97, 246402 (2006).
  * Bettelheim and Wiegmann [2011] E. Bettelheim and P. B. Wiegmann, Phys. Rev. B 84, 085102 (2011).
  * Bettelheim and Glazman [2012] E. Bettelheim and L. Glazman, Phys. Rev. Lett. 109, 260602 (2012).
  * Note [2] $\alpha_{d}$ is a geometric factor arising form angular integration. $\alpha_{d}\,=\,1$ for $d$ odd, $\alpha_{d}\,=\,\sqrt{\pi/2}$ for $d$ even.
  * Fürst  _et al._ [2013] M. L. Fürst, J. Lukkarinen, P. Mei, and H. Spohn, J. Phys. A: Math. 46, 485002 (2013).
  * Fürst  _et al._ [2013] M. L. R. Fürst, C. B. Mendl, and H. Spohn, Phys. Rev. E 88, 012108 (2013).
  * Zechmann  _et al._ [2022] P. Zechmann, A. Bastianello, and M. Knap, Phys. Rev. B 106, 075115 (2022).

SUPPLEMENTAL MATERIAL

Large-scale universality in Quantum Reaction-Diffusion from Keldysh field
theory

Federico Gerbino,1,2 Igor Lesanovsky,1,3 and Gabriele Perfetto1

1Institut für Theoretische Physik, Universität Tübingen,
Auf der Morgenstelle 14, 72076 Tübingen, Germany
2Physics Department, University of Trento, via Sommarive, 14 I-38123 Trento,
Italy
3School of Physics and Astronomy and Centre for the Mathematics
and Theoretical Physics of Quantum Non-Equilibrium Systems,
The University of Nottingham, Nottingham, NG7 2RD, United Kingdom

This Supplemental Material provides additional information on the model and
the calculations at the basis of the results presented in the main text. In
Sec. I, we derive the continuum-space Lindblad dynamics, Eqs. (1)-(4) of the
main text, by taking the continuum limit of the corresponding lattice model.
In Sec. II, we briefly recall the basic aspects of the Keldysh field theory
for open quantum systems, which are needed for the understanding of the action
formulation in Eqs. (5) and (6) of the main text. In Sec. III, we briefly
report the main steps underlying the diagrammatic calculation of the self-
energy Keldysh matrix at first order in $\Gamma$ and the associated derivation
of the Boltzmann equation in Eq. (8) of the main text. In Sec. IV, we compute
the effective exponent for the deconfining dynamics ensuing from the trap-
release quench of the trapping potential (see Fig. 3 of the main text).

##  I Lindblad dynamics in the space-continuum limit

We start by introducing, for the sake of illustrative purposes, a one-
dimensional lattice model, which is pictorially represented in Fig. S1(a). Let
us consider an infinite chain whose sites labelled by index $j\in\mathbb{Z}$
can be either occupied by one fermion
$n_{j}\ket{...\bullet_{j}...}=\ket{...\bullet_{j}...}$ or empty
$n_{j}\ket{...\circ_{j}...}=0$. The number operator
$n_{j}=c_{j}^{\dagger}c_{j}$ is written in terms of fermionic destruction
($c_{j}$) and creation ($c_{j}^{\dagger}$) operators satisfying the
anticommutation relation $\\{c_{i},c_{j}^{\dagger}\\}=\delta_{i,j}$. We also
introduce the lattice spacing $a$, i.e., the shortest length scale of the
problem. The hopping Hamiltonian reads

$H\,=\,-\,\frac{D}{a^{2}}\,\sum_{j}\,(\,c_{ja}^{\dagger}c_{(j+1)a}+c^{\dagger}_{(j+1)a}c_{ja}\,)\,+\,\sum_{j}\,c^{\dagger}_{ja}V_{ja}c_{ja}\,.$
(S1)

We remark that we here follow the convention of Ref. [21] for the hopping
amplitude notation so that $D/\hbar$ has units
[$\mbox{length}^{2}/\mbox{time}$] of a diffusion constant, while
$D/(a^{2}\hbar)$ has units [$\mbox{time}^{-1}$] of a rate. We further include
the effect of a position-dependent trapping potential $V_{ja}$. In the absence
of the latter, the fermionic hopping Hamiltonian (S1) has been studied in
Refs. [24, 48] for quantum RD models on a chain. The reaction part is encoded
into the jump operators of the Lindbladian and it is given by two-body
annihilation process $A+A\to\emptyset$. This destroys two neighbouring
fermions at a rate $\Gamma/a^{3}$ ($\Gamma$ has units
[$\mbox{length}^{3}/\mbox{time}$])

$\dot{\rho}(t)=-i[H,\rho(t)]+\sum_{j}\left[L_{ja}\rho{L_{ja}}^{\dagger}-\frac{1}{2}\left\\{{L_{ja}}^{\dagger}L_{ja},\rho\right\\}\right],\quad\mbox{with}\quad
L_{ja}=\sqrt{\frac{\Gamma}{a^{3}}}\,c_{ja}c_{(j+1)a}.$ (S2)

Clearly, for fermions, one cannot define reactions occurring on the same site
as this would violate the exclusion principle. The annihilation process
between neighbouring sites is sketched in Fig. S1(a).

The generalization of Eqs. (S1) and (S2) to generic space dimension $d$ can be
simply accomplished by introducing a $d-$dimensional hypercubic lattice.
Lattice points on this lattice are identified by a $d-$dimensional vector
$\textbf{j}=(j_{1},\dots j_{\alpha},\dots j_{d})$ of integer numbers
$j_{\alpha}\in\mathbb{Z}$ ($\alpha=1,2\dots d$). Then one replaces $c_{ja}$
and $c_{(j+1)a}$ (the same applying for the corresponding creation operators)
with $c_{\textbf{j}a}$ and $c_{\textbf{j}a+a\textbf{e}_{\alpha}}$,
$\textbf{e}_{\alpha}$ being the unit-vector pointing to direction
$\alpha\,=\,1,...,d$. The Hamiltonian and the Lindblad master equation then
read

$H=-\,\frac{D}{a^{2}}\,\sum_{\alpha=1}^{d}\,\sum_{\textbf{j}}\,(\,c_{\textbf{j}a}^{\dagger}c_{\textbf{j}a+a\textbf{e}_{\alpha}}+c_{\textbf{j}a+a\textbf{e}_{\alpha}}^{\dagger}c_{\textbf{j}a})\,+\,\sum_{\alpha=1}^{d}\,\sum_{\textbf{j}}\,c_{\textbf{j}a}^{\dagger}V_{\textbf{j}a}c_{\textbf{j}a},$
(S3)

and

$\dot{\rho}(t)=-i[H,\rho(t)]+\sum_{\alpha=1}^{d}\sum_{\textbf{j}}\left[L_{\textbf{j}a}^{(\alpha)}\rho{L_{\textbf{j}a}^{(\alpha)}}^{\dagger}-\frac{1}{2}\left\\{{L_{\textbf{j}a}^{(\alpha)}}^{\dagger}L_{\textbf{j}a}^{(\alpha)},\rho\right\\}\right],\quad\mbox{with}\quad
L_{\textbf{j}\alpha}^{(\alpha)}=\sqrt{\frac{\Gamma}{a^{d+2}}}\,c_{\textbf{j}a}c_{\textbf{j}a+a\textbf{e}_{\alpha}}\,.$
(S4)

Equation (S3) shows that particles can hop isotropically (with the same
amplitude) to one of the neighbouring sites in any direction $\alpha$.
Annihilation reactions, according to Eq. (S4), also take place isotropically.

Let us take the continuum limit of this model defined by the substitutions
$\textbf{j}a\,\rightarrow\,\vec{x}$,
$\sum_{\textbf{j}}\,\rightarrow\,\int\,d^{d}x/a^{d}$ and by considering
$a\,\rightarrow\,0$. We may then write:

$\lim_{a\rightarrow
0}\,\frac{L_{\textbf{j}a}^{(\alpha)}}{a^{d/2}}=\lim_{a\rightarrow
0}\,\sqrt{\frac{\Gamma}{a^{d}}}\,\frac{c_{\textbf{j}a}c_{\textbf{j}a+a\textbf{e}_{\alpha}}}{a^{d/2+1}}=\lim_{a\rightarrow
0}\,\sqrt{\frac{\Gamma}{a^{d}}}\,\frac{c_{\textbf{j}a}}{a^{d/2}}\,\frac{c_{\textbf{j}a+a\textbf{e}_{\alpha}}-c_{\textbf{j}a}}{a}.$
(S5)

Introducing the field operators

$\psi(\vec{x})=\frac{c_{\textbf{j}a}}{a^{d/2}},\quad\mbox{with}\quad\\{\psi(\vec{x}),\psi^{\dagger}(\vec{x}^{\prime})\\}=\delta(\vec{x}-\vec{x}^{\prime}),$
(S6)

one recognises that Eq. (S5) is the definition of a partial derivative in
direction $x_{\alpha}$, namely:

$\lim_{a\rightarrow
0}\,\frac{L_{\textbf{j}a}^{(\alpha)}}{a^{d/2}}=\sqrt{\Gamma}\,\psi(\vec{x})\,\frac{\partial\psi}{\partial
x_{\alpha}}=L_{\alpha}(\vec{x})\,.$ (S7)

Repeating the same procedure for the Hamiltonian (S3), the definition of a
second partial derivative $\frac{\partial^{2}}{\partial x_{\alpha}^{2}}$
appears

$H=\int\,d^{d}x\,\psi^{\dagger}(\vec{x})\,[\,-D\nabla^{2}+V(\vec{x})\,]\,\psi(\vec{x}).$
(S8)

This equation is readily recognized as Eq. (2) of the main text. For the
dissipator in Eq. (S4), we also have

$\dot{\rho}(t)=-i[H,\rho(t)]+\Gamma\,\int\,d^{d}x\,\Big{[}\,\psi\vec{\nabla}\psi\,\rho\,\cdot\,\vec{\nabla}\psi^{\dagger}\psi^{\dagger}\,-\,\\{\,\vec{\nabla}\psi^{\dagger}\psi^{\dagger}\,\cdot\,\psi\vec{\nabla}\psi,\,\rho\,\\}\,\Big{]}.$
(S9)

Writing the dot product $\cdot$ explicitly in terms of the associated
$\alpha=1,2\dots d$ components, Eq. (S9) is recognized as Eq. (3) of the main
text. In Fig. S1(b), we pictorially show the continuum formulation of the
model (in one dimension for the sake of simplicity in the illustration).

Figure S1: Space-continuum limit. (a) Quantum RD dynamics on a one-dimensional
lattice with lattice spacing $a$. Fermions coherently hop (blue-wiggly lines)
between neighbouring sites of an optical lattice (black wavy line) at rate
$D/(a^{2}\hbar)$. Binary annihilation $A+A\to\emptyset$ reactions take place
between neighbouring sites at rate $\Gamma/a^{3}$, as sketched in yellow. (b)
Taking the continuum limit $a\,\rightarrow\,0$ of the discrete-lattice model,
defines the continuous quantum RD dynamics under examination in Eqs. (S8) and
(S9). Here, $D/\hbar$ has units [$\mbox{length}^{2}/\mbox{time}$] and it
controls coherent motion, while binary annihilation reactions occur via
derivative couplings of the fields with coupling constant $\Gamma$, with units
[$\mbox{length}^{3}/\mbox{time}$].

##  II Keldysh path integral formulation of open systems

We report here the main steps underlying the construction of Keldysh field
theory for dissipative many-body quantum systems [61, 62, 63, 60, 64]. We then
report the main definitions concerning the associated Green functions, which
are necessary to understand the derivation of the quantum Boltzmann equation
in the main text.

Let us consider the formal solution to a generic quantum Master equation of
the form, e.g., in Eq. (S4). We start by considering a zero-dimensional system
with no spatial structure (no sum over lattice sites $\mathbf{j}$). We can
Trotter-decompose the time-evolution operator in infinitely many time slices,
thus writing:

$\rho(t)=e^{\mathcal{L}t}[\rho(0)]=\lim_{N\rightarrow\infty}\,(\,1+\mathcal{L}\delta
t\,)^{N}[\rho(0)],$ (S10)

with $\delta t\,=\,t/N$. Hence, at any time instant $t_{n}\,=\,n\delta t$,
with $n=1,2\dots N$, the subsequent density matrix operator evolves via:

$\rho_{n+1}=(\,1+\mathcal{L}\delta t\,)[\rho_{n}]+\,O(\delta t^{2}).$ (S11)

It is now possible to introduce two sets of fermionic coherent states, labeled
by Grassmann fields $\psi_{n}^{\pm}$ which are eigenvalues of the fermionic
destruction operator. From the resolution for the identity operator one has:

$\rho_{n}=\int\,d\psi_{n}^{+}\,d\bar{\psi}_{n}^{+}d\psi_{n}^{-}\,d\bar{\psi}_{n}^{-}\,e^{-\bar{\psi}^{+}_{n}\psi^{+}_{n}-\bar{\psi}^{-}_{n}\psi^{-}_{n}}\,\bra{\psi_{n}^{+}}\rho_{n}\ket{-\psi_{n}^{-}}\,\ket{\psi_{n}^{+}}\bra{-\psi_{n}^{-}}.$
(S12)

We label $\ket{\psi_{n}^{+}}$ the states acting at time slice $n$ on the
density operator from the left, hence evolving the latter in the forward
direction along the so-called forward (+) contour. Conversely, states
$\ket{\psi_{n}^{-}}$ act on the density operator from the right and thus
evolve backward in time along a backward (-) contour (cf. the central panel of
Fig. (1) of the main text). Doubling the field variables from two to four
fields evolving along anti-parallel time axes is a peculiarity of out-of-
equilibrium quantum field theories [61, 62, 63, 60, 64]. Besides, an
additional minus sign in front of $\ket{-\psi_{n}^{-}}$ and
$\bra{-\psi_{n}^{-}}$ coherent states has been added with no effect on the
result of the calculation. Inserting Eq. (S12) into (S11), we arrive to:

$\begin{split}\bra{\psi_{n+1}^{+}}\rho_{n+1}\ket{-\psi_{n+1}^{-}}=\int\,d\psi_{n}^{+}d\bar{\psi}_{n}^{+}\,&d\psi_{n}^{-}d\bar{\psi}_{n}^{-}\,e^{-\bar{\psi}^{+}_{n}\psi^{+}_{n}-\bar{\psi}^{-}_{n}\psi^{-}_{n}}\\\
&\bra{\psi_{n}^{+}}\rho_{n}\ket{-\psi_{n}^{-}}\,\bra{\psi_{n+1}^{+}}\big{\\{}(1+\mathcal{L}\delta
t)\left[\ket{\psi_{n}^{+}}\bra{-\psi_{n}^{-}}\right]\big{\\}}\ket{-\psi_{n+1}^{-}}\,+\,O(\delta
t^{2})\,.\end{split}$ (S13)

Evaluation of the matrix element in squared brackets, together with the
normalisation factors of coherent states completeness relations
$e^{-\bar{\psi}^{+}_{n}\psi^{+}_{n}-\bar{\psi}^{-}_{n}\psi^{-}_{n}}$, leads
to:

$\begin{split}\bra{\psi_{n+1}^{+}}\rho_{n+1}\ket{-\psi_{n+1}^{-}}=\int\,d\psi_{n}^{+}d\bar{\psi}_{n}^{+}\,&d\psi^{-}_{n}d\bar{\psi}^{-}_{n}\,\bra{\psi_{n}^{+}}\rho_{n}\ket{-\psi_{n}^{-}}\\\
&\mathrm{exp}\Big{\\{}i\delta
t\Big{[}-i\partial_{t}\bar{\psi}^{+}_{n}\psi^{+}_{n}-i\bar{\psi}^{-}_{n}\partial_{t}\psi^{-}_{n}-i\mathcal{L}[\psi_{n}^{+},\bar{\psi}_{n+1}^{+},\psi_{n+1}^{-},\bar{\psi}_{n}^{-}]\Big{]}\Big{\\}}\,+\,O(\delta
t^{2}).\end{split}$ (S14)

In the previous equation, the function
$\mathcal{L}[\psi_{n}^{+},\bar{\psi}_{n+1}^{+},\psi_{n+1}^{-},\bar{\psi}_{n}^{-}]$
is the Lindbladian evaluated in terms of Grassmann fields after the
Hamiltonian and the jump operators have been normally ordered so that creation
operators always lie to the left of destruction ones. We also introduced the
notation $\partial_{t}\psi_{n}\,=\,\frac{\psi_{n+1}-\psi_{n}}{\delta t}$.

Repeated iteration of the previous equation from time slice $n\,=\,1$ to
$n\,=\,N$ connects the initial time to the final time density matrix. One can
then take the limit $N\,\rightarrow\,\infty$, thus neglecting the $O(\delta
t^{2})$ terms and performing the double substitution $\sum_{n}\delta
t\,\rightarrow\,\int dt$, $\prod_{n}d\psi_{n}\,\rightarrow\,\mathbf{D}\psi$.
The trace operation joins the forward and backward contour leading to the
Keldysh partition function:

$Z(t)=\int\mathbf{D}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]\,\bra{\psi_{+}(0)}\rho(0)\ket{-\psi_{-}(0)}\,e^{iS[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}\,.$
(S15)

In the following, we will neglect the boundary term
$\bra{\psi_{+}(0)}\rho(0)\ket{-\psi_{-}(0)}$ as it does not affect the large-
scale Boltzmann equation, which is solely determined by the bulk part of the
action. We therefore focus on the bulk part of the action only and we require
the associated normalisation constraint

$Z(t)=\mathrm{tr}[\rho(t)]=\int\mathbf{D}[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]\,e^{iS[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]}=1.$
(S16)

The bulk part of the Keldysh action is denoted as $S$ and it reads

$S[\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}]=\int_{0}^{t}\,dt^{\prime}\,(\,\bar{\psi}_{+}i\partial_{t^{\prime}}\psi_{+}-\bar{\psi}_{-}i\partial_{t^{\prime}}\psi_{-}-i\mathcal{L}(\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-}),$
(S17)

with the Lindbladian:

$\mathcal{L}(\psi_{+},\bar{\psi}_{+},\psi_{-},\bar{\psi}_{-})\,=\,-i(H_{+}-H_{-})\,+\,\sum_{\alpha}\,\big{[}\,\bar{L}_{\alpha,-}L_{\alpha,+}-\frac{1}{2}(\bar{L}_{\alpha,+}L_{\alpha,+}+\bar{L}_{\alpha,-}L_{\alpha,-})\,\big{]}.$
(S18)

Functions $H_{\pm}=H(\bar{\psi}_{\pm},\psi_{\pm})$,
$L_{\alpha,\pm}=L_{\alpha}(\bar{\psi}_{\pm},\psi_{\pm})$ and
$\bar{L}_{\alpha,\pm}L_{\alpha,\pm}=L^{\dagger}_{\alpha}L_{\alpha}(\bar{\psi}_{\pm},\psi_{\pm})$
are evaluated on the forward (+) or backward (-) contour after normal
ordering, as explained also after Eq. (S14).

Clearly, the field theory described by the preceding Keldysh action applies to
zero-dimensional systems without a space structure. However, one can
generalise the derivation to spatially-extended systems, with the
representation of the identity via a set of coherent states
$\ket{\\{\psi_{n}\\}}=\ket{\psi_{n,\mathbf{1}}}\ket{\psi_{n,\mathbf{2}}}\dots\ket{\psi_{n,\mathbf{j}}}\dots$
labelled on a time slice $n$ by the lattice site coordinate
$\mathbf{1},\mathbf{2}\dots\mathbf{j}\dots$ (the bold face notation
representing a vector of integer coordinates with the notation introduced
after Eq. (S2)):

$\mathbb{I}=\int\prod_{\textbf{j}}d\psi_{n,\textbf{j}}d\bar{\psi}_{n,\textbf{j}}e^{-\sum_{\textbf{j}}\bar{\psi}_{n,\textbf{j}}\psi_{n,\textbf{j}}}\ket{\\{\psi_{n}\\}}\bra{\\{\psi_{n}\\}}.$
(S19)

In this case, the continuum-space limit turns the sum over lattice sites j
into a space integral, while it includes in the functional measure
$\mathbf{D}[\psi,\bar{\psi}]$ an infinite product over the spatial field
configurations evaluated at each time slice. Accordingly, the Keldysh action
for the fermionic binary annihilation dynamics, described by the Lindblad
equation in Eqs. (S8) and (S9), in the $\pm$ basis reads:

$\displaystyle S=S_{\mathrm{diff}}+S_{\mathcal{D}}$ $\displaystyle=\int
d^{d}xdt\,\frac{1}{\hbar}\,\Big{[}\,\bar{\psi}_{+}(i\hbar\partial_{t}+D\nabla^{2}-V(\vec{x}))\psi_{+}-\bar{\psi}_{-}(i\hbar\partial_{t}+D\nabla^{2}-V(\vec{x}))\psi_{-}\,\Big{]}+$
$\displaystyle+i\Gamma\int
d^{d}x\,dt\,\Big{[}\,\frac{1}{2}\vec{\nabla}\bar{\psi}_{+}\bar{\psi}_{+}\cdot\psi_{+}\vec{\nabla}\psi_{+}+\frac{1}{2}\vec{\nabla}\bar{\psi}_{-}\bar{\psi}_{-}\cdot\psi_{-}\vec{\nabla}\psi_{-}-\vec{\nabla}\bar{\psi}_{-}\bar{\psi}_{-}\cdot\psi_{+}\vec{\nabla}\psi_{+}\,\Big{]}\,,$
(S20)

where the r.h.s. of the second equality on the first line is the quadratic
action $S_{\mathrm{diff}}$ describing coherent motion in a continuous
$d$-dimensional space (S8). The quartic action $S_{\mathcal{D}}$ on the second
lines represents, instead, the interaction vertices due to nearest-neighbour
reactions (S7). By means of the Keldysh rotation of Eq. (S23) which we shall
later introduce, the former line gives rise to the diffusion action in Eq.
(5), while the latter generates the dissipative sector discussed in Eq. (6).
The two-point Green’s functions associated to the Keldysh action are of
central relevance. We list them here as they will be needed in the derivation
of the quantum Boltzmann equation in Sec. III. We use henceforth in the
Supplemental Material, the compact 4-vector notation $x=(\vec{x},t)$:

$\displaystyle
iG^{T}(x_{1},x_{2})=\langle\psi_{+}(x_{1})\bar{\psi}_{+}(x_{2})\rangle=\braket{T\left[c(x_{1})c^{\dagger}(x_{2})\right]},$
(S21a) $\displaystyle
iG^{<}(x_{1},x_{2})=\langle\psi_{+}(x_{1})\bar{\psi}_{-}(x_{2})\rangle=-\langle
c^{\dagger}(x_{2})c(x_{1})\rangle,$ (S21b) $\displaystyle
iG^{>}(x_{1},x_{2})=\langle\psi_{-}(x_{1})\bar{\psi}_{+}(x_{2})\rangle=\braket{c(x_{1})c^{\dagger}(x_{2})},$
(S21c) $\displaystyle
iG^{\tilde{T}}(x_{1},x_{2})=\langle\psi_{-}(x_{1})\bar{\psi}_{-}(x_{2})\rangle=\braket{\widetilde{T}\left[c(x_{1})c^{\dagger}(x_{2})\right]}.$
(S21d)

In Eqs. (S21a) and (S21d)) the symbols $T$ and $\widetilde{T}$ denote time
ordering and anti-time ordering along the Keldysh contour, respectively (the
operator evaluated at the larger/smaller time goes to the left, respectively,
possibly taking a minus sign if a permutation is needed). Equation (S21) also
shows the direct connection between the Green functions and dynamical
correlation functions of creation and destruction operators. We refer the
reader to Section $5.2$ of Ref. [68] for a detailed discussion of correlation
functions in dissipative systems. Because of probability conservation, the
four functions are not independent, i.e., they satisfy the identity:

$G^{T}(x_{1},x_{2})+G^{\tilde{T}}(x_{1},x_{2})-G^{<}(x_{1},x_{2})-G^{>}(x_{1},x_{2})\,=\,0\,.$
(S22)

It is then customary to introduce the Keldysh rotation of the fields, which
automatically deletes the spurious degrees of freedom, namely, one correlation
function is identically vanishing after the rotation. Following the convention
by Larkin and Ovchinnikov [70], we define:

$\phi_{1}=\frac{\psi_{+}+\psi_{-}}{\sqrt{2}},\;\;\;\phi_{2}=\frac{\psi_{+}-\psi_{-}}{\sqrt{2}},\;\;\;\bar{\phi}_{1}=\frac{\bar{\psi}_{+}-\bar{\psi}_{-}}{\sqrt{2}},\;\;\;\bar{\phi_{2}}=\frac{\bar{\psi}_{+}+\bar{\psi}_{-}}{\sqrt{2}}.$
(S23)

Applying this linear transformation of the fields $\psi_{\pm}$ to the Keldysh
action in Eq. (S20) gives the action reported in Eqs. (5) and (6) of the main
text. The rotation also clarifies the role of the remaining Green’s functions

$\displaystyle iG^{R}(x_{1},x_{2})$
$\displaystyle=\langle\phi_{1}(x_{1})\bar{\phi}_{1}(x_{2})\rangle=\Theta(t-t^{\prime})\langle\\{\,c(x_{1}),c^{\dagger}(x_{2})\\}\rangle,$
(S24a) $\displaystyle iG^{A}(x_{1},x_{2})$
$\displaystyle=\langle\phi_{2}(x_{1})\bar{\phi}_{2}(x_{2})\rangle=-\Theta(t^{\prime}-t)\langle\\{\,c(x_{1}),c^{\dagger}(x_{2})\\}\rangle,$
(S24b) $\displaystyle iG^{K}(x_{1},x_{2})$
$\displaystyle=\langle\phi_{1}(x_{1})\bar{\phi}_{2}(x_{2})\rangle=\braket{[c(x_{1}),c^{\dagger}(x_{2})]}=-2\langle
c^{\dagger}(x_{2})c(x_{1})\rangle+\langle\\{\,c(x_{1}),c^{\dagger}(x_{2})\\}\rangle,$
(S24c)

with $\braket{\phi_{2}(x_{1})\bar{\phi}_{1}(x_{2})}=0$ being zero. In fact,
the retarded and the advanced Green’s functions $G^{R}$ (S24a) and $G^{A}$
(S24b), respectively, give information about the quasi-particle spectrum of
the system, while the Keldsyh Green’s function $G^{K}$ (S24c) gives
information on the statistical occupation of the quasi-particle spectrum. The
retarded and advanced Green functions are related by Hermitian adjoint
$G^{R}(x_{1},x_{2})=[G^{A}(x_{1},x_{2})]^{\dagger}$ (complex conjugate and
swap of the space-time variables $x_{1}$ and $x_{2}$), while $G^{K}$ is anti-
hermitian $[G^{K}(x_{1},x_{2})]^{\dagger}=-G^{K}(x_{1},x_{2})$. The three
Green functions (S24) can be organised in the following 2x2 Keldysh-space
matrix:

$\hat{G}(x_{1},x_{2})=\begin{pmatrix}G^{R}(x_{1},x_{2})&G^{K}(x_{1},x_{2})\\\
0&G^{A}(x_{1},x_{2})\end{pmatrix}.$ (S25)

It is eventually customary to parametrise the Keldysh Green’s function in
terms of the hermitian distribution function
$F(x_{1},x_{2})=[F(x_{1},x_{2})]^{\dagger}$, which is defined from the
relation

$G^{K}=G^{R}\circ F-F\circ G^{A}.$ (S26)

In the previous equation, and in Sec. III, the symbol $\circ$ denotes the
convolution product in the space-time variables. Whenever the symbol $\circ$
appears together with a propagator $\hat{G}$, as in Eq. (S25), contraction of
the Keldysh 2x2 indices is also implied. The distribution function $F$ is the
object of the kinetic Boltzmann equation, whose derivation we detail in the
next section.

##  III Self-energy and Boltzmann equation

The evolution equation for the distribution function $F$ follows from the
Dyson equation $[\hat{G}^{-1}_{0}-\hat{\Sigma}]\circ\hat{G}\,=\,\mathds{1}$.
Here $\hat{G}_{0}$ denotes the bare propagators of the quadratic action
$S_{\mathrm{diff}}$ in Eq. (5) of the main text (we do not report their
expressions as they are not needed in the derivation), while $\hat{G}$ the
propagators dressed by the interaction via the self-energy $\hat{\Sigma}$. The
latter has the same structure in Keldysh indices as the propagator $\hat{G}$
in Eq. (S25)

$\hat{\Sigma}(x_{1},x_{2})=\begin{pmatrix}\Sigma^{R}(x_{1},x_{2})&\Sigma^{K}(x_{1},x_{2})\\\
0&\Sigma^{A}(x_{1},x_{2})\end{pmatrix}.$ (S27)

Taking the Keldysh component of the Dyson equation and parametrizing $G^{K}$
in terms of $F$ as in Eq. (S26), one has

$\left[-i(\partial_{t_{1}}+\partial_{t_{2}})-D(\nabla_{\vec{x_{1}}}^{2}-\nabla_{\vec{x_{2}}}^{2})+(V(\vec{x_{1}})-V(\vec{x_{2}}))\right]F(x_{1},x_{2})=I[F]=\Sigma^{K}\circ\mathds{1}-(\,\Sigma^{R}\circ
F-F\circ\Sigma^{A}\,).$ (S28)

The latter is named kinetic equation, and the right-hand side $I[F]$ is called
the collision integral. The self-energy matrix $\hat{\Sigma}$ encodes the
effect of the interaction. Within this description, the computation of $F$ and
of correlation functions necessarily amounts to calculating $\hat{\Sigma}$.
The latter cannot be evaluated exactly and therefore one needs to develop an
approximation scheme. In this manuscript, we consider the so-called Born
perturbation theory approximation [61, 62, 60]. Within this approximation,
$\hat{\Sigma}$ is evaluated by perturbatively expanding the interaction
vertices around the quadratic part and thereby expressing the resulting
Feynman diagrams in terms of the bare Green functions $\hat{G}_{0}$.

As summarised in the main text, we consider, in particular, first-order
contributions in $\Gamma$ to $\hat{\Sigma}\approx\hat{\Sigma}_{(1)}$. Such
contributions are encoded into tadpole diagrams, where internal loops are
connected to the external legs via a four-fielded vertex. There are four
possible choices of the external legs, i.e., two possible fields for each of
the two legs. Hence, we can sum up all internal contributions, and organise
them in the 2x2 Keldysh-space matrix (S27). Besides, at first-order in
perturbation theory
$\hat{\Sigma}(x_{1},x_{2})\approx\hat{\Sigma}_{(1)}(x_{1},x_{2})=\hat{\Sigma}_{(1)}(x)$
as tadpole graphs are defined at one space-time point $x_{1}=x_{2}=x$ only.
Its components read:

$\hat{\Sigma}_{(1)}(x)\,=\,\begin{pmatrix}\Sigma^{R}_{(1)}(x)&\Sigma^{K}_{(1)}(x)\\\
0&\Sigma^{A}_{(1)}(x)\end{pmatrix}\,=\,\frac{\Gamma}{4}\,\begin{pmatrix}[f(G^{A}_{0})+f(G^{K}_{0})]&[2f(G^{K}_{0})+f(G^{A}_{0})-f(G^{R}_{0})]\\\
0&[f(G^{R}_{0})-f(G^{K}_{0})]\end{pmatrix}\,.$ (S29)

Since the interaction vertices of the theory contain gradients of the fields
$\phi_{1,2}$, the self-energy $\hat{\Sigma}$ is a differential operator acting
on Feynman diagrams’ external legs. In the previous equation, this is
expressed by the appearance of the operator $f(G_{0})$, which is defined as:

$f(G_{0})=[\,G_{0}(x,x)\overleftarrow{\nabla}\,]\cdot\overrightarrow{\nabla}-[\,\overrightarrow{\nabla}G_{0}(x,x)\cdot\overleftarrow{\nabla}\,]+\overleftarrow{\nabla}\cdot[\,\overrightarrow{\nabla}G_{0}(x,x)\,]-\overleftarrow{\nabla}\cdot
G_{0}(x,x)\overrightarrow{\nabla}.$ (S30)

In Eq. (S30), the over-script arrows indicate the direction of the gradient
operator. In particular, the notation
$(G_{0}(x_{1},x_{2})\overleftarrow{\nabla})$ is used to denote differentiation
of the Green’s function $G_{0}(x_{1},x_{2})$ with respect to their second
variable, i.e.,
$(G_{0}(x_{1},x_{2})\overleftarrow{\nabla})=(\vec{\nabla}_{x_{2}}G_{0}(x_{1},x_{2}))$.
Conversely, the notation $(\overrightarrow{\nabla}G_{0}(x_{1},x_{2}))$
indicates differentiation with respect to the first variable of $G_{0}$, i.e.,
$(\overrightarrow{\nabla}G_{0}(x_{1},x_{2}))=(\vec{\nabla}_{x_{1}}G_{0}(x_{1},x_{2}))$.
The differential operators in Eq. (S30) not enclosed into square brackets act
on the diagrams’ external legs, with $\overleftarrow{\nabla}$ differentiating
with respect to the second variable of the external leg and
$\overrightarrow{\nabla}$ with respect to the first one. From Eq. (S29), using
the Hermitian adjoint properties of $G_{0}^{R,A,K}$ recalled after Eq. (S24),
it follows that $\Sigma^{R}_{(1)}=[\Sigma^{A}_{(1)}]^{\dagger}$ and
$[\Sigma^{K}_{(1)}]^{\dagger}=-\Sigma^{K}_{(1)}$. The self-energy therefore
shares the same structure as $\hat{G}$, as anticipated before Eq. (S27) (this
is expected to hold at any order in perturbation theory in $\Gamma$). We now
derive a more useful expression for (S29) and (S30), which can be better used
into the kinetic equation (S28). Since the self-energy matrix elements
$\Sigma^{K}$, $\Sigma^{R}$ are convoluted with the identity and with $F$,
respectively, we now want to turn operator $f(G_{0})$ into a standard
differential operator which acts on functions in subsequent positions on the
right. We manage to do so by introducing a space-time variable $x_{2}$ which
is constrained to take the same value of $x_{1}$ via a delta
$\delta(x_{1}-x_{2})$ (due to the tadpole structure of the diagrams). With
this trick, we can write:

$\Sigma_{(1)}^{R}(x_{1},x_{2})=\frac{\Gamma}{4}\,\Big{[}\,\tilde{f}(G^{A}_{0})+\tilde{f}(G^{K}_{0})\,\Big{]},\quad\mbox{and}\quad\Sigma_{(1)}^{K}(x_{1},x_{2})=\frac{\Gamma}{4}\,\Big{[}\,2\,\tilde{f}(G^{K}_{0})+\tilde{f}(G^{A}_{0})-\tilde{f}(G^{R}_{0})\,\Big{]},$
(S31)

with

$\displaystyle\tilde{f}(G_{0})$
$\displaystyle=(\vec{\nabla}_{x_{1}}\delta(x_{1}-x_{2}))\cdot[\,G_{0}(x_{1},x_{2})\vec{\nabla}_{x_{2}}-(\vec{\nabla}_{x_{1}}G_{0}(x_{1},x_{2}))\,]+$
$\displaystyle+\delta(x_{1}-x_{2})\,[\,(\vec{\nabla}_{x_{1}}G_{0}(x_{1},x_{2}))\cdot\vec{\nabla}_{x_{2}}+(\vec{\nabla}_{x_{2}}G_{0}(x_{1},x_{2}))\cdot\vec{\nabla}_{x_{2}}-(\nabla^{2}_{x_{1}}G_{0}(x_{1},x_{2}))-(\vec{\nabla}_{x_{1}}\cdot\vec{\nabla}_{x_{2}}G_{0}(x_{1},x_{2}))\,].$
(S32)

Derivatives of the Dirac deltas are understood in the sense of distribution,
i.e., they have to be evaluated via integration by parts when inserted in the
convolution products of the kinetic equation. We can also derive an analogous
expression for $\Sigma^{A}$, which acts on the distribution function from the
right, c.f. Eq. (S28). where gradient operators act on the distribution
function from the right only:

$\Sigma^{A}_{(1)}(x_{1},x_{2})=\frac{\Gamma}{4}\,\Big{[}\,\tilde{\tilde{f}}(G^{R}_{0})-\tilde{\tilde{f}}(G^{K}_{0})\,\Big{]},$
(S33)

with

$\displaystyle\tilde{\tilde{f}}(G_{0})$
$\displaystyle=\delta(x_{1}-x_{2})\,[\,(\nabla^{2}_{x_{1}}G_{0}(x_{1},x_{2}))+(\nabla^{2}_{x_{2}}G_{0}(x_{1},x_{2}))\,]\,-\,\overleftarrow{\nabla}_{x_{1}}\delta(x_{1}-x_{2})\cdot[\,(\vec{\nabla}_{x_{1}}G_{0}(x_{1},x_{2}))+(\vec{\nabla}_{x_{2}}G_{0}(x_{1},x_{2}))\,]-$
$\displaystyle+\,[\,(\vec{\nabla}_{x_{2}}G_{0}(x_{1},x_{2}))-\overleftarrow{\nabla}_{x_{1}}G_{0}(x_{1},x_{2})\,]\cdot(\vec{\nabla}_{x_{2}}\delta(x_{1}-x_{2})).$
(S34)

In order to extract the large-scale universal limit of the kinetic equation
(S28) and explicitly evaluate the convolution products appearing in the
collision term, we resort to the Wigner transformation (WT). We first move to
the set of Wigner coordinates $x=(x_{1}+x_{2})/2$, $x^{\prime}=x_{1}-x_{2}$,
again using the 4-vector notation $x^{\prime}=(\vec{x}^{\prime},t^{\prime})$
and $x=(\vec{x},t)$. These set of variables clearly distinguishes between the
slow, macroscopic, evolution in $x$ and the fast dynamics in $x^{\prime}$. The
chain rule yields
$\vec{\nabla}_{x_{1}}\,=\,\vec{\nabla}_{x}/2+\vec{\nabla}_{x^{\prime}}$. The
Wigner transform $A(x,k)$ corresponds to the Fourier transform of
$A(x_{1},x_{2})$ with respect to the set of relative coordinates $x^{\prime}$,
which introduces the conjugated 4-momentum $k=(\vec{k},\epsilon)$:

$A(x,k)=\int
d^{d}x^{\prime}\,dt^{\prime}\,A(x+x^{\prime}/2,x-x^{\prime}/2)\,e^{-i\vec{k}\cdot\vec{x}^{\prime}+i\epsilon
t^{\prime}},\,\,\,A(x_{1},x_{2})=\int\,\frac{d^{d}k}{(2\pi)^{d}}\,\frac{d\epsilon}{2\pi}\,A\left(\frac{x_{1}+x_{2}}{2},k\right)\,e^{i\vec{k}\cdot(\vec{x_{1}}-\vec{x_{2}})-i\epsilon_{k}(t_{1}-t_{2})}.$
(S35)

The main advantage of this procedure is to turn convolution products of two
functions $A$ and $B$ in Eq. (S28) into Moyal products [77, 78, 79, 80, 81] of
their Wigner transforms

$C(x_{1},x_{2})=(A\circ B)(x_{1},x_{2})=\int
d^{d}x_{3}A(x_{1},x_{3})B(x_{3},x_{2})\,\xlongrightarrow{\text{WT}}\,C(x,p)=A(x,k)\,\mathrm{exp}\left\\{\frac{i}{2}[\overleftarrow{\partial}_{x}\overrightarrow{\partial}_{k}-\overleftarrow{\partial}_{k}\overrightarrow{\partial}_{x}]\right\\}\,B(x,k).$
(S36)

Crucially, in the reaction-limited regime $\hbar n\Gamma/D\ll 1$, time
variations due to reactions are slow. Space modulations, induced by trapping
potential $V$, can be also considered weak as long as $V$ varies on
macroscopic length scales. In this limit, one can assume that the dependency
of $G^{R,A,K}(x,x^{\prime})$ on the macroscopic variable $x$ is much slower
than the one on the relative one $x^{\prime}$, entailing that

$\vec{\nabla}_{x_{1}}G^{R,A,K}(x_{1},x_{2})=\frac{1}{2}\vec{\nabla}_{x}G^{R,A,K}(x,x^{\prime})+\vec{\nabla}_{x^{\prime}}G^{R,A,K}(x,x^{\prime})\simeq\vec{\nabla}_{x^{\prime}}G^{R,A,K}(x,x^{\prime})\xlongrightarrow{WT}i\vec{k}G^{R,A,K}(x,k).$
(S37)

Moreover, for slow space-time variations, one can simplify the Moyal product
in Eq. (S36) by expanding the complex exponential. This leads to a derivative
expansion, which can be truncated at the lowest order. From this expansion,
the left-hand-side of the kinetic equation (S28) becomes

$\left[-i(\partial_{t_{1}}+\partial_{t_{2}})-D(\nabla_{\vec{x_{1}}}^{2}-\nabla_{\vec{x_{2}}}^{2})+(V(\vec{x_{1}})-V(\vec{x_{2}}))\right]F(x_{1},x_{2})\xlongrightarrow{\text{WT}}\,i\,\Big{[}\,\partial_{t}+v_{g}(\vec{k})\cdot\vec{\nabla}_{x}-\frac{1}{\hbar}\,\vec{\nabla}_{x}V(\vec{x})\cdot\vec{\nabla}_{k}\,\Big{]}F(x,k).$
(S38)

We remark that Eq. (S38) is exact for harmonic potentials. For anharmonic
potentials, the first correction to (S38) is proportional to the cubic-space
derivative of the potential [81]. To deal with the WT of the collision
integral, we first write $G_{0}^{K}=F(G_{0}^{R}-G_{0}^{A})$ in Eqs.
(S31)-(S34). The previous equation follows from the WT of Eq. (S26) at lowest
order in derivatives. Here, albeit $G_{0}^{R,A,K}$ are bare Green functions,
$F$ does not obey the free/quadratic-action dynamics, but it must be
determined self-consistently from the kinetic equation. Truncating the Moyal
expansion (S36) of convolution products at lowest order in derivatives and
using Eq. (S37), we eventually obtain for the collision integral
[$q=(\vec{q},\omega)$]

$I[F]=-\frac{i\Gamma}{2}\int\frac{d^{d}q}{(2\pi)^{d}}\frac{d\omega}{2\pi}(\vec{k}-\vec{q})^{2}[1-F(x,k)-F(x,q)+F(x,k)F(x,q)]A_{0}(x,q),$
(S39)

where we introduced the bare spectral function
$A_{0}(x,q)=-2\mbox{Im}(G^{R}_{0}(x,q))=i[G^{R}_{0}(x,q)-G^{A}_{0}(x,q)]$. One
can also replace the bare green functions $G_{0}^{R,A}$ in (S39) with the
dressed ones $G^{R,A}$, and therefore $A_{0}(x,q)\to A(x,q)$, within the self-
consistent Born approximation. One then writes an equation for $G^{R,A}$ which
has to be solved self-consistently together with the equation for $F$. The
self-consistent Born approximation is necessary when the perturbative Born
series is divergent, as shown in Ref. [71]. In our case, however, the
perturbative expression (of order $\Gamma$) (S39) is finite and the self-
consistent Born approximation therefore only gives sub-leading corrections.
For this reason, we do not consider it here.

From Eq. (S39), one can see that the collision integral $I[F]$ is a purely
imaginary quantity: $F(y,k)$ is real since it is the WT of the hermitian
distribution function $F(x_{1},x_{2})$. This result implies that the self-
energy does not generate any renormalization of the quasi-particle energy
spectrum $\epsilon_{k}(\vec{x})=D\vec{k}^{2}+V(\vec{x})$, i.e., the kinetic
term in Eq. (S38) is left unchanged. On the contrary, the imaginary part of
the collision integral represents the finite particle lifetime
$\sim(n\Gamma)^{-1}$ introduced by the dissipative interactions. In the
reaction-limited regime of weak dissipation, this lifetime
$\sim(n\Gamma)^{-1}$ is much longer than the coherent time-scale set by
$(D/\hbar)^{-1}$ and therefore quasi-particles remain well defined during time
evolution. The spectral function $A(x,k)$ is accordingly sharply peaked in
$\epsilon$ around $\epsilon_{k}(\vec{x})$, with a width $\sim\Gamma$. In the
non-interacting limit $\Gamma=0$, one exactly has
$A_{0}(x,k)=2\pi\delta(\epsilon-\epsilon_{k}(\vec{x}))$. This allows us to
consider the on-shell distribution function
$\widetilde{F}(\vec{x},\vec{k},t)\equiv
F(\vec{x},\vec{k},t,\epsilon=\epsilon_{k}(\vec{x}))$, which is given by

$\widetilde{F}(\vec{x},\vec{k},t)=\int\frac{d\epsilon}{2\pi}\,F(x,k)A(x,k)=i\int\frac{d\epsilon}{2\pi}G^{K}(x,k)=iG^{K}(\vec{x},\vec{k},t,t)=1-2n(\vec{x},\vec{k},t).$
(S40)

In the first equality, we used the peaked structure of the spectral function
$A(x,k)$, in the second the expression of $G^{K}(x,k)$, in the third the
inverse WT definition in $\epsilon$ (S35), and in the fourth equality the WT
in space of Eq. (S24c) (evaluated at equal times). The evolution equation for
the one-body Wigner function $n(\vec{x},\vec{k},t)$ is thus equivalent to that
for $\widetilde{F}(\vec{x},\vec{k},t)$. The latter is obtained by multiplying
(S38) and (S39) by $A(x,k)$ and integrating in $\epsilon$. The integral in
$\omega$ in Eq. (S39) is similarly handled by exploiting the Dirac delta form
of $A_{0}(x,q)=2\pi\delta(\omega-\epsilon_{k}(\vec{x}))$. Applying the
substitution $\widetilde{F}(\vec{x},\vec{k},t)=1-2n(\vec{x},\vec{k},t)$, the
Boltzmann equation (8) of the main text eventually follows.

##  IV Effective exponent for the trap-release protocol

Figure S2:  Effective exponent for the trap-release quench. Effective exponent
$\xi_{\mathrm{eff}}$ for increasing values of $\Omega$ (from top to bottom),
for the decay of the rescaled particle number $\tilde{N}(\tilde{t})$ as a
function of the rescaled time interval $\tilde{t}\in[0,10]$. Particle number
and time are rescaled as explained in the main text and the values of $\Omega$
reported are identical to those of Fig. 3(c) of the main text. The plot shows
the existence of an approximate power-law decay
$\tilde{N}(\tilde{t})\sim\tilde{t}^{-\xi_{M}}$ at intermediate times, where
$\xi_{\mathrm{eff}}$ reaches a maximum $\xi_{M}$. For longer times, a non-
algebraic decay sets in and $\xi_{\mathrm{eff}}$ decreases with respect to
$\xi_{M}$. Such slow non power-law decay is present for every finite $\Omega$
value, and it takes place earlier in time as $\Omega$ is increased. All
effective exponents have been computed with $b=1.1$ in Eq. (S41).

We give here additional details concerning the decay of the particle number in
the trap-release protocol discussed in Fig. 3(c)-(d) of the main text. As we
presented in the main text, the behavior of the re-scaled particle number
$\tilde{N}(\tilde{t})=N(\tilde{t})/N_{0}$ at short rescaled times $\tilde{t}$
approaches the power-law $\tilde{N}(\tilde{t})\sim\tilde{t}^{-\xi}$, where
$\xi$ is an exponent depending on parameter $\Omega$. In order to quantify
$\xi$ and to see how the asymptotic decay changes as $\Omega$ is tuned, we
compute the effective decay exponent $\xi_{\mathrm{eff}}$ [2]

$\xi_{\mathrm{eff}}=-\,\frac{\log\,\Big{[}\tilde{N}(b\tilde{t})/\tilde{N}(\tilde{t})\Big{]}}{\log(b)},$
(S41)

with $b$ an arbitrary adimensional positive parameter. One can check that if
$\tilde{N}(\tilde{t})$ asymptotically approaches a power law, i.e.,
$\tilde{N}(\tilde{t})\sim\tilde{t}^{-\xi}$ at long times, then
$\lim_{t\rightarrow\infty}\xi_{\mathrm{eff}}=\xi$, $\forall b$. In order to
evaluate the effective decay exponent in a neighborhood of a given time
coordinate, we choose $b=1.1$.

The plot in Fig. S2 shows the effective exponent
$\xi_{\mathrm{exp}}(\tilde{t})$ in the interval $\tilde{t}\in[0,10]$ for
various values of $\Omega$ (the same used in Fig. 3(c) of the main text). For
very small $\Omega=0.01$ value, one can identify a power-law behavior with
effective exponent $\xi_{\mathrm{eff}}(\tilde{t})=\xi\approx 1/2$ throughout
the considered domain, since $\xi_{\mathrm{eff}}$ converges to that value.
Conversely, by increasing $\Omega$, an approximate algebraic decay
$\tilde{N}(\tilde{t})\sim\tilde{t}^{-\xi_{M}}$ is observed, where the
effective exponent reaches a maximum value $\xi_{M}$ and remains constant
within an intermediate rescaled time window. For instance, the selected values
$\Omega=0.05$, $\Omega=0.1$ are characterized by $\xi_{M}\approx 0.45$,
$\xi_{M}\approx 0.4$ within the rescaled time interval $\tilde{t}\sim[3,6]$,
$\tilde{t}\sim[2,4]$, respectively. For longer times the effective exponent
$\xi_{\mathrm{eff}}$ drifts away from the aforementioned values and slowly
decreases. This signals the onset of a non-algebraic (since
$\xi_{\mathrm{eff}}$ is not constant) slow decay. As $\Omega$ is increased,
the approximate decay exponent $\xi_{M}$ therefore decreases, the time window
where it is observed is anticipated and shrinks. For $\Omega=0.01$, the non-
algebraic behavior takes place for longer times than those shown in the plot
($\tilde{t}<10$). Conversely, for large values $\Omega=1$, the intermediate
power law decay is completely wiped out, no exponent can be identified and the
non-algebraic decay immediately sets in.