Multiplicity Scaling of Light Nuclei Production in Relativistic Heavy-Ion Collisions

Using the nucleon coalescence model based on kinetic freeze-out nucleons from the 3D MUSIC+UrQMD and the 2D VISHNU hybrid model with a crossover equation of state, we study the multiplicity dependence of deuteron ($d$) and triton ($t$) production from central to peripheral Au+Au collisions at $\sqrt{s_\mathrm{NN}}=$ 7.7, 14.5, 19.6, 27, 39, 62.4 and 200 GeV and Pb+Pb at $\sqrt{s_\mathrm{NN}}=2.76$ TeV, respectively. It is found that the ratio $N_t N_p/N_d^2$ of the proton yield $N_p$, deuteron yield $N_d$ and triton yield $N_t$ exhibits a scaling behavior in its multiplicity dependence, i.e., decreasing monotonically with increasing charged-particle multiplicity. A similar multiplicity scaling of this ratio is also found in the nucleon coalescence calculation based on kinetic freeze-out nucleons from a multiphase transport (AMPT) model. The scaling behavior of $N_t N_p/N_d^2$ can be naturally explained by the interplay between the sizes of light nuclei and the nucleon emission source. We further argue that the multiplicity scaling of $N_t N_p/N_d^2$ can be used to validate the production mechanism of light nuclei, and the collision energy dependence of this yield ratio can further serve as a baseline in the search for the QCD critical point in relativistic heavy-ion collisions.


I. INTRODUCTION
Light nuclei production in high energy heavy-ion collisions have been extensively studied both experimentally and theoretically [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Besides its intrinsic interest, studying light nuclei production provides the possibility to probe the local baryon density and the space-time structure of the emission source in relativistic heavy-ion collisions [16][17][18]. However, the production mechanism for light nuclei in heavy-ion collisions is still under debate. On the one hand, the statistical model, which assumes that light nuclei, like hadrons, are produced from a thermally and chemically equilibrated source, provides a good description of measured yields of light (anti-)nuclei in central Pb+Pb collisions at the Large Hadron Collider (LHC) [19][20][21]. On the other hand, the coalescence model, which assumes that light nuclei are formed from the recombination of kinetically freeze-out protons and neutrons [22][23][24][25], can also successfully describe the transverse momentum spectra and the collective flow of light nuclei in heavy-ion collisions [26][27][28][29][30][31][32][33]. In between these two extreme scenarios for light nuclei production in heavy-ion collisions is the transport model, which aims to study how light nuclei evolve during the hadronic evolution by including their production and annihilation, such as the processes π+p+n ↔ d+π for the deuteron [34,35]. In this approach, it is assumed that the deuteron can exist in hot dense hadronic matter, although its temperature of more than 100 MeV is significantly higher than the deuteron binding energy of 2.2 MeV.
Recently, more experimental data on light nuclei production in relativistic heavy-ion collisions have become available. For example, the STAR Collaboration at Relativistic Heavy-Ion Collider (RHIC) and the ALICE Col-laboration at the LHC have collected a wealth of data on light nuclei, such as (anti-)deuteron (d, d), (anti-)triton (t, t) and (anti-)helium-3 ( 3H e, 3 He), in Au+Au collisions from 7.7 GeV to 200 GeV [36][37][38] and Pb+Pb collisions at 2.76 TeV [10,39,40], respectively. A noteworthy result from these experiments is the observation that the yield ratio of triton, deuteron and proton, N t N p /N 2 d , in central Au+Au collisions shows a non-monotonic energy dependence with a peak around √ s NN =20 GeV [37,38]. A physically interesting explanation for this non-monotonic behavior is that it is due to the density fluctuations indeuced by the change in the quark-gluon plasma (QGP) to hadronic matter phase transition from a crossover at small baryon chemical potential to a first-order transition at large baryon chemical potential and the associated critical point [41][42][43][44][45][46][47]. Another interesting observation from these experiments is the suppressed production of light nuclei in peripheral Au+Au collisions or in collisions of small systems. A possible explanation for this suppression is due to the non-negligible sizes of light nuclei compared to the size of the nucleon emission source in such collisions [17,48]. Alternatively, introducing the canonical suppression in the thermal model can also lead to the suppression of light nuclei production in collisions of small systems [20,21].
In this paper, we investigate the multiplicity dependence of deuteron and triton production from peripheral to central Au+Au collisions at √ s NN = 7.7 -200 GeV from the Beam Energy Scan (BES) program at RHIC by using the nucleon coalescence model with the needed phase-space distribution of nucleons generated by the 3D MUSIC+UrQMD hybrid model [49] and also the AMPT model [50]. We further use the 2D VISHNU hybrid model [51] to study deuteron and triton production in Pb+Pb collisions at 2.76 TeV. Specifically, nucleons at the kinetic freeze-out of the expanding hadronic matter in heavy-ion collisions are obtained from these models with a smooth crossover transition from the QGP or the partonic matter to the hadronic matter without the presence of any density fluctuations. Since the MUSIC+UrQMD hybrid model used in the present study has been shown to nicely reproduce the measured yield and collective flow of various hadrons in heavy-ion collisions at both RHIC and LHC energies [49,[51][52][53][54][55][56][57], although with a different initial state, it provides a robust baseline for the study of light nuclei production in the RHIC BES program. We find from our study that the yield ratio N t N p /N 2 d in collisions from √ s NN =7.7 GeV to 2.76 TeV is essentially independent of the collision energy and decreases with increasing charged-particle multiplicity. Such a multiplicity scaling of this ratio contradicts to the prediction of the thermal model, which shows an opposite multiplicity dependence, i.e., an increasing with increasing charged-particle multiplicity [20,21].
This paper is organized as follows: Section II briefly introduces the nucleon coalescence model for deuteron and triton production, and the MUSIC+UrQMD hybrid model for the description of heavy-ion collision dynamics at RHIC. Section III presents and discusses results from the MUSIC+UrQMD hybrid model on the multiplicity dependence of the transverse momentum spectra and the rapidity distributions of protons, deuterons and tritons, the coalescence parameters for deuteron and triton production, and the yield ratio N t N p /N 2 d in Au+Au collisions at RHIC BES energies. The obtained charged particle multiplicity dependence of the yield ratio N t N p /N 2 d is then compared with that from the AMPT+coalescence model and also extended to Pb+Pb collisions at LHC by using the VISHNU+coalescence model. Section IV concludes this paper with a summary.

II. THE THEORETICAL FRAMEWORK
A. The coalescence model for light nuclei production In the coalescence model [22][23][24][25], the number of a light nucleus of atomic number A and consisting of Z protons and N neutrons (A = Z+N ) is given by the overlap of the Wigner function f A of the nucleus with the phase-space distributions f p (x i , p i , t) of protons and f n (x j , p j , t) of neutrons [24,25], where g A = (2J A + 1)/2 A is the statistical factor for A spin 1/2 nucleons to form a nucleus of angular momentum J A . The p µ d 3 σ µ denotes the differential hyper-surface of emitted nucleons such that p µ d 3 σ µ . The x i and p i are the coordinate and momentum of i-th nucleon in the frame of the nucleon emission source, while x i and p i are those in the rest frame of the produced nucleus, which can be obtained from the coordinate x i and momentum p i by the Lorentz transformation. For the proton (neutron) phase space distributions f p/n (x, p, t), they are obtained in the present study from the positions, momenta and times of protons and neutrons at their last scatterings, i.e., their kinetic freeze out, in the UrQMD or AMPT model. Following Ref. [24], we take the Wigner functions of deuteron and triton to be Gaussian and use the same proton and neutron masses. For the deuteron, its Wigner function is then with the relative coordinate ρ and the relative momentum p ρ defined as For the triton, its Wigner function is with the additional relative coordinate λ and relative momentum p λ defined as The width parameter σ d in Eq. (2) is related to the root-mean-square matter radius r d of deuteron by σ d = 2 √ 3 r d . Similarly, the width parameter σ t in Eq. (4) can be related to the root-mean-square radius r t of triton by σ t = r t . For the values of r d and r t , they are taken to be 1.96 fm and 1.59 fm, respectively, from experimental measurements [58]. We note that because of the very small yields of deuterons and tritons in heavy-ion collisions at the RHIC BES energies, protons used in the coalescence processes have negligible effects on the final proton spectra.
It is worth to point out that there is an ambiguity of a prefactor in Eq.(1) for the three-body coalescence production of triton via p+n+n → t. In Refs. [23][24][25], where nucleons are treated as classical particles, a prefactor of 1/2 is introduced to avoid double counting the neutrons in the production of tritons. This factor is, however, absent in Ref. [16]. According to Refs. [41,45], without this factor the triton yield in the coalescence model using nucleons from a thermally and chemically equilibrated emission source is identical to that in the thermal model when the triton binding energy is neglected. In this work, we only consider the three-body coalescence channel for triton production and do not include a prefactor of 1/2 in Eq.(1) as in Refs. [16,41,45]. Compared to our previous study on triton production from the three-body channel in Ref. [32], which includes the prefactor of 1/2, both the triton yield and the ratio N t N p /N 2 d at high multiplicity in this work will thus be about a factor of two larger.

B. The MUSIC +UrQMD hybrid model for heavy-ion collisions
For the description of heavy-ion collision dynamics, we employ the MUSIC+UrQMD hybrid model [49,54,[61][62][63]. The (3+1)D viscous hydrodynamics in MUSIC, which conserves both the energy-momentum and baryon number of produced QGP, is developed for describing the collective dynamics of QGP and the soft hadrons produced from the hadonization hypersurface. At the RHIC BES energies, this hybrid model uses a smooth initial condition with the net baryon density taken from an extended 3D Glauber model as the input for the subsequent hydrodynamic evolution [49]. For the present study, we use a crossover type of equation of state (EoS) (NEOS-BQS) with the strangeness neutrality condition of vanishing net strangeness density, n s = 0, and the net electric charge-to-baryon density ratio n Q = 0.4n B [64]. We use such an EoS without a QCD critical point because the aim of our study is to provide a reliable baseline calculation without any critical effects on light nuclei production in heavy-ion collisions at energies from the RHIC BES program. Following Ref. [49], we include in the hydrodynamic evolution a temperature and baryon chemical potential dependent specific shear viscosity η/s. After the hydrodynamic evolution, we convert fluid cells to hadrons when their energy densities drop below the switching energy density e sw =0.26 GeV/fm 3 . For the evolution and decoupling of the resulting hadronic matter, the MUSIC code switches to the microscopic UrQMD hadron cascade model [65][66][67]. Finally, we obtain the phase-space distributions of kinetic freeze-out nucleons for the coalescence model calculations.
We would like to emphasize that the MUSIC+UrQMD hybrid model employed in this paper has achieved a quantitative description of soft particle production, including the p T -spectra and collective flow, in the central and peripheral Au+Au collisions at GeV, as demonstrated in Ref. [49]. This indicates that the MUSIC+UrQMD hybrid model can provide, for the coalescence model calculations, the proper phase-space distributions of nucleons at the kinetic freeze-out, i.e., their last scatterings in the center-of-mass frame of the hadronic matter or the fireball frame. In the coalescence calculations, the coalescing nucleons are boosted to their center of mass frame with those of earlier times further propagated to the time of the latest one. These nucleons then have the probability to form a light nucleus that is given by the product of their spin statistical factor to form the nucleus and the Wigner function of the nucleus evaluated at their equal-time spatial and momentum separation distances. The position and momentum of formed nucleus are finally Lorentz transformed back to the fireball frame.

III. RESULTS
In this section, we study the transverse momentum spectra, the rapidity distribution dN/dy, the mean transverse momentum p T as well as the multiplicity dependence of the coalescence parameters GeV. It is seen that the MUSIC+UrQMD hybrid model can quantitatively describe measured proton spectra below 2.5 GeV but slightly underestimates the data at higher p T , where the quark recombination contribution has been shown to become increasingly important [69][70][71][72][73][74][75] 1 With the phase-space distributions of protons and neutrons at kinetic freeze-out from the MUSIC+UrQMD hy-      [37] and from the STAR and PHENIX Collaborations for protons [59,60]. The lower panels give the ratio of the model results to the experimental data.
brid model, we calculate the spectra of deuterons and tritons within the framework of nucleon coalescence model. As shown by blue lines in Fig. 1, without any free parameters, the coalescence model calculations can reasonably describe the p T -spectra of deuterons measured by the STAR Collaboration at 0-10%, 10-20%, 20-40% and 40-60% centrality bins in Au+Au collisions at √ s NN = 7.7 − 200 GeV. The ratio of the theoretical result to the experimental data, denoted by Theory/Exp in Fig. 1, is mostly between 0.5 and 1.5. Figure 2 shows the rapidity distribution dN/dy of protons, deuterons, and tritons in 0-10% Au+Au collisions at √ s NN = 7.7−200 GeV. It is seen that our model calculations give an excellent description of the STAR data at mid-rapidity. For the proton yield at mid-rapidity, it is determined in the (3+1)D MUSIC model by the interplay between the initial baryon stopping and thermal production at chemical freeze-out. At low collision energies, the initial baryon stopping and baryon current evolution in the hydrodynamic phase are more important than ther-mal production at particalization, resulting in a larger proton yield as observed in experiments. The calculated dN/dy of deuterons and tritons also show a similar trend with respect to the change in the collision energy, which are again consistent with the STAR data. For the predicted shapes of the dN/dy of protons, deuterons and tritons at √ s NN =7.7-200 GeV, they become wider with increasing collision energy, which can be tested in upcoming experimental measurements.
In the upper panel (a) of Fig. 3, we show the atomic mass (m A , A = 1, 2, 3 for proton, deuteron and triton) dependence of the rapidity distribution per degree of freedom, dN/dy 2J+1 , in 0-10% Au+Au collisions at √ s N N =7.  GeV. An exponentially decreasing dN/dy 2J+1 with the atomic mass is observed as expected from the nucleon coalescence model, in which the binding energies of produced light nuclei are ignored. The lower panel (b) shows the dependence of the mean transverse momentum p T on the atomic mass m A . The p T increases with increas-  ing m A because light nuclei are produced in the coalescence model from nucleons close in phase space, thus enhancing their momenta. Besides, larger values of p T are observed at higher collision energies due to the stronger radial flow, which pushes particles to larger p T .

B. Coalescence parameters
In interpreting light nuclei production from nuclear reactions, one usually expresses the yield d 3 N A /d 3 p A of a nucleus with the mass number A = Z + N and momentum p A in terms of the yield d 3 N p /d 3 p p of protons with momentum p p and the yield d 3 N n /d 3 p n of neutrons with momentum p n in terms of the coalescence parameter B A as follows: where E p,n are the proton and neutron energies. The coalescence parameter B A (A=2,3) contains information on the freeze-out properties of the nucleon emission source, i.e., B A ∝ V 1−A eff with V eff being the effective volume of the nucleon emission source [16][17][18]76]. Figure 4 shows the multiplicity dependence of the coalescence parameters B 2 (d) and and B 3 (t) from our model calculations decrease monotonically with increasing multiplicity, and this is because the size of the nucleon emission source increases monotonically with increasing multiplicity. It is seen that the B 3 (t) decreases more rapidly than B 2 (d) as a function of multiplicity. Specifically, the values of B 3 (t) are very close to B 2 (d) at low multiplicities, but they start to deviate when dN ch /dη > 200, which can be attributed to the different deuteron and triton sizes [76]. This is also consistent with the multiplicity dependence of the yield ratio N t N p /N 2 d shown in the next subsection. The slightly smaller coalescence parameters at p T /A=0.65 GeV/c than those at p T /A=1.05 GeV/c are consistent with the decreasing correlation lengths (HBT radii) with mean pair transverse momentum extracted in the STAR expeiments [77]. For a fixed multiplicity, the coalescence parameter is systematically larger at a higher collision energy. For example, at dN ch /dη ∼ 250, the B 2 (d) in collisions at √ s NN =200 GeV is clearly greater than that at √ s NN =14.5 GeV. This is attributed to the stronger collective expansion at the higher collision energy, which leads to a decrease of the effective volume at kinetic freeze-out because of the focusing effect due to the larger transverse flow, although the total volume is larger at higher collision energy [77].
Although the coalescence parameters B 2 (d) and B 3 (t) are expected to depend inversely on the effective volume of the nucleon emission source, their dependence on the charged particle multiplicity, which is proportional to the true volume of the emission source, is nontrivial because of the complicated relation between the effective and true volumes of the emission source. For the coalescence parameters B 2 (d) and B 3 (t) from the MUSIC+UrQMD hybrid model shown in Fig. 4, it is found that their scaled values B 2 (d)(dN ch /dη) 0.765 and B 3 (t)(dN ch /dη) 0.765 then depend only weakly on the charged particle multiplicity.
C. The yield ratio NtNp/N 2 d We have also studied the charged particle multiplicity dependence of the yield ratio N t N p /N 2 d in Au+Au at √ s NN =7.7-200 GeV from the MUSIC+UrQMD hybrid model, and it is shown in Fig. 5 by black open symbols. The yield ratio N t N p /N 2 d is seen to decrease monotonically with increasing charged particle multiplicity, and the same behavior is seen for different collision centralities and energies as they almost lie on the same curve. This scaling behavior in the multiplicity dependence indicates that the yield ratio N t N p /N 2 d in relativistic heavyion collisions is insensitive to the baryon density and collective flow effects.
The above scaling behavior can be qualitatively understood by considering a thermally equilibrated spherical Gaussian nucleon emission source with a width parameter or radius R at certain temperature and nucleon chemical potentials, i.e., ∝ e −r 2 /(2R 2 ) . The yield ratio N t N p /N 2 d in the coalescence model can then be obtained analytically [48], and it is given by Since (2/ √ 3)r d = 2.26 fm is greater than r t = 1.59 fm, the yield ratio thus decreases with increasing source size or the particle multiplicity in heavy-ion collisions as one expects R ∝ (dN ch /dη) 1/3 for a static source. Eq. (7) also shows that this behavior is independent of the temperature and nucleon chemical potentials of the emission source. This result is similar to the multiplicity scaling behavior of the yield ratio N t N p /N 2 d seen in the coalescence model calculations based on kinetic freeze-out nucleons from the MUSIC+UrQMD hybrid model. Eq. (7) further shows that the ratio N t N p /N 2 d has an asymptotic value of 4/9 when R r d , r t , which is again similar to the value shown in Fig. 5. However, fitting the ratio N t N p /N 2 d according to Eq.(7) using the relation R = α(dN ch /dη) β fm by treating α and β as parameters leads to an unrealistic large radius for a given charged particle multiplicity. This indicates the inadequacy of using a static Gaussian form to describe the emission source from the hybrid and AMPT models, which includes dynamic effects such as the space and momentum correlations. As a possible improvement to the modeling of the emission source, we multiply Eq.(7) by a factor p 0 . A good description of the MUSIC+UrQMD hybrid model results in Fig. 5 can then be obtained with p 0 = 0.683 and a more realistic relation of R = 0.547(dN ch /dη) 0.331 as shown by the blue solid line in the figure. We note that if we had taken the nucleon emission source to be uniform in space as assumed in Refs. [41,42], the yield ratio N t N p /N 2 d would have a smaller asymptotic value of (4/9) × (3/4) 3/2 = 1 2 √ 3 . Therefore, including the surface diffuseness of the nucleon emission source, which is expected in realistic heavy-ion collisions, enhances this from the MUSIC+UrQMD, VISHNU and AMPT models, as well as the statistical model that includes only stable nuclei [21]. Also shown by the blue solid line is the result from fitting the results of MUSIC+UrQMD hybrid model by Eq. (7) with the mutiplication factor p0 = 0.683 and using the relation R = 0.547(dN ch /dη) 0.331 fm between the radius R and the charged particle multiplicity dN ch /dη.
value. Also, the size effects discussed here will be even stronger for hypernuclei due to their larger sizes [48].
Similar results are obtained from the nucleon coalescence model using kinetic freeze-out nucleons from the AMPT model as shown by pink open symbols in Fig. 5. The AMPT model is a multiphase transport model that uses initial conditions from the HIJING model [78,79], the ZPC model [80] for the partonic cascade, and the ART model [81] for the hadronic transport as well as a spatial quark coalescence model to convert kinetic freezeout quarks and anti-quarks to initial hadrons.
Also shown by the filled square in Fig. 5 [10], which is much larger than the theoretical value of about 0.33 from Refs. [31,51]. Although the uncertainty in the ALICE data is still large, the almost factor of three difference between the coalescence model calculation and the experimental data at LHC needs to be understood in the future.
We note that the thermal model gives an opposite behavior in the multiplicity dependence of N t N p /N 2 d , i.e., it increases with increasing multiplicity as shown by the red solid line in Fig. 5 [20,21]. According to Eq.(7), the larger N t N p /N 2 d at lower multiplicity is due to the larger deuteron than triton radius. This enhancement in the value of N t N p /N 2 d would disappear if the deuteron and triton radii were similar. Indeed, changing the triton radius to that of deuteron in the coalescence model calculations using nucleons from the AMPT model, we have obtained an almost constant value of about 0.33 for all particle multiplicities. This result strongly indicates the importance of the size effect on coalescence production of light nuclei in heavy-ion collisions. The upcoming measurements of N t N p /N 2 d as a function of multiplicity would thus help discriminate different production mechanisms of light nuclei in heavy-ion collisions. We would like to point out that neither the hydrodynamic model nor the AMPT model used in our study includes any dynamical density fluctuations related to the critical point or first-order phase transition of the QCD matter. Therefore, our results can serve as the baseline predictions for the yield ratio of light nuclei in the search for the possible QCD critical point from experimental beam energy scan of heavy-ion collisions.

IV. SUMMARY
In this paper, we have studied the multiplicity dependence of light nuclei production in Au+Au collisions at √ s NN = 7.7 -200 GeV by using the nucleon coalescence model. With the phase-space distributions of protons and neutrons at kinetic freeze-out taken from the MU-SIC+UrQMD model with a smooth crossover transition from the QGP to hadronic matter, the theoretical results nicely reproduce the measured p T -spectra of protons and deuterons in 0-10%, 10-20%, 20-40% and 40-60% Au+Au collisions at √ s NN = 7.7 − 200 GeV. The coalescence parameters B 2 (d) and √ B 3 (t) for deuteron and triton production, respectively, are found to decrease monotonically with increasing multiplicity. Using also kinetic freeze-out nucleons from the VISHNU hybrid model with a crossover transition for the coalescence production of deuterons and tritons in Pb+Pb collision at √ s NN = 2.76 TeV, the yield ratio N t N p /N 2 d in heavy-ion collisions from 7.7 GeV to 2.76 TeV is found to exhibit a scaling behavior in its multiplicity dependence, i.e., decreasing monotonically with increasing charged-particle multiplicity, but not much on the collision energy and centrality. This scaling behavior is further seen in calculations that use the kinetic freeze-out nucleons from the AMPT model for deuteron and triton production in heavy-ion collisions at RHIC energies. The predicted multiplicity scal-ing of the yield ratio N t N p /N 2 d can be used to validate the production mechanism of light nuclei and to extract their sizes in relativistic heavy-ion collisions. As to the collision energy dependence of this yield ratio predicted by the MUSIC+UrQMD hybrid model and the AMPT model with a crossover QGP to hadronic matter transition as used in the present study, both give essentially a constant value [32,45]. Although these models cannot describe the non-monotonic collision energy dependence in the preliminary STAR data [38], which has been argued to be related to dynamical density fluctuations from the critical point and first-order phase transitions of QCD [41,42,83], results from these models can serve as a baseline against which one can compare the experimental data from the beam energy scan of heavy-ion collisions to search for the critical point in the QCD phase diagram.