Light nuclei production as a probe of the QCD phase diagram

It is generally believed that the quark-hadron transition at small values of baryon chemical potentials $\mu_B$ is a crossover but changes to a first-order phase transition with an associated critical endpoint (CEP) as $\mu_B$ increases. Such a $\mu_B$-dependent quark-hadron transition is expected to result in a double-peak structure in the collision energy dependence of the baryon density fluctuation in heavy-ion collisions with one at lower energy due to the spinodal instability during the first-order phase transition and another at higher energy due to the critical fluctuations in the vicinity of the CEP. By analyzing the data on the $p$, d and $^3$H yields in central heavy-ion collisions within the coalescence model for light nuclei production, we find that the relative neutron density fluctuation $\Delta \rho_n=\langle(\delta \rho_n)^2\rangle/\langle \rho_n\rangle^2$ at kinetic freeze-out indeed displays a clear peak at $\sqrt{s_{NN}}=8.8$ GeV and a possible strong re-enhancement at $\sqrt{s_{NN}}=4.86$ GeV. Our findings thus provide a strong support for the existence of a first-order phase transition at large $\mu_B$ and its critical endpoint at a smaller $\mu_B$ in the temperature versus baryon chemical potential plane of the QCD phase diagram.

1. Introduction.-Understanding the phase diagram of strongly interacting matter is of fundamental importance in nuclear physics, astrophysics and cosmology. Lattice quantum chromodynamics (LQCD) calculations [1] and various effective models [2][3][4] have suggested that the transition between quark-gluon plasma (QGP) and hadronic matter is a smooth crossover at vanishing baryon chemical potential (µ B ), but likely changes to a first-order phase transition at large µ B , with an associated critical endpoint (CEP) or a tricritical endpoint [5]. Experimentally, heavy-ion collisions provide a unique tool to probe the QCD phase diagram [6][7][8][9][10]. In particular, to search for the CEP and locate the phase boundary in the QCD phase diagram is the main motivation for heavy-ion collision experiments being carried out in the Beam Energy Scan (BES) program at the Relativistic Heavy Ion collider (RHIC) as well as those planned at the future Facility for Antiproton and Ion Research (FAIR), Nuclotron-based Ion Collider facility (NICA) and the SPS Heavy Ion and Neutrino Experiment (NA61/SHINE).
In heavy-ion collisions, the created matter could develop strong baryon density fluctuations when its evolution trajectory in the QCD phase diagram passes across * sunkaijia@sjtu.edu.cn † Corresponding author: lwchen@sjtu.edu.cn ‡ ko@comp.tamu.edu § pujiephy@sjtu.edu.cn ¶ xzb@bnl.gov the first-order phase transition line (due to the spinodal instability [11][12][13][14][15][16]) or approaches the CEP (because of a rapid increase of correlation length in the critical region [3,17]). In particular, for heavy-ion collisions at lower energies when the system enters the region of firstorder phase transition, it is expected that the density fluctuation should reach a maximum at a collision energy √ s M at which the system stays the longest time inside the unstable spinodal region [13]. With increasing collision energy, a second maximum in the density fluctuation is expected at a collision energy √ s E when the first-order phase transition ends and the CEP is reached [3,17]. This double-peak structure of the density fluctuation as a function of the collision energy √ s is plotted schematically in Fig. 1 where the corresponding phase regions in the QCD phase diagram are also indicated.
To extract the density fluctuations in heavy-ion collisions from experimental observables is a challenging task as only the particle momentum distributions are generally measured. Because of the rapid expansion of the fireball formed in heavy-ion collisions, the enhanced density fluctuations caused by the spinodal instability or critical fluctuation could, however, survive the finalstate interactions and affect observables that are sensitive to nucleon density fluctuations and correlations at kinetic freeze-out (at which particles cease interacting). Besides the Hanbury-Brown-Twiss (HBT) interferometry of identical particles, which can provide information on the space-time structure of the particle emission source [18,19] and thus the effect of density fluctuations at hadronization [7], these observables also include light nuclei production via nucleon coalescence [20][21][22][23][24][25][26][27][28]. Indeed, we have recently shown [26,27] that the yield ratio O p-d-t = N p N3 H /N 2 d of produced proton (N p ), deuteron (N d ), and triton (N3 H ) depends on the relative neutron density fluctuation ∆n as well as the neutron and proton density correlation C np at kinetic freeze-out. Assuming the ratio of C np to ∆n to be independent of the collision energy, we have found [27] from comparing with the measured yield ratio O p-d-t in Pb+Pb collisions by the NA49 Collaboration [29] at the CERN Super Proton Synchrotron (SPS) that the ∆n as a function of √ s shows a possible peak structure.
In this Letter, we argue that the phase-space volume of the fireball produced in heavy-ion collisions at kinetic freeze-out is related to the entropy per nucleon and can be determined from the temperature and volume at chemical freeze-out (at which chemical equilibrium is reached and particles ratios are fixed). With the well determined chemical freeze-out temperature and volume from the statistical model fit to available experimental data, we can then uniquely determine the collision energy dependence of C np and ∆n as well as their ratio in heavy-ion collisions from the yields of p, d and 3 H ( 3 He), instead of assuming that this ratio is independent of the collision energy as in Ref. [27]. By analyzing the data in central collisions of Au+Au [30,31] measured at the Brookhaven Alternating Gradient Synchrotron (AGS) and Pb+Pb measured by the NA49 Collaboration [29] at SPS, we find that the ∆n displays a clear peak at 2. Coalescence production of light nuclei and density fluctuations and correlations.-We start by calculat-ing the deuteron yield (dN/dy at midrapidity) through a quantum mechanical treatment in full phase space. Assuming nucleons (with mass m, proton number N p and neutron number N n ) are emitted from an isotropic thermalized fireball with temperature T and volume V and uniformly distributed in space, one can obtain the deuteron yield as which is exactly the same as the formula in Ref. [27] obtained by assuming Bjorken boost invariance for the expanding fireball. For non-uniform distributions, the deuteron yield becomes where n(x) and n p (x) denote the neutron and proton number densities in space, respectively. Similarly, the 3 H yield is given by To account for the density fluctuations, the neutron and proton densities in the emission source can be expressed as where · = 1 V dx denotes the average over space and δn(x) (δn p (x)) with δn = 0 ( δn p = 0) denotes the neutron (proton) density fluctuation from its average value n ( n p ). By neglecting the term (δn) 2 δn p , we can rewrite Eqs. (2) and (3) as where C np = δnδn p /( n n p ) and ∆n = (δn) 2 / n 2 characterize the neutron and proton density correlation and the relative neutron density fluctuation, respectively. We would like to emphasize that both C np and ∆n are defined to be dimensionless to eliminate possible effects due to the collision energy dependence of average nucleon densities. It is seen from Eqs. (6) and (7) that the deuteron yield depends on C np but not ∆n, while the triton yield depends on both. In addition, one sees from Eqs. (6) and (7) that C np and ∆n can be uniquely determined to be as defined earlier in the Introduction. The ratio R np = N p /N n = n p / n denotes the isospin factor and can be estimated from the pion yield ratio as R np = (π + /π − ) 1/2 , and V ph = (2πmT ) 3/2 V is the effective phase-space volume of nucleons in the fireball at kinetic freeze-out [26].
The phase-space volume V ph is directly related to the entropy per nucleon (S/N ), which is given by the Sackur-Tetrode equation [32] S/N = 5/2 + ln(V ph /N ) in the non-relativistic Boltzmann approximation. If we assume T 3/2 V = λT 3/2 ch V ch with the chemical freeze-out temperature T ch and volume V ch , λ is then unity if a gas of constant number of nucleons expands isentropically after chemical freeze-out. A recent microscopic transport model study [33] has shown, however, that it is the entropy per particle that remains a constant after chemical freeze-out in heavy-ion collisions, which is dominated by pion production if the collision energy is high, and that all particle ratios remain essentially unchanged from the chemical to the kinetic freeze-out in these collisions as assumed in the statistical model for particle production. Based on the time evolution of T and V after chemical freeze-out from this study [33], we find that the λ value ranges from 1.5 at 7.7 GeV to 1.7 at 200 GeV, suggesting that although T 3/2 V increases as fireball expands, λ has a very weak dependence on the collision energy. Therefore, we can uniquely determine the values of C np and ∆n by using a constant λ = 1.6.
It should be mentioned that light nuclei (like d and 3 H) are formed from nucleons in a very restricted phasespace volume of ∆x ∼2 fm and ∆p ∼0.1 GeV, and their production in heavy-ion collisions can thus be used to probe the local nucleon density fluctuations. In contrast to the studies [9,[34][35][36] that are mainly focusing on the event-by-event fluctuations of conserved quantities within a specific window in momentum space, the ∆n in Eq. (9) has the advantage that it directly measures the spacial density fluctuation. Mapping the fluctuations in momentum space to those in coordinate space is highly non-trivial, especially at SPS and AGS energies where the longitudinal boost invariance is not well satisfied.
With C np and ∆n, we can further investigate the fluctuation of neutron and proton density difference ( (δn − δn p ) 2 ) which is closely related to the isospin density fluctuation as neutron and proton have opposite isospin quantum numbers. While the fluctuations of baryonic (B), electric (Q) and strange (S) charges have been extensively investigated [9,34], the fluctuation of isospin is less studied in heavy-ion collisions. Recently it has been shown [37] that the isospin effect at BES energies can dramatically change the critical behavior of baryon and charge number fluctuations, suggesting that the isospin effects would be strong in heavy-ion collisions at SPS and AGS energies. Defining the isospin density fluctuation ∆n I as (10) one sees that a negative C np can increase ∆n I . Although the isospin in high energy heavy-ion collisions is mostly carried by pions, the ∆n I defined in terms of nucleons may still carry important information on the isospin density fluctuation due to the frequent interactions between pions and nucleons.
3. Results and discussions.-Shown in Fig. 2 is the collision energy dependence of C np , ∆n and ∆n I . The results are obtained from Eqs. (8)-(10) using the data on the p, d and 3 H yields together with the yield ratio π + /π − from central collisions of Pb+Pb at √ s N N = 6.3 GeV, 7.6 GeV, 8.8 GeV, 12.3 GeV and 17.3 GeV measured by the NA49 Collaboration at the CERN SPS [29,38,39]. It is seen from Fig. 2 (a) that the C np shows a non-monotonic behavior with a valley located at √ s N N =8.8 GeV. Besides, the extracted C np at SPS energies are all negative, indicating a strong negative correlation between neutron and proton densities. From Fig. 2 (b), one sees that ∆n and ∆n I show very similar non-monotonic behaviors with peaks also located at √ s N N = 8.8 GeV. The obtained peak structure of ∆n is consistent with that obtained in Ref. [27], where the ratio α = C np /∆n is fixed at some constant values (e.g., Although the α value only has a weakly non-monotonic dependence on the collision energy, justifying its constancy assumed in Ref. [27], its magnitude is significantly larger, leading to a much more pronounced peak of ∆n at √ s N N = 8.8 GeV than that found in Ref. [27]. The similar behavior of ∆n and ∆n I with both peaked at √ s N N = 8.8 GeV could be due to the same underlying physics, namely, the critical fluctuation in the vicinity of CEP. According to the universality of critical behavior, the singular parts of ∆n and ∆n I in the second-order phase transition both scale with the correlation length l as l 2−η , where η denotes the critical exponent and is zero in the mean-field approximation, and diverge at the CEP in the QCD phase diagram. Due to the effects of critical slowing down [17] and dynamical expansion in heavy-ion collisions, only modest but simi- From the parametrization in Ref. [40] for the chemical freeze-out conditions based on the statistical model fit to available experimental data of hadron yields, the temperature at √ s N N = 8.8 GeV is estimated to be T CEP ∼ 144 MeV with a corresponding baryon chemical potential µ CEP B ∼ 385 MeV, which is close to the predicted CEP from the LQCD [1], the Dyson-Schwinger equation (DSE) [41] and the hadronic bootstrap approach [42], but is much larger than that (∼ 95 MeV) inferred from a finite size scaling (FSS) analysis of twopion correlations [43]. It is very likely that the critical region is reached within the energy region √ s N N = 6.3 ∼ 17.3 GeV. With the temperature and baryon chemical potential of about 131 ∼ 159 MeV and 481 ∼ 229 MeV [40], respectively, as determined from the statistical model, one can estimate the size of critical region to be ∆T /T CEP ≈ 0.1 and ∆µ B ≈ 0.1 GeV, which is consistent with the effective model calculations [4]. The critical exponents, however, can not be determined from present data.
In the scenario depicted in Fig. 1, signals from the firstorder phase transition and the CEP are connected, because the CEP stays well on the top of the spinodal unstable region in the T -µ B plane of the QCD phase diagram.
Indeed, it has been estimated [12] that the temperature (T M ) at point 'M' is related to T CEP by T M /T CEP ≈ 1 3 − 1 2 . From the parametrization in Ref. [40] for the chemical freeze-out conditions, one can roughly estimate that T M ∼ 50 − 70 MeV and √ s M ∼ 2 − 3 GeV, which is much smaller than the energies available at SPS, indicating that the signals from the first-order phase transition and the CEP are well separated from each other. Hence, the scenario depicted in Fig. 1 allows us to unambiguously identify the signals from both the first-order phase transition and the CEP. In fact, according to Ref. [13], the enhanced density fluctuation due to spinodal instability happens at around 2-4 AGeV ( and is in nice agreement with our estimate. Of course, the definitive value of √ s M should depend on the specific equation of state (EOS) adopted in the calculations, which is still largely unknown.
With decreasing collision energy to √ s N N ∼ 2−6 GeV, the ∆n is expected to rise again and reach a second maximum as a result of the spinodal instability as shown in Fig. 1 GeV are large, and more precise measurements are extremely important to confirm the present conclusion. Some effects can potentially affect the above numerical results. These include the uncertainties in the value of λ associated with the entropy of the expanding fireball as well as the chemical freeze-out temperature T ch and volume V ch . Also, the finite size of the fireball and the use of non-relativistic approximation for the nucleon momentum distributions [26,27] can affect the extracted values of C np , ∆n and ∆n I . However, these effects on the above dimensionless quantities are expected to have a weak dependence on the collision energy, and the non-monotonic behaviors shown in Fig. 2 should remain qualitatively similar.
4. Conclusion.-We have proposed that there should exist two separated maxima in the collision energy dependence of the baryon density fluctuation in heavy-ion collisions, with the lower energy one from the spinodal instability due to the first-order phase transition and the higher energy one induced by the second-order phase transition at the CEP. This double-peak structure seems to be supported by the collision energy dependence of the relative neutron density fluctuation ∆n = (δn) 2 / n 2 at kinetic freeze-out extracted from analyzing the measured yields of p, d and 3 H in central heavy-ion collisions at AGS and SPS energies within the coalescence model. In particular, we have found the ∆n to display a clear peak at It is particularly interesting to verify the present conclusion by using microscopic transport model simulations as well as hydrodynamics calculations with the proper treatment of the EOS (with and without the CEP) and critical fluctuations. Future measurements of light nuclei production at BES/RHIC, FAIR, NICA and NA61/SHINE are extremely useful to confirm the present observations and to more precisely determine the structure of the QCD phase diagram.