Deuteron production in AuAu collisions at $\sqrt{s_{NN}} = 7-200$ GeV via pion catalysis

We study deuteron production using no-coalescence hydrodynamic + transport simulations of central AuAu collisions at $\sqrt{s_{NN}} = 7 - 200$ GeV. Deuterons are sampled thermally at the transition from hydrodynamics to transport, and interact in transport dominantly via $\pi p n \leftrightarrow \pi d$ reactions. The measured proton, Lambda, and deuteron transverse momentum spectra and yields are reproduced well for all collision energies considered. We further provide a possible explanation for the measured minimum in the energy dependence of the coalescence parameter, $B_2(\sqrt{s_{NN}})$ as well as for the difference between $B_2(d)$ for deuterons and that for anti-deuterons, $B_2(\bar{d})$.


I. INTRODUCTION
Heavy ion collisions are often called "Little Bang" due to a rapid expansion, cooling, and a sequence of freezeouts reminiscent of the evolution of the early Universe. Another common feature of the Little and Big Bangs is nucleosynthesis, or production of light nuclei. The Big Bang nucleosynthesis for deuterons primarily occurred via pn ↔ dγ reaction. In relativistic heavy ion collisions this reaction does not have sufficient time to create the observed amount of deuterons, which follows from its small cross section and the typical time of collision ∼ 10 −23 − 10 −22 s. Here other reactions are at work, depending on collision energy and rapidity region. In particular, we have previously suggested that the pion catalysis reaction πpn ↔ πd plays the dominant role in deuteron production at √ s N N = 2760 GeV in the midrapidity region [1,2]. In the present paper we shall argue that the same reaction is still the most important one down to collision energies of √ s 7 GeV for deuteron production at mid-rapidity.
We note, however, that the two most popular models of deuteron production -thermal [3][4][5][6][7] and coalescence [8][9][10][11][12][13][14][15][16][17][18][19] models -do not need to explicitly involve any particular reactions. The thermal model postulates that light nuclei are created from a fireball in chemical equilibrium with hadrons. At the chemical freeze-out the reactions that change hadron yields cease and hadrons only continue to collide (quasi-)elastically. These collisions change the momentum distributions, but do not change the yields. Thus, for hadrons the chemical freeze-out temperature, T CF O , which is determined from hadron yields, is larger than the kinetic temperature T KF O which is extracted from the momentum spectra, T CF O > T KF O . This picture is supported by the fact that the yieldchanging reactions typically have smaller cross sections, so they cease earlier during the expansion of the fireball. Deuteron yields and spectra are consistent with the same T CF O and T KF O for nuclei as for hadrons [20]. This means that they have to be colliding with other particles between chemical and kinetic freeze-out. However the 2.2 MeV binding energy of deuterons is much smaller than T CF O 150 MeV or T KF O 110 MeV. Simple intuition tells that a deuteron must be easily destroyed at such temperatures. Due to this intuition light nuclei in heavy ion collisions were called "snowballs in hell" [21], where light nuclei would be "snowballs" and the fireball of the heavy ion collisions is referred to as "hell". However, this simple intuition fails in two ways. Firstly, even despite the small binding energy of a deuteron, the elastic cross section of d + π → d + π reaches as high as 70 mb at the kinetic energies of pion and proton corresponding to temperatures of 100-150 MeV (see Fig. 1 of [1]). One assumes that a thermal pion at T 150 MeV should easily break up a deuteron, but in σ elastic πd /σ total πd 1/4 of all π+d collisions this does not happen: instead, the pion excites one of the nucleons, which subsequently de-excites emitting a pion back while leaving the deuteron intact. Secondly, the inelastic reactions that destroy deuterons (d+X ↔ p+n+X, where X is an arbitrary hadron) have backreactions that can also create deuterons. We have shown in [1,2] that for Pb+Pb collisions at √ s N N = 2.76 TeV deuteron creation and destruction occur at approximately equal rates between T CF O and T KF O . This mechanism, thus, resolves the "snowballs in hell" paradox. It justifies calculating the deuteron yield in the hadron resonance gas model at T CF O while determining the deuteron momentum spectrum in a blast wave model at T KF O .
In contrast to the thermal model, coalescence models postulate that deuterons are produced from nucleons at the kinetic freeze-out. Nucleons coalesce into a deuteron in case they are close in the phase space. By energy and momentum conservation proton and neutron cannot just form a deuteron without a particle carrying away the energy that binds a deuteron. Coalescence implies that there is a reaction like pn → dγ, or pn → dπ, or pn → Xd, or Xpn → Xd, where X is some particle carrying away the energy. The concrete species of the particle is not important for coalescence models, because this energy is always neglected in the actual coalescence calculation, although in principle a method was developed to effectively account for it by putting nucleons off-shell [19]. Coalescence models sometimes consider not only pn → d coalescence, but also coalescence of two protons or two neutrons to a deuteron, see for example [22]. This implies charge exchange reactions like π − pp ↔ π 0 d, π 0 pp ↔ π + d, π 0 nn ↔ π − d, π + nn ↔ π 0 d. However, coalescence models never consider coalescence of particles other than nucleons to deuteron, although the reactions like ∆p ↔ πd, N * p ↔ πd, ∆∆ ↔ dππ, or Λp ↔ dK are possible.
Thermal and coalescence models are considered to be opposite and conflicting scenarios for deuteron production. Indeed, the thermal model postulates deuteron production at chemical freeze-out and coalescence postulates production at kinetic freeze-out. In the thermal model nucleons from resonances do not contribute to deuterons, while in coalescence all nucleons, including the ones from resonance decays, are able to produce deuterons. In coalescence the spatial extend of deuteron wavefunction relative to the size of colliding system matters, while in the thermal model it does not. Despite these differences, we are able to accommodate the core ideas of both thermal and coalescence models in our dynamical simulation in the following way: We use relativistic hydrodynamics to simulate the locally equilibrated fireball evolution until the chemical freeze-out. Then we switch to particles and allow them to rescatter using a hadronic cascade.
To this end we sample deuterons and all other hadrons according to the local temperature and chemical potential of the switching hypersurface, which in our work is controlled by a certain value of the energy density, sw . This transition from hydrodynamic to transport is often referred to as particlization. The deuteron yield at this moment is the yield one would obtain in the thermal model, where the volume V is determined by the particlization hypersurface. The majority of these initial (thermal) deuterons is destroyed in the subsequent kinetic evolution, but at the same time the new ones are created, so that the average yield remains approximately the same. Like in the coalescence model deuterons that finally survive are mostly, although not exclusively, those which are created very late in the hadronic evolution, and thus do not experience any more collisions. Also our rate of deuteron production has a large peak at low relative momentum between nucleons, which means they have to be close in the phase space to make a deuteron, as the coalescence model postulates.
Conceptually our approach here is the same as in our previous study at 2760 GeV [1,2]. The key difference is that at lower energies one has to account for the evolution of the net-baryon current, which does not vanish anymore. This requires additional equations for baryon current conservation in the hydrodynamic simulation and the specification of initial condition for the baryon density. These extensions are presently not available in the CLVisc code we used in previous study. Therefore, for the present study the MUSIC code for the 3D hydrodynamical evolution of the background medium together with a geometric-based initial condition, which had already been tuned to reproduce hadron spectra at the considered energy range [23]. In our previous study we found that due to reactions like πpn ↔ πd the deuteron yields and spectra and intimately related to those of the pro-tons. Therefore, our primary concern is a good description of the measured proton spectra. For the considered collision energies ( √ s N N ≥ 7 GeV) light nuclei production is merely a small perturbation over the space-time evolution of baryon density. For example, the d/p ratio at 7.7 GeV at midrapidity is around 0.03, and at higher energies this ratio is even smaller. Therefore, one may view any dynamical model of light nuclei production as a combination of a "background" of expanding fireball with nucleon density evolving in space-time, and of a "perturbation" that acts on this background and creates (and possibly disintegrates) nuclei. No matter how detailed and realistic the nuclei production model is, the overall precision cannot be better than that of the background model. That is why we pay attention to fitting proton yields and spectra well. It turned out that a precise account of weak decays is important, so we also make sure that we reproduce the yields of Λ baryon. As we shall show, applying our model of deuteron production to Au+Au collisions at energies from 7 to 200 GeV we observe that we are able to describe the measured deuteron spectra and yields using the same reactions πpn ↔ πd, N pn ↔ N d, N N ↔ πd with the same cross sections that we used to describe the yields, spectra and flow at 2760 GeV [1,2]. The work is organized as follows: in Section II we explain the details of our simulation, and in Section III we present and discuss the resulting proton, Λ, and deuteron spectra, yields, and reaction rates relevant for deuteron production. Finally, we explore the role of the correction from the weak decays feeddown to protons in the context of various observables involving deuterons.

II. HYDRODYNAMICS + TRANSPORT SIMULATION METHODOLOGY
To simulate the full evolution of a system created in heavy ion collisions we employ a hybrid relativistic hydrodynamics + hadronic transport approach. Hydrodynamics is applied at the earlier stage of collision, where the density is high and, therefore, hadrons cannot be treated as individual particles. The transition from the quarkgluon plasma to a hadron gas is handled implicitly via the equation of state used in the hydrodynamics. Hadronic transport is applied at the later stage of collision, when the fireball is dilute enough so that mean free paths of the particles are larger than their Compton wavelengths.

A. Initial state
The initial conditions and hydrodynamic simulations are based on the work in Ref. [23]. In our study the main role of this initial state is that it allows to reproduce the measured proton phase-space distributions and midrapidity yields at different collision energies. Here we mention some features of the initial state, and the full de-scription can be found in Ref. [23]. We start simulations by generating event-averaged initial energy density and net baryon density profiles for hydrodynamics at a constant proper time τ = √ t 2 − z 2 = τ 0 , where z is the coordinate along the collision axis. The proper time τ 0 is a function of collision energy, which is chosen to be slightly longer than the overlapping time that would take nuclei to pass through each other in absence of interactions The values of τ 0 were chosen such that the the hydrodynamics + hadronic transport simulations can reasonably reproduce the measured mean transverse momentum of identified particles [23]. The initial energy-momentum tensor is assumed to have a diagonal ideal-fluid form T µν = ( +p)u µ u ν −pg µν . The initial baryon current is also assumed to have an ideal-fluid form, j µ = n B u µ . At τ = τ 0 , Bjorken flow is assumed: u µ = (cosh η s , 0, 0, sinh η s ). Based on the local energy and momentum conservation, the local rest frame energy-density (x, y, η s ) profile in case of our symmetric Au+Au collision system is parametrized as described in [23]: where the normalization factor N e (x, y) is, and T A (x, y) and T B (x, y) are the nuclear thickness functions for the incoming projectile and target nucleus. Although our colliding system is symmetric, the local nuclear thickness functions T A (x, y) = T B (x, y) at a nonzero impact parameter. For the initial baryon density profile n B (x, y, η s ), we use parametrization as was in Refs. [23,24], Here the longitudinal profiles are parametrized with asymmetric Gaussian profiles, and The normalization factor N n B is chosen such that n B (η s ) = 1. All the model parameters are specified in Table 1 of Ref. [23].

B. Hydrodynamics
To calculate the hydrodynamic evolution numerically we use the open-source hydrodynamic code, MUSIC v3.0 [24][25][26][27][28]. The hydrodynamics equations include energymomentum and baryon number conservation, equation of state p = p( , n B ), and Israel-Stewart-type relaxation equations for the viscous stress tensor. In this work we include only shear viscous corrections, while bulk viscous corrections and baryon diffusion are neglected. We used a lattice QCD based "NEOS-BSQ" equation of state p = p( , n B ) described in Ref. [29]. It smoothly interpolates between the equation of state at high temperature from lattice QCD [30] and hadron resonance gas at low temperature. Higher-order susceptibilities were used to extend this EoS to finite baryon chemical potential as a Taylor expansion. This equation of state explicitly imposes strangeness neutrality, n S = 0, and constrains the ratio of the local net charge density to net baryon density to that of the colliding nuclei, n Q /n B = 0.4. A temperature and chemical potential dependent specific shear viscosity η/s(T, µ B ) was included. The detailed parameterization in given in Fig.4 in Ref. [23]. The hydrodynamic equations are evolved in τ until all computational cells reach energy density below sw . A particlization hypersurface of constant "switching" energy density sw is constructed and its normal four-vectors dσ µ are computed as described in [31]. The value of sw = 0.26 GeV/fm 3 is adjusted to fit bulk observables in [23]. This value results in a good simultaneous fit of pion, kaon, and net proton observables across the range of energies that we study. However, at 39 GeV and above the proton yields at midrapidity are somewhat underestimated. To fine-tune proton yields at midrapidity we take higher sw at highest STAR energies: sw = 0.35 GeV/fm 3 at √ s N N = 39 GeV, sw = 0.45 GeV/fm 3 at √ s N N = 62.4 GeV, and sw = 0.5 GeV/fm 3 at √ s N N = 200 GeV.
In principle the fine-tuning of proton yields can be performed by adjusting other parameters as well, but tuning sw was the simplest solution, because pion and kaon yields are known to be rather insensitive to sw , in contrast to proton and Lambda yields (see Fig. 8 in [29]). Fine-tuning of sw allowed us to reproduce proton and Lambda yields slightly better, at the cost of reproducing anti-proton yields slightly worse.
Performing particlization on the obtained hypersurfaces and allowing the generated hadrons to subsequently re-scatter via a hadronic transport approach, one obtains a good description of the measured pion, kaon, and proton yields, transverse momentum and rapidity spectra, and flow v 2 [23]. The particlization is a standard grandcanonical Cooper-Frye particlization, conducted by the iSS sampler v1.0, which was described and tested in [32] and is available publicly at [33].

C. Transport simulation
The hydrodynamic evolution and particlization are followed by a hadronic afterburner, where particles are allowed to scatter and decay. For this purpose we use the relativistic transport code SMASH [34], version 1.7, in the cascade mode (= no mean field potentials). Transport simulation is initialized from particles at particlization. Then the whole system is propagated from action (collision or decay) to action, using a list of actions sorted by time in the computational frame. A collision is added to the list by the geometrical criterion, πd 2 tr < σ, where d tr is the transverse distance in the center of mass frame of colliding particles, and σ is the total cross section. The collision time is defined as a time of the closest approach, the decay time is randomly drawn from the exponential distribution, which takes time dilation into account. For new particles produced by actions, we search for their subsequent collisions and decays, and merge the found ones into the sorted list of actions. This timestep-less collision finding in SMASH 1.7 is an improvement compared to the timestep-based one in SMASH 1.0, described in the original publication [34], see the comparison and thorough testing in [35]. The end time of our transport simulation is set to 100 fm/c, when the system is already too dilute to sustain even the reactions with very large cross sections, such as deuteron formation, as evident from Fig. 2.
Possible reactions in SMASH include: elastic collisions, resonance formation and decay, 2 → 2 inelastic reactions such as N N → N ∆, N N → N N * , N N → N ∆ * (N * and ∆ * denote all nucleon-and delta-resonances), and strangeness exchange reactions; string formation and immediate decay into multiple hadrons. The SMASH resonance list comprises most of the hadron resonances listed in the Particle Data Group collection [36] with pole mass below 2.6 GeV. The main update relevant for this study since the publication [34] is the high-energy hadronic interactions via string formation, in particular baryonantibaryon annihilations. All the reactions, except the ones with strings, obey the detailed balance principle. The implementation of hadronic interactions in SMASH is described in detail in [34], while [37] is devoted specifically to reactions involving strangeness. Soft string formation and fragmentation are similar to the the UrQMD code [38] and described in detail in [39]. The main difference to the UrQMD implementation is that the Lund fragmentation functions from Pythia 8 [40] are employed for string fragmentation.
The deuteron treatment is the same as in [1], with the same reactions and the same cross sections following the detailed balance principle: πd ↔ πnp, N d ↔ N np, N d ↔N np, πd ↔ N N reactions, elastic πd, N d andN d and all of their CPT-conjugates. The latter produce and destroy anti-deuterons, and one can observe their role in Fig. 3, where anti-deuteron yields are reproduced rather well. The deuteron is modelled as a pointlike particle, as it is done in [1,[41][42][43]. Treating deuterons as pointlike particles is only justified, when the mean free path is at least twice larger than the deuteron size. In our simulation this condition is fulfilled only after t ≈ 10 − 20 fm/c depending on the collision energy. At earlier time our deuterons are not defined as particles and should be understood as correlated nucleon pairs. The reactions πd ↔ πnp, N d ↔ N np,N d ↔N np, are implemented in two steps using a fake d resonance: where X can be a pion, a nucleon, or an anti-nucleon. The d pole mass is taken to be m d = m d + 10 MeV and width is Γ d = 100 MeV, the spin of d is assumed to be 1. The motivation for this width is to have the d lifetime close to the time that proton and neutron spend flying past each other. The deuteron disintegration cross sections reach 200 mb, see Figs. 1-2 of [1]. As a consequence of the detailed balance relations and sharp d spectral function, the πd and N d cross sections are even larger, reaching up to 1500 mb peak values. To find these reactions correctly, the default collision finding cutoff in SMASH is increased to 2000 mb -the same value that was used in [1]. This cutoff is the only change to the SMASH publically available code [44] necessary to reproduce our results related to deuterons.
To reduce the effects of finite range interaction due to the geometric collision criterion, the testparticle method is used in SMASH. Specifically, at particlization the amount of particles is oversampled by factor N test , and, at the same time, all cross sections in SMASH are reduced by factor N test . We have shown previously (see Appendix of [1]) for N test = 10 that this helps to maintain the detailed balance for deuteron reactions with better than 1% precision in an equilibrated box simulation. The effect of N test on our results is discussed further and is shown in Fig. 1.

III. DEUTERON PRODUCTION
Before comparing the deuteron production to experimental data, let us first explore its general features in our simulation. For this let us consider deuterons at our lowest energy (7.7 GeV) at midrapidity, |y| < 0.5. Similarly to [1], we try two scenarios: (i) deuterons are sampled at the particlization (ii) deuterons are not sampled at the 0-10% AuAu, 7.7 GeV, |y| < 0 particlization. In the second scenario all deuterons are created in the hadronic afterburner. The deuteron yield in this case is around 30% smaller than in case (i), and the data favors sampling of deuterons at particlization. In Fig. 1 we also show the effect of the number of testparticles, N test used in the simulation. Larger N test reduces the non-locality of our geometric collision criterion, to which deuterons seem to be rather sensitive because of the large production and destruction cross sections. The deuteron yield increases by almost 30%, when N test is increased from 1 to 10. Further increase of N test does not change the deuteron yield significantly, as shown in Fig. 1, therefore in our simulations we use N test = 10.
In Fig. 1 one can see that the deuteron yield does not change significantly over time in case when deuterons are sampled at particlization. At the particlization hypersurface (which in our case can also be called hadronic freeze-out surface, because stable hadron yields including resonance decay contributions are changing at most by 10% in the hadronic afterburner) the deuteron yield is already close to the measured yield, however its transverse momentum spectra at this point correspond to a hadronic chemical freeze-out temperature. Later the deuteron momentum spectrum changes, but the yield stays approximately constant. We understand it as a result of deuteron being in relative equilibrium with nucleons: the amount of deuterons is determined by the amount of nucleons, which stays approximately constant. The relative equi-librium is kept mostly by πd ↔ πpn reactions. These are the same features of deuteron production that we observed in [1] for PbPb collisions at a much higher 2.76 TeV energy. In Fig. 1 there is a small initial dip in the deuteron yield as a function of time. We observed a similar dip in [1] at 2.76 TeV. We attribute it to the fact, that we do not sample d : although deuterons start in relative equilibrium with protons, it takes time to equilibrate d and d together. The dip is, therefore, an unwanted, but luckily small, artifact of the fake d resonance.
To further illustrate the picture described above, in Fig. 2 we show the reaction rates at midrapidity (rapidity of the reaction was computed from the total momentum of the incoming particles) at 7.7 GeV. One can see that the rates of forward and reverse πd ↔ πpn almost coincide, the differences not exceeding 5%. The N d ↔ N pn reactions occur at several times lower rate than πd ↔ πpn, as one can see in Fig. 2. It means that the πd ↔ πpn reaction is dominant even at the energy as low as 7.7 GeV. In fact, we have observed in the separate pure transport simulation that at 7.7 GeV πd ↔ πpn reactions alone are fast enough to drive deuteron into relative equilibrium with nucleons. Only below 4-5 GeV the N d ↔ N pn reactions become more important than πd ↔ πpn, because for lower energies nucleon abundance increases, while pion abundance decreases. The πd ↔ N N reactions are out of equilibrium, with deuteron destruction dominating at late time. However, their rate is negligible compared to πd ↔ πnp and the integrated rate over time is too small to influence the deuteron yield.
Altogether, above we have established that our simulation behaves in a similar way from 7 up to 200 GeV, as at 2.76 TeV. Apriori it is not obvious that this behaviour should still be consistent with the experimental observables. It is not excluded that some new physical phenomena become important at 7-200 GeV, that did not play role at 2.76 TeV; it could be for example contributions from excited states of 4 He [47] (expected to be small for deuteron, here we just use them as an example) or a vicinity of the critical point.
Already in [1] we have noticed that a reasonable proton description is crucial for meaningful deuteron studies. Therefore, our first step in comparison of our simulation results to experiment is to test, if the hybrid MUSIC + SMASH approach is able to reproduce proton yields and transverse momentum spectra. One caveat in such test is that proton yields measured by NA49 collaboration are corrected for weak decays [46,48], while those measured by many other collaborations are not. Specifically, in STAR [45,49,50], PHENIX [51] (although a correction for Λ decays is available [52]), E895 [53], E802 [54,55], and preliminary HADES data proton yields and spectra [56] are not corrected for the weak decay feeddown. This causes an apparent disagreement between NA49 and STAR data shown in Fig. 3. However, from the left panel of Fig. 3 one can see, that when the weak decays are included, the experimental results both from NA49 and In contrast, πd ↔ pn (right panel) has a much lower rate, it is not equilibrated, and destroys more deuterons than produces. However, the latter reaction contributes negligibly to the deuteron yield because of its low rate.
STAR are described well by our approach. The proton yields with weak decays in Fig. 3 include all possible weak decays into protons, which comprise contributions from Λ, Σ 0 , Σ + , Ξ − , and Ω. This allows a fair comparison to the STAR data, because protons from STAR are truly inclusive with respect to weak decays [45]. To make sure that we correctly account for weak decays, we check the midrapidity yield of the Λ-hyperon. In Fig. 3 we show Λ + Σ 0 yields for a fair comparison, because Σ 0 has a very short lifetime of 7.4 · 10 −20 s and decays with almost 100% branching ratio as Σ 0 → Λγ [36]. This makes Σ 0 experimentally indistinguishable from Λ in heavy ion collisions. The Λ yields in Fig. 3 do not include weak decay contributions from Ξ and Ω baryons, both in our model and in experiment [57]. We also reproduce the proton p T spectra rather well, as one can see in Fig. 4. The p T spectra are characterized comprehensively by the integrated yield dN/dy and mean transverse momentum p T . In Figs. 3 and 4 one can see that they are reproduced in our calculation for protons. A small cusp in proton and deuteron p T at 19.6 GeV originates from the fact that the starting time of hydrodynamics τ 0 is tuned individually at each collision energy. While the p T of Λ is not shown, we have checked that the m T is the same as that for protons within error bars, both in our model and in experiment. To sum up, as demonstrated in Figs. 3 and 4 proton and Λ yields, spectra, and p T are described very well by our approach. Furthermore, one can see in Fig. 3, that the deuteron yields from different experiments [54,56,[63][64][65], as well as spectra and p T are in good agreement with the MUSIC + SMASH simulations. We notice that wherever the proton spectrum in our model deviates from experiment, the deuteron spectrum qualitatively deviates in the same way. For example, at 7.7 GeV, where our description of proton spectra is the least accurate (to improve it the initial longitudinal baryon density profile in the hy-drodynamics has to be considered more carefully) and over(under)shoots the data, the deuteron spectrum also over(under)shoots. Therefore we conjecture that if we tune the model to reproduce proton observables even better, the deuteron description will also improve.
The reactions involving anti-deuterons in SMASH are the CPT-conjugated deuteron reactions with the same cross sections. Consequently, proton and deuteron yields and spectra are connected in the same way as anti-proton and anti-deuteron yields and spectra. Just like for deuterons, pion catalysis πd ↔ πnp is the most important reaction for anti-deuteron production, and it leads to a reasonable description of the anti-deuteron yields, as shown in Fig. 3. At 62.4 and 200 GeV one can see in Fig. 3 that the anti-deuteron yield overshoots in our model. This is mainly because the anti-proton yield overshoots. A better description of protons and anti-protons at 62.4 and 200 GeV will require simultaneous fine-tuning of the initial hydrodynamical baryon density profile (aka baryon stopping) together with the switching energy density sw .
As we have shown, the quality of the model description of proton and deuteron spectra are strongly related. This suggests that ratios of these spectra should be described even better than the yields. Therefore we construct the so called B 2 ratio of the spectra, which plays an important role in coalescence models: In coalescence models this ratio is inversely proportional to an emission volume. The B 2 ratio is known both theoretically and experimentally to grow with p T , consistent with the volume from femtoscopic measurements [66] decreasing with p T . The B 2 ratio is also known to decrease with increasing collision energy, again consistent with the increase of the volume from femtoscopic measurements with the energy. However, two non-trivial   [45,65]. Solid lines correspond to the hydrodynamics + transport simulation, stars are experimental data points. Thin dotted lines in panel (b) are deuteron spectra at particlization. The difference between dotted and solid lines demonstrates the effect of the afterburner for deuterons.
features are present in the STAR measurement of B 2 [65]: (i) the experimentally measured B 2 ( √ s, p T /A = 0.65 GeV) has a broad minimum at 20-60 GeV; (ii) the anti-deuteron Bd 2 ( √ s, p T /A = 0.65 GeV) is smaller than the deuteron B 2 . In our previous work [1] we speculated that the minimum of B 2 might be connected to a switch of the dominant deuteron production mechanism from πpn ↔ πd to N np ↔ N d reactions. However, as shown in this work this conjecture is not supported by our calculations, because the πpn ↔ πd reaction is dominating all the way down to √ s = 7.7 GeV, which is well below the location of the minimum in B 2 . Therefore, let us inspect the behaviour of B 2 closer and suggest a possible explanation, why it exhibits the aforementioned minimum.
We find that the shape of the transverse momentum dependence of B 2 (p T ) is similar for all considered energies and matches the experiment rather well. However, comparing the magnitude of B 2 ( √ s, p T /A = 0.65 GeV) our simulation significantly overestimates the experimental values in the energy range of 20 GeV < ∼ √ s < ∼ 60 GeV, where the minimum is located (see Fig. 5). This discrepancy is surprising given that we reproduce proton, Λ, and deuteron spectra rather well; after all B 2 is nothing but the ratio of the spectra. Investigating this closer we find that it is the weak decays that play a crucial role in this discrepancy. Indeed, if we compute B 2 by dividing the deuteron spectrum from STAR over the weak-decayinclusive proton spectrum from STAR, we reproduce this weak-decay inclusive B 2 rather well, see Fig. 5b. This shows that the contribution from weak decays is much larger in our model than the weak decay correction in the  , and deuteron to proton ratio (panel c) to our simulations. We tend to describe weak-decay-inclusive observables well, in contrast to poorly described weak-decay-corrected observables.
STAR data. This is exacerbated by the fact that in the definition of B 2 the proton spectrum, and thus the weak decay correction, enters in square. Furthermore, comparing the B 2 ratio from STAR with and without weak decays we find that the B 2 ratio with weak-decay-inclusive protons does not exhibit a minimum, whereas that with weak-decay-corrected protons does show minimum (see Fig. 5b). This suggests that the minimum structure in the energy dependence of Bd 2 ( √ s, p T /A = 0.65 GeV) may originate from the weak decay corrections. To further explore the effect of weak decays we consider the the d/p ratio (see Fig. 5c). Again, our model calculation reproduces the weak-decay-inclusive d/p ratio rather well, but not the preliminary weak-decay corrected d/p ratio from STAR [67].
Since our model describes both proton and Λ yields well, one is led to the conjecture that the weak decay corrections to the measured proton yields might be underestimated. Our conjecture is inspired by our model, but it also has some model-independent support from experimental data. First, we notice in Fig. 5 that at the energies where STAR and NA49 data intersect, the weakdecay-corrected B 2 from NA49 is always higher, even though the Pb+Pb collision system is slightly larger than the Au+Au system of the STAR measurement. Since B 2 scales with the inverse size of the system, one would expect the B 2 ratio obtained by NA49 to be below that measured by STAR. If, on the other hand, the weak decay correction to proton yields were larger in the STAR data, it would improve the agreement between STAR and NA49 results for B 2 . Second, one can estimate the weak decay correction in a data-driven way from the recent STAR measurements of strange particle production. Let us consider such an estimate at 39 GeV. The measured yield of Λ + Σ 0 is around dN Λ /dy 10, see Fig. 3c. Let us assume that the Σ + yield, which is not measured, is approximately equal to the Σ 0 and Λ yields. This assumption is mainly motivated by the ther-mal model, where the yields are determined by hadron masses, which are close for Σ + , Σ 0 , and Λ. Additional contribution from Ξ and Ω decays to Λ constitutes around 10-20% of the Λ yield. Taking into account the branching ratios BR(Λ → pπ − ) 0.63 and BR(Σ + → pπ 0 ) 0.52, we obtain the yield of protons from weak decays dN p−weak /dy 9 − 11. Therefore, at 39 GeV at midrapidity around 20 protons per event are prompt and approximately 10 originate from weak decays. The weak decay correction coefficient is thus 30/20 = 1.5. The STAR estimate is roughly 1.15-1.25, both for the B 2 and d/p ratio, see Fig. 6. We note, that the weak decay correction estimate in [65] is not data-driven, but involves the UrQMD model, and may possibly be model dependent. Needles to say, that a data-driven weak decay correction would be beneficial to understand the B 2 ( √ s) behaviour. To sum up, our results as well as data-driven analysis suggest that the observed minimum in B 2 ( √ s) may originate from the weak decay corrections of proton spectra. Moreover, there are both theoretical and experimental indications that these weak decay corrections are underestimated in [65]. If these indications are true, then it has intriguing consequences, which we discuss next.
This work was largely inspired by the study [18], which relates the yields of light nuclei to spatial fluctuations of nucleon densities. Specifically, spatial nucleon density fluctuations are connected to the NtNp N 2 d . This ratio has been measured by STAR recently [67] and it exhibits a peak, which might be a signal of the enhanced nucleon density fluctuations and therefore potentially a critical point. One can see the preliminary STAR data in Fig. 6. The proton yields in this measurement are corrected for weak decays. Suppose our conjecture about weak decay corrections turns out to be correct and the corrections have to be re-evaluated. What will the corrected NtNp N 2 d ratio be? To answer this question quantitatively we extract the ratio of total to non-weak-decay (prompt) pro- tons from STAR data in two ways: from the published B 2 and from preliminary d/p ratio, see Fig. 6. These two ways do not have to give necessarily identical results. Indeed, in case of d/p ratio the relevant observable is the p T -integrated proton yield, while in B 2 it is p Tdifferential yield at proton p T = 0.65 GeV. However, one can clearly see in Fig. 6 that the two ways are in agreement. For convenience we parameterize the STAR weak decay correction x -total proton yield over primordial proton yield without weak decays -as where √ s is the collision energy in GeV. Our correction estimated from the hydro + transport approach (which reproduces the STAR experimental yields) can be approximately parametrized as We have checked that the weak decay correction in our hydro + transport simulation closely resembles that of the thermal model using Thermal-FIST package [68], see Fig. 6. Similarly to our model, the thermal model describes both proton and Λ yields rather well at STAR BES energies . This supports our argument and allows to test it with a much simpler thermal model setup. Indeed, our calculation of weak decay correction for proton yields is in agreement with the thermal model calculation presented by the STAR collaboration in Ref. [45]: we obtain (N p primordial + N p weak )/N p primordial = 1.23 at 7.7 GeV (Fig. 6a) or N p weak /(N p primordial + N p weak ) = 18.7%, consistent with 18% in Ref. [45]; at 39 GeV we obtain (N p primordial + N p weak )/N p primordial = 1.44 or N p weak /(N p primordial + N p weak ) = 30.6 %, again consistently with the 29% STAR has obtained in their estimate, Ref. [45].
After correcting the preliminary STAR data for the NtNp N 2 d ratio by the factor x model /x ST AR , we observe that the peak in the NtNp N 2 d ( √ s) dependence becomes much less pronounced. In fact, the ratio scaled in this way becomes more consistent with a constant value predicted by the coalescence models in absence of non-trivial effects, such as enhanced nucleon density fluctuations. Finally we would like to underline that our analysis is by no means conclusive. We simply want to draw attention to the "technical" issue of weak decay corrections pointing out that it may influence the interpretation of the data in a profound way.
We have already mentioned that the B 2 ratio for antideuterons measured by STAR is smaller than that for deuterons. In Fig. 5a one can see that our simulation produces the opposite trend: the B 2 of anti-deuterons is larger than the B 2 of deuterons. It appears that baryon annihilation, BB → mesons plays a prominent role in that difference. In SMASH this reaction is not balanced: the annihilation BB → mesons is allowed, but the reverse process is not possible. After switching off the BB → mesons annihilation reaction we obtain a lower B 2 (d) (see Fig. 5a), while B 2 (d) remains the same except for 62.4 and 200 GeV, where B 2 (d) is also slightly reduced. As one can see in Fig. 5a, without baryon annihilation B 2 (d) ≈ B 2 (d) within statistical error bars. Thus it seems that the experimentally measured difference of B 2 for deuterons and anti-deuterons may be due to baryon-antibaryon annihilations. However, other possibilities are not excluded, for example, the effect of a weak decay correction, which is larger for anti-protons than for protons.

IV. SUMMARY AND DISCUSSION
In summary, the results of this work and [1], based on hydrodynamics plus transport simulations show that it is possible to reproduce deuteron yields and spectra at energies 7-2760 GeV using pion catalysis reactions πd ↔ πpn with large cross sections obeying the detailed balance principle. One important detail is that the underlying hydrodynamical simulation has to be tuned to reproduce protons and Lambdas well, which we successfully accomplish. The conclusions from [1] regarding deuteron staying in relative equilibrium with nucleons and its yield being almost constant starting from hadronic chemical freeze-out are still valid down to a collision energy of √ s 7 GeV. This also explains the apparent puzzle why the deuteron yield is determined at chemical freeze-out while their spectra correspond to the kinetic freeze-out. At lower energies the deuteron production mechanism is expected to change: N d ↔ N pn reactions will start to dominate.
Analysing the B 2 ratio, which is a ratio of the deuteron spectrum over square of the proton spectrum (Eq. 10), we realized that weak decay corrections to the proton spectrum play a significant role. In particular we found that the observed minimum in B 2 ( √ s) is most likely related to the weak decay corrections. Furthermore, we noticed that weak decay corrections to the proton yield also affect the interpretation of the NtNp N 2 d ratio, which exhibits an unexplained peak as a function of the collision energy. We also demonstrated that the difference between the measured B 2 (d) of anti-deuterons and B 2 (d) of deuterons might be related to BB annihilations.
It seems that we have reached a satisfactory understanding of proton and deuteron production across the STAR beam energy scan energies. The next step is to consider the production of A = 3 nuclei: triton, helium-3, and possibly hypertriton. Unfortunately, the extension of our method to A = 3 nuclei requires an additional fake resonance, t . It can be avoided, however, through implementing 2 ↔ 3 reactions via stochastic rates method. This work is currently in progress.
Another possible extension of the present work is to consider the role of the mean field nuclear potentials on the light nuclei production. Since the light nuclei are mostly formed at the late stage of collision, their yields may be sensitive to the nuclear mean fields at few normal nuclear densities and below. Therefore, it would be interesting to verify to which extent (if any) the nuclear matter liquid-gas phase transition influences the light nuclei yields.