"Smashing more than two": Deuteron production in relativistic heavy ion collisions via stochastic multi-particle reactions

We study the deuteron production via the deuteron pion and nucleon catalysis reactions, $\pi p n \leftrightarrow \pi d$ and $N p n \leftrightarrow N d$, by employing stochastic multi-particle reactions in the hadronic transport approach SMASH for the first time. This is an improvement compared to previous studies, which introduced an artificial fake resonance $d'$ to simulate these $3 \leftrightarrow 2$ reactions as a chain of $2\leftrightarrow 2$ reactions. The derivation of the stochastic criterion for multi-particle reactions is presented in a comprehensive fashion and its implementation is tested against an analytic expression for the scattering rate and the equilibrating particle yields in box calculations. We then study Au + Au collisions at $\sqrt{s_{\mathrm{NN}}} = 7.7$ GeV, where we find that multi-particle collisions substantially reduce the time required for deuterons to reach partial chemical equilibrium with nucleons. Subsequently, the final yield of $d$ is practically independent from the number of $d$ at particlization, confirming the results of previous studies. The mean transverse momentum and the integrated elliptic flow as a function of centrality are rather insensitive to the exact realization of the $2\leftrightarrow 3$ reactions.

We study the deuteron production via the deuteron pion and nucleon catalysis reactions, πpn ↔ πd and N pn ↔ N d, by employing stochastic multi-particle reactions in the hadronic transport approach SMASH for the first time. This is an improvement compared to previous studies, which introduced an artificial fake resonance d to simulate these 3 ↔ 2 reactions as a chain of 2 ↔ 2 reactions. The derivation of the stochastic criterion for multi-particle reactions is presented in a comprehensive fashion and its implementation is tested against an analytic expression for the scattering rate and the equilibrating particle yields in box calculations. We then study Au + Au collisions at √ sNN = 7.7 GeV, where we find that multi-particle collisions substantially reduce the time required for deuterons to reach partial chemical equilibrium with nucleons. Subsequently, the final yield of d is practically independent from the number of d at particlization, confirming the results of previous studies. The mean transverse momentum and the integrated elliptic flow as a function of centrality are rather insensitive to the exact realization of the 2 ↔ 3 reactions.

I. INTRODUCTION
In transport simulations of heavy ion collisions the most commonly implemented reactions are 2 → 2 elastic and inelastic scattering, 1 → 2 (and sometimes 1 → many) resonance decays and 2 → 1 resonance formations. These reactions are the most likely ones to occur in the dilute limit, where transport approaches are applicable and are easier to implement. However, in certain cases multi-particle reactions, where more than two particles are incoming or outgoing become important. Such cases can be identified by employing the detailed balance principle.
The detailed balance principle states that in an equilibrated system the rate of each elementary process is equal to the rate of the reverse reaction. A sufficient condition for this is the equality of the squared matrix elements for any forward and reverse reaction. This, in turn, follows from time reversal invariance, which is strictly fulfilled for the strong interaction. One consequence of the detailed balance principle is that if a 1 → n decay or 2 → n reaction is frequent, then in a system close to equilibrium the reverse reaction is frequent too. Moreover, the equality of matrix elements allows to compute the rate of the reverse reaction. Practical cases, where multi-particle reactions * Electronic address: staudenmaier@fias.uni-frankfurt.de † Electronic address: oliiny@uw.edu ‡ Electronic address: torres-rincon@itp.uni-frankfurt.de § Electronic address: h.elfner@gsi.de have been demonstrated to play a role are the following: • Deuteron production by N pn ↔ N d reactions in heavy ion collisions at low energies (Nb+Nb collisions at projectile energy of 650 MeV/nucleon in case of [1]).
• Reactions reverse to many-body decays, such as 3π → ω need to be present to fulfill the detailed balance principle [6].
In this work we focus on the three-body deuteron catalysis reactions involving π and N : πpn ↔ πd and N pn ↔ N d. At lower collision energies baryons dominate at mid-rapidity, therefore nucleons are the most frequent catalysts of deuteron production. At energies above √ s NN ≈ 5 GeV pion catalysis, i.e. πpn ↔ πd reactions, start to dominate. In contrast to this work, deuteron observables in heavy ion collisions are usually computed either by final-state coalescence from nucleons or by a thermal approach assuming chemical equilibrium of deuterons with hadrons, see [7] for an overview of the available models. Microscopic deuteron production by pion catalysis is a recent idea introduced in [8] for central Pb+Pb collisions at LHC energies and further tested for non-central collisions and at lower energies down to √ s NN = 7.7 GeV [9,10].

arXiv:2106.14287v3 [hep-ph] 26 Nov 2021
In all these cases the catalysis reactions proceed rapidly enough to keep deuterons in relative equilibrium with nucleons. This circumstance explains the "snowballs in hell" paradox -the apparent survival of the light nuclei bound by just few MeV at temperatures of more than hundred MeV. The light nuclei do not survive, instead they are destroyed and created at similar rates. While the pion catalysis approach seems to be successful, it also received criticism [11], because the implementation of the πpn ↔ πd reaction involved a non-existent intermediate resonance d : the reaction was split into two stages pn ↔ d , πd ↔ πd. This purely technical simplification is alleviated in the present work -the catalysis reactions are implemented directly as a 3 ↔ 2 process, without involving an artificial intermediate resonance.
The difficulty of implementing multi-particle reactions arises from the collision criterion. In transport approaches, the most common criterion whether particles collide is based on the comparison of the distance between the particle and a geometric interpretation of the cross-section [12,13]. So far no generalization of this criterion from two particle to multi-particle reactions is available. Another possible criterion that is easily generalized to multi-particle reactions is directly derived from the collision integral of the Boltzmann equation and formulates a probability (stochastic rates) for a reaction [2]. Apart from the advantage of the straight-forward generalization to multi-particle reactions, it is inherently boost-invariant.
This work introduces the stochastic criterion in order to treat multi-particle reactions and its application to the deuteron catalysis reactions for the recently established transport approach SMASH [14]. The new treatment of multi-particle collisions allows to assess the previously found conclusions for the deuteron production on a more solid basis. It also allows to investigate the difference arising from modelling multi-particle reactions as a chain of two-body reactions (e.g. πππ → ρπ → ω). This helper construct with intermediate resonances is employed in SMASH in several places to adhere to geometric criterion and maintain detailed balance.
Employing a probabilistic collision criterion to enable direct multi-particle reactions is already explored in the literature, starting with [1,15]. While [15] focused on the optimisation of computing time and only applied the method for two-body reactions, the authors in [1] also discussed the production of deuterons. The approach presented here and the one from [1] both include the deuteron catalysis reaction involving nucleons, N pn → N d. For the low energies discussed in [1] this reaction is likely sufficient, since the system is dominated by N , but for the higher (intermediate) beam energy discussed here the pion catalysis reaction, πpn → πd, has to be included as well [10]. This is similar to a recent study for high beam energies [16]. Other authors focus on the BB annihilation reactions that produce multiple mesons [2][3][4] with the PHSD (Parton-Hadron String Dynamics) approach [17]. On the parton level, stochas-tic rates are studied with the parton transport approach BAMPS (Boltzmann Approach of Multi-Parton Scatterings) including the ggg ↔ gg gluon bremsstrahlung reactions [5]. The hadronic transport approach GiBUU also allows to employ a probabilistic collision criterion [18].
The structure of this article is as follows: First the employed transport approach and the treatment of deuterons is described in Sec. II. After the description of the geometric collision criterion so far employed in SMASH, the stochastic collision criterion is introduced by providing a sketch of its derivation. The stochastic criterion is validated by applying it to two-body reactions. Then, the treatment of multi-particle deuteron reactions is discussed. The stochastic multi-particle reactions are again validated by comparing the scattering rates to analytic results obtained from rate equations in Sec. III. The main results for the production of deuterons in Au+Au collision are presented in Sec. IV and Sec. V closes with a summary and conclusions.

II. MODEL DESCRIPTION
A. Transport Approach: SMASH The hadronic transport approach employed for the present study is SMASH [14], which is extended to include a stochastic collision criterion (see Sec. II C) and multiparticle reactions (see Sec. II D) that obey the detailed balance principle.
SMASH includes most hadronic degrees of freedom listed by the Particle Data Group [19] up to a mass of 2.3 GeV. See [14] for the list of included degrees of freedom and [20] for an update focused on the strange sector. The transport approach constitutes an effective solution of the Boltzmann equation by mapping the collision term to binary elastic and inelastic scatterings as well as the formation and decay of excited resonances. Additionally, string fragmentation to describe the high-energy cross sections and baryon-antibaryon annihilation together with the mentioned newly-introduced multi-particle reactions (see Sec. II D) are included. The reader is referred to [21] for details on the string fragmentation treatment. Resonances have vacuum properties and are tuned to reproduce the elementary cross sections. The Breit-Wigner spectral functions include a mass-dependent width, which follows the idea from Manley and Saleski in [22] (but with updated resonance properties). The approach is able to perform infinite-matter (box ), expanding sphere, collider and afterburner (list) calculations.
The version used throughout this work is SMASH-2.0.1 [23].
This new major version introduces the here discussed multi-particle reactions together with an integration of a hydrodynamics phase and improvements for high-energy AA collisions, when Pythia is required for hard scatterings. Also, the distance definition for the geometric criterion is updated to a fully covariant formulation [24].

B. Deuterons in SMASH
Deuterons in SMASH [8,10] are treated like on-shell point-particles as in [1,25,26]. The applicability of this approach is not strictly justified for the whole time of the evolution: the deuteron mean free path at the start of the afterburner evolution is comparable to twice the geometric size of deuteron wavefunction, although at later stages the mean free path increases making the approach more justifiable. Because of this we should strictly speaking call our "deuterons" correlated nucleon pairs. However, it seems that the inclusion of the finite deuteron size is not important in Pb+Pb collisions -a recent study [16] on the inclusion of the finite deuteron size in a similar approach only found a significant effect for a much smaller fireball created in pp collisions, which are not discussed in this work. Our approach with on-shell point-like d describes the data for heavy-ion collisions reasonably well, and is able to capture the fact that deuterons (or "correlated nucleon pairs") are created and destroyed. In our model the formation and disintegration of deuteron is catalyzed by pions or nucleons in the following reactions and their CPT-conjugates: πd ↔ πnp, N d ↔ N np, N d ↔N np and πd ↔ N N . Additionally, elastic reactions for πd, N d andN d are included.
The challenge in a microscopic transport approach is to treat the above mentioned Xd ↔ Xnp (X = {π, N,N }) catalysis reactions: Usually a geometric collision criterion is employed and it is not clear how the 3 → 2 reaction should be treated, since no generalization of a geometric collision criterion for more than 2 particles is available. Therefore, the 3 ↔ 2 is broken down into two steps with a fake decaying resonance (called d ) in the intermediate step, so that the reaction (chain) only contains 2-body collisions: Xd ↔ Xd ↔ Xnp. The properties of the d are chosen in order to reproduce the πd, N d andN d cross section and to achieve a lifetime that lasts the time n and p spend flying by each other. This treatment of the deuteron reactions is introduced in detail in [8] and will be labeled in the following as the 2-to-2 stochastic or 2-to-2 geometric treatment depending on the employed collision criterion.
The presented work introduces an alternative to the above described approach in the following, by employing a stochastic collision criterion, which allows for a direct (one-step) treatment of 3 ↔ 2 reactions (labeled as 3-to-2 stochastic, see Sec. II D).

C. Collision Criteria
At the core of microscopic (transport) approaches like SMASH, where point-like particles propagate and interact, the decision if and when they collide has to be made. Two different categories are known for this decision: Geometric and stochastic collision criteria.

Geometric Criterion
The most common approach to decide when two particles collide is based on a geometric interpretation of the cross section. The transverse distance of closest approach d T has to smaller than the interaction distance d int given by the cross section: There are several definitions available for the transverse distance of closest approach d T of two particles [13,24,27]. SMASH followed the approach introduced by UrQMD [13] up to version SMASH-2.0, where the approach is changed to the fully covariant formulation given in [24]. One main disadvantage of the criterion is the lack of a straight-forward generalization of d T to more than two particles.

Stochastic Criterion
In contrast to the geometric criterion, the stochastic criterion defines a probability for a reaction of a given particle set. The probability is defined as the number of reactions ∆N reactions over the number of all possible particle combinations inside a sub-volume ∆ 3 x and time interval ∆t, The final expression of P n→m is worked out in App. A. The stochastic criterion is inherently boost-invariant. This idea was already explored in the literature before in [1-5, 15, 18], also often referred to as the local-ensemble method.
Several special cases (for n, m = 1, 2, 3) are presented in App. A. In particular for a decay process n = 1, where m is the number of final particles after the decay, and M is the total energy in the center-of-mass of the system (or the mass of the initial resonance in its rest frame). This decay probability is the same already used in SMASH (cf. Eq. (38) in [14]). For n = 2 one can write the probability in terms of the total cross section, where the relative velocity is defined in Eq. (A14). For m = 1 the cross section is given in Eq. (A15), and coincides with σ 2→1 used in [14]. In this paper we use the case m = 3, and σ 2→3 is fitted to available experimental data in [8] as described above.
Other cases can be worked out along the lines given in App. A. Even though, this straight-forward generalization to arbitrary n → m reactions is a large advantage compared to geometric criteria, Eq. (A2) cannot be applied directly, since the matrix element is generally unknown.
When employing the test-particle method the particle number is scaled with the test-particle number ∆N i → ∆N i N test in the here presented notation. The number of collision scales similar: ∆N n→m reactions → ∆N n→m reactions N test . Therefore the scaled probability P n→m when employing the test-particle method is For the numerical realization of the stochastic criterion, the space is divided into equally sized cells (∆ 3 x) and the probability is calculated for all 2-particle and 3particle combinations within each time-step, so only particles within a cell are able to interact. The calculated probability is used for a Monte-Carlo decision i.e. if a generated random number between 0 and 1 is smaller or equal than the probability the reaction is accepted. The collision time is randomly chosen within the given time-step ∆t.
In contrast to the geometric criterion, the stochastic criterion is a strictly time-step based method. The timesteps, therefore, have to be chosen small enough that the assumption that each particle only interacts once per time-step is justified and the defined probability is not exceeding 1 for a given ∆ 3 x. Also the cells ∆ 3 x have to be chosen sufficiently small, since only in the limit of ∆ 3 x → 0 (and ∆t → 0) the numerical solution matches the exact solution of the Boltzmann equation. At the same time, the cells still have to be sufficiently filled with particles.
To verify the stochastic collision criterion, first the collision rate is studied in Fig. 1. The expected collision rate for the case of a box filled with π 0 that interact via a constant elastic cross section is given by ρσ. Fig. 1 shows that the criterion performs well, if the time-step is chosen small enough as mentioned above. Otherwise it exhibits the usual limitation of time-step based approaches: When the density is large and the time-step is also large more than one interaction per time-step is expected and therefore the collision rate is underestimated. It is furthermore checked that the scattering rate is consistently reproducing the analytic expectation also when varying other parameters like the box size, the test-particle number, the cross section and the temperature of the box.
Complementing the result for the collision rate in the box, produced particle multiplicities from heavy-ion re- actions are shown in Fig. 2 for different beam energies.
Since the collision criterion is in the end an implementation detail, it is reassuring that the particle yield does not depend on the chosen collision criterion as Fig. 2 exemplifies. The stochastic collision criterion matches the previous results with the geometric criterion for a set of abundant hadronic species as a function of beam energy.

D. Multi-particle Reactions
The main class of multi-particle reactions realized with the stochastic collision criterion are 3 ↔ 2 reactions. These are needed for the description of creation and annihilation of deuterons in the hadron gas. Expressions for P 2→3 and P 3→2 are derived in Eqs. (A22) and (A24), respectively. For completeness, similar calculations for 3 ↔ 1 reactions are given in Appendix B.
Since the scattering matrix element is commonly unknown, Eq. (A2) cannot be directly used to calculate the required probability. Assuming that the scattering matrix element is independent of the final momenta, the probability for the reverse process of a decay (n → 1) or a 2-body scatterings (n → 2) can be expressed in terms of the known decay width or cross section of the {1, 2} → n process. This is the case, because the scattering amplitude is the same for the forward and backward process (time-reversal invariance).
Let us focus here on the P 3→2 probability applied to the deuteron case. The stochastic criterion allows for a direct (one-step) treatment of 3 → 2 reactions, which is applied to the deuteron formation catalysis reactions πpn ↔ πd, N pn ↔ N d andN pn ↔N d (and labeled as the 3-to-2 stochastic treatment in the results below.).
The probability for a 3-to-2 process is given as where Φ 2 , Φ 3 are the 2-and 3-body phase space corresponding to final and initial states, respectively (see also (A24)). The symmetry factors for the reactions πpn → πd, N pn →N d and their charge conjugated are S! S ! = 1, while for N pn → N d and it charge conjugated, as there are 2 identical particles in the initial state, are S! S ! = 2 . The spin degeneracies read Note that the derivation presented in the Appendix A resembles the one presented in [2] and [5]. The absence of a known scattering matrix element makes additional steps necessary in the probability derivation compared to [5]. A similar idea as in [2] is therefore followed, where the scattering matrix element is assumed to be only dependent on the initial center-of-mass energy. This approach differs from the idea in [1], where the matrix element is factorized into two terms as an approximation (also discussed in [5]). This factorization results in a different formula for the 3-to-2 collision probability containing the two-body N N cross section and a momentum dependent volume [1]. The approach presented here, while also making an assumption about the matrix element, is more general in the sense that it allows to treat all processes, where a decay width or cross section for the reverse reaction is known.

E. Hybrid Approach: Hydrodynamics and Cooper-Frye
In this work we focus on the deuteron production in the hadronic afterburner. Therefore, we describe mod-elling of the earlier stages of the heavy-ion collisions only briefly. The expansion and cooling of the fireball during the dense stage is simulated by the relativistic 3+1-dimensional open-source hydrodynamic code MUSIC v3.0 [28][29][30][31]. The initial condition for MUSIC is a smooth parametrized energy density and baryon density as functions of spatial coordinates. The parametrization is described in [32], it is tuned to reproduce charged particle rapidity distribution, pion midrapidity yields, and net-proton midrapidity yields. The initial energymomentum tensor is assumed to have an ideal fluid form T µν = ( + p)u µ u ν − pg µν . The flow at initial eigentime τ 0 (τ 0 = 3.6 fm/c at √ s NN =7.7 GeV) is assumed to be only longitudinal and have Bjorken form u µ = (cosh η s , 0, 0, sinh η s ). The equation of state combined with hydrodynamic equations is a lattice QCD based "NEOS-BSQ" equation of state p = p( , n B ) described in Ref. [33]. Shear viscous corrections are included with a specific shear viscosity ηT /(e + P ) = 0.1, while bulk viscous corrections and baryon number diffusion are neglected. Particlization is performed at a constant energy-density hypersurface, (τ, x, y, η s ) = 0.26 GeV/fm 3 . This corresponds to a line in (T, µ) plane; the distribution of hypersurface cells number at |η s | < 0.5 peaks at (T, µ) = (141±4, 340±50) MeV, where the numbers after ± denote an approximate width of the peak. This is slightly lower energy density and baryon chemical potential, but almost the same temperature compared to the chemical freeze-out in the thermal model fit of the hadron yields in central AuAu collisions at √ s NN = 7.7 GeV, where (T, µ) = (143.8 ± 2.7, 399.8 ± 13.3) MeV [34]. The particlization procedure is a standard grandcanonical Cooper-Frye particlization with a Grad 14moments ansatz for shear viscous corrections, see [28,29] for details. Particles from particlization are inserted into hadronic afterburner and rolled back to equal time, but are forbidden to interact until their actual emergence time. Deuterons may be sampled or not at the particlization, this is mentioned specifically in the text. The d is never sampled. Overall, the most important features of the hydrodynamic part of our simulation is that it provides a reasonable space-time distribution of nucleons and pions at particlization. The nucleon and pion rapidity and transverse momentum spectra reproduce experimental data rather well. Based on our previous work, these are sufficient prerequisites to describe the measured deuteron spectra by introducing deuteron-creating reactions into the afterburner.

III. VALIDATION OF STOCHASTIC MULTI-PARTICLE RATES IN THE BOX
Before showing the physics results, let us start by demonstrating that the implementation of the stochastic criterion and the newly introduced 3-to-2 multi-particle reactions work as expected. For this purpose, a particle configuration is initialized in a box with periodic boundary conditions and compared to the analytic expectation. The goal is to ensure that the content of the box equilibrates, that it equilibrates to the correct state, and that particle multiplicities change at the expected rate in the process of equilibration. This method of testing has proven very useful, because it checks detailed balance and reaction rates, and is sensitive even to minor errors in implementation.

A. Rate equations
The analytic expectation is provided by rate equations of the same form as introduced in [35] with only one exception -we take the spectral function into account, when we compute thermally averaged resonance widths Γ . The necessary derivation and notation for writing the rate equations for the specific reaction systems is given in Appendix C. Below, only the resulting systems of rate equations are given.
The first case is the d, d , π, N system with 2 ↔ 2 reactions, where in addition to elastic collisions the following reactions are allowed: Denoting time derivative dλ dt asλ the system of rate equations for this case is given as, Here the amount of protons and neutrons is assumed equal. The initial conditions are determined by the content of the box at time t = 0. The thermally averaged width of d , Γ d , was artificially divided by factor 2 to agree with the simulation in Fig. 3. The need in factor 2 might emerge from the fact that the d spectral function has a long high-mass tail, which takes more time to equilibrate than expected.
The second case is a system of d, π, N with 3 ↔ 2 reactions, where the following reactions are allowed: Resulting in the corresponding rate equations as follows, In this second case, no additional correction factors need to be applied, since it avoids the fictitious d resonance in the first place.

B. Comparing analytic results to simulations
Let us compare the solutions of the rate equations above to corresponding simulations in a box, which allows to probe the equilibrium properties of the employed reaction treatment and collision criteria. The verification of the approach in a static box scenario is the basis for further exploration of dynamic non-equilibrium systems like nucleus-nucleus collisions in the remainder of this work. The box size is set to be V = (10 fm ) 3 , particles are initialized uniformly in coordinate space and according to a Boltzmann distribution in momentum space with temperature T = 0.155 GeV. The initial multiplicities are chosen to be 30 for each pion species and 60 for each nucleon species. The cross sections of πd and N d are the ones described in [8] and taken further for Au+Au simulations. The πd ↔ N N reaction is switched off for the box test. The only allowed reactions are elastic collisions and the reactions present in Eq. (9). In Eq. (9) we assume that the temperature stays constant over time. In general, this does not have to be true in the box simulation, but we check by fitting the momentum distributions that the temperature in fact stays constant. The comparison between the particle yields of the analytic calculation and the simulations is presented in Figure 3. Since the box is only filled with π and N in the beginning, the production of d over time is observed. The simulation agrees for the two different reaction treatments for d production (2 ↔ 2 with resonance and direct 2 ↔ 3) as well as employing either of the two criteria (geometric and stochastic criterion). The agreement allows not only to validate equilibration to the correct yields, but also the equilibration process itself. It is separately verified that detailed balance principle is fulfilled for all allowed reactions once the yields are equilibrated.
Physically more interesting than the verification of the different reaction treatments, is the difference observed between the now possible 3-to-2 treatment and the modeling of the same reaction via 2-to-2 reaction. Employing direct multi-particle reactions leads to a significantly faster rise of the deuteron yield. Consequently the equilibrated yields for d (and p) are reached significantly earlier, as predicted by the rate equations. The impact of this observed faster equilibration of the system, when employing direct 3-to-2 reactions in an expanding medium created in heavy-ion collisions, is discussed throughout the remainder of this work.
A similar result is presented for 3-to-1 reactions as part of Appendix B.

IV. DEUTERON PRODUCTION IN AU+AU COLLISIONS
In the following, the results for the deuteron production in Au+Au collisions at √ s NN = 7.7 GeV are discussed employing the newly introduced 3-to-2 reactions. All presented results display the deuteron evolution in the afterburner stage of the hybrid approach. Special emphasis is placed on the difference to the calculations modeling the same process with a 2-to-2 reaction chain including the d resonance. Furthermore, two scenarios at particlization are distinguished in the results. In one case deuterons are assumed to be produced in the hydrodynamic stage of the collision (with d at particlization), in the other case no deuterons are present at the start of the afterburner calculation (without d at particlization).
The distinction allows to compare two different pictures for deuteron production. The thermal model-like picture, where d are produced early at high temperatures and the coalescence-like picture, where d are assumed to be formed at later times. Figure 4 shows the number of deuterons propagated in the afterburner stage over time. The d production is enhanced when employing direct 3-to-2 reactions. Es- pecially in the case with no d at particlization, a more rapid increase of the d number is observed, which drives the number close to the case with d at particlization. The final number of deuterons is almost identical for the two particlization scenarios. The remaining difference is on the order of experimental errors. The difference is smaller for the calculation with the 3-to-2 reactions in comparison to the 2-to-2 approach. Those findings are understandable considering the above observed faster equilibration when employing multi-particle reactions. The 3to-2 reactions drive the system faster to statistical equilibrium before it freezes out due to its expansion. The expansion is also the reason why the yields without d at particlization are not in full agreement with d at particlization, the d reactions seize too quickly due to the cooling before enough d can be produced (cf. Figure 5  and 6). Both particlization scenarios are also in agreement with the experimental values for 0-10% centrality, which shows that d yield is possible to understand in terms of multi-particle catalysis reaction being the main production mechanism.
Comparing the two presented centrality classes in Fig-Centrality N   Going from central to more peripheral collisions, the equilibrated yield without d at particlization (N without d, equil ) in comparison to the final yield with d at particlization (N with d, equil ) is less in agreement as Table I shows. The smaller medium created in peripheral collision seems to suppress the full statistical equilibration of the system before freeze-out when all d are produced in the late (af-terburner) stages.
Employing the stochastic criterion for the 2-to-2 reaction chain (not shown), yields the same results as obtained with the geometric criterion. The evolution of the d yield is the result of the competing 3 → 2 formation and 2 → 3 break-up reaction rates shown in Figure 5 for the case with d at particlization and in Figure 6 for without d at particlization. Forward and backward (direct) 3-to-2 reactions rates are close for central collisions in Figure 5, but as also seen for the yield some time is necessary before they are close to being equilibrated. The collision rates also clearly indicate the dominance of the π catalysis reactions. This underlines the necessity to include the π in addition to the N reactions for this beam energy, which is the main extension compared to [1]. The π reactions are also closer to being equilibrated than the N reactions. For more peripheral reactions and with this a smaller medium, the 2 → 3 rate is dominating over the formation reactions. This relatively lower 3-to-2 reaction rate again hints at incomplete statistical equilibration in the smaller system. Considering the calculation without d in Figure 6, the 3 → 2 reaction dominates for all rates, as expected by the rapid rise of d numbers at the beginning of the evolution. The reaction rate figures also allow to pinpoint a chemical freeze-out at least for the shown 3 ↔ 2 d reactions at around 50 fm for central collisions (0-10%).
Since the system is smaller for more peripheral collisions (30-40%), the freeze-out is earlier at around 30 fm. Note that while the N N ↔ πd reaction is also included in the calculation, its contribution is only sub-leading, even compared to the N catalysis reaction, and therefore the rate is not shown. Experimental data from [34,38].
In addition to the yields, first the average transverse momentum is presented for protons and deuterons in Figure 7 for 6 centrality classes. The mean-p T slightly declines towards more peripheral collisions. Here and in the following, three treatments for the deuteron catalysis reactions are compared: calculating with the 2-to-2 reaction chain for the geometric (blue round points) and the stochastic (orange triangles) criterion as well as direct 3-to-2 reactions (red squares) that are only possible with the stochastic criterion. The mean-p T results are unaffected by the different approaches and all agree with the available experimental data [36,37] within errors for both, p and d, validating the transverse dynamics of the calculations. Note that this confirms the previous findings in [10], where the here employed hybrid approach was carefully constrained by an extensive experimental data set.  Figure 8 presents the results for the elliptic flow (v 2 ) for p and d for the same 6 centrality classes, as a more sensitive probe of the momentum distribution. Due to limited statistics of the calculation only the integrated v 2 is presented, which still allows to contrast the different approach for the d reactions. Even though no experimental data is available for this observable, the order of magnitude of around 0.1 is comparable to the p T dependent v 2 reported in [38]. Regarding the integrated v 2 of d in the bottom panel of the figure, the different reaction treatments are found to have an effect. While for central collision they agree, for more peripheral collisions the elliptic flow is decreased for the 3-to-2 reactions and by employing the stochastic criterion. The latter is also found for the protons in the upper panel. The ad-ditional effect for 2-to-2 reactions for the two different collision criteria might hint at limitations of the stochastic criterion for the small systems in peripheral collisions. However, a more uniform, thermalized medium and a subsequent lower v 2 might also be expected when employing multi-particle reactions. In Figure 9 the impact of d solely being produced in the afterburner on the elliptic flow is studied. Observing a clear difference here would potentially enable to disentangle the two pictures for d production times. Even though a clear difference is not found, a small decrease is seen for the case without d at particlization. The difference is not significant for all centrality classes, the effect for the average over all centralities, however, is significant.

V. CONCLUSION
This work presents the first application of multiparticle reactions in the hadronic transport approach SMASH by employing a stochastic collision criterion. This is of major importance since it allows to treat multiparticle reactions while fulfilling detailed balance.
A comprehensive derivation of the stochastic collision criterion is shown for different reaction classes, from an one-and two-body initial state to the multi-particle reactions with 3 particles. The derivation for 3-body reactions focuses on the inverse reaction for meson Dalitz decays (3 → 1) and the deuteron catalysis reactions (3 ↔ 2): πpn ↔ πd and N pn ↔ N d. The stochastic collision criterion allows to avoid the introduction of an artificial resonance by treating the deuteron 3 → 2 catalysis reactions in one-step while simultaneously fulfilling detailed balance.
The newly-introduced stochastic criterion is validated by an agreement with the analytic expectation for the two-body collision rate in a box setup with elastic collisions. In addition, the particle multiplicities for heavy-ion reactions at different beam energies are reproduced, when including the same reactions as in the default SMASH version. An analytic expression for the chemical equilibration is derived with rate equations. The stochastic multi-particle reactions for the deuteron catalysis are validated by agreeing with the analytic results, exhibiting the correct chemical equilibration and detailed balance. As predicted by the rate equations, the time required to equilibrate is significantly reduced for the direct 3 ↔ 2 treatment compared to the previous 2 ↔ 2 multi-step process. The finding of a faster equilibration employing multi-particle reactions is also confirmed in 3-to-1 reactions ( Figure 10).
Studying the deuteron production in gold-gold collision at a beam energy of 7.7 GeV in a hybrid approach, the faster equilibration process of the multi-particle reactions leads to a more rapid increase in the deuteron yield before the system freezes out chemically due to the expansion of the system. The d yield is consequently enhanced. The difference in the final number of d, when comparing the scenarios of d being produced in the hydrodynamic stage or just in the hadronic afterburner, is greatly reduced when employing multi-particle reactions due to three-body reaction driving the d faster to statistical equilibrium. The yield agrees with the experimental data if the nuclei are produced at time of particilization or not. Hereby confirming the previous findings when employing the slower equilibrating two-body reaction chain involving the fake d resonance [8,10]. In addition, a decrease in the elliptic flow is found when employing the stochastic criterion as well as multi-particle reactions for more peripheral collisions. Similarly a small decline in flow for all centralities is reported, if all d are produced during the late (afterburner) stages of the collision. No dependency is found for mean transverse momentum, where an agreement with experimental data for all centrality classes is seen.
The stochastic criterion for multi-particle reactions also opens the possibility to investigate other reactions in the future. Of particular interest is the theoretically straightforward extension of this work to include other light nuclei like (hyper-) triton and 4 He. While the conclusions employing the 2 ↔ 2 reaction treatment involving the artificial d are confirmed by this work, introducing additional fake resonances like t is tedious and would certainly enhance the found differences to the theoretical preferable direct 3 ↔ 2 treatment. In addition, the stochastic reaction approach is directly extendable to study larger n-particle reactions (n > 3). Here, the study of the back-reaction of pp annihilation (5π → pp) will be especially interesting. Appendix A: Reaction probability for a n → m process The stochastic collision criterion uses a reaction probability which depends on the number and type (stable particle or resonances) of incoming and outgoing states. For a generic reaction from n to m particles (we will denote with primed indices the states in the final state) the reaction probability reads Using the Boltzmann equation to quantify the number of reactions in a phase-space element during a time ∆t we finally obtain, where S is the number of identical particles in the final state (to avoid double counting when integrating the phase-space variables), e.g. if the final state contains N pn, then S ! = 2. |T n→m | 2 is the scattering matrix squared of the process, summed over final internal states and averaged over the initial ones, where g j is the internal degeneracy of the state j. Since we distinguish isospin states, it will be used to account for spin degeneracy g j = 2s j + 1, where s j is the spin of the state. In Eq. (A2) the factor dΦ m is the m−particle phasespace element, where P is the total 4-momentum of the reaction. For stable particles we will use (keeping implicit the particle subindex k for simplicity) with E 2 = p 2 + m 2 , while for resonances where M = √ s = p 0,2 − p 2 , E 2 M = p 2 + M 2 and A(M ) is the spectral function of the state, normalized to In the limit of very narrow resonance one recovers the stable particle case by doing A(M ) = 2M δ(M 2 − m 2 ). Usually one assumes for simplicity that the scattering amplitude (when it is unknown) only depends on the initial center-of-mass energy. In such cases one can pull it out of the integral over the final phase-space. When this is used, the integrated m−body phase space reads for m = 1, 2, 3, For special cases n, m = 1, 2 one can write the probabilities in terms of the decay width of interaction cross sections. Let us consider the particular cases considered in this work.
• n = 1, m = 2 In this case we can introduce the decay width of a resonance (with the degeneracy and symmetry factors incorporated) to arrive to • n = 2, m = 1 The inverse reaction reads For practical purposes we introduce the relative velocity and the 2 → 1 cross section (cf. Eq. (39) in [14]) The P 2→1 reads finally • n = 1, m = 3 The decay to three particles is similar, where we introduced the decay width • n = 3, m = 1 This case can be written as where the factor S counts the number of identical particles in the initial state of the 3 → 1 reaction. It appears here because the same factor is included in the definition of Γ 1→3 [called S in Eq. (A18)], and needs to be canceled.
• n = m = 2 This case represents a binary collision where • n = 2, m = 3 The next two cases are needed for the deuteron formation/annihilation.
where the 2 → 3 cross section reads . (A23) • n = 3, m = 2 It is important to remember that the factor S (S ) always refers to the number of identical particles in the initial (final) state of the particular reaction which is considered. For example the factor S appearing in Eq. (A23) is, in general, different from the one appearing in Eq. (A24). While the 3-to-1 back-reaction are found to be rare in nucleus-nucleus collisions, they still allow to test and study the approach and effects for multi-particle reactions in general. Therefore, the equilibration in a box employing the reaction is checked in a similar fashion as presented above for the deuteron reactions in Figure 10 for the πππ → ω reaction as an example. The (10 fm) 3sized box is filled with 300 ω mesons at initialization. The calculation including the 3-to-1 multi-particle reactions is displayed as the solid lines and clearly shows to equilibrate chemically.
To verify the equilibration of the system, the ratio of equilibrated particle numbers of R ωπ = Nω N 3 π is compared to the thermal i.e. grand-canonical ideal gas expectation for this ratio (R th ωπ ). For this the temperature of the box is extracted (T = 160 MeV) by fitting the π energy spectrum assuming an exponential shape in equilibrium. Taking the ratio cancels the for the employed reactions specific chemical potential.
The calculated and thermal ratio are consistent with each other: Rωπ−R th ωπ R th ωπ = 0.04, thereby verifying the here presented approach for 3-to-1 reactions. That detailed balance, in particular for the 3-to-1 reaction is fulfilled, is also checked. Figure 10 furthermore shows the results for modeling the same 3-to-1 reaction with two-body reactions (πππ → ρπ → ω) labeled 2-to-2. Note that the ρ (and the the ρ ↔ ππ) reaction were only included in the 3-to-1 calculation to have the same degrees of freedom as in the 2-to-2 case. Its inclusion is not necessary to employ the 3-to-1 reactions of interest. The equilibrated yield is matching for both, multi-particle and multi-step, treatment. Interestingly, comparing the two further, a very similar trend as for the same comparison for d reactions is observed (cf. Figure 3). The direct treatment of the 3to-1 reaction leads to a faster equilibration of the yields. Hereby, confirming the findings made for deuterons that multi-particle reactions drive the medium faster to equilibrium.
Appendix C: Rate equations derivation As in Ref. [35] we start by assuming that all species in the box are in kinetic, but not chemical equilibrium. This means that their momentum distribution is a Boltzmann distribution at temperature T , but its normalization, the yield, is arbitrary. Therefore, every species i is assigned a fugacity λ i and the multiplicity N i is written as where V is the box volume, g i is the degeneracy of the species, A(M ) its spectral function as defined after Eq. (A6), and K 2 (x) is the modified Bessel function of the second kind. Finally, E 2 M = p 2 + M 2 . We assume that the temperature is constant in time, which indeed turns out to be the case for temperature from momentum spectra in the simulations. Therefore, our rate equations are going to be equations for fugacities. The reaction rates (number of reactions per unit volume per unit time) for 2 ↔ n reactions indexed as 1 + 2 ↔ 1 + 2 + · · · + n are dN 2→n with A ≡ g 1 d 3 p 1 (2π ) 3 g 2 d 3 p 2 (2π ) 3 σ 2→n v rel e −(E1+E2)/T = σ 2→n v rel n th 1 (T )n th 2 (T ) .
In our case the cross section σ 2→n depends only on the center of mass energy of the reaction √ s. In this case we find that σ 2→n v rel expression simplifies to [41] σ 2→n v rel = This expression is equivalent to Eq. (30) from Ref. [42]. The rate of 1 ↔ n decays and formations is As mentioned in Section II C 2, the stochastic collision criterion divides the space into a grid with equally sized cells. The cell size (∆ 3 x) is used for the collision probability calculation. Only particles within a cell interact. The general requirement is that the cells are small, but still sufficiently filled with (test) particles. An appropriate grid choice balances these requirements with the calculation runtime. For the presented Au+Au afterburner calculations (Section IV), this is achieved by choosing a fixed (minimal) cell length of l min = 1.78 fm together with a test particle number of N test = 20. The grid is chosen to be the same size as the medium i.e. all particles are inside of grid cells and updated at every timestep. The number of cells is determined by dividing l min by the total grid length in each dimension and rounding down to an integer value. If necessary, all cell lengths are increased from their minimal value to fill the total grid length. This increase is negligible for the presented results, since the number of cells is large. The timestep for all calculations is ∆t = 0.1 fm. It was separately verified that the results are stable when employing a smaller grid length, more test particles and a smaller timestep.