Quantum reaction-limited reaction–diffusion dynamics of noninteracting Bose gases

We investigate quantum reaction–diffusion (RD) systems in one-dimension with bosonic particles that coherently hop in a lattice, and when brought in range react dissipatively. Such reactions involve binary annihilation ( A+A→∅ ) and coagulation ( A+A→A ) of particles at distance d. We consider the reaction-limited regime, where dissipative reactions take place at a rate that is small compared to that of coherent hopping. In classical RD systems, this regime is correctly captured by the mean-field approximation. In quantum RD systems, for noninteracting fermionic systems, the reaction-limited regime recently attracted considerable attention because it has been shown to give universal power law decay beyond mean field for the density of particles as a function of time. Here, we address the question whether such universal behavior is present also in the case of the noninteracting Bose gas. We show that beyond mean-field density decay for bosons is possible only for reactions that allow for destructive interference of different decay channels. Furthermore, we study an absorbing-state phase transition induced by the competition between branching A→A+A , decay A→∅ and coagulation A+A→A . We find a stationary phase-diagram, where a first and a second-order transition line meet at a bicritical point which is described by tricritical directed percolation. These results show that quantum statistics significantly impact on both the stationary and the dynamical universal behavior of quantum RD systems.


I. INTRODUCTION
The systematic classification of universal properties in nonequilibrium many-body systems is a timely and challenging research direction in statistical mechanics.In this field, reaction-diffusion (RD) models are an important class of classical many-body nonequilibrium systems where dynamical critical behavior has been thoroughly studied [1][2][3][4][5][6].These systems describe particles which diffuse at rate Ω through random motion in space and, when brought in contact, can react with each other at rate Γ.The particle density n(t) of RD systems at long times t decays as a power law.The latter is described in the "reaction-limited regime" of strong diffusion Γ ≪ Ω by mean-field approximation [1,3,4,[7][8][9].Fast diffusion, indeed, quickly erases spatial structures rendering the density profile homogeneous in space.Under these conditions, the dynamics is described by law of mass action rate equations, which for reactions as mA → lA (l < m) predict n(t) ∼ (Γt) −1/(m− 1) . ( In particular, annihilation 2A → ∅ at rate Γ α (m = 2, l = 1) and coagulation 2A → A (m = 2, l = 1) at rate Γ γ within mean field behave as n(t) ∼ (Γ α,γ t) −1 .In the opposite "diffusion-limited regime" Γ/Ω ∼ 1, instead, spatial fluctuations are relevant.Annihilation and coagulation, for instance, belong the same universality class and in one dimension, D = 1, in both the cases one finds algebraic decay [5,[10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] beyond mean field For D > 1, instead, the mean-field prediction ( 1) is retrieved as diffusion is effective in filling the whole available space thus eventually rendering the density homogeneous.Reaction-diffusion systems display also stationary-state universal properties.This happens when branching processes A → A + A, at rate Γ β , injecting particles into the system are considered.Namely, in the contact process (CP) [2,3,5,6], the competition between branching and one-body decay A → ∅, at rate Γ δ , induces a stationarystate phase transition.In the thermodynamic limit and at long times, tuning the relative strength Γ β /Γ δ between the branching and the decay rate, one has either an active phase, where the stationary density is nonzero, or an inactive phase with zero density.In the latter phase, the stationary state is an absorbing-state void of particles and no fluctuations are possible.The absorbing-state phase transition of the contact process is of second-order and it belongs to the directed percolation universality class [1,3,4,7,8].For restricted occupancy of the lattice, the stationary density in the active phase is finite, while in the unrestricted case one needs to include also annihilation and/or coagulation in order to have a finite stationary density [6,19,25].
Quantum RD systems have been treated analytically, in the thermodynamic limit, only in the reaction-limited regime in Refs.[33][34][35][36][37][38][39]41] for fermionic systems both in the lattice [33,34,36,38,39] and in the continuum [35,37,41].Spin systems subject to m-body losses have been also considered in Ref. [40].The studies [34][35][36][37][38][39][40][41] are all based on the time-dependent generalized Gibbs ensemble (TGGE) method [82][83][84][85][86], while [41] employs Keldysh field-theory diagrammatics to show the emergence of the TGGE dynamical equation.It has been found in these references that universal decay beyond mean field is present for initial states displaying quantum coherences in real space.The law of mass action prediction (1) is, instead, recovered when the dynamics starts from incoherent initial states.This is specifically true for annihilation processes of m neighbouring particles mA → ∅ [38][39][40] (on site multibody losses are not possible for fermions).For 2A → ∅, for coherent initial states, one specifically finds n(t) ∼ (Γ α t) −1/2 .Coagulation reactions 2A → A, on the contrary, decay as n(t) ∼ (Γ γ t) −1 with the mean-field exponent is found both for coherent and incoherent initial states.This, in particular, shows that for quantum systems binary annihilation and coagulation do not belong to the same universality class, in stark contrast to the classical case.In Refs.[34,38], annihilation reactions that create quantum superpositions between two nearest-neighbor decay channels were considered.Annihilation jump operators with this feature introduce inherent quantum interference effects and allows for a non-mean-field decay n(t) ∼ (Γ α t) −1/2 for both coherent and incoherent initial states.For the case where the CP is considered with branching and decay, the associated absorbing-state phase transition is of second order and it is found to belong to the mean-field directed percolation universality class.
The reaction-limited dynamics of the interacting Bose gas has been studied in Ref. [37].In the limiting case of no interaction, it is therein found that the density for onsite losses mA → ∅ decays according to the law of mass action (1).For other kinds of reactions additional analysis is, however, needed.In particular, it is interesting to understand whether reactions coupling bosons at different sites possibly creating superpositions between two, or more, decay channels induce beyond mean-field decay, as in the aforementioned fermionic cases.Furthermore, the quantum reaction-limited dynamics of the Bose gas in the presence of reactions, such as branching, creating particles has not been analyzed so far.In this perspective, for bosons, the universality class and the order of the absorbing-state phase transition taking place in the presence of branching is also not understood.
The goal of this manuscript is to address the impact of quantum bosonic statistics on the universal behavior of the quantum RD dynamics.In particular, we study the impact of the bosonic statistics both on the relaxation dynamics of the noninteracting Bose gas and on the universality class of its stationary phase diagram obtained in the presence of branching reactions.In both the relaxation and the stationary case, our goal is then to show that the universal properties of quantum RD dynamics of the Bose gas significantly differ both from the fermionic formulation and from their classical counterpart.
Our approach is based on analytically studying the quantum RD dynamics of the noninteracting Bose gas in the reaction-limited regime and in the thermodynamic limit.We employ the time-dependent generalised Gibbs ensemble (TGGE) method [82][83][84][85][86].We first consider distance-selective binary losses, which amount to a loss of two bosons at a distance d.We further allow for the possibility of interference between two or three decay channels of the distance selective loss.We then move to the case of onsite coagulation reactions and eventually consider the absorbing-state phase transition induced by the competition between coagulation, branching and onebody decay.
We find that distance selective losses eventually always render, independently of the initial state, the law of mass action prediction for any value of the distance d.This is in contrast to the case of fermionic systems, where non-mean-field behavior is found for nearest neighbor annihilation (d = 1).The case where interference between two annihilation decay channels is present is then discussed.When the channels are unequally weighted, we retrieve mean-field behaviour as in Eq. ( 3).For a fine-tuned equal balance between the two channels, instead, we find algebraic decay beyond mean-field for any value d ≥ 1. Namely, we observe the same asymptotic decay as in Refs.[34,38] for fermions, (4) The main difference with respect to the fermionic case is that in the latter case Eq. ( 4) does not require the two decay channels to be equally weighted.We eventually consider the case where the interference takes place between three neighbouring decay channels (rate Γ ᾱ).In this case, we also find beyond mean-field decay n(t) ∼ (Γ ᾱt) −0.28 , three-channels interference, 2A → ∅. ( for any d ≥ 1.These results show that for bosons interference effects within the annihilation decay channels are necessary in order to get beyond mean-field algebraic decays.For onsite coagulation, we find mean-field decay when coagulation is considered alone.This decay is valid independently of the initial state considered.Equation (6) implies that for bosons annihilation (3) and coagulation share the same decay exponent and they belong, at least in the reaction-limited regime, to the same universality class.This again contrasts the fermionic case of Ref. [38], where the same reaction-limited decay exponent is observed only for incoherent initial states.The competition between coagulation, one-body decay and branching in the CP reveals gives rise to an interesting phase diagram.It is of mean-field nature but considerably richer than the phase diagram obtained within the classical and fermionic reaction-limited cases.In particular, we find a change of the absorbing-state phase transition from first to second order as the relative strength Γ β /Γ δ between classical branching and decay is tuned.These first and second-order transition lines meet at a bicritical point, which is characterized by the mean-field exponents of tricritical directed percolation [87][88][89].Mean-field exponents are here observed already in one dimensions since the reaction-limited regime with fast hopping mixing is considered.A similar phase diagram has been previously observed in Refs.[57,58] for a different spin 1/2 model, where classical branching/coagulation competes with quantum branching/coagulation.Therein the order of the transition is changed as the relative strength between coherent and incoherent branching is tuned.The phase transitions observed in our work, however, solely rely on classical reactions, with the only quantum effect stemming from the hopping of the bosons.
This manuscript is structured as follows.In Sec.II, we discuss the formulation of quantum RD dynamics in terms of the quantum master equation.We further present all the reactions processes analysed in the text.In Sec.III, we briefly review known results of classical RD dynamics and then we discuss the quantum RD dynamics.In particular, we focus on the reaction-limited regime, which is the main subject of the text.We briefly introduce the TGGE method and we write the associated dynamical equation for the bosonic occupation function in momentum space.In Sec.IV, we specialize this equation to the various processes introduced in Sec.II.In Sec.V, we summarize and discuss our results.The Appendices A-D contain technical aspects and details regarding the calculations at the basis of the results presented in the main text.

II. THE SYSTEM
We consider bosonic particles on a quantum chain of length L with periodic boundary conditions.The j-th lattice site of the chain can therefore either be occupied, j bj the number operator at site j.Here, bj and b † j are bosonic destruction and creation operators, respectively, obeying the bosonic commutation rule [ bi , b † j ] = δ i,j .Each site can be occupied by n j bosons with n ∈ N, so that the total number of particles at each lattice site is unrestricted.The quantum reaction-diffusion dynamics of the system is ruled by the Lindblad master equation [42][43][44] (we set ℏ = 1 in the whole manuscript) where ρ is the density matrix.In the quantum formulation of the RD dynamics, classical-incoherent hopping (diffusion in the continuum limit) is replaced by coherent hopping of particles between neighbouring sites according to the Hamiltonian H: with the hopping rate Ω.The Hamiltonian, thus, conserves the number of particles ([H, N (t)] = 0, with N = j nj the total particle number).We note that H describes the noninteracting Bose gas (vanishing interaction strength) on a optical lattice.The dissipator D[ρ] accounts for the irreversible reaction processes and it is of Lindblad form where L ν j are dubbed jump operators.The superscript ν labels the different reaction types.We consider various reaction types, which are pictorially represented in Fig. 1.All these dissipative processes violate particle number conservation.Distance-selective binary annihilation A + A → ∅ (rate Γ α ) destroys pairs of particles at a distance d and is represented by the jump operator Here, the parameter θ ∈ [0, π) allows for interferences between two decay channels.At θ = 0 or θ = π/2, the annihilation operator reduces to the classical-incoherent binary annihilation of particles at distance d.Such incoherent distance-selective losses have been analyzed in Refs.[30,31] for hard-core bosons, where they were implemented via laser-exciting highly excited electronic Rydberg states in the facilitation regime.For θ ̸ = 0, π/2, on the other hand, the jump operators (10) effectuate transitions into quantum superposition states, thereby creating quantum coherence.Such setting, whose microscopic origin is the quantum interference of two annihilation channels, was previously studied in Refs.[34,38] for fermions.
To further study the impact of interferences onto the quantum bosonic RD dynamics, we consider a binary annihilation process (rate Γ ᾱ) which couple three decay channels between particles at fixed distance d.We refer to the jump operator (11) henceforth as second-order annihilation reaction as L ᾱ j is constructed such that in the continuum space limit, the zeroth and the first order of the expansion the jump operator in the lattice spacing vanish (see Appendix A for the details).We anticipate (see the discussion in Subsecs.IV B for further details) that both the jump operator (10), with θ = π/4, and ( 11) have a zero-momentum Bose-Einstein condensate state (BEC) as an exact dark state.Similar jump operators having the BEC as a many-body dark state, albeit still conserving the total particle number, have been studied, e.g., in Refs.[46,49].
Other types of reactions we are considering is onsite coagulation and one-body decay A → ∅ (rate Γ δ ) In all the cases (10)-( 13), the density of particles n(t) = ⟨ N (t)⟩ /L is a monotonically decreasing function of time t.The universal behavior of the dynamics lies in the asymptotic late-time approach of the system towards the steady state devoid of particles.
In order to have a nontrivial steady state, supporting a nonzero density of particles, we need to include reactions that increase the particle number (see Fig. 1).We do this by introducing onsite branching reactions which lead to the creation of an offspring from a lattice site already occupied by at least one particle.The form of the jump operators ( 12) and ( 14) is dictated by bosonic statistics, allowing reactions to take place at the same lattice site.The competition between branching (14) and decay (13) in systems with restricted occupation of lattice sites (spin or fermionic systems) leads to a second-order (continuous) stationary-state phase transitions separating an active (nonzero density) phase from an inactive-aborbing one (with zero density) [25].For bosons the maximum occupation number per site is unrestricted.Therefore, in order to have a finite density in the active phase also coagulation (12) (or annihilation (10)) has to be considered [6,19,25].Here, universal behavior emerges near the critical point of the ensuing second-order phase transition, as we recall in the next Sec.III.10) -( 14) of the Lindblad dynamics.Coagulation, branching and decay are sketched.Note that coagulation requires at least two particles on the lattice site (transition from darker-grey to brighter-grey circle), while decay can act also on single occupied site (transition from grey to white circle).In the middle line, the annihilation decay (10) at distance d = 1, for the sake of illustrative purposes, is shown.

Hopping
For θ ̸ = 0, π/2 coherent superposition between different quantum mechanical states is introduced into the dynamics.In the bottom line, the action of the jump operator (11) is sketched (d = 1 again).In this case, superposition of three different states is introduced.

III. DIFFUSION AND REACTION-LIMITED DYNAMICS
The dynamics of the density n(t) as a function time is the central quantity characterizing the emergent universal behavior of RD systems.Two different timescales rule the RD dynamics: the reaction time ∼ Γ −1 and the diffusion time ∼ Ω −1 .The former quantifies the time (on average) two nearby particles take to react, while the latter the time (on average) two particles take to get in contact.The dynamics changes qualitatively depending on the ratio Γ/Ω.The limit Γ/Ω ≪ 1 is named "reaction-limited regime", while the limit where Γ/Ω is at least Γ/Ω ∼ 1, is dubbed "diffusion-limited regime".The exponent δ of the asymptotic power-law decay n(t) ∼ t −δ is different in the two regimes.In Subsec.III A, we first briefly recall and discuss previous results concerning the bosonic RD dynamics for classical systems.In Subsec.III B, we then recall aspects of the quantum RD dynamics that are relevant for the treatment of the present manuscript.
In particular, we focus on the quantum-reaction limited regime, which is the main focus of this manuscript.

A. Classical RD Dynamics
Within the classical RD realm, the reaction-limited regime is well-described by the so-called "law of mass action" rate equation [1,3,4,7,8].This equation states that the rate of reactions equals the product of the concentrations of reactants.For the annihilation reaction, where m particles react to l particles (with m > l ∈ N), the law of mass action equation yields, Equation ( 15) is valid whenever reactions occur everywhere in space with the same probability, i.e., it is valid in a well-stirred regime where the density n(t) is homogeneous.This is precisely the reaction-limited regime Γ/Ω ≪ 1, where the fast hopping (diffusion in continuum space) mixing wipes out spatial patterns in the density profile coming from local depletion of the particle number.For m ≥ 2, the long-time behaviour associated to (15) is The exponent of the power-law decay depends on the number of reactants m involved in the reaction.Note that in this regime, the density is a function of the rescaled time τ = Γt according to the reaction rate Γ.This further shows that in the reaction-limited regime the limiting factor for the asymptotic decay of the density is the reaction time (controlled by Γ −1 ).The diffusion-limited regime, Γ/Ω ∼ 1, is qualitatively different.In this regime the hopping mixing is finite and therefore spatial fluctuations in the density profile, due to local depletion of the particle number, are relevant in the kinetics.For classical systems with unrestricted occupation, this problem has been analytically tackled in Refs.[5,[17][18][19][20][22][23][24] by mapping the classical microscopic master equation into a bosonic second-quantization problem via the Doi-Peliti formalism [20,21].From the latter, a field-theory description can be eventually derived.The associated renormalization group scaling analysis yields for 2A → ∅ in one dimension D = 1: Here, the density depends on the rescaled time τ = Ωt according to the diffusion rate.The decay (17) applies classically also for coagulation 2A → A, which belongs, indeed, to the same universality class of annihilation [90][91][92][93][94].The result (17) shows that the limit factor for the decay of the density of particles is, in this case, the time needed for two far apart particles to meet.This is a universal property of the random walk, which is recurrent only for D ≤ 2. This explain both the universal character of Eq. ( 17) and the emergence of the upper critical dimension D c = 2.For D ≥ D c , diffusion is effective in filling the whole space and one recover the mean-field prediction ( 16) (up to a logarithmic correction at Reaction-diffusion model can also host stationary state phase transitions between an active phase and an absorbing one.This happens when branching reactions A → A + A are included.The CP model, namely, describes the stationary phase transition ensuing from the competition between branching ( 14) and ( 13).In the case where site occupancy is unrestricted, the so-called bosonic contact process [6,25], one further needs to include coagulation (12) (or annihilation (10) at θ = 0 and d = 0) in order to have a finite stationary density.This can be heuristically understood by writing the law of mass action (15) in the presence of branching and decay This equation trivially admits the stationary solution n stat = 0 corresponding to the absorbing phase.Furthermore, one has the additional solution In the case of restricted CP, one does not need coagulation to have a finite stationary density since a term proportional to n 2 is produced on the right hand side of ( 18) due to the restricted occupation constraint [n(t)(1 − n(t)) from the branching].Equation ( 18) immediately locates the mean-field critical point Γ c β = Γ δ of the phase transition.This transition is classically of second order and it belongs to the directed percolation universality class [1][2][3][4][5][6].Universal behavior lies in the algebraic behavior of the order parameter n stat ∼ (Γ β − Γ c β ) β close to the transition point.The associated critical exponent is β = 1 within the reaction-limited regime.Away from this limit, the exponent β ̸ = 1 differs from the mean-field prediction and it depends on the space-dimensionality for D ≤ D c , and D c = 4 [1][2][3][4][5][6].We do not report these values here as they are not central for the understanding of the results of this manuscript, which is focused on the quantum analogue of Eqs. ( 18) and (19).We introduce the quantum reaction-limited regime in the next Subsection.We will see in Sec.IV that the ensuing quantum reaction-limited equation for the CP displays a much richer phase diagram than that of Eqs. ( 18) and (19).

B. Quantum RD Dynamics
For quantum RD systems little is currently known.The quantum master equation in Eqs. ( 7)-( 14) cannot be solved analytically since it is not quadratic as a consequence of the structure of the jump operators.At the same time, large-scale numerics, involving large system sizes and long times, is severely hindered since the simulation of the quantum trajectories associated to Eq. ( 7) requires the knowledge of the whole many-body wavefunction.Due to the exponential scaling of the Hilbert space with the system size L such simulations become rapidly infeasible.For example, in Ref. [26], the quantum XX spin-chain Hamiltonian (which maps to free fermions via Jordan-Wigner transformation) subject to binary annihilation 2A → ∅ or coagulation 2A → A has been studied.Therein the diffusion-limited regime Ω = Γ = 1 is considered, the initial state is the fully filled product state |• • • • • •⟩ and a maximal system size of L = 22 is taken.Under such conditions, a power-law decay n(t) ∼ t −b , with 1/2 < b < 1 is obtained for annihilation.The dynamics is therefore slower than the mean-field prediction (15) but faster than the classical diffusion-limited analogue (17).This might be related to the faster ballistic spreading of quantum particles, compared to classical diffusion, albeit the extrapolation of these results to the thermodynamic limit is difficult to assess.No analytical prediction for the quantum-diffusion limited decay exponent is far present.A field-theory description for 2A → ∅ has been proposed in Ref. [41] using the Keldysh field-theory representation of the quantum master equation.A systematic renormalization group scaling analysis of this action is, however, still missing.
The quantum reaction-limited regime (Γ/Ω ≪ 1) is in this perspective unique since it allows exact analytical calculations in the thermodynamic limit and at long times.This regime has been recently studied in a number of works for fermionic and spin systems [34][35][36][38][39][40], where occupation restrictions are therefore present.It has been therein shown that for fermionic particles the binary annihilation (2A → ∅) decays as n(t) ∼ t −1/2 for quantum coherent (in real space) initial states.The dynamics is therefore not always well-described by the mean-field approach despite the system being in the reaction-limited regime.Quantum effects due to coherences therefore allow for novel form of emergent behavior, which do not have a classical counterpart.The interacting Bose gas, modelled by the Lieb-Liniger Hamiltonian, subject to m-body annihilation, mA → ∅, has been considered in Ref. [37].In the noninteracting limit, mean-field decay (15) is found for every m.In the presence of nonzero interaction, the numerical evaluation of the differential equation for the density turns out to be cumbersome and therefore no estimate for the density decay exponent is therein made.The analysis of all the works [34][35][36][37][38][39][40] is based on the analytical approach of the time-dependent Gibbs ensemble (TGGE) method [82][83][84][85].We briefly recall here this method as it will be used for all the results presented in Sec.IV.
The reaction-limited regime Γ/Ω ≪ 1 corresponds to weak dissipation regime of the Lindblad dynamics.In this limit, the reaction time ∼ 1/Γ is much larger than the hopping time ∼ 1/Ω.Because of this separation of time scales, the state of the system ρ(t) quickly relaxes to a stationary state ρ SS (t) of the Hamiltonian [H, ρ SS (t)] = 0.The stable quasi-particles characterizing the integrable Hamiltonian H (8) (it is trivially integrable since it is quadratic) can still be defined but they acquire a finite lifetime ∼ 1/Γ.This causes the stationary state of the system becoming time-dependent ρ SS (t) and changing over the long time scale 1/Γ.The TGGE method makes an ansatz for ρ SS (t) in the form of a GGE [95,96].This ansatz is motivated by the idea of local generalized thermalization of the Hamiltonian dynamics in the limit of slow reactions.In this limit, the density matrix quickly relaxes towards a maximal-entropy generalized Gibbs state, which keeps into account all the conserved charges of the Hamiltonian.In the case of the Hamiltonian (8), all the conserved charges are linear combination of the occupation number nk in momentum space [95,96].For the Hamiltonian (8), the time-dependent GGE (TGGE) then reads as [82][83][84][85] with is the quasi-momentum and nk = b † k bk the number operator in Fourier space with bk , b † k the bosonic Fouriertransformed operators (see Appendix B).The GGE state describes averages ⟨O⟩ GGE of local observables O in the thermodynamic limit.The time-dependence of the TGGE is contained in the Lagrange multipliers (or generalised inverse temperatures) λ k → λ k (τ ).This implies that conserved charges can still be defined, but they slowly (on a time scale ∼ 1/Γ) drift in time as consequence of the dissipation.Since [n k , H] = 0 for all k, [ρ GGE (t), H] = 0 and the master equation ( 7) takes the form: The state in Eq. ( 20) is diagonal in momentum space (because of translation invariance) and Gaussian.The entire dynamics ensuing from (20) is consequently encoded in the bosonic occupation function in momentum space B q (t) ≡ ⟨ b † q bk ⟩ GGE (t) = δ q,k /(exp(λ q (t)) − 1).The time evolution for the Lagrange parameters λ q (t) is therefore in one-to-one correspondence with that of B q (t).The dynamical differential equation for the latter immediately follows from Eq. ( 21) exploiting that [ρ GGE , nq ] = 0 [34][35][36]38]: In the Appendices B-D, we report the expression of the jump operators L ν j in Fourier space and we detail all the calculations based on Eq. (22).
It is also clear from the previous equation that B q (t) = B q (τ = Γt), i.e., the bosonic occupation function is a function of the rescaled time τ = Γt according to the reaction rate only.This is consistent with the reaction-limited description of the classical RD dynamics ( 15) and (16) where time obeys the very same rescaling.As a matter of fact, the TGGE ansatz of Eqs. ( 20) and ( 22) is valid in the scaling limit t → ∞, Γ → 0 with τ = Γt fixed, as shown in Refs.[82][83][84][85] We also mention that in the present case B q (τ ) does not depend on the space coordinate x since we only deal in this manuscript with homogeneous Lindbladians evolving from homogeneous initial states.This is also the reason why the hopping rate Ω does not appear in the equation (22).In homogeneous systems, indeed, there is no transport of particles.Consequently, in the limit Γ/Ω ≪ 1, the Hamiltonian contribution is effectively integrated out assuming local relaxation to the GGE state (20).Whenever spatial inhomogeneities are included, the Hamiltonian additionally contributes to Eq. ( 21) via a convective term describing ballistic transport of particles.Equation (21) takes in this case the form of a Boltzmann equation [35,37,40,41].
The right hand side of Eq. ( 22) can be exactly computed from the Wick's theorem since it amounts to computing higher-point bosonic correlation functions over the Gaussian state (20).This allows to derive exact rate equations for B q (τ ).Within this perspective, the TGGE method presented here bears similarities with the Hartree-Fock decoupling of four (and higher) body terms in the dynamical equation for the two-point function ⟨ b † n bm ⟩ (with m, n lattice-site indices).This decoupling is often used in dissipative many-body systems, see, e.g., Ref. [97], in order to truncate the infinite hierarchy of equations for the correlation functions to a closed and calculable form.We, however, remark that the Hartree-Fock decoupling is an uncontrolled approximation, while the validity of the TGGE ansatz is controlled by the parameter Γ → 0 (reaction-limited regime) according to the scaling limit τ = Γt fixed explained above.Furthermore, the similarity between the two methods applies only in the present case of the noninteracting Bose gas with quadratic Hamiltonian (8).As discussed in Ref. [37], for the interacting Bose gas described by the Lieb-Liniger Hamiltonian, the GGE is not Gaussian and higher point function in the GGE state cannot be reduced to the two-point function via Wick theorem.In the next Section, we specialize (22) for the various reaction processes introduced before in Eqs. ( 10)- (14).

IV. RESULTS
In this Section, we present our results for the quantum reaction-limited RD dynamics of bosonic particles.In Subsec.IV A, we first discuss the distance-selective annihilation process (10) (θ = 0).In Subsec.IV B, we move to the discussion of annihilation processes with interferences between two decay channels, (10) with θ ̸ = 0, π/2.In Subsec.IV C, we discuss the second-order annihilation process with three interfering decay channels (11).In Subsec.IV D, coagulation decay (12) is studied.In Subsec.IV E, we consider the branching process (14).We then study the corresponding absorbing-state phase transition arising from the competition between branching (14), coagulation (12) and one-body decay (13).The associated stationary phase diagram is reported and discussed.
For all the aforementioned cases, we solve the TGGE rate equation in Eq. ( 22) for the bosonic momentum occupation function B q (τ ).Namely, we consider three different initial conditions: the Bose Einstein condensate (BEC) , Eq. ( 23), the flat mode filling, Eq. ( 24), and the Gaussian distributed mode filling, Eq. ( 25).The implemented initial conditions are the following: • Bose Condensate: • Flat Filling: • Gaussian State: In the BEC, the quasi-momentum q = 0 is macroscopically populated by the total number N of bosons initially in the system.In one spatial dimension, this state corresponds to a quasi-condensate pure state [98,99].In the flat filling, every mode in momentum space is equally occupied by the same number n 0 of bosons.The parameter n 0 is therefore the initial density of particles.In the state (25), the occupation of the modes is Gaussian distributed.The density is computed from the occupation function B q (τ ) as the latter equality being valid as L → ∞.Henceforth, we set N = L in Eq. ( 23), and n 0 = 1 in Eq. ( 24), so that for all the three initial conditions we have an initial density n 0 = 11 .The BEC state (23) has quantum coherences in real space since the associated density matrix has non-zero offdiagonal matrix elements in the bosonic Fock-space basis.The latter being spanned by states of the form |{n j }⟩ = j∈Λ bnj j |0⟩, with Λ denoting an arbitrary set of lattice sites.Similarly, the Gaussian state (25) corresponds to an initial state of the GGE form (20) identified by the occupation function B GS,q inhomogeneous in momentum space.The latter is associated to a bosonic two-point correlation function in real space ⟨ b † n bm ⟩ GGE not diagonal.The latter fully characterizes the Gaussian initial state, which is therefore also in this case quantum coherent in the real space basis.On the contrary, for the flat initial occupation in momentum space (24), the two-point bosonic correlation function is diagonal.The associated initial density matrix ρ 0 = exp(−λN )/Z is diagonal in the classical basis of the Fock space introduced above.For these reasons, we will refer henceforth to initial states (23) and (25) as quantum coherent, in contrast with (24) which is incoherent.
The asymptotic decay exponent ⟨n⟩ GGE (τ ) ∼ τ −δ is evaluated by computing the effective exponent [3]: Here, b > 0 is a scaling parameter.In all the calculations, we set b = 1.5.In the case of a power-law decay, the effective decay exponent δ eff (τ ) approaches asymptotically in time the exponent δ.We plot in the next Subsections δ eff as a function of τ for the three different initial states ( 23)-( 25) for the various reactions considered ( 10)- (12).

A. Distance-selective loss
In this Subsection, we consider the binary annihilation jump operator in Eq. ( 10) without interferences (that is θ = 0).We consider generic values of the distance d between the two bosons.The rate equation Eq. ( 22) then yields (see Appendix B) with τ = Γ α t the re-scaled time.The results for different values of the distance d are shown in Fig. 2. In particular, we see that for d = 0 Eq. ( 28) does not couple q with other momenta k ̸ = q and it therefore gives a closed equation for the density ⟨n⟩ GGE (t): This equation coincides with the law of mass action (15) up to a factor 4 (instead of 2).This factor only affects the amplitude of the asymptotic decay, but not the exponent ⟨n⟩ GGE (τ ) ∼ τ −1 .This result is consistent with that of Ref. [37] for onsite binary losses mA → ∅ and m = 2.For m > 2, in Ref. [37], the asymptotic mean-field exponent n(τ ) ∼ τ −1/(m−1) is also recovered.
For d ̸ = 0, Eq. ( 29) still holds true for the BEC initial state (23).In this case, indeed, binary losses do not populate momenta q ̸ = 0, so that the momentum distribution is concentrated on q = 0 B q (τ ) = N (τ )δ q,0 at any time τ .In this case, Eq. ( 28) reduces to (29) for any d value.For the flat filling initial state Eq. ( 24), instead, the occupation function B q (τ ) = n(τ ) remains flat in q at any time τ and therefore (28) exactly reduces, at any time τ , for d ̸ = 0 to the law of mass action This result is analogous to the one valid for fermions starting from incoherent initial states identified by a flat  29), while for d = 1, 20 Eq. ( 30) is recovered.In all the cases, the asymptotic exponent (right plot) follows the mean-field prediction ⟨n⟩ GGE (τ ) ∼ τ −1 .b) Dynamics from the initial Gaussian occupation function (25) for the same values of d as above.For d = 0, the density dynamics again follows (29), while for d = 1, 20, Eq. ( 30) is asymptotically obtained.Also in this case the asymptotic decay exponent follows from mean field, as the convergence of δ eff (τ ) → 1 shows (right plot).
in q occupation function [38][39][40].Also in those cases the fermionic quantum reaction-limited RD dynamics reduces to its classical counterpart.Meanwhile, for the initial Gaussian occupation function, studied in Refs.(25), Eq. ( 28) reduces to Eq. ( 30) only asymptotically for long τ values.The time needed to asymptotically approach (30) depends on the distance d > 0, the larger the latter, the faster the approach to the law of mass action dynamical behavior.This behavior holds generically for any initial state identified by an occupation function B q (0) inhomogeneous in momentum q.Therefore, all the initial states considered show the same mean-field algebraic decay: This applies to both quantum coherent (23), (25) and incoherent (24) initial conditions.This is a fundamental difference from the RD dynamics of fermions.For fermionic systems, indeed, nearest neighbour binary annihilation (d = 1 with the notation of this manuscript) and quantum coherent initial states give rise to the decay law ⟨n⟩ GGE (τ ) ∼ τ −1/2 [38], which deviates from the mean-field result.In the next Subsection, we study whether such bosonic mean-field decay is robust against the introduction of interference effects (θ ̸ = 0, π/2) in the annihilation reaction channels.

B. Annihilation reaction with interferences
The TGGE rate equation (22) for the jump operators (10), i.e. generic values of θ, reads as (see Appendix B for the details of the calculations) with τ = tΓ α the re-scaled time and This differential equation is similar to the corresponding one for fermions of Ref. [38] (cf.Eq. (S10) therein), the only difference being in the opposite sign of the terms cos(d(k − q)) and cos(d(k + q)).In Fig. 3, we plot the solution of the differential Eq. ( 32) for the density ⟨n⟩ GGE (τ ) as a function of τ .We find that the long-time dynamics at θ = π 4 is qualitatively different than that obtained for all other values, i.e., θ ̸ = π/4.Namely, for all the considered initial states, we find that for θ ̸ = π/4, the density decay exponent is gain the mean-field one ⟨n⟩ GGE (τ ) ∼ τ −1 .
For θ = π/4, however, the density is constant in time ⟨n⟩ GGE (τ ) = ⟨n⟩ GGE (0).The BEC state is, indeed, a dark state of the jump operator in Eq. ( 10) at θ = π 4 , which amounts to saying that |BEC⟩ is annihilated by all the annihilation (10) jump operators Moreover, the BEC is an eigenstate (the ground state in this case) of the Hamiltonian (8).This observation is consistent with the results of Refs.[46,49], where the |BEC⟩ state was similarly identified as an exact dark state of the quantum master equation of the noninteracting Bose gas.However, in Refs.[46,49], the jump operators do not involve particle losses, but still they annihilate the antisymmetric part of the wavefunction due to terms bi − bi+d , thus rendering the BEC state dark.For the flat filling (24) and the Gaussian occupation states (25), instead, we find a slower decay ⟨n⟩ GGE (τ ) ∼ τ −1/2 compared to mean field, as shown in Fig. 3.The asymptotic 1/2 decay exponent (and the amplitude) can be exactly computed from the asymptotic of Eq. ( 32) (see again Appendix B for the details of the calculations).In particular, Eq. ( 32) admits the following implicit solution for the mode occupation function B q (τ ) Here, for the sake of brevity, we denoted ⟨n⟩ GGE (τ ) = n(τ ).The previous equation shows that q = 0 is the slowest decaying mode.Furthermore, for B q (0) = n(0)δ(q) (the BEC occupation in the limit L → ∞), Eq. (35) gives n(τ ) = n(0) consistently with the previous discussion concerning the dark-state condition for the BEC.Using Eq. ( 35) into Eq.( 26), the q integral for the density n(τ ) can be asymptotically computed using the saddle-point approximation (which is justified since τ 0 dt n(t) → ∞ as τ → ∞).The density asymptotic one obtains with the saddle-points q * n = nπ/d, n = 1, 2 . . .d − 1, determined by the d − 1 zeros of the sin 2 (dq) function in the interval q ∈ (0, π).Therefore, when θ = π/4 the decay exponent, 1/2, does not depend on the distance d between the pair of bosons involved in the loss.The amplitude of the decay, instead, does generically depend on d.Only in the case of the flat-filling initial occupation (24), the amplitude of the decay takes the d-independent value 1/2 n 0 /π.
To sum up, for the decay exponent of the particle density in the presence of binary annihilation with interference between two decay channels (10), we find both for incoherent initial states (24) and for coherent ones, e.g., with a Gaussian initial occupation function B q (0) (25).The same behavior (37) applies generically for initial quantum coherent states identified by a non-flat occupation function B q (0) of generic form.The |BEC⟩ state (23) for θ = π/4 does not decay since it is dark with respect to the dissipation.The result (37) significantly deviates from the one valid for fermions discussed in Ref. [38].In the case of fermions, the decay ⟨n⟩ GGE (τ ) ∼ τ −1/2 is valid for any value of θ ̸ = 0, π/2, not just θ = π/4.The different behavior between fermions and bosons subject to binary annihilation with interference, and at the same time, the emergence of the beyond mean-field decay (37) can be understood by looking at the space continuum limit of Eq. (10).In the continuum limit, the lattice spacing a (so far set to 1) is reintroduced, so that lattice point are identified as x j = ja and ℓ = La is the dimensionful length of the system.The continuum limit is obtained by sending a → 0, L → ∞ with ℓ fixed, as we detail in Appendix A. The jump operator (10) L α ja reduces, in this limit, to onsite pair annihilation L α ja → L α (x) ∼ b2 (x) as long as θ ̸ = π/4, with b(x) = bja / √ a being the continuum bosonic destruction field operator.This operator yields mean-field decay (31), following the same steps as those done in Sec.IV A for the lattice case.Only for θ = π/4 (see again Appendix A for the details), the binary annihilation operator (10) in the continuum limit takes a different This effect eventually causes the decay law (37), which is beyond mean-field.In the fermionic case, the leading term L α (x) ∼ ĉ2 (x) of the continuum limit expansion, ĉ(x) being the fermionic field destruction operator, is always zero because of the fermionic statistics.The continuum limit of fermionic nearest-neighbour annihilation is therefore always of the form L α (x) ∼ ĉ(x)∂ x ĉ(x) and the decay ⟨n⟩ GGE (τ ) ∼ τ −1/2 applies for any θ ̸ = 0, π/2.
In order to further understand the emergence of nonmean-field decay exponents due to interfering decay channels, we proceed by discussing in the following Subsection annihilation reactions (11) with interferences among three different decay channels.

C. Second-order annihilation
The TGGE rate equation (22) can be specialized to the jump operator (11) following similar steps as those performed in Secs.IV A and IV B. We report the main steps in Appendix B, while here we give the final result: with the re-scaled time τ = Γ ᾱt.The results of the numerical solution of the previous equation are reported in Fig. 4. The BEC initial state ( 23) is also dark, as in Eq. ( 34), with respect to the jump operator L ᾱ j (for any j) and therefore it is a many-body dark state of the quantum master equation.The density ⟨n⟩ GGE (τ ) = ⟨n⟩ GGE (0) = n 0 consequently remains constant in time.38).In the right plot (log scale on the horizontal axis only), the effective exponent δ eff (τ ) as a function of τ (Eq.( 27)) is shown.Both the initial state (24) with a flat in q occupation function, and ( 25) with an initial occupation function of Gaussian form are studied.In both cases, the same algebraic decay ⟨n⟩ GGE (τ ) ∼ τ −0.28 is observed.The exponent of the algebraic decay is quantified by plotting the effective exponent, which at long times converges to δ eff (τ ) ≃ 0.28 for both initial states (24) and (25).
For both the flat filling (24) and the Gaussian occupation function (25) initial states, we, instead, observe algebraic decay asymptotically in time with a non-meanfield exponent This decay exponent is numerically computed by plotting the effective exponent δ eff (τ ) as a function of τ , as shown in the right panel of Fig. 4. The effective exponent (27) converges at long times to the value δ eff (τ ) ≃ 0.28 for both the initial states considered.The choice of the initial state solely affects how fast the asymptotic value of δ eff is met, but not the value itself.We also verified that the result in Eq. ( 39) holds generically for initial states identified by a different q-dependent initial occupation function B q (0).
In the case of Eq. ( 39), we remark, however, that it is not possible to derive analytically the asymptotic exponent, as we did in Eqs. ( 35) and ( 36) of Subsec.IV B, due to the more complicated structure of Eq. ( 38).
The result ( 39) is a non-mean-field prediction for the algebraic decay.Following the same reasoning used at the end of the previous Subsec.IV B, we can understand this behavior by looking at the continuum limit of the annihilation process (11) coupling three different decay channels.In this case, the continuum limit (see again Appendix A for the details) leads to a jump operator of the form L ᾱ(x) ∼ b(x)∂ 2 x b(x).The jump operator L ᾱ(x) couples now to second-order spatial derivatives (this is why we also name the process as "second-order annihilation") and it induces spatial fluctuations over larger regions and therefore deviations from the mean field are more sensible.In particular, the exponent decreases compared to Eq. ( 36) for the first-order interference drifting further from the mean-field value (δ = 1).
Apparently, interference effects generally slow down the decay of the density compared to the classical reactionlimited case.In both cases analysed in Eqs.(36) and (39), we see that quantum interferences are necessary in order to observe beyond mean-field decay in the noninteracting Bose gas.

D. Coagulation
The onsite coagulation decay is modelled by the jump operator (12) introduced in Sec.II.It describes the reaction 2A → A where two (or more) bosons meet on the same site j leading to the destruction of one of the particles at the same lattice site.The resulting rate equation (cf.Appendix C for the details of the calculations) for the bosonic occupation function is with τ = Γ γ t the re-scaled time according to the coagulation rate.We note that for one-site coagulation, similarly to the case of onsite binary annihilation (29), the TGGE equation for the density is closed, and it reads This equation contains both a two-body term, ⟨n⟩ 2 GGE , as in the corresponding classical law of mass action (15), and a three-body term, ⟨n⟩ 3 GGE .The latter is notably absent in the classical law of mass action description.It is also absent in the fermionic quantum reaction-limited dynamics [38] and therefore a genuine property of the quantum bosonic reaction-limited dynamics.In the presence of coagulation only (and more generally in the presence of reactions depleting the system) this three-body term does not, however, impact the late time asymptotics.The density ⟨n⟩ GGE , as a matter of fact, decays asymptotically in time to zero and therefore the three-body term ⟨n⟩ 3 GGE can be neglected with respect to the two-body term ⟨n⟩ 2 GGE .The latter immediately gives the mean-field algebraic decay In Fig. 5, we report the numerical solution of Eq. ( 40) together with the plot of the associated effective exponent δ eff (τ ) versus the re-scaled time τ .The effective exponent converges to the mean-field exponent δ eff (τ ) ≃ 1 at long times.This mean-field behaviour is, importantly, valid for any initial state ( 23)-( 25), quantum coherent or not.
A similar result has been derived in Ref. [38] for fermionic systems, where the coagulation decay similarly shows ⟨n⟩ GGE (τ = Γ γ t) = τ −1 decay both for coherent and incoherent initial states.In the fermionic case, for a Fermisea initial state, quantum coherences are only observed in Ref. [38] to rescale time by a factor dependent on the initial density n 0 , leaving the asymptotic decay exponent unchanged.In the bosonic case, instead, Eq. ( 41) gives the very same dynamics for all kind of initial conditions.We also note that the coefficient 2 in front of the twobody term on the right-hand side of Eq. ( 41) is half of the coefficient 4 of the onsite annihilation two-body term (29).This implies that the density ⟨n⟩ coag/ann GGE (τ, n 0 ) (n 0 being the initial density) as a function of time in the annihilation and coagulation processes are asymptotically (as long as the ⟨n⟩ 3 GGE term for coagulation can be neglected) related by In the context of classical RD [90][91][92][93][94] this relation is taken as the hallmark that signals that both coagulation and annihilation dynamics belong to the same universality class.In particular, this implies that the two processes share the same asymptotic decay.Apparently in the bosonic quantum RD system the two processes still belong to the same universality class, at least in the reactionlimited regime.This statement applies to both quantum coherent and incoherent initial states, in contrast with the fermionic case of Ref. [38], where the equivalence (43) applies only for incoherent initial states.24), but the very same curve is obtained for the BEC (23) and the Gaussian occupation state (25).In all the cases, the density decays as in the mean-field approximation ⟨n⟩ GGE ∼ τ −1 .In the right plot, the effective exponent δ eff (τ ) (27) is plotted (log scale on the horizontal axis only) and it converges to mean-field value 1 at long times.

E. Absorbing-state phase transition in quantum bosonic RD
In order to obtain a process with a stationary state featuring a non-zero density of particles, we need to include reactions increasing the particle number, thereby competing with the loss processes ( 10)- (12).Therefore, we consider onsite branching reactions A → A + A (14) at rate Γ β .For the branching reaction process, we found the following rate equation (see Appendix D for the calculation): In the following, we denote ⟨n⟩ GGE (t) = n, dropping the time dependence for the sake of brevity.As in the case of onsite coagulation (41), the particle density obeys the closed equation We note that, as well as for coagulation, branching produces the three-body term, n 3 , which is absent both in the classical and in the fermionic reaction-limited RD dynamics.The branching process leads to the creation of particles and therefore the density increases in time in Eq. ( 45).This persistent growth can be countered by introducing annihilation reactions, and in the following we will study an absorbing-state phase transition that results from such competing processes.We consequently consider a model where (onsite) branching (14), at rate Γ β = Γβ, is competing with (onsite) coagulation (12), at rate Γ γ = Γγ, and one-body decay (13), at rate Γ δ = Γδ.This parametrization of the rates allows to identify Γ as the overall dissipation (inverse) time scale, while β, γ, δ encode the relative strength of the three considered processes.Considering Eq. ( 45) together with (41) one then finds the rate equation where τ = Γt.Note, that the linear term, −δn, stems from the one-body decay process.We already note that Eq. ( 46) is qualitatively different from the classical mean-field description (18) of the unrestricted contact process due to the presence of cubic and quartic terms in the density.However, Eq. ( 46) is still of mean-field type since it implies that the GGE is factorized also in real space: ρ GGE (t) ∼ exp(−β(t) N ).In this sense, for bosons, the quantum mean-field description is different and richer than the classical mean-field one.Eq. ( 46) can be equivalently written as where has the meaning of a potential.We remark that a similar fourth order potential W (n) was found in Refs.[57,58], which studied the so-called quantum contact process.Microscopically this process is, however, different to the one discussed here: it is formulated in terms of spin 1/2 particles (not bosons) and the Hamiltonian embodies coherent branching and coagulation, i.e., spin flips conditional on the excitation of a neighbouring spin.In our case, instead, the Hamiltonian (8) leads to free, unconstrained, hopping of particles.The dissipation is described by classical incoherent processes in both cases.Despite their apparent differences both models yield the same potential W (n).
It is noteworthy that in Refs.[57,58] the potential W (n) describes the stationary phases of the quantum contact process only in the absence of spatial fluctuations, i.e., within the mean-field approximation.The latter is valid only in high spatial dimensions, while it fails in one dimension.In the present work, however, W (n) emerges in one dimension in the reaction limited regime Γ/Ω ≪ 1 as an exact result.
The stationary-state phase diagram of the process (47) is reported in Fig. 6, where we plot the stationary state density n SS (represent by the color code) as a function of the ratios β/γ and δ/γ.The stationary density n SS is one of the (possibly multiple) minima of W (n). The number of minima of W (n) and their position depends on β/γ and δ/γ.For β/γ > 1, the potential W (n) is unbounded from below and the active phase supports an infinite density.For this reason, we take β/γ < 1, in the following.In the Fig. 6, we also sketch the potential W (n) so that the stationary state density n SS value in the various regions of the phase diagram can be readily understood.We can identify three different phases a) − c): a) If δ/γ < 2β/γ, an active, i.e. finite density solution, n act ̸ = 0, of Eq. ( 46) is present: (49) The density n act is the only minimum of W (n), while n abs = 0 is a maximum.In this case, the system is therefore in the active phase n SS = n act .b) For β/γ > 1/4 we define the curve (δ/γ) c as For δ/γ < (δ/γ) c and δ/γ > 2β/γ, the potential W (n) has two minima n abs = 0 and n act (49).For this choice of parameters the stationary state n SS can be either n abs = 0 or n act depending on the value of the initial density n(τ = 0) = n 0 .If n 0 is small enough so as it falls in the basin of attraction of n abs = 0, then n SS = 0. On the contrary, if n 0 is large enough so as to fall in the basin of attraction of n act , then the system is active n SS = n act .For this reason, we identify this region of the phase diagram as a bistable phase.
For β/γ < 1/4, the boundary between a) and c) is given by the line β/δ = 1/2, which is shown in red solid in Fig. 6.Along this line the order parameter n SS of the transition varies continuously.This is shown in the bottom-right panel of Fig. 6, where the stationary density n SS is displayed as a function β/γ for δ/γ = 0.4 (horizontal red-dotted line in the density plot in Fig. 6).In particular, expanding Eq. ( 49) for β − 2δ → 0 and β < γ/4 one obtains At the same time, setting β/δ = 1/2 and β/γ < 1/4, one finds that the density n(t) decays asymptotically as Equations ( 51) and ( 52) identify the critical exponents of the classical mean-field directed percolation universality class [1-3, 5, 6], as also recalled after Eq. (19).Along the second-order transition line, the cubic nonlinearity n 3 of ( 46) is therefore negligible and the transition is driven by the second-order n 2 nonlinearity.This renders (19) effectively analogous to the classical mean-field equation (19).The classical mean-field directed percolation exponents are therefore obtained.For β/γ > 1/4, the boundary between b) and c) is given by the curve (δ/γ) c in Eq. ( 50) and it reveals, instead, a richer physics compared to the mean-field description (18) of the classical contact process.In particular, along the line (δ/γ) c given by (50), plotted in solid green in Fig. 6, we find a first-order transition since the density jumps discontinuously from the value (49) to zero.This is shown in the middle panel in Fig. 6.Therein we plot the stationary density n SS as a function of β/γ for fixed δ/γ = 0.6 (horizontal green-dotted line).
Equation (50), for the first-order line, is obtained by calculating when the potential W (n) has an inflection point d 2 W (n)/dn 2 = 0.For δ/γ < (δ/γ) c and β/γ > 1/4, the inflection point morphs into a first minimum n act (49), while a second minimum is given by n abs = 0. We emphasize that the relative height of the two minima W (n act )/W (n abs ) does not fix the stationary density since in this nonequilibrium case there is no principle of free-energy minimization as in equilibrium systems.In the latter case, indeed, W (n) would be interpreted as a free energy and the equilibrium solution would necessarily be the absolute minimum of W (n) in the thermodynamic limit.In the present nonequilibrium scenario, instead, the stationary density can take either the value n SS = n act or the value n SS = n abs = 0 depending on the initial density n 0 .This phenomenon is referred to as bistability and it happens in proximity of the first-order branch (δ/γ) c , where the absorbing and the active phase coexist due to the presence of two minima in W (n). The bistable phase (plotted in grey blurred scale in Fig. 6) is delimited by the coexistence curves (δ/γ) c (green-solid line) and δ = 2β (red-dashed line).Coexistence of difference phases in the presence of a first-order phase transition has been also observed in a different context, both in equilibrium and in nonequilibrium, for Dicke-Bose-Hubbard systems [100,101].
The two coexistence lines meet the second-order line at the bicritical point (β/γ, δ/γ) = (1/4, 1/2) (represented in orange in Fig. 6).At this point, the transition is continuous, but it belongs to a different universality class than that of Eqs. ( 51) and (52).Namely, for β/γ = 1/4 and 2β − δ → 0, we find from (49) the power law for the stationary density as a function of β/δ.This is shown in the central-right panel of Fig. 6, where n SS is plotted as a function of β/γ for fixed δ/γ = 0.5 (horizontal orange-dashed line in the density plot in Fig. 6).The jump discontinuity in n SS across the first-order branch goes continuously to zero as the bicritical point is approached.Similarly, at β/δ = 1/2 and β = γ = 1/4, one has that the density n(t) decays algebraically in time as The same critical exponents as in Eqs. ( 53) and ( 54) have been found also in Refs.[57,58] for the mean-field regime of the quantum contact process involving both classicalincoherent and quantum-coherent branching/coagulation.Here, a similar bicritical point to that of Fig. 6 was found in the plane spanned by the classical and quantum branching rates (note the difference with Fig. 6, where both the axes β/γ and δ/γ are classical rates).The bicritical point is argued in Refs.[57,58] to be described by the tricritical directed percolation universality class [87][88][89].For spatial dimension D ≥ 3, the mean-field exponents (53)  In this table, we summarise our results from Eq. ( 20) for the different annihilation reaction types ( 10)-( 12) treated in Subsecs.IV A-IV D in the reaction-limited regime Γ/Ω ≪ 1.
In all the listed cases, the density decays as a power-law ⟨n⟩ GGE (τ ) = τ −δ with the listed exponent δ.The exponent is given for each reaction and for the considered initial states.The flat filling ( 24) is incoherent, while the BEC ( 23) and the Gaussian state ( 25) are quantum coherent.In the case of the BEC, the value 0 of the exponent, on the second and third line, refers to the fact that the BEC is dark to these jump operators and therefore it does not decay.

V. DISCUSSION
In this manuscript, we have studied the quantum reaction-limited reaction-diffusion dynamics of the noninteracting Bose gas.The main finding of the manuscript is that bosonic statistics considerably impacts on the universal properties of the quantum RD dynamics.In particular, we find that the combination of bosonic statistics with quantum effects, such as coherent hopping and quantum interferences, leads to universal nonequilibrium behavior different from both fermionic and classical RD systems.
The dynamics for QRD systems is formulated in terms of the quantum master equation ( 7)- (9).Within this formulation, we considered various types of reactions introduced in Sec.II.Our analytical study is based on the time-dependent generalized Gibbs ensemble method, briefly recalled in Sec.III.This method gives direct access to the quantum reaction-diffusion dynamics in the thermodynamic limit in the so-called reaction-limited regime Γ/Ω ≪ 1 of weak dissipation.Our results for the density algebraic decay ⟨n⟩ GGE (τ ) ∼ τ −δ in the presence of only reactions depleting the system are summarized in Table I.We see that the law of mass action (16) decay exponent is generically recovered for all classical-incoherent reactions, such as distance selective losses (10) (θ ̸ = 0) and coagulation.This result applies both to quantum coherent and incoherent initial states.Algebraic decay beyond mean field is possible for bosons only when interference between different decay channels is allowed, as shown in Figs. 3  and 4. In the case of first-order interference (10) with θ = π/4, in Subsec.IV B, we find, namely, a power-law decay with exponent 1/2 (cf.Eq. ( 36)).For all the other values θ ̸ = π/4, mean-field decay is observed.For the second-order interference (11) in Subsec.IV C, the decay exponent is, instead, approximately 0.28.Incoherent coagulation (12) decay also yields mean-field decay (42).These results show markedly different behavior between bosons and fermions.In the latter case, power-law decay beyond mean field is, indeed, attained under broader conditions, such as also for classical annihilation reactions from incoherent initial states, and for generic interference decay θ ̸ = 0, π/2 [38][39][40].Furthermore, for bosons, coagulation and binary annihilation belong to the same universality class, at least as long as the reaction-limited regime is concerned.For fermions, instead, this is not true [38].
In table II, we summarize the main critical exponents associated to the absorbing-state phase transition observed in bosonic systems in the presence of branching (14), coagulation and decay.contact process (β, γ, δ) This transition can be either of first or second order.In the table, we summarize the exponents (ξ, χ) associated to the latter.In the second column, the exponents associated to the second-order line β = δ/2 are reported.These exponents are those of the mean-field directed percolation universality class.We further give (third column) the exponents associated with the bicritical point of Fig. 6.These exponents are those of the mean-field tricritical directed percolation universality class.
The corresponding stationary-state phase diagram in Fig. 6 is of mean-field nature.However, it is qualitatively different from both the classical mean-field (18) and the fermionic analogue [38].In the bosonic case, we find that both a first-order and a second-order line are present.The first-order transition line has no classical counterpart.The second-order line is again characterized by the critical exponents of the mean-field directed percolation universality class (middle column of Table II) as in the classical reaction-limited description (18).Interestingly, the second-order line and the first-order one meet at a bicritical point.This point is identified by the mean-field exponents of the tricritical directed percolation universality class (rightmost column of Table II).The universality class of the bicritical point is consequently different from directed percolation, which characterizes the classical dynamics.For bosons, the quantum meanfield description of the contact process phase diagram is therefore richer than the classical mean-field description in Eq. (18).
Our work paves the way for many future research directions concerning the Bose gas subject to dissipative loss processes.It would be interesting to study the quantum RD dynamics in the presence of reactions conserving the particle number as in Refs.[46,49].Therein, dissipation allows for the dynamical preparation of a pure |BEC⟩ state.This state is a many-body dark state of the dynamics and it is therefore attained at long times.It is important and experimentally relevant to study the stability of this dynamical preparation protocol of the BEC state against particle losses.
It is also interesting to look at different observables than the particle density.The two-point bosonic correlation function ⟨ b † n bm ⟩ GGE (τ ) can be readily computed within the TGGE method and its temporal decay would shed light on the decay of quantum coherences as the stationary state is approached.In addition, it is interesting to study Bose superfluids with weak dissipative losses.These models can be also addressed with the TGGE ansatz of this manuscript.We expect that for Bose superfluids the dynamics of the superfluid order parameter gives additional information on the decoherence in time of the system due to dissipation, as it has been shown in Ref. [33] for lossy fermionic superfluids.
Ultimately, it is desirable to go beyond the quantum reaction-limited regime Γ/Ω ≪ 1 discussed in this manuscript.The quantum analogue of the diffusionlimited decay (2) is, indeed, not known.Quantifying the anticipated quantum diffusion-limited algebraic decay might be possible by exploiting the field theory formulation of the quantum master equation via the Keldysh path integral method [102,103].In the diffusion-limited regime, the significance of spatial fluctuations due to the finite hopping is inherently increased.A systematic renormalization group scaling analysis is therefore needed in order to quantify the space dependence and the impact of spatial fluctuations on the decay exponents.
With respect to the absorbing-state phase transition, it is important to understand the reason behind the similarity between the phase diagram of Fig. 6 and the one found for the quantum contact process in Refs.[57,58].The model studied therein is different from the quantum RD model here proposed since coherent effects are introduced via a different Hamiltonian giving coherent branching/coagulation.In our case the Hamiltonian (8) simply gives free hopping.The field-theory scaling analysis of the quantum RD model of Subsec.IV E could, in this way, shed light on the, largely not understood, universality class of the quantum contact process away from the reaction-limited, mean-field, limit.
Equations (B8) and (B9) coincide with Eqs.(32) and (33) of the main text.Upon setting θ = 0, Eq. (B9) immediately renders (28).For θ = π/4, the expression in Eq. (B9) simplifies to g π/4,d (k, q) = 2[sin(kd) + sin(qd)] 2 = 2[sin 2 (kd) + sin 2 (qd) + 2 sin(kd) sin(qd)]. (B10) This function is very similar to the one reported in Refs.[34,38] for the case of fermions (setting d = 1 therein).The only difference with respect to the fermionic case lies in the sign change [sin(kd) − sin(qd)] 2 , which in turn causes the last term on the right hand side of the second equality of Eq. (B10) to have the opposite sign −2 sin(kd) sin(qd).This difference, however, is not important as long as one considers initial states whose distribution is even in q: B q (0) = B −q (0).This is true for all the initial states ( 23)- (25) we considered in this work.If the initial state is, indeed, invariant under quasi-momenta reversal q → −q, then this symmetry is kept at all times B q (τ ) = B −q (τ ).This, in turn, implies that the last term on the right hand side the second equality of (B10) is zero since k sin(kd)B k (τ ) = 0 (the sin(kd) function is odd in k).In this case, Eq. ( 32) reduces to the very same form discussed in Refs.[34,38] for fermions: From this equation, we can derive the implicit solution (35) for B q (τ ), following similar derivations as those performed in Refs.[34,35].We report here this derivation for the sake of completeness.The equation for the density ⟨n⟩ GGE (τ ) ≡ n(t) from (B8) reads as dn(τ We remark that (B12) is not a closed equation for the density n(t) since the occupation function B k (τ ) appears explicitly on the right hand side.The evolution for the occupation function B q of the mode q is, indeed, coupled to that of all the other modes B k through the function g π/4,d (k, q).This function arises from the structure of the jump operator (10) involving the superposition of two decay channels between nearest neighbouring sites.In the case of onsite processes, e.g., coagulation and branching discussed in Appendix C and D below, respectively, different Fourier modes are not coupled and one can write closed equations for the density of particles.Substituting Eq. (B12) into (B11) yields dB q (τ ) dτ = B q ṅ(τ ) 2n(τ ) − 2 sin 2 (dq)B q (τ )n(τ ), (B13) with ṅ(τ ) = dn(τ )/dτ .The time integration of Eq. (B13) gives ln (B q (τ )/B q (0)) = 1 2 ln (n(t)/n(0)) − 2 sin 2 (qd) The exponentiation of (B14) eventually yields Eq. ( 35) of the main text.The derivation of Eq. ( 38) for the jump operator L ᾱ j (11) describing second-order annihilation decay follows the very same steps as those leading to (B8).In particular, Eq. (B4) is still valid upon changing the definition of the function f θ,d (k) → f d (k).For the jump operator (11)  The calculation of the six-point functions follows to very same lines of the calculation of Eqs.(C4) and (C5) and we therefore do not report it here for the sake of brevity.The final result for the evolution equation for the occupation function B q in momentum space reads dB q (t) dt = Γ β [2(2n − B q ) + 8n 2 + 6B q n 2 ], (D4) which coincides with Eq. ( 44) of the main text.As in the case of onsite coagulation, for onsite branching the evolution of a mode q is not coupled to the evolution of the other modes k ̸ = q.This allows to derive the closed equations ( 45) and (46) for the density of particles.These equations are akin to the classical law of mass action equation (18), but, at the same time, they are fundamentally different from those equations due to the three body terms.These terms generate the phase diagram discussed in Subsec.IV E of the main text, which shows a richer behavior than the classical mean-field contact process in Eq. ( 18).

Figure 2 .
Figure 2. Distance-selective binary annihilation quantum RD dynamics.(Left) Log-log plot of the density ⟨n⟩ GGE (τ ) as a function of τ = Γαt from the solution of Eq. (28).(Right) Plot of the associated effective exponent δ eff (τ ) (27) as a function of τ (log scale on the horizontal axis only).a) Dynamics from the flat filling initial state (24) for three different values of the distance d = 0, 1, 20.For d = 0, the density follows Eq. (29), while for d = 1, 20 Eq. (30) is recovered.In all the cases, the asymptotic exponent (right plot) follows the mean-field prediction ⟨n⟩ GGE (τ ) ∼ τ −1 .b) Dynamics from the initial Gaussian occupation function(25) for the same values of d as above.For d = 0, the density dynamics again follows(29), while for d = 1, 20, Eq. (30) is asymptotically obtained.Also in this case the asymptotic decay exponent follows from mean field, as the convergence of δ eff (τ ) → 1 shows (right plot).

Figure 4 .
Figure 4.Quantum RD dynamics of second-order annihilation reaction.Log-log plot (left) of the density ⟨n⟩ GGE (τ ) as a function of τ = Γᾱt from Eq. (38).In the right plot (log scale on the horizontal axis only), the effective exponent δ eff (τ ) as a function of τ (Eq.(27)) is shown.Both the initial state(24) with a flat in q occupation function, and (25) with an initial occupation function of Gaussian form are studied.In both cases, the same algebraic decay ⟨n⟩ GGE (τ ) ∼ τ −0.28 is observed.The exponent of the algebraic decay is quantified by plotting the effective exponent, which at long times converges to δ eff (τ ) ≃ 0.28 for both initial states(24) and(25).

Figure 5 .
Figure 5. Coagulation quantum RD dynamics.Log-log plot (left) of the density ⟨n⟩ GGE (τ ) as a function of time τ .The initial condition in the given case is the flat filling state in Eq. (24), but the very same curve is obtained for the BEC(23) and the Gaussian occupation state(25).In all the cases, the density decays as in the mean-field approximation ⟨n⟩ GGE ∼ τ −1 .In the right plot, the effective exponent δ eff (τ )(27) is plotted (log scale on the horizontal axis only) and it converges to mean-field value 1 at long times.

Figure 6 .
Figure 6.Stationary phase diagram of quantum RD with incoherent branching and coagulation.The stationary state density nSS, represented by the color code, is plotted as a function of the ratios β/γ, incoherent branching over incoherent coagulation, and δ/γ, incoherent decay over incoherent coagulation.We identify three different regions a) − c) in the phase diagram.In each region, we sketch the potential W (n), Eq. (48).A second-order continuous phase transition separates the absorbing c) and the active a) phases and it is plotted with the red-solid line β = δ/2.The stationary density nSS vanishes continuously as a function of β/γ, as shown in the bottom-right panel for a fixed value of δ/γ = 0.4.This transition belongs to the mean-field classical directed percolation universality class (51).A first-order discontinuous phase transition takes place across the green-solid line (50).This line separates the absorbing phase c) from a bistable phase b) (shown in grey-blurred scale).The bistable phase is delimited by the coexistence lines (δ/γ)c (green solid) and δ = 2β (red dashed).Across the first-order line, the density discontinuously jumps from zero to the value nSS, as shown in the top-right panel for a fixed value of δ/γ = 0.6 (greendotted horizontal cut).The two coexistence lines meet at the bicritical point (orange dot) (β/γ, δ/γ) = (1/4, 1/2).At this point, the transition is still second order, as shown in the central-right panel, but it belongs to the mean-field tricritical directed percolation universality class (53) and (54).

Table I .
Summary of quantum reaction-limited density decay exponents for the noninteracting Bose gas.
(54)the exponents(53)and(54)are then a consequence of considering the reaction-limited regime Γ/Ω ≪ 1, where spatial fluctuations can be neglected.Consequently, in the reaction-limited regime here considered, such exponents are exact already in one spatial dimension D = 1.

Table II .
Summary of quantum reaction-limited exponents for the absorbing-state phase transition in the noninteracting Bose gas.In this table, we summarise our results for the quantum RD model with absorbing-state phase transition in the reaction-limited regime (see Subsec.IV E).