font=small  ††thanks<EMAIL_ADDRESS>
# Quantum reaction-limited reaction-diffusion dynamics of annihilation
processes

Gabriele Perfetto  Institut für Theoretische Physik, Universität Tübingen, Auf
der Morgenstelle 14, 72076 Tübingen, Germany  Federico Carollo  Institut für
Theoretische Physik, Universität Tübingen, Auf der Morgenstelle 14, 72076
Tübingen, Germany  Juan P. Garrahan  School of Physics, Astronomy, University
of Nottingham, Nottingham, NG7 2RD, UK.  Centre for the Mathematics,
Theoretical Physics of Quantum Non-Equilibrium Systems, University of
Nottingham, Nottingham, NG7 2RD, UK  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

###### Abstract

We investigate the quantum reaction-diffusion dynamics of fermionic particles
which coherently hop in a one-dimensional lattice and undergo annihilation
reactions. The latter are modelled as dissipative processes which involve
losses of pairs $2A\to\emptyset$, triplets $3A\to\emptyset$ and quadruplets
$4A\to\emptyset$ of neighbouring particles. When considering classical
particles, the corresponding decay of their density in time follows an
asymptotic power-law behavior. The associated exponent in one dimension is
different from the mean-field prediction whenever diffusive mixing is not too
strong and spatial correlations are relevant. This specifically applies to
$2A\to\emptyset$, while the mean-field prediction just acquires a log
correction for $3A\to\emptyset$ and is exact for $4A\to\emptyset$. A mean-
field approach is also valid, for all the three processes, when the diffusive
mixing is strong, i.e., in the so-called reaction-limited regime. Here, we
show that the picture is quite different for quantum systems. We consider the
quantum reaction-limited regime and we show that for all the three processes
power-law behavior beyond mean-field is present as a consequence of quantum
coherences, which are not related to space dimensionality. In the case of
$3A\to\emptyset$, we further show that power-law behavior only appears within
an intermediate time window, while at long times the density decay is not
power-law. Our results show the markedly different origin of the richer
universal behavior in quantum reaction-diffusion dynamics as compared to its
classical formulation.

##  I Introduction

The investigation and classification of universal behavior in non-equilibrium
many-body systems is a timely and challenging research area. Far from
equilibrium the Boltzmann-Gibbs measure does not correctly describe the system
and therefore the emergence of universal behavior cannot be clearly pinpointed
as in equilibrium-critical systems Wilson and Kogut (1974); Kogut (1979).
Within this perspective, reaction-diffusion (RD) systems, where classical
particles hop on a lattice (diffusion in the continuum limit) and react upon
meeting, are genuine non-equilibrium systems where universal quantities can be
identified and characterized Henkel  _et al._ (2008); Hinrichsen (2000);
Vladimir (1997); Tauber ; Täuber  _et al._ (2005); Täuber (2014). At long
times the density of particles decays as a power-law, which can be obtained,
for instance, from a mean-field treatment whenever diffusive mixing is
effective. This is valid when the diffusive hopping is strong, i.e., in the
reaction-limited regime (sometimes also called “well stirred mixture” regime)
Hinrichsen (2000); Vladimir (1997); Kang and Redner (1985); Privman and
Grynberg (1992), or in the opposite diffusion-limited regime of weak hopping
in dimensions larger than the upper critical dimension $d_{c}$. This dimension
can be identified both via exact lattice calculations Toussaint and Wilczek
(1983); Spouge (1988); Privman (1994); Torney and McConnell (1983); Kang and
Redner (1984a, b); Kang  _et al._ (1984); Rácz (1985) and/or renormalization
group methods Tauber ; Täuber  _et al._ (2005); Täuber (2014); Doi (1976);
Peliti (1985, 1986); Mattis and Glasser (1998). In the case of binary
annihilation, $2A\to\emptyset$, it is $d_{c}=2$, so that in one dimension
spatial fluctuations are relevant and they are responsible for universal
power-law decay that is different from mean-field Toussaint and Wilczek
(1983); Spouge (1988); Privman (1994); Torney and McConnell (1983); Kang and
Redner (1984a, b); Tauber ; Täuber  _et al._ (2005); Täuber (2014); Peliti
(1985, 1986). The upper critical dimension depends, however, on the number of
particles involved in the reaction. For triplet annihilation $3A\to\emptyset$
it is $d_{c}=1$, and thus in one dimension spatial correlations are marginal
and the mean-field power-law decay only acquires a logarithmic correction
Tauber ; Täuber  _et al._ (2005); Täuber (2014). For quadruplet annihilation,
$4A\to\emptyset$, or higher, $d_{c}<1$, so that spatial fluctuations are
irrelevant and the mean-field decay is observed in all dimensions both for the
diffusion-limited and reaction-limited regimes.

Quantum RD systems, where particles move via coherent hopping while subject to
dissipative reactions, represent a class of dynamical processes which are
currently under intense investigation. Annihilation reactions have, indeed, a
direct connection to cold-atomic experiments with two Syassen  _et al._
(2008); Traverso  _et al._ (2009); Yamaguchi  _et al._ (2008); Yan  _et al._
(2013); Zhu  _et al._ (2014); Sponselee  _et al._ (2018); Kinoshita  _et al._
(2005); Tomita  _et al._ (2017), three Bouchoule and Schemmer (2020); Söding
_et al._ (1999); Tolra  _et al._ (2004) and four Ferlaino  _et al._ (2009);
Gurian  _et al._ (2012) body losses, which recently received significant
attention also at a theoretical level Diehl  _et al._ (2008); Bouchoule and
Dubail (2021, 2022); García-Ripoll  _et al._ (2009); Ates  _et al._ (2012);
Everest  _et al._ (2014); Bouchoule  _et al._ (2021); Rossini  _et al._
(2021); Rosso  _et al._ (2022, 2023); Bouchoule  _et al._ (2020). Analyzing
their universal properties is particularly challenging since these
investigations require the consideration of dissipative quantum dynamics of
large systems at long times. Such dynamics, and that of similar dissipative
models, such as those with kinetic constraints Lesanovsky and Garrahan (2013);
Olmos  _et al._ (2014); Everest  _et al._ (2016); Marcuzzi  _et al._ (2016);
Buchhold  _et al._ (2017); Gutiérrez  _et al._ (2017); Roscher  _et al._
(2018); Wintermantel  _et al._ (2020); Helmrich  _et al._ (2020); Nigmatullin
_et al._ (2021); Kazemi and Weimer (2021); Carollo  _et al._ (2022), serve as
benchmark problems for numerics van Horssen and Garrahan (2015); Carollo  _et
al._ (2019); Gillman  _et al._ (2019, 2020); Jo  _et al._ (2021) and for
quantum simulators Zeiher  _et al._ (2016); Kim  _et al._ (2018); Ebadi  _et
al._ (2021); Jo and Kim (2022).

In recent work, Ref. Perfetto  _et al._ (2022), we presented an analytical
study of quantum RD systems in their reaction-limited regime. We showed there
that power-law behavior distinct from mean-field can occur for binary
reactions, in contrast to the classical case. In this paper we address the
same question for three and four-body reactions, where spatial fluctuations
are expected to be irrelevant already in one dimension at the classical level.

We do this by considering the quantum RD dynamics of fermionic chains, where
particles coherently hop on a lattice and are subject to dissipative
annihilation reactions, as sketched in Fig. 1. We analytically study this
dynamics in the thermodynamic limit in the reaction-limited regime by using
the time-dependent generalized Gibbs ensemble approach (tGGE) Lange  _et al._
(2018); Mallayya  _et al._ (2019); Lange  _et al._ (2017); Lenarčič  _et al._
(2018). In order to introduce the problem, we first review the case of binary
annihilation, $2A\to\emptyset$, which had been previously studied with the
tGGE method in Refs. Rossini  _et al._ (2021); Rosso  _et al._ (2022, 2023).
In this case, the exponent of the power-law decay deviates from the mean-field
one as a consequence of quantum coherences in the initial state Perfetto  _et
al._ (2022). For the case of three-body annihilation, $3A\to\emptyset$, and
four-body annihilation, $4A\to\emptyset$ annihilation, we find that the effect
of quantum coherences is even more dramatic. For $3A\to\emptyset$, we find
that power-law behavior (also distinct from mean-field) is only transient,
with an asymptotic decay of the density which is not a power-law. This might
be related to the fact that the associated classical (and quantum) diffusion-
limited dynamics acquires a logarithmic correction. Our analysis, however,
shows that the quantum decay is richer and the correction to the algebraic
time scaling is slower than just a logarithm. For $4A\to\emptyset$, we find a
non mean-field power-law decay for initial states that possess quantum
coherences, similarly to $2A\to\emptyset$ case Perfetto  _et al._ (2022). For
$4A\to\emptyset$, this result is even more surprising given that spatial
fluctuations are irrelevant classically in one dimension. Our results for
$3A\to\emptyset$ and $4A\to\emptyset$ show that the universal quantum RD
behavior is not solely determined by the relevance of spatial fluctuations as
in the classical case, but also depends on quantum effects which are present
regardless of space dimensionality (such as coherences in the initial state).

The remainder of the manuscript is organized as follows. In Sec. II, we
formulate quantum RD dynamics with annihilation reactions by means of a
quantum master equation approach. In Sec. III, we first briefly discuss known
results about classical RD dynamics. We then move to the reaction-limited
regime of quantum RD dynamics and the associated analysis using the tGGE
method. This method is used in Sec. IV to address the annihilation reactions
$2A\to\emptyset$, $3A\to\emptyset$ and $4A\to\emptyset$. In Sec. V, we report
our conclusions. The appendices A and B contain technical aspects concerning
the tGGE calculations for three and four-body annihilation, respectively.

##  II The system

We consider a one-dimensional lattice with $L$ sites and periodic boundary
conditions. Each site $j$ can be either occupied by a fermion
$n_{j}\ket{\cdots\bullet_{j}\cdots}=\ket{\cdots\bullet_{j}\cdots}$ or be empty
$n_{j}\ket{\cdots\circ_{j}\cdots}=0$. Here, $n_{j}=c_{j}^{\dagger}c_{j}$ is
the number operator and the operators $c_{j},c_{j}^{\dagger}$ obey the
fermionic anticommutation relations
$\\{c_{j},c_{j^{\prime}}^{\dagger}\\}=\delta_{j,j^{\prime}}$. Particles
occupying adjacent sites may be lost into the environment through an
annihilation reaction. Therefore, the ensuing dynamics is not unitary and we
assume it to be governed by the quantum master equation Gorini  _et al._
(1976); Lindblad (1976); Breuer and Petruccione (2002) ($\hbar=1$ henceforth)

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

Here, $\rho$ is the density matrix, $H$ is the quantum Hamiltonian, and the
irreversible reaction processes are encoded in the dissipator $\mathcal{D}$,
which we take in Lindblad form Gorini  _et al._ (1976); Lindblad (1976);
Breuer and Petruccione (2002)

$\mathcal{D}[\rho]=\sum_{j}\left[L_{j}^{\nu}\rho{L_{j}^{\nu}}^{\dagger}-\frac{1}{2}\left\\{{L_{j}^{\nu}}^{\dagger}L_{j}^{\nu},\rho\right\\}\right].$
(2)

The $L_{j}^{\nu}$ are local jump operators. We consider here three different
types of reaction processes, namely, binary $2A\to\emptyset$ ($\nu=2\alpha$),
triplet $3A\to\emptyset$ ($\nu=3\alpha$) and quadruplet $4A\to\emptyset$
annihilation ($\nu=4\alpha$). These processes are pictorially represented in
Fig. 1. In particular, for binary annihilation, $2A\to\emptyset$, of a pair of
neighboring particles at rate $\Gamma_{2\alpha}$, we have

$L_{j}^{2\alpha}=\sqrt{\Gamma_{2\alpha}}\,c_{j}c_{j+1}.$ (3)

This process corresponds to two-body losses, which can be implemented in cold
atomic gases, e.g., via inelastic scattering interactions Syassen  _et al._
(2008); Traverso  _et al._ (2009); Yamaguchi  _et al._ (2008); Yan  _et al._
(2013); Zhu  _et al._ (2014); Sponselee  _et al._ (2018), or photoassociation
into excited compounds Kinoshita  _et al._ (2005); Tomita  _et al._ (2017).

We also consider annihilation of triplets of neighbouring particles
$3A\to\emptyset$ at rate $\Gamma_{3\alpha}$

$L_{j}^{3\alpha}=\sqrt{\Gamma_{3\alpha}}\,c_{j}c_{j+1}c_{j+2}.$ (4)

This triplet annihilation is also present in cold atomic gases, where it is
caused by recombination of atoms into molecules Bouchoule and Schemmer (2020);
Söding  _et al._ (1999); Tolra  _et al._ (2004). Finally, we consider
annihilation of quadruplets of neighbouring particles $4A\to\emptyset$ at rate
$\Gamma_{4\alpha}$, which is described by the jump operator

$L_{j}^{4\alpha}=\sqrt{\Gamma_{4\alpha}}\,c_{j}c_{j+1}c_{j+2}c_{j+3}.$ (5)

Such four-body losses have also been experimentally detected in cold atomic
gases Ferlaino  _et al._ (2009); Gurian  _et al._ (2012).

For all the three losses mechanisms in Eqs. (3)-(5), the Lindblad master
equation (1)-(2) [with the Hamiltonian introduced below in Eq. (6)] we employ
should be understood as an effective description of the dynamics. The
derivation of the Lindblad equation from a microscopic weak system-bath
coupling procedure (see e.g., Ref. Breuer and Petruccione (2002)) for the jump
operators (3)-(5) is an open problem. This derivation has been so far only
pursued in Ref. Bouchoule  _et al._ (2021) for one-body decay
$L_{j}^{\alpha}=\sqrt{\Gamma_{\alpha}}c_{j}$.

We also mention that the form of the jump operators (3)-(5) is dictated by the
fermionic statistics, which forces annihilation reactions to occur between
neighbouring particles since double (or higher) occupancy of the same lattice
site is not possible. This is different from the case of bosonic systems,
where annihilation reactions are typically defined between particles residing
on the same site. These processes are therefore modelled Bouchoule and Dubail
(2021, 2022); García-Ripoll  _et al._ (2009); Rossini  _et al._ (2021); Rosso
_et al._ (2022); Bouchoule  _et al._ (2020) for bosons through jump operators
of the form $L_{j}^{k\alpha}=\sqrt{\Gamma_{k\alpha}}\,b_{j}^{k}$, with
$k=2,3,4$ and $b_{j}$ a bosonic destruction operator.

Figure 1: Sketch of quantum RD processes. We consider a fermionic chain where
each lattice site is either empty $n_{j}\ket{\cdots\circ_{j}\cdots}=0$ or
filled $n_{j}\ket{\cdots\bullet_{j}\cdots}=\ket{\cdots\bullet_{j}\cdots}$ with
a fermion, with $n_{j}$ the number operator at site $j$. Reactions are
irreversible and they are modelled through the jump operators of the Lindblad
dynamics in Eqs. (1) and (2). In particular, we consider binary
$2A\to\emptyset$ (3) (rate $\Gamma_{2\alpha}$), triplet $3A\to\emptyset$ (4)
(rate $\Gamma_{3\alpha}$) and quadruplet $4A\to\emptyset$ (5) (rate
$\Gamma_{4\alpha}$) annihilation of neighbouring particles. Coherent-
Hamiltonian hopping (6) between neighbouring sites at rate $\Omega$ replaces
classical diffusion. The particle density $\braket{n(t)}$ decays at long times
towards the vacuum $\ket{\circ\circ\dots\circ\circ}$ in a universal way. This
decay is different from the one obtained in the classical setup as a
consequence of quantum coherences in the initial state.

For the coherent dynamics in Eq. (1), we take quantum hopping between adjacent
lattice sites at a rate $\Omega$ (cf. Fig. 1):

$H=-\Omega\sum_{j=1}^{L}(c_{j}^{\dagger}c_{j+1}+c_{j+1}^{\dagger}c_{j})\,.$
(6)

This Hamiltonian provides a “natural” quantum generalization of the diffusive
motion considered in classical RD dynamics Henkel  _et al._ (2008); Hinrichsen
(2000); Vladimir (1997); Tauber ; Täuber  _et al._ (2005); Täuber (2014). The
Hamiltonian (6) is quadratic and it can be exactly diagonalized by Fourier
transform. The whole Lindblad dynamics (1) and (2), however, is not quadratic
for the jump operators in Eqs. (3)-(5) and therefore it cannot be solved
exactly.

Quantum hopping introduces coherence in the dynamics, which is generically
expected to affect emergent universal dynamical behavior. For the purpose of
this work we quantify the latter through the long-time asymptotic decay of the
particle density $\braket{n(t)}=\braket{N(t)}/L$, where
$N(t)=\sum_{j}n_{j}(t)$ the total particle number at time $t$. The idea is as
follows: while $N$ is conserved by the Hamiltonian, $[H,N]=0$, it is not
conserved by the Lindblad dynamics due to the reactions (3)-(5). They deplete
the system, taking it towards a trivial vacuum stationary state. The non-
trivial universal behavior lies in the way this stationary state is approached
at long times, e.g., through a density decaying via a power-law. The
functional form of this asymptotic behavior is controlled by the relative
strength, $\Gamma/\Omega$, of incoherent dissipation with respect to coherent
Hamiltonian hopping, as we explain in the next Section.

##  III Diffusion and reaction-limited dynamics

The RD dynamics is characterized by two fundamental time scales. On the one
hand, the reaction time $\sim\Gamma^{-1}$, gives the time needed for two
nearby particles to react. On the other hand, the diffusion [hopping in the
quantum case of Eq. (6)] time $\sim\Omega^{-1}$ yields the characteristic time
needed by two separated particles to get close to each other. The dynamics is
diffusion-limited when the ratio $\Gamma/\Omega$ is at least
$\Gamma/\Omega\sim 1$, while it is reaction-limited when $\Gamma/\Omega\ll 1$.
The universal long-time decay of the particle density $\braket{n(t)}$ is
markedly different in the two regimes. In Subsec. III.1 and III.2, we briefly
recall some previous results concerning the RD dynamics in the diffusion and
reaction-limited regimes for classical and quantum systems, respectively.

###  III.1 Classical RD dynamics

For classical RD systems, the reaction-limited regime is described by mean-
field Hinrichsen (2000); Vladimir (1997); Kang and Redner (1985); Privman and
Grynberg (1992). This is due to the fast diffusion mixing, which renders the
particle density homogeneous in space. Spatial fluctuations in particle
concentration are rapidly smoothed out and reactions can therefore take place
everywhere with the same probability. The ensuing mean-field description, for
the case of annihilation processes of $k$ particles $kA\to\emptyset$, is given
by the law of mass action equation

$\frac{\mbox{d}\braket{n}}{\mbox{d}t}=-k\,\Gamma_{k\alpha}\braket{n}^{k},\,\,\braket{n}\sim(\Gamma_{k\alpha}t)^{-1/(k-1)}.$
(7)

In this limit the density depends therefore only on the rescaled time
$\tau=\Gamma_{k\alpha}t$.

In the diffusion-limited regime, local density fluctuations are relevant for
the dynamics. For this reason, the exponent of the density decay generically
deviates from the mean-field prediction (7) Toussaint and Wilczek (1983);
Spouge (1988); Privman (1994); Tauber ; Täuber  _et al._ (2005); Täuber
(2014); Doi (1976); Peliti (1985, 1986). This non mean-field universal
behavior simply comes from the fact that at long times the dynamics relies on
few remaining particles separated by large distances. In the classical case,
the decay is, consequently, controlled by the large-distance properties of a
random-walk, i.e., by the probability of far apart particles to meet. This is
a universal quantity which depends only on the diffusion constant $\Omega$, a
macroscopic quantity, and on the space dimensionality $d$, which sets the
dimensionality of the underlying diffusion process.

The random walk is recurrent in $d\leq 2$ Feller (1968), which implies that in
continuum space two particles can surely meet. This identifies the upper-
critical dimension $d_{c}=2$ for binary annihilation $2A\to\emptyset$. The
associated decay in $d=1<d_{c}$ is

$\braket{n(t)}\sim(\Omega t)^{-1/2},\,\,2A\to\emptyset,\,\,\,d=1.$ (8)

For $d>2$, diffusive mixing is not effective, as particles far apart in the
continuum do not meet, universality is lost and the mean-field description (7)
is recovered.

For reactions $kA\to\emptyset$ involving $k>2$ particles, the upper critical
dimension is $d_{c}=2/(k-1)$ Tauber ; Täuber  _et al._ (2005); Täuber (2014),
as space dimensionality must be further lowered in order to constrain a larger
number of particles to meet. For triplet annihilation $3A\to\emptyset$,
$d_{c}=1$. In $d=d_{c}=1$, the mean-field power-law (7) is valid up to a
logarithmic correction

$\braket{n(t)}\sim\left(\frac{\ln\Omega t}{\Omega
t}\right)^{1/2},\,\,3A\to\emptyset,\,\,\,d=1.$ (9)

Both in Eq. (8) and (9), the density depends on the rescaled time $\tau=\Omega
t$, which differs to the scaling valid in the reaction-limited regime (7).

For quadruplet annihilation $4A\to\emptyset$, instead, $d_{c}<1$ and the mean-
field description applies also in the one-dimensional $d=1>d_{c}$ diffusion-
limited dynamics

$\braket{n(t)}\sim(\Gamma_{4\alpha}t)^{-1/3},\,\,4A\to\emptyset,\,\,\,d=1.$
(10)

The mean-field decay is similarly valid also for five-body, and higher-order,
annihilation.

We therefore see that in classical RD systems universal non mean-field
behavior can only arise in the diffusion-limited regime as a consequence of
spatial fluctuations in low dimensions. In the quantum RD dynamics, this is
not the case: universal behavior that departs from mean-field is not only due
to spatial fluctuations, but can also be present in the the reaction-limited
regime where spatial fluctuations are smoothed out Perfetto  _et al._ (2022).
We will explore this in the next section for general multi-body annihilation
reactions.

###  III.2 Quantum RD dynamics

Currently, very little is known about quantum RD dynamics. The Eqs. (1)-(2)
cannot be analytically solved since the Lindbladian is not quadratic due to
the structure of the reaction jump operators (3)-(5). Likewise, numerical
simulations are hard due to the exponential scaling of the Hilbert space
dimension with the system size: in contrast to the classical case where one
can generate trajectories of configurations, a quantum trajectory unfolding of
(1) requires the propagation of whole state, making large scale numerics
unfeasible.

The diffusion-limited regime, with $\Gamma/\Omega=1$, for binary annihilation
$2A\to\emptyset$ [cf. Eqs. (3) and (6)] has been addressed numerically via
exact diagonalization up to $L=22$ sites in Ref. van Horssen and Garrahan
(2015). The system has been therein initialized in product state
$\ket{\bullet\bullet\dots\bullet\bullet}$ with unitary filling. The density is
found to decay algebraically in time $\braket{n(t)}\sim t^{-b}$, with
$1/2<b<1$. The decay exponent is different from the mean-field prediction,
which for quantum systems is analogous to the classical one (7), showing that
the dynamics is controlled in one-dimension by spatial fluctuations. Quantum
effects seem in this case to affect the decay exponent by making it presumably
larger than the corresponding classical value — $1/2$ — in Eq. (8), though
extrapolation of the finite-size results to the thermodynamic limit is not
trivial.

The reaction-limited regime $\Gamma/\Omega\ll 1$ has been only recently
studied in Ref. Perfetto  _et al._ (2022). Using an analytical approach, this
study could show that for binary reactions, such as $2A\to\emptyset$, the
quantum reaction-limited dynamics is not always described by the mean-field
approximation (7). This turns out to be the case whenever quantum coherences
in the initial state are present. The analysis of Ref. Perfetto  _et al._
(2022) is based on the time-dependent generalized Gibbs ensemble method (tGGE)
introduced in Refs. Lange  _et al._ (2018); Mallayya  _et al._ (2019); Lange
_et al._ (2017); Lenarčič  _et al._ (2018). In the following, we briefly
recall the main aspects of the tGGE, as this method will be employed in the
next sections.

The tGGE approach is based on a separation of time scales between reactions
and hopping dynamics, which is possible in the limit $\Gamma/\Omega\ll 1$.
Here, reactions are slow and long time intervals, on average, elapse between
consecutive reaction events. The much faster hopping dynamics (6) can thus be
integrated out by considering the state $\rho(t)$ in between consecutive
reactions as being relaxed with respect to the Hamiltonian $H$:
$[H,\rho(t)]=0$. The time dependence of the state $\rho(t)$ accounts for the
remaining slow evolution, taking place on the time scale $\Gamma^{-1}$, due to
the reactions. The tGGE method then moves forward by making an ansatz for the
state $\rho(t)=\rho_{\mathrm{GGE}}(t)$ of the form of a generalized Gibbs
ensemble, see, e.g., the reviews Essler and Fagotti (2016); Vidmar and Rigol
(2016). In the case of the Hamiltonian (6), the GGE takes the form

$\rho_{\mathrm{GGE}}(t)=\frac{1}{\mathcal{Z}(t)}\mbox{exp}\left(-\sum_{k}\lambda_{k}(t)\hat{n}_{k}\right),$
(11)

where $\mathcal{Z}(t)=\prod_{k}[1+e^{-\lambda_{k}(t)}]$. In the previous
equation, $\lambda_{k}(t)$ are dubbed Lagrange mulipliers (or generalized
inverse temperatures), $k\in(-\pi,\pi)$ is the quasi-momentum and
$\hat{n}_{k}=\hat{c}_{k}^{\dagger}\hat{c}_{k}$ is the number operator in
Fourier space, with $\hat{c}_{k}$ ($\hat{c}_{k}^{\dagger}$) fermionic
desctruction (creation) operators (see Appendix A). We note that the GGE state
(11) gives direct access to the RD dynamics in the thermodynamic limit, as it
describes the local relaxation of the system, i.e., the expectation
$\braket{\dots}_{\mathrm{GGE}}(t)$ provides the exact average behavior of
local observables in the thermodynamic limit. The state (11) is Gaussian and
diagonal in momentum space. Its dynamics is therefore fully characterized by
the momentum occupation functions
$\braket{\hat{n}_{q}}_{\mathrm{GGE}}(t)=C_{q}(t)=1/(\mbox{exp}(\lambda_{q}(t))+1)$,
which obey the equations Rossini  _et al._ (2021); Rosso  _et al._ (2022);
Bouchoule  _et al._ (2020); Rosso  _et al._ (2023); Perfetto  _et al._ (2022)

$\frac{\mbox{d}C_{q}(t)}{\mbox{d}t}=\sum_{j}\braket{{L_{j}^{\nu}}^{\dagger}[\hat{n}_{q},L_{j}^{\nu}]}_{\mathrm{GGE}}(t),\,\quad\forall
q.$ (12)

The structure of the previous equation clearly shows that $C_{q}(\tau)$ is a
function of the rescaled time $\tau=\Gamma_{\nu}t$ according to the reaction
rate, as in Eq. (7) for the classical reaction-limited regime. This aspect
further shows that the tGGE method is naturally apt to describe the reaction-
limited dynamics.

##  IV Quantum reaction-limited annihilation dynamics

In this Section, we present our results for the quantum annihilation process
in the reaction-limited regime. In order to set the stage, we first discuss in
Subsec. IV.1 the binary annihilation $2A\to\emptyset$ case of Eq. (3). We then
consider in Subsec. IV.2, triplet annihilation $3A\to\emptyset$ in Eq. (4). In
Subsec. IV.3, quadruplet annihilation $4A\to\emptyset$ of Eq. (5) is
eventually discussed.

In all the three cases, we solve Eq. (12) considering two different classes of
initial states. First, the Fermi sea (FS), i.e., the ground state of the
Hamiltonian (6) with an initial filling $n_{0}$, which is uniquely identified
by the occupation functions

$C_{q}(t=0)=\left\\{\begin{array}[]{lr}1\;\;\;\;\mbox{if}\,\,q\in[-\pi
n_{0},\pi n_{0}],\\\ 0\;\;\;\;\mbox{otherwise}.\end{array}\right.$ (13)

This initial state displays quantum coherences in real space, i.e., the
associated density matrix possesses off-diagonal elements in the Fock-space
basis $\prod_{j\in\Lambda}c_{j}^{\dagger}\ket{\circ\circ\circ\dots\circ}$,
with $\Lambda$ an arbitrary set of lattice sites. In the case $n_{0}=1$,
$\Lambda=\\{1,2,\dots L\\}$ and the state corresponds to the simple product
state $\ket{\bullet\bullet\bullet\dots\bullet}$, where every lattice site is
filled. Thus, as long as $n_{0}\neq 1$ the FS initial state is, instead,
quantum coherent.

Second, we consider also initial incoherent states of the form
$\rho_{0}=\mbox{exp}(-\lambda N)/\mathcal{Z}_{0}$. This initial state,
differently from the FS (13), is diagonal in the (classical) basis introduced
above and it is associated with a momentum-independent occupation function
equal to the initial filling $n_{0}$:

$C_{q}(0)=n_{0}.$ (14)

In both the cases (13) and (14), the density of particles
$\braket{n(\tau)}_{\mathrm{GGE}}$ is computed from the occupation function
$C_{q}(\tau)$ as

$\braket{n(\tau)}_{\mathrm{GGE}}=\frac{1}{L}\sum_{q}C_{q}(\tau).$ (15)

The decay exponent of the particle density as a function of time is quantified
by computing the effective exponent $\delta_{\mathrm{eff}}(\tau)$, which is
defined as Hinrichsen (2000)

$\delta_{\mathrm{eff}}(\tau)=-\frac{\log\left(\braket{n(b\tau)}_{\mathrm{GGE}}/\braket{n(\tau)}_{\mathrm{GGE}}\right)}{\log
b},$ (16)

with $b$ a scaling parameter. In the case of an asymptotic in time power-law
behavior $\braket{n(\tau)}_{\mathrm{GGE}}\sim a\tau^{-\delta}$, the effective
exponent $\delta_{\mathrm{eff}}(\tau)$ converges at long times to the exponent
$\delta$ of the power-law. In all the calculations of
$\delta_{\mathrm{eff}}(\tau)$ in this manuscript, we use $b=2$.

###  IV.1 Binary annihilation

The application of Eq. (12) to binary annihilation $2A\to\emptyset$ (3) leads
to the evolution equation Perfetto  _et al._ (2022)

$\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}=-\frac{1}{L}\sum_{k}f_{2}(k,q)C_{k}(\tau)C_{q}(\tau),$
(17)

with $\tau=\Gamma_{2\alpha}t$ and the function $f_{2}(k,q)$ given by

$f_{2}(k,q)=2(1-\cos(k-q)).$ (18)

Figure 2: Quantum binary annihilation dynamics in the reaction limited regime.
(a) Log-log plot of the density $\braket{n(\tau)}_{\mathrm{GGE}}$ as a
function of the rescaled time $\tau=\Gamma_{2\alpha}t$. The blue solid line
refers to the initial FS state (13) at filling $n_{0}=0.6$, while the red-
dashed line corresponds to the initial incoherent state (14) at the same
filling $n_{0}=0.6$. In the latter case, the tGGE rate equation reduces to the
law of mass action (7) and the density decays asymptotically as
$\braket{n(\tau)}_{\mathrm{GGE}}=\braket{n}_{\mathrm{MF}}(\tau)\sim\tau^{-1}$.
In the case of the FS state, on the contrary, the power-law exponent changes
$\braket{n(\tau)}_{\mathrm{GGE}}\sim\tau^{-1/2}$. (b) Plot of the effective
exponent $\delta_{\mathrm{eff}}(\tau)$ (16) as a function of $\tau$ (log scale
only on the horizontal axis) for the FS initial state with $n_{0}=0.6$. The
effective exponent converges to the value $0.5$.

In Fig. 2(a), we plot the density of particles as a function of time for the
two different initial conditions (13) and (14). In the case of the incoherent
initial state (14) (red-dashed line), we observe that the momentum
distribution function $C_{q}(\tau)$ remains flat in $q$ at all times so that
Eq. (17) exactly reduces to the classical reaction-limited equation (7) (with
$k=2$ therein). The asymptotic decay exponent is therefore the mean-field one
$\braket{n(\tau)}_{\mathrm{GGE}}\sim\tau^{-1}$. The coherences in the initial
state (13), on the other hand, deeply affect the asymptotic decay (blue-solid
line), whose exponent

$\braket{n(\tau)}_{\mathrm{GGE}}\sim\tau^{-1/2},$ (19)

differs from the mean-field one. The algebraic decay exponent $1/2$ is in
agreement with previous results concerning fermionic gases subject to two-body
losses in a lattice Rossini  _et al._ (2021); Perfetto  _et al._ (2022) and in
continuum space Rosso  _et al._ (2022). In Fig. 2(b), one sees that
$\delta_{\mathrm{eff}}(\tau)$ (16) neatly converges to $0.5$. This non mean-
field asymptotic decay goes beyond the classical reaction-limited description,
captured by the mean-field equation (7), and it is determined by quantum
coherences in the initial state. Coherences in the initial state are, indeed,
fundamental as they allow for an inhomogeneous in $q$ initial occupation
function $C_{q}(0)$. For this kind of initial conditions, the second term in
Eq. (18) for $f_{2}(k,q)$ does not vanish in the sum over $k$ in Eq. (17) and
it is responsible for the deviation of the power-law exponent from the mean-
field description. The nontrivial function $f_{2}(k,q)$ is determined by the
fermionic statistics, which hinders double (or higher) occupancy of lattice
sites and it therefore determines the form of the considered jump operators
(3). The joint effect of fermionic statistics and quantum coherences in the
initial state is thus responsible for the emergence of universal behavior in
the quantum RD dynamics beyond mean-field. In the next Sections, we show that
this behavior is even richer in the case of three and four-body annihilation
reactions.

###  IV.2 Three-body annihilation

In the case of three-body annihilation, $3A\to\emptyset$ (4), the rate
equation (12) takes the form

$\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}=-\frac{C_{q}(\tau)}{L^{2}}\sum_{k,k^{\prime}}f_{3}(k,k^{\prime},q)C_{k}(\tau)C_{k^{\prime}}(\tau),$
(20)

with the rescaled time $\tau=\Gamma_{3\alpha}t$ and the function
$f_{3}(k,k^{\prime},q)$ given by

$f_{3}(k,k^{\prime},q)=2\left[\sin(k-q)-\sin(k-k^{\prime})+\sin(q-k^{\prime})\right]^{2}.$
(21)

Equations (20) and (21) are obtained by first expressing the fermionic
operators $c_{j}$ in $L_{j}^{3\alpha}$ (4) in terms of the Fourier-space
operators $\hat{c}_{k}$ and then computing the right hand side of Eq. (12) via
Wick’s theorem, exploiting the fact that the tGGE state (11) is Gaussian and
diagonal in momentum space. The derivation is reported in Appendix A.

In Fig. 3(a), we show the density $\braket{n(\tau)}_{\mathrm{GGE}}$ of
particles as a function of time $\tau$ both for the FS (13) and the incoherent
(14) initial state. In the latter case, similarly as in the case of
$2A\to\emptyset$, elementary manipulations on the function
$f_{3}(k,k^{\prime},q)$ show that Eq. (20) reduces to

$\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}=-3C_{q}(\tau)\braket{n}^{2}_{\mathrm{GGE}}(\tau).$
(22)

Taking the sum over the quasi-momenta $q$ on both sides of the previous
equation, the law of mass action (7) with $k=3$ is retrieved. In this case the
quantum reaction-limited regime therefore exactly coincides with its classical
analogue and the asymptotic decay is
$\braket{n(\tau)}_{\mathrm{GGE}}=\braket{n}_{\mathrm{MF}}(\tau)\sim\tau^{-1/2}$.

Figure 3: Quantum three-body annihilation dynamics in the reaction limited
regime. (a) Log-log plot of the density $\braket{n(\tau)}_{\mathrm{GGE}}$ as a
function of the rescaled time $\tau=\Gamma_{3\alpha}t$. The initial states
chosen are identical to those considered in Fig. 2 for $2A\to\emptyset$. In
particular, the blue solid line refers to the initial FS state (13), while the
red-dashed line corresponds to the initial incoherent state (14). In both the
cases the same initial filling value $n_{0}=0.6$ is taken. In the case of the
initial state (14), the tGGE rate equation reduces to the law of mass action
(7) and (22) and the density decays asymptotically as
$\braket{n}_{\mathrm{GGE}}=\braket{n}_{\mathrm{MF}}\sim\tau^{-1/2}$. In the
case of the FS state, power-law decay is valid only up to $\tau\lesssim
10^{5}$ with a non mean-field exponent
$\braket{n}_{\mathrm{GGE}}\sim\tau^{-0.25}$. (b) Plot of the effective
exponent $\delta_{\mathrm{eff}}(\tau)$ (16) as a function of $\tau$ (log scale
only on the horizontal axis). The effective exponent converges only for times
$\tau\lesssim 10^{5}$. For longer times, plotted in the inset,
$\delta_{\mathrm{eff}}(\tau)$ slowly drifts in time pointing out that the
decay acquires a non-algebraic correction. This correction is slower than the
logarithmic correction (9) observed for the classical analogue of the dynamics
in $d=1$ in the different diffusion-limited regime.

The case of the FS initial state (13) yields instead a much richer physics
with a density dynamics which is markedly different from the mean-field
prediction. At first inspection, from Fig. 3(a), the density
$\braket{n(\tau)}_{\mathrm{GGE}}$ seems to follow a power-law in time $\tau$
with an exponent different from the mean-field value, similarly as in the case
of $2A\to\emptyset$ in Fig. 2. The calculation of the effective exponent
$\delta_{\mathrm{eff}}(\tau)$, in the main panel of Fig. 3(b), shows that the
density follows the power-law behavior for times of the order $\tau\lesssim
10^{5}$

$\braket{n(\tau)}_{\mathrm{GGE}}\sim\tau^{-0.25},\,\,\mbox{for}\,\,\tau\lesssim
10^{5}.$ (23)

Remarkably, in the quantum reaction-limited regime, the power-law decay
exponent is different from the mean-field prediction ($1/2$) even in
$d=d_{c}=1$. This contrasts the classical RD dynamics, where the power-law
decay is always determined by the mean-field exponent $1/2$, both in the
reaction-limited regime (7) and in the diffusion-limited one (9). In the
latter case (9), deviations from the mean-field are solely given by the
logarithmic correction $\sim(\ln\Omega t)^{1/2}$, the exponent of the
algebraic decay $(\Omega t)^{-1/2}$ still being equal to $1/2$. The appearance
of the non mean-field exponent at $d=d_{c}$ is only possible in the quantum RD
dynamics where quantum effects, not determined by space dimensionality, induce
correlations beyond mean-field. In order to observe the decay dynamics (23),
one, indeed, needs the simultaneous presence of quantum coherences in the
initial state and the fermionic statistics. The latter forces the annihilation
reaction to take place only among adjacent particles (4), which is reflected
in the non-trivial function $f_{3}(k,k^{\prime},q)$ in Eq. (20). Ultimately,
this function renders the dynamics different from the law of mass action
prediction (22) for quantum coherent initial conditions.

A closer inspection of the effective exponent $\delta_{\mathrm{eff}}(\tau)$,
in the inset of Fig. 3(b) for $\tau>10^{5}$, shows, however, that the power-
law (23) observed for triplet annihilation, $3A\to\emptyset$, remains valid
only up to times $\tau\lesssim 10^{5}$. For longer times, the effective
exponent $\delta_{\mathrm{eff}}(\tau)$ slowly increases as a function of time,
indicating that the behavior in Eq. (23) acquires a non-trivial non power-law
correction. It is then natural to attempt to link this behavior to the
logarithmic correction in Eq. (9) for the classical triplet annihilation in
the diffusion-limited regime. Our analysis, however, shows that a logarithmic
correction to the power-law (23) as $\braket{n(\tau)}_{\mathrm{GGE}}\sim
a/\tau^{0.25}+b\,(\log(t)/t)^{\gamma}$, with $a$ the amplitude of the decay in
Eq. (23) while $\gamma$ and $b$ are fitting parameters, is not compatible with
the behavior of $\delta_{\mathrm{eff}}(\tau)$ in Fig. 3(b). In particular, we
find that the asymptotic decay is slower than the one predicted by the
aforementioned logarithmic correction.

This highly non-trivial non-power-law behavior makes the dynamics of three-
body annihilation particularly distinct from their classical counterpart: in
the classical reaction-limited regime, indeed, only power-law behavior (7) is
possible since spatial fluctuations are absent due to the rapid mixing through
coherent hopping. In the diffusion-limited regime, instead, non-algebraic
asymptotic is possible only in the logarithmic form of Eq. (9) at $d=d_{c}$.
The latter is, in turn, determined by the spatial fluctuations induced by the
small diffusive mixing, as recalled in Subsec. III.1. The non-algebraic
behavior in Fig. 3 is non logarithmic and it has a different origin determined
not only by space dimensionality, but also by the interplay of the latter with
quantum effects due to coherences in the initial states.

###  IV.3 Four-body annihilation

In the case of quadruplet annihilation $4A\to\emptyset$ (5), equation (12)
reads as

$\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}\\!=\\!-\frac{C_{q}(\tau)}{L^{3}}\sum_{k,k^{\prime},k^{\prime\prime}}\\!\\!\\!f_{4}(k,k^{\prime},k^{\prime\prime},q)C_{k}(\tau)C_{k^{\prime}}(\tau)C_{k^{\prime\prime}}(\tau),$
(24)

with the rescaled time $\tau=\Gamma_{4\alpha}t$. The function
$f_{4}(k,k^{\prime},k^{\prime\prime},q)$ takes a rather cumbersome form. It
results from a long calculation involving the application of Wick’s theorem
onto the eight-point function of the fermionic operators $\hat{c}_{k}$
deriving from the Fourier expression of the quadruplet annihilation jump
operator (5). The expression for $f_{4}$ is therefore reported in Eq. (46) of
Appendix B, where the calculations leading to (24) are summarized.

In Fig. 4(a), the density $\braket{n(\tau)}_{\mathrm{GGE}}$ is shown as a
function of $\tau$ for the FS (13) initial condition as well as for an initial
incoherent state (14). The dynamics from the incoherent initial state reduces,
as in the case of binary $2A\to\emptyset$ and triplet $3A\to\emptyset$
annihilation, to the law of mass action (7) with $k=4$

$\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}=-4C_{q}(\tau)\braket{n(\tau)}^{3}_{\mathrm{GGE}}.$
(25)

The decay of the density is therefore ruled by
$\braket{n}_{\mathrm{GGE}}=\braket{n}_{\mathrm{MF}}\sim\tau^{-1/3}$. This is
power-law decay is plotted with the dashed-red line in Fig. 4(a). We
emphasize, as recalled in Subsec. III.1, that the decay (25) applies for
classical RD both in the diffusion-limited and in the reaction-limited regime
as for this process spatial fluctuations are irrelevant already in
$d=1>d_{c}$.

In the case of the FS initial state, however, the decay of the density,
depicted with the blue-solid line of Fig. 4(a),
$\braket{n(\tau)}_{\mathrm{GGE}}$ does not follow the mean-field prediction.
In particular, we observe the algebraic decay

$\braket{n(\tau)}_{\mathrm{GGE}}\sim\tau^{-0.1},$ (26)

obtained by computing the effective exponent $\delta_{\mathrm{eff}}(\tau)$
(16), which is plotted in Fig. 4(b). We observe that the convergence of
$\delta_{\mathrm{eff}}(\tau)$ to the asymptotic value is slower compared to
the case of binary annihilation $2A\to\emptyset$ in Fig. 2. In the latter case
$\delta_{\mathrm{eff}}(\tau)\simeq 0.5$ for $\tau\gtrsim 10^{5}$, while, in
Fig. 4(b), $\delta_{\mathrm{eff}}\simeq 0.1$ for $\tau\gtrsim 10^{11}$. In all
the cases discussed, $2A\to\emptyset$ (19), $3A\to\emptyset$ (23) and
$4A\to\emptyset$ (26), we observe that quantum coherences slow down the
density decay, which is reflected in an exponent that is smaller than the one
of the corresponding classical reaction-limited process (7).

Figure 4: Quantum four-body annihilation dynamics in the reaction limited
regime. (a) Log-log plot of the density $\braket{n(\tau)}_{\mathrm{GGE}}$ as a
function of the rescaled time $\tau=\Gamma_{4\alpha}t$. The blue solid line
represents the dynamics from the FS state (13), while the red dashed one the
dynamics from the state (14). The initial filling value is $n_{0}=0.6$, as in
the case of Figs. 2 and 3. For the incoherent initial state (14), the quantum-
reaction limited dynamics exactly coincides with the mean-field law of mass
action prediction (25). The associated asymptotic density decay is
$\braket{n}_{\mathrm{GGE}}=\braket{n}_{\mathrm{MF}}=(\tau)\sim\tau^{-1/3}$.
For the coherent Fermi-sea initial state, the decay is slower and it is given
by the power-law $\braket{n(\tau)}_{\mathrm{GGE}}\sim\tau^{-0.1}$. (b) Plot of
the effective exponent $\delta_{\mathrm{eff}}(\tau)$ (16) as a function of
$\tau$ (log scale only on the horizontal axis). The effective exponent
converges in a slower way than in the case of $2A\to\emptyset$ of Fig. 2. In
particular, observes that $\delta_{\mathrm{eff}}(\tau)\simeq 0.1$ at long
times $\tau\gtrsim 10^{11}$, as shown in the inset.

The non mean-field algebraic decay (26) is valid in $d=1$ and, therefore,
above the upper-critical dimension $d_{c}$ of the quadruplet annihilation
reaction $4A\to\emptyset$. This might look surprising as one expects for
$d>d_{c}$ the system dynamics to be captured by the mean-field approximation,
as reported in Sec. III. One should note that the upper critical dimension
$d_{c}$ characterizes diffusion-limited RD dynamics. In the diffusion-limited
regime, correlations leading to non mean-field behavior are caused by spatial
fluctuations in the density profile due to the diffusive motion of particles
far apart from each other. Spatial fluctuations are relevant in low
dimensions, while in higher dimensions diffusive mixing is effective in
filling the whole space. This explains the emergence of an upper critical
dimension $d_{c}$, beyond which critical exponents are given by the mean-field
prediction. In the reaction-limited regime, where Eq. (26) applies, spatial
fluctuations in the density are by construction absent as one assumes the
system to be relaxed to the homogeneous GGE state. Fluctuations beyond mean-
field are of strict quantum origin due quantum coherences in the initial
state. These effects are valid in any space dimensionality. Consequently one
has non mean-field behavior even for $d>d_{c}$. This is the hallmark of the
different origin of universal dynamical behavior in quantum reaction-limited
RD systems as compared to their classical counterparts.

On the basis of these results we expect also that in quantum RD systems with
five, $5A\to\emptyset$, and higher-body annihilation reactions will have
algebraic decay with a non mean-field exponent due to to the same quantum
effects [though possibly with a numerical value different from the one in Eq.
(26) for $4A\to\emptyset$]. The calculation of the decay exponents for such
higher order systems is similar to that presented here (and in Appendix B),
only more cumbersome due to combinatorics, since higher fermionic correlation
functions are involved (for instance a ten-point function for
$5A\to\emptyset$). In any case, such behavior is again expected to be distinct
to that of the corresponding classical RD dynamics, where for the reaction-
limited regime all processes are described by the mean-field law of mass
action, cf. (7).

##  V Discussion

In this manuscript, we analyzed the quantum RD dynamics of fermionic quantum
gases on a one-dimensional lattice subject to annihilation processes,
formulated in terms of a Lindblad master equation (1)-(2). Here classical
diffusive motion is replaced by quantum coherent hopping (6), while
annihilation reactions are irreversible and they are encoded into the jump
operators (3)-(5). We considered binary annihilation $2A\to\emptyset$ of pairs
of neighbouring particles, triplet annihilation $3A\to\emptyset$ and
quadruplet $4A\to\emptyset$ annihilation. We studied the dynamics from both
coherent initial states, such as the FS (13), and from incoherent initial
states as (14). Quantum effects enter through the coherent Hamiltonian hopping
and, possibly, through coherences within the initial state.

We solved the problem analytically in the thermodynamic limit by exploiting
the approximation of the dynamics obtained via the time-dependent generalized
Gibbs ensemble method (11) Lange  _et al._ (2018); Mallayya  _et al._ (2019);
Lange  _et al._ (2017); Lenarčič  _et al._ (2018). This method describes the
reaction-limited, or weak-dissipation, regime $\Gamma_{\nu}/\Omega\ll 1$,
where the irreversible reaction rate is much smaller than the coherent hopping
rate.

In all the three cases considered, reactions $2A\to\emptyset$,
$3A\to\emptyset$ and $4A\to\emptyset$, we observed that for incoherent initial
states the quantum reaction-limited dynamics reduces to classical mean-field
description in terms of the law of mass action (7). Quantum coherences in the
initial state are therefore necessary in order to have a non mean-field
universal decay of the particle density at long times. The fermionic
statistics is also important as it determines the structure (3)-(5) of the
jump operators, which, consequently, determines the non-trivial functions
$f_{2}$ [(18) in Eq. (17)], $f_{3}$ [(21) in Eq. (20)] and $f_{4}$ [(46) in
Eq. (24)]. These functions eventually render the tGGE rate equations different
from the law of mass action for quantum coherent initial conditions. In
particular, in the case of binary annihilation $2A\to\emptyset$, the density
decays algebraically in time as in Eq. (19) with exponent $1/2$ (cf. Fig. 2).
For triplet $3A\to\emptyset$ and quadruplet $4A\to\emptyset$ annihilation, the
impact of the simultaneous presence of quantum coherences and fermionic
statistics onto the asymptotic decay is even richer. In particular, for
$3A\to\emptyset$, we find algebraic decay as in Eq. (23) only in an
intermediate time regime $\tau\lesssim 10^{5}$ with exponent approximately
$0.25$. For later times, this decay acquires, however, a non power-law
correction. In the classical RD dynamics non-algebraic corrections are only
possible in the diffusion-limited regime with a logarithmic form (9), which
comes from the fact that the upper-critical dimension of the process is
$d_{c}=d=1$ Tauber ; Täuber  _et al._ (2005); Täuber (2014). Spatial density
fluctuations are therefore marginal in one dimension for the reaction
$3A\to\emptyset$. The decay we observe in the quantum reaction-limited regime
in Fig. 3 is, instead, slower than that predicted by a logarithmic correction.
For quadruplet annihilation $4A\to\emptyset$, we find the algebraic decay in
Eq. (26) with exponent approximately equal to $0.1$ shown in Fig. 4. This
result is in contrast with the classical description of the process which is
always, both in the diffusion and in the reaction-limited regime, in agreement
with the mean-field approximation (10). This is a consequence of the fact that
for $4A\to\emptyset$ spatial fluctuations are irrelevant in any physical
dimension since $d_{c}<1$. Our results for $3A\to\emptyset$ and
$4A\to\emptyset$ show that in quantum RD dynamics correlations beyond mean-
field are not only determined by spatial fluctuations, but also by inherently
quantum effects. These effects are present in any space dimensionality and
therefore universal behavior, without any classical correspondence, is
possible even at $d\geq d_{c}$.

As a future direction, one may extend the present analysis to bosonic systems.
In this case, we expect the absence of the exclusion principle to lead, in the
reaction-limited regime, to mean-field results. In the case of annihilation
channels exhibiting interference effects (see, e.g., Perfetto  _et al._
(2022); García-Ripoll  _et al._ (2009); Syassen  _et al._ (2008); Diehl  _et
al._ (2008)), one may, however, still see a universal non mean-field decay. At
the same time, it would be interesting to investigate the RD dynamics of
fermionic and bosonic gases in the continuum. The time-dependent GGE
description of the reaction-limited regime can be carried out in analogy to
the lattice gases discussed here. The extension of the results here presented
to spatial dimensions larger than one is important, as well. In particular,
since the universal behavior in the quantum reaction-limited regime is
dictated by quantum coherent effects, we expect the exponent for the algebraic
decay to be non-trivial even for $d>d_{c}$. However, one expects the impact of
these effects to depend on space dimensionality $d$ itself, and it would thus
be interesting to understand the dependence of the decay exponent on $d$.

###### Acknowledgements.

G.P. acknowledges support from the Alexander von Humboldt Foundation through a
Humboldt research fellowship for postdoctoral researchers. We are grateful for
financing from the Baden-Württemberg Stiftung through Project No.
BWST_ISF2019-23. We also acknowledge funding from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Project No.
435696605, through the Research Unit FOR 5413/1, Grant No. 465199066 and
through the Research Unit FOR 5522/1, Grant No. 499180199. F.C. is indebted to
the Baden-Württemberg Stiftung for the financial support by the Eliteprogramme
for Postdocs. We also acknowledge financial support from EPSRC Grant no.
EP/V031201/1.

##  Appendix A Three-body annihilation dynamics

In this Appendix we provide details regarding the derivation of Eqs. (20) and
(21) for triplet annihilation $3A\to\emptyset$. We consider the case of
periodic boundary conditions $c_{j+L}=c_{j}$ for the Hamiltonian (6). As the
analysis based on the tGGE method applies in the thermodynamic limit
$L\to\infty$, the choice of boundary conditions does not impact on the final
result in Eqs. (20) and (21). We first introduce the Fourier transform
$\hat{c}_{k}$ of the lattice fermionic operators $c_{j}$ as

$\hat{c}_{k_{n}}=\frac{1}{\sqrt{L}}\sum_{j=1}^{L}e^{-ik_{n}j}c_{j},\qquad\mbox{with
inverse}\qquad
c_{j}=\frac{1}{\sqrt{L}}\sum_{k_{n}}e^{ik_{n}j}\hat{c}_{k_{n}},$ (27)

where $k_{n}=2\pi n/L$ are the quasi-momenta and they are parametrized on the
lattice in terms of the integer number $n=1,2\dots L$. In the previous
summation the sum $\sum_{k_{n}}$ denotes a summation over the integer number
$n$. In the main text and in the rest of the Supplemental material we use for
brevity the shorter notation $\sum_{k_{n}}\to\sum_{k}$. When multiple
summations over the quasi-momenta are present, e.g.,
$\sum_{k_{1}}\sum_{k_{2}}\dots\sum_{k_{n}}$, we also use the compact notation
$\sum_{k_{1},k_{2}\dots k_{n}}$ (see Eqs. (20) and (24)). We write the triplet
annihilation jump operator $L_{j}^{3\alpha}$ in Fourier space as

$L_{j}^{3\alpha}=\sqrt{\Gamma_{3\alpha}}\,c_{j}c_{j+1}c_{j+2}=\frac{\sqrt{\Gamma_{3\alpha}}}{L^{3/2}}\sum_{k,k^{\prime},k^{\prime\prime}}e^{ij(k+k^{\prime}+k^{\prime\prime})}g(k^{\prime},k^{\prime\prime})\hat{c}_{k}\hat{c}_{k^{\prime}}\hat{c}_{k^{\prime\prime}},$
(28)

with

$g(k^{\prime},k^{\prime\prime})=\mbox{exp}(i(k^{\prime}+2k^{\prime\prime})).$
(29)

In order to simplify the commutator appearing in Eq. (12), we use the
identities

$[\hat{n}_{q},\hat{c}_{k}\hat{c}_{k^{\prime}}\hat{c}_{k^{\prime\prime}}]=-\hat{c}_{k}\hat{c}_{k^{\prime}}\hat{c}_{q}\delta_{k^{\prime\prime},q}-\hat{c}_{k}\hat{c}_{q}\hat{c}_{k^{\prime\prime}}\delta_{k^{\prime},q}-\hat{c}_{q}\hat{c}_{k^{\prime}}\hat{c}_{k^{\prime\prime}}\delta_{k,q},\quad\mbox{from}\quad[\hat{n}_{q},\hat{c}_{k}]=-\delta_{k,q}\hat{c}_{q}.$
(30)

Using Eq. (30), we then get for $[\hat{n}_{q},L_{j}^{3\alpha}]$

$\displaystyle[\hat{n}_{q},L_{j}^{3\alpha}]=-\frac{\sqrt{\Gamma_{3\alpha}}}{L^{3/2}}\left[\sum_{k,k^{\prime}}\hat{c}_{k}\hat{c}_{k^{\prime}}\hat{c}_{q}e^{i(k+k^{\prime}+q)j}g(k^{\prime},q)+\sum_{k,k^{\prime}}\hat{c}_{k}\hat{c}_{q}\hat{c}_{k^{\prime}}e^{i(k+k^{\prime}+q)}g(q,k^{\prime})+\sum_{k,k^{\prime}}\hat{c}_{q}\hat{c}_{k}\hat{c}_{k^{\prime}}e^{i(k+k^{\prime}+q)j}g(k,k^{\prime})\right].$
(31)

Upon inserting Eq. (31) into Eq. (12) one obtains three terms

$\displaystyle\sum_{j}(L_{j}^{3\alpha})^{\dagger}[\hat{n}_{q},L_{j}^{3\alpha}]$
$\displaystyle=-\frac{\Gamma_{3\alpha}}{L^{3}}\sum_{j}\sum_{k_{1},k_{2},k_{3}}e^{-ij(k_{1}+k_{2}+k_{3})}\hat{c}^{\dagger}_{k_{3}}\hat{c}^{\dagger}_{k_{2}}\hat{c}^{\dagger}_{k_{1}}g^{\ast}(k_{2},k_{3})\left[\sum_{k,k^{\prime}}g(k^{\prime},q)\hat{c}_{k}\hat{c}_{k^{\prime}}\hat{c}_{q}e^{ij(k+k^{\prime}+q)}\right.$
$\displaystyle\left.+\sum_{k,k^{\prime}}g(q,k^{\prime})\hat{c}_{k}\hat{c}_{q}\hat{c}_{k^{\prime}}e^{ij(k+k^{\prime}+q)}+\sum_{k,k^{\prime}}g(k,k^{\prime})\hat{c}_{q}\hat{c}_{k}\hat{c}_{k^{\prime}}e^{ij(k+k^{\prime}+q)}\right]$
$\displaystyle=-\frac{\Gamma_{3\alpha}}{L^{2}}\left[\sum_{k_{1},k_{2},k_{3},k^{\prime}}g^{\ast}(k_{2},k_{3})g(k^{\prime},q)\hat{c}_{k_{3}}^{\dagger}\hat{c}_{k_{2}}^{\dagger}\hat{c}_{k_{1}}^{\dagger}\hat{c}_{k_{1}+k_{2}+k_{3}-k^{\prime}-q}\hat{c}_{k^{\prime}}\hat{c}_{q}\right.$
$\displaystyle\quad\quad\quad\,\,+\sum_{k_{1},k_{2},k_{3},k^{\prime}}\left.g^{\ast}(k_{2},k_{3})g(q,k^{\prime})\hat{c}_{k_{3}}^{\dagger}\hat{c}_{k_{2}}^{\dagger}\hat{c}_{k_{1}}^{\dagger}\hat{c}_{k_{1}+k_{2}+k_{3}-k^{\prime}-q}\hat{c}_{q}\hat{c}_{k^{\prime}}\right.$
$\displaystyle\quad\quad\quad\,\,\left.+\sum_{k_{1},k_{2},k_{3},k^{\prime}}g^{\ast}(k_{2},k_{3})g(k_{1},k^{\prime})\hat{c}_{k_{3}}^{\dagger}\hat{c}_{k_{2}}^{\dagger}\hat{c}_{k_{1}+k^{\prime}+q-k_{2}-k_{3}}^{\dagger}\hat{c}_{q}\hat{c}_{k_{1}}\hat{c}_{k^{\prime}}\right],$
(32)

where in the last equality the Fourier representation of the Kronecker delta
has been used. From the previous equation, one sees that the Lindblad dynamics
of the two point function $\braket{\hat{n}_{q}}$ for the triplet annihilation
$3A\to\emptyset$ is coupled to the dynamics of six-point functions. The
evolution equation for $\braket{\hat{n}_{q}}$ is therefore not closed and one
has a hierarchy of equations coupling the dynamics of correlation functions to
higher-order correlation functions. In order to break this hierarchy, the tGGE
assumption (11) is fundamental. The time-dependent GGE state (11) describes
the dynamics in the reaction-limited, or, equivalently, weak-dissipation,
limit $\Gamma_{3\alpha}/\Omega\ll 1$ and it amounts to replace
$\braket{\hat{n}_{q}}\to\braket{\hat{n}_{q}}_{\mathrm{GGE}}$ (and analogously
for other expectation values in the previous equation). The GGE state is
Gaussian for the free-fermionic Hamiltonian (6) and diagonal in momentum
space. Higher-point correlation functions of fermionic operators $\hat{c}_{k}$
can therefore be computed in the GGE solely on the basis of the two-point
function $\braket{\hat{c}^{\dagger}_{q}\hat{c}_{k}}_{\mathrm{GGE}}\equiv
C_{q}\delta_{k,q}$ via Wick’s theorem. In the case of the six-point function
appearing on the right hand side of the first line of the second equality in
Eq. (32), one has, applying Wick’s theorem,

$\displaystyle\braket{\hat{c}_{k_{3}}^{\dagger}\hat{c}_{k_{2}}^{\dagger}\hat{c}_{k_{1}}^{\dagger}\hat{c}_{k_{1}+k_{2}+k_{3}-k^{\prime}-q}\hat{c}_{k^{\prime}}\hat{c}_{q}}_{\mathrm{GGE}}(\tau)=$
$\displaystyle-
C_{k^{\prime}}C_{q}C_{k_{3}}\delta_{k_{2},k^{\prime}}\delta_{q,k_{1}}+C_{k_{3}}C_{q}C_{k^{\prime}}\delta_{k_{1},k^{\prime}}\delta_{k_{2},q}+C_{k^{\prime}}C_{k_{2}}C_{k_{1}}\delta_{k_{3},k^{\prime}}\delta_{k_{1},q}$
$\displaystyle-
C_{k^{\prime}}C_{k_{2}}C_{k_{1}}\delta_{k_{3},k^{\prime}}\delta_{k_{2},q}+C_{k_{3}}C_{k^{\prime}}C_{k_{1}}\delta_{k_{3},q}\delta_{k_{2},k^{\prime}}-C_{k_{3}}C_{k_{2}}C_{k_{1}}\delta_{k_{3},q}\delta_{k_{1},k^{\prime}}.$
(33)

The other two six-point correlation functions in Eq. (32) can treated in an
analogous way. In the previous equation, all momentum occupation functions
$C_{k}(\tau)$ are function of the rescaled time $\tau=\Gamma_{3\alpha}t$ (not
reported explicitly for brevity). We also note that $[\hat{n}_{q},H]=0$ and
therefore the right hand side of the evolution equation for $\braket{n}_{q}$
is only given by Eq. (32). Upon using therein the tGGE assumption and the
expression in Eq. (33), one obtains Eqs. (20) and (21) of the main text (after
some lengthy but straightforward algebraic manipulation). In the numerical
solution of the equation, plotted in Fig. 3 of the main text, we take a value
of $L$ large and we check the stability of the obtained curves upon further
increasing $L$. In Fig. 3, we used $L=150000$ since for longer times,
$\tau\simeq 10^{7},10^{8}$ as in the inset of Fig. 3(b), larger values of $L$
are needed to get stable numbers.

We remark that Eqs. (20) and (21) can be written in the equivalent form

$\displaystyle\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}=-3C_{q}\braket{n}^{2}_{\mathrm{GGE}}(\tau)+\frac{C_{q}}{L^{2}}\sum_{k,k^{\prime}}C_{k}C_{k^{\prime}}$
$\displaystyle\left[\cos(2k+2k^{\prime})+2\cos(k-k^{\prime})+2\cos(2k-2q)-2\cos(k+k^{\prime}-2q)\right.$
$\displaystyle\left.+4\cos(k-q)-4\cos(2k-k^{\prime}-q)\right].$ (34)

This form of the equation makes more transparent the connection with the
classical reaction-limited dynamics ruled by the law of mass action (7). In
particular, the first term on the right hand side of (34) corresponds to the
law of mass action for $3A\to\emptyset$, upon summing on both sides over all
the possible values of $q$. The second term on the right hand side of the
previous equation, instead, couples the quasi-momentum $q$ to all the other
quasi-momenta through a three-body term $C_{q}C_{k}C_{k^{\prime}}$. The weight
of this term is, however, non trivial and it is determined by the fermionic
nature of the particles and the exclusion principle. Remarkably, this term
contributes to the dynamics only if the initial momentum occupation function
$C_{q}(0)$ is not flat in $q$. This amounts to considering initial states of
the GGE form (11) with $\lambda_{q}(0)$ not flat in $q$ and, therefore,
initial GGE states with a non diagonal, in real space, two-point fermionic
correlation matrix. These kind of states display coherences in real space. The
simulataneous presence of fermionic statistics and quantum coherences in the
initial state makes the second term of (34) relevant and it causes the non
mean-field universal behavior of Fig. 3 discussed in the main text.

##  Appendix B Four-body annihilation dynamics

In this Appendix we report the main steps of the derivation of Eq. (24) for
quadruplet annihilation $4A\to\emptyset$. We also give the expression of the
function $f_{4}(k,k^{\prime},k^{\prime\prime},q)$. We write the jump operator
$L_{j}^{4\alpha}$ in Eq. (5) in Fourier space (27) as

$L_{j}^{4\alpha}=\sqrt{\Gamma_{4\alpha}}\,c_{j}c_{j+1}c_{j+2}c_{j+3}=\frac{\sqrt{\Gamma_{4\alpha}}}{L^{2}}\sum_{k_{1},k_{2},k_{3},k_{4}}e^{ij(k_{1}+k_{2}+k_{3}+k_{4})}g(k_{2},k_{3},k_{4})\hat{c}_{k_{1}}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}},$
(35)

with

$g(k_{2},k_{3},k_{4})=\mbox{exp}(i(k_{2}+2k_{3}+3k_{4})).$ (36)

The commutator in Eq. (12) can be written by exploiting the chain rule for
commutators and (30) as

$\displaystyle[\hat{n}_{q},L_{j}^{3\alpha}]=-\frac{\sqrt{\Gamma_{4\alpha}}}{L^{2}}\Big{[}\sum_{k_{2},k_{3},k_{4}}e^{i(k_{2}+k_{3}+k_{4}+q)j}$
$\displaystyle[g(k_{2},k_{3},k_{4})\hat{c}_{q}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}}+g(q,k_{3},k_{4})\hat{c}_{k_{2}}\hat{c}_{q}\hat{c}_{k_{3}}\hat{c}_{k_{4}}$
$\displaystyle+g(k_{2},q,k_{4})\hat{c}_{k_{3}}\hat{c}_{k_{2}}\hat{c}_{q}\hat{c}_{k_{4}}+g(k_{2},k_{3},q)\hat{c}_{k_{4}}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{q}]\Big{]}.$
(37)

From Eq. (37) into (12) one has

$\displaystyle\sum_{j}(L_{j}^{4\alpha})^{\dagger}[\hat{n}_{q},L_{j}^{4\alpha}]=-\frac{\Gamma_{4\alpha}}{L^{4}}\sum_{j}\sum_{k,p,m,n}\sum_{k_{2},k_{3},k_{4}}$
$\displaystyle e^{ij(q+k_{2}+k_{3}+k_{4}-k-p-
m-n)}g^{\ast}(p,m,n)\hat{c}^{\dagger}_{n}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}$
$\displaystyle\qquad\qquad\qquad\qquad\Big{[}[g(k_{2},k_{3},k_{4})\hat{c}_{q}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}}+g(q,k_{3},k_{4})\hat{c}_{k_{2}}\hat{c}_{q}\hat{c}_{k_{3}}\hat{c}_{k_{4}}$
$\displaystyle\qquad\qquad\qquad\qquad+g(k_{2},q,k_{4})\hat{c}_{k_{3}}\hat{c}_{k_{2}}\hat{c}_{q}\hat{c}_{k_{4}}+g(k_{2},k_{3},q)\hat{c}_{k_{4}}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{q}]\Big{]}.$
(38)

The previous equation can be split into four terms:

$\displaystyle\frac{1}{L^{4}}\sum_{j}\sum_{k,p,m,n}\sum_{k_{2},k_{3},k_{4}}e^{ij(q+k_{2}+k_{3}+k_{4}-k-p-
m-n)}g^{\ast}(p,m,n)g(k_{2},k_{3},k_{4})\hat{c}^{\dagger}_{n}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{q}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}}$
$\displaystyle=\frac{1}{L^{3}}\sum_{k,p,m}\sum_{k_{2},k_{3},k_{4}}g^{\ast}(p,m,q+k_{2}+k_{3}+k_{4}-k-p-m)g(k_{2},k_{3},k_{4})\hat{c}^{\dagger}_{q+k_{2}+k_{3}+k_{4}-k-p-m}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{q}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}},$
(39)

$\displaystyle\frac{1}{L^{4}}\sum_{j}\sum_{k,p,m,n}\sum_{k_{2},k_{3},k_{4}}e^{ij(q+k_{2}+k_{3}+k_{4}-k-p-
m-n)}g^{\ast}(p,m,n)g(q,k_{3},k_{4})\hat{c}^{\dagger}_{n}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{q}\hat{c}_{k_{3}}\hat{c}_{k_{4}}$
$\displaystyle=\frac{1}{L^{3}}\sum_{k,p,m}\sum_{k_{2},k_{3},k_{4}}g^{\ast}(p,m,q+k_{2}+k_{3}+k_{4}-k-p-m)g(q,k_{3},k_{4})\hat{c}^{\dagger}_{q+k_{2}+k_{3}+k_{4}-k-p-m}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{q}\hat{c}_{k_{3}}\hat{c}_{k_{4}},$
(40)

$\displaystyle\frac{1}{L^{4}}\sum_{j}\sum_{k,p,m,n}\sum_{k_{2},k_{3},k_{4}}e^{ij(q+k_{2}+k_{3}+k_{4}-k-p-
m-n)}g^{\ast}(p,m,n)g(k_{3},q,k_{4})\hat{c}^{\dagger}_{n}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{q}\hat{c}_{k_{4}}$
$\displaystyle=\frac{1}{L^{3}}\sum_{k,p,m}\sum_{k_{2},k_{3},k_{4}}g^{\ast}(p,m,q+k_{2}+k_{3}+k_{4}-k-p-m)g(k_{3},q,k_{4})\hat{c}^{\dagger}_{q+k_{2}+k_{3}+k_{4}-k-p-m}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{q}\hat{c}_{k_{4}},$
(41)

$\displaystyle\frac{1}{L^{4}}\sum_{j}\sum_{k,p,m,n}\sum_{k_{2},k_{3},k_{4}}e^{ij(q+k_{2}+k_{3}+k_{4}-k-p-
m-n)}g^{\ast}(p,m,n)g(k_{3},k_{4},q)\hat{c}^{\dagger}_{n}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}}\hat{c}_{q}$
$\displaystyle=\frac{1}{L^{3}}\sum_{k,p,m}\sum_{k_{2},k_{3},k_{4}}g^{\ast}(p,m,q+k_{2}+k_{3}+k_{4}-k-p-m)g(k_{3},k_{4},q)\hat{c}^{\dagger}_{q+k_{2}+k_{3}+k_{4}-k-p-m}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}}\hat{c}_{q}.$
(42)

The four summations in Eqs. (39)-(42) can be grouped into a single sum
exploiting the fermionic anticommutation relations, so that Eq. (38) reads as

$\displaystyle\sum_{j}$
$\displaystyle(L_{j}^{4\alpha})^{\dagger}[\hat{n}_{q},L_{j}^{4\alpha}]=$
$\displaystyle=\frac{1}{L^{3}}\sum_{k,p,m}\sum_{k_{2},k_{3},k_{4}}g^{\ast}(p,m,q+k_{2}+k_{3}+k_{4}-k-p-m)f(k_{2},k_{3},k_{4},q)\hat{c}^{\dagger}_{q+k_{2}+k_{3}+k_{4}-k-p-m}\hat{c}^{\dagger}_{m}\hat{c}^{\dagger}_{p}\hat{c}^{\dagger}_{k}\hat{c}_{k_{2}}\hat{c}_{k_{3}}\hat{c}_{k_{4}}\hat{c}_{q},$
(43)

with the function $f(k_{2},k_{3},k_{4},q)$ given by

$f(k_{2},k_{3},k_{4},q)=g(k_{2},k_{3},k_{4})-g(q,k_{3},k_{4})+g(k_{3},q,k_{4})-g(k_{3},k_{4},q).$
(44)

From the previous equation the calculation proceeds similarly as in Appendix A
for $3A\to\emptyset$. In particular, one resorts to the tGGE approximation in
order to decompose the eight-point fermionic correlation function in Eq. (43)
in terms of the two-point function $C_{k}(t)$ via Wick’s theorem. A lengthy
calculation leads after some algebraic manipulations to Eq. (24)

$\frac{\mbox{d}C_{q}(\tau)}{\mbox{d}\tau}\\!=\\!-\frac{C_{q}(\tau)}{L^{3}}\sum_{k,k^{\prime},k^{\prime\prime}}\\!\\!\\!f_{4}(k,k^{\prime},k^{\prime\prime},q)C_{k}(\tau)C_{k^{\prime}}(\tau)C_{k^{\prime\prime}}(\tau),$
(45)

with $f_{4}(k,k^{\prime},k^{\prime\prime},q)$ written as

$\displaystyle f_{4}(k,k^{\prime},k^{\prime\prime},q)$
$\displaystyle=4-6\cos(k-k^{\prime})-6\cos(k-q)-4\cos(2k-2k^{\prime})-4\cos(2k-2q)-2\cos(3k-3k^{\prime})-2\cos(3k-3q)$
$\displaystyle+2\cos(2k-k^{\prime}-k^{\prime\prime})+2\cos(k-2k^{\prime}+k^{\prime\prime})+8\cos(2k^{\prime}-k-q)+4\cos(k+k^{\prime}-2q)$
$\displaystyle+4\cos(3k-2k^{\prime}-k^{\prime\prime})+4\cos(3k-2k^{\prime}-q)+4\cos(3k-k^{\prime}-2q)+4\cos(2k+k^{\prime}-3q)$
$\displaystyle-6\cos(3k-k^{\prime}-k^{\prime\prime}-q)-2\cos(k+k^{\prime}+k^{\prime\prime}-3q)+4\cos(k+k^{\prime}-k^{\prime\prime}-q)+2\cos(3k-3k^{\prime}+k^{\prime\prime}-q)$
$\displaystyle+2\cos(3k+k^{\prime}-k^{\prime\prime}-3q)-4\cos(2k+k^{\prime}-2k^{\prime\prime}-q)-4\cos(2k+k^{\prime}-k^{\prime\prime}-2q)$
$\displaystyle-2\cos(3k-2k^{\prime}-2k^{\prime\prime}+q)-4\cos(3k-2k^{\prime}+k^{\prime\prime}-2q)-2\cos(2k-k^{\prime}+2k^{\prime\prime}-3q)$
$\displaystyle+4\cos(2k-2k^{\prime}+2k^{\prime\prime}-2q).$ (46)

Notice that in passing from Eqs. (43) and (44) to Eqs. (45) and (46) we
renamed the dummy summation variables $m,p$ with $k^{\prime},k^{\prime\prime}$
in order to match the notation used in the main text in Eq. (24). We also
introduced the rescaled time $\tau=\Gamma_{4\alpha}t$. Equations (45) and (46)
have been numerically solved in order to the get the results in Fig. 4. At the
technical level, we remark that the calculation of the effective exponent
$\delta_{\mathrm{eff}}(\tau)$ at long times, $\tau\simeq 10^{11}$ as in the
inset in Fig. 4(b), requires a large value of $L$. We used $L=180000$ and we
checked that the obtained numbers for $\delta_{\mathrm{eff}}(\tau)$ are stable
upon further increasing $L$ to $L=240000$.

In the case of Eqs. (45) and (46) considerations similar to those done in
Appendix A for Eq. (34) apply. In particular, the factor $4$ on the first line
of the right hand side of (46) represents the law of mass action term (25).
All the other remaining terms couple the quasi-momentum $q$ to all the other
quasi momenta through the non-trivial form of the function $f_{4}$ determined
by the fermionic statistics. These terms contribute to the dynamics when the
initial occupation function $C_{q}(0)$ is not flat in $q$ and the initial
state therefore displays quantum coherences. This determines the universal
behavior of Fig. 4. We remark that this collective non mean-field behavior is
valid already in one dimension, above the upper-critical dimension of
$4A\to\emptyset$, as it is originates from quantum effects, due to quantum
coherences in the initial state, that are present in any spatial dimension.

## References

  * Wilson and Kogut (1974) K. G. Wilson and J. Kogut, Phys. Rep. 12, 75 (1974).
  * Kogut (1979) J. B. Kogut, Rev. Mod. Phys. 51, 659 (1979).
  * 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).
  * (6) 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).
  * Kang and Redner (1985) K. Kang and S. Redner, Phys. Rev. A 32, 435 (1985).
  * Privman and Grynberg (1992) V. Privman and M. D. Grynberg, J. Phys. A. Math. Gen. 25, 6567 (1992).
  * Toussaint and Wilczek (1983) D. Toussaint and F. Wilczek, J. Chem. Phys. 78, 2642 (1983).
  * Spouge (1988) J. L. Spouge, Phys. Rev. Lett. 60, 871 (1988).
  * Privman (1994) V. Privman, Phys. Rev. E 50, 50 (1994).
  * Torney and McConnell (1983) D. C. Torney and H. M. McConnell, J. Phys. Chem. 87, 1941 (1983).
  * Kang and Redner (1984a) K. Kang and S. Redner, Phys. Rev. A 30, 2833 (1984a).
  * Kang and Redner (1984b) K. Kang and S. Redner, Phys. Rev. Lett. 52, 955 (1984b).
  * Kang  _et al._ (1984) K. Kang, P. Meakin, J. Oh, and S. Redner, J. Phys. A: Math. Gen. 17, L665 (1984).
  * Rácz (1985) Z. Rácz, Phys. Rev. Lett. 55, 1707 (1985).
  * Doi (1976) M. Doi, J. Phys. A: Math. Gen. 9, 1479 (1976).
  * Peliti (1985) L. Peliti, J. Phys. France 46, 1469 (1985).
  * Peliti (1986) L. Peliti, J. Phys. A: Math. Gen. 19, L365 (1986).
  * Mattis and Glasser (1998) D. C. Mattis and M. L. Glasser, Rev. Mod. Phys. 70, 979 (1998).
  * 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, Science 320, 1329 (2008).
  * 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).
  * 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).
  * Kinoshita  _et al._ (2005) T. Kinoshita, T. Wenger, and D. S. Weiss, Phys. Rev. Lett. 95, 190406 (2005).
  * Tomita  _et al._ (2017) T. Tomita, S. Nakajima, I. Danshita, Y. Takasu, and Y. Takahashi, Sci. Adv. 3, e1701513 (2017).
  * Bouchoule and Schemmer (2020) I. Bouchoule and M. Schemmer, SciPost Phys. 8, 60 (2020).
  * 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).
  * Ferlaino  _et al._ (2009) F. Ferlaino, S. Knoop, M. Berninger, W. Harm, J. P. D’Incao, H.-C. Nägerl, and R. Grimm, Phys. Rev. Lett. 102, 140401 (2009).
  * Gurian  _et al._ (2012) J. H. Gurian, P. Cheinet, P. Huillery, A. Fioretti, J. Zhao, P. L. Gould, D. Comparat, and P. Pillet, Phys. Rev. Lett. 108, 023005 (2012).
  * Diehl  _et al._ (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. Büchler, and P. Zoller, Nat. Phys. 4, 878 (2008).
  * 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).
  * 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).
  * Ates  _et al._ (2012) C. Ates, B. Olmos, W. Li, and I. Lesanovsky, Phys. Rev. Lett. 109, 233003 (2012).
  * Everest  _et al._ (2014) B. Everest, M. R. Hush, and I. Lesanovsky, Phys. Rev. B 90, 134306 (2014).
  * Bouchoule  _et al._ (2021) I. Bouchoule, L. Dubois, and L.-P. Barbier, Phys. Rev. A 104, L031304 (2021).
  * 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).
  * Rosso  _et al._ (2023) L. Rosso, A. Biella, J. De Nardis, and L. Mazza, Phys. Rev. A 107, 013303 (2023).
  * Bouchoule  _et al._ (2020) I. Bouchoule, B. Doyon, and J. Dubail, SciPost Phys. 9, 44 (2020).
  * Lesanovsky and Garrahan (2013) I. Lesanovsky and J. P. Garrahan, Phys. Rev. Lett. 111, 215305 (2013).
  * Olmos  _et al._ (2014) B. Olmos, I. Lesanovsky, and J. P. Garrahan, Phys. Rev. E 90, 042147 (2014).
  * Everest  _et al._ (2016) B. Everest, M. Marcuzzi, J. P. Garrahan, and I. Lesanovsky, Phys. Rev. E 94, 052108 (2016).
  * Marcuzzi  _et al._ (2016) M. Marcuzzi, M. Buchhold, S. Diehl, and I. Lesanovsky, Phys. Rev. Lett. 116, 245701 (2016).
  * Buchhold  _et al._ (2017) M. Buchhold, B. Everest, M. Marcuzzi, I. Lesanovsky, and S. Diehl, Phys. Rev. B 95, 014308 (2017).
  * Gutiérrez  _et al._ (2017) R. Gutiérrez, C. Simonelli, M. Archimi, F. Castellucci, E. Arimondo, D. Ciampini, M. Marcuzzi, I. Lesanovsky, and O. Morsch, Phys. Rev. A 96, 041602 (2017).
  * Roscher  _et al._ (2018) D. Roscher, S. Diehl, and M. Buchhold, Phys. Rev. A 98, 062117 (2018).
  * Wintermantel  _et al._ (2020) T. M. Wintermantel, Y. Wang, G. Lochead, S. Shevate, G. K. Brennen, and S. Whitlock, Phys. Rev. Lett. 124, 070503 (2020).
  * Helmrich  _et al._ (2020) S. Helmrich, A. Arias, G. Lochead, T. M. Wintermantel, M. Buchhold, S. Diehl, and S. Whitlock, Nature 577, 481 (2020).
  * Nigmatullin  _et al._ (2021) R. Nigmatullin, E. Wagner, and G. K. Brennen, Phys. Rev. Research 3, 043167 (2021).
  * Kazemi and Weimer (2021) J. Kazemi and H. Weimer, arXiv:211.05352 (2021).
  * Carollo  _et al._ (2022) F. Carollo, M. Gnann, G. Perfetto, and I. Lesanovsky, Phys. Rev. B 106, 094315 (2022).
  * van Horssen and Garrahan (2015) M. van Horssen and J. P. Garrahan, Phys. Rev. E 91, 032132 (2015).
  * Carollo  _et al._ (2019) F. Carollo, E. Gillman, H. Weimer, and I. Lesanovsky, Phys. Rev. Lett. 123, 100604 (2019).
  * 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).
  * Zeiher  _et al._ (2016) J. Zeiher, R. Van Bijnen, P. Schauß, S. Hild, J.-y. Choi, T. Pohl, I. Bloch, and C. Gross, Nat. Phys. 12, 1095 (2016).
  * Kim  _et al._ (2018) H. Kim, Y. Park, K. Kim, H.-S. Sim, and J. Ahn, Phys. Rev. Lett. 120, 180502 (2018).
  * Ebadi  _et al._ (2021) S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, _et al._ , Nature 595, 227 (2021).
  * Jo and Kim (2022) M. Jo and M. Kim, arXiv:2203.14295 (2022).
  * Perfetto  _et al._ (2022) G. Perfetto, F. Carollo, J. P. Garrahan, and I. Lesanovsky, arXiv:2209.09784 (2022).
  * 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).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * Lindblad (1976) G. Lindblad, Comm. Math. Phys. 48, 119 (1976).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Feller (1968) W. Feller,  _An introduction to probability theory and its applications_ (Wiley, 1968).
  * Essler and Fagotti (2016) F. H. Essler and M. Fagotti, J. Stat. Mech.: Theory Exp. 2016, 064002 (2016).
  * Vidmar and Rigol (2016) L. Vidmar and M. Rigol, J. Stat. Mech.: Theory Exp. 2016, 064007 (2016).