Probing QCD critical fluctuations from light nuclei production in relativistic heavy-ion collisions

Based on the coalescence model for light nuclei production, we show that the yield ratio $\mathcal{O}_\text{p-d-t} = N_{^3\text{H}} N_p / N_\text{d}^2$ of $p$, d, and $^3$H in heavy-ion collisions is sensitive to the neutron relative density fluctuation $\Delta n= \langle (\delta n)^2\rangle/\langle n\rangle^2$ at kinetic freeze-out. From recent experimental data in central Pb+Pb collisions at $\sqrt{s_{NN}}=6.3$~GeV, $7.6$~GeV, $8.8$~GeV, $12.3$~GeV and $17.3$~GeV measured by the NA49 Collaboration at the CERN Super Proton Synchrotron (SPS), we find a possible non-monotonic behavior of $\Delta n$ as a function of the collision energy with a peak at $\sqrt{s_{NN}}=8.8$~GeV, indicating that the density fluctuations become the largest in collisions at this energy. With the known chemical freeze-out conditions determined from the statistical model fit to experimental data, we obtain a chemical freeze-out temperature of $\sim 144~$MeV and baryon chemical potential of $\sim 385~$MeV at this collision energy, which are close to the critical endpoint in the QCD phase diagram predicted by various theoretical studies. Our results thus suggest the potential usefulness of the yield ratio of light nuclei in relativistic heavy-ion collisions as a direct probe of the large density fluctuations associated with the QCD critical phenomena.

Understanding the properties of strongly interacting matter under extreme conditions, particularly the phase transition between the quark-gluon plasma (QGP) and the hadronic matter, is a topic of great current interest [1][2][3][4]. Results from lattice quantum chromodynamics (LQCD) calculations [5][6][7][8][9][10] and effective model studies [11][12][13][14][15][16] have indicated that the QGP to hadronic matter transition is likely first-order phase transition if the system has a large baryon chemical potential but changes to a crossover if its baryon chemical potential is small. This suggests the existence of a critical endpoint (CEP), where the first-order phase transition ends, in the temperature versus baryon chemical potential (T, µ B ) plane of QCD phase diagram. To search for the CEP and locate its position in the QCD phase diagram, experiments have been carried out through the Beam Energy Scan (BES) program and will also be performed at the future Facility for Antiproton and Ion Research (FAIR) in Germany.
It has been argued that the enhanced long-wavelength fluctuations near the CEP can lead to singularities in all thermodynamic observables [12]. The resulting eventby-event fluctuations of conserved quantities in relativistic heavy-ion collisions have thus been extensively studied both theoretically and experimentally. For exam-ple, the energy dependence of the fourth-order fluctuation (κσ 2 ) of net-proton distribution measured in the BES program by the STAR Collaboration is found to exhibit the largest deviation from unity in Au+Au collisions at √ s N N = 7.7 GeV [27]. Also, owing to the different features between a first-order phase transition and a rapid crossover, one expects a non-monotonic behavior in the collision energy and centrality dependence of certain properties of the produced matter in heavyion collisions as it approaches the CEP, such as the ratio of its shear viscosity to entropy density [17,18], expansion speed [19,20] and the slope of direct flow of light cluster [21,22]. Furthermore, a non-monotonic excitation function for the Gaussian emission source radii difference (R 2 out − R 2 side ) extracted from two-pion interferometry measurements [23][24][25] in Au+Au ( √ s N N =7.7-200 GeV) and Pb+Pb ( √ s N N = 2.76 TeV) collisions has recently been observed with a maximum value located around √ s N N = 40 GeV [26].
In analogy to the phenomenon of critical opalescence observed in the liquid-gas phase transition [28,29], the matter created in relativistic heavy-ion collisions could develop large baryon density fluctuations when its evolution trajectory in the (T, µ B ) plane of QCD phase diagram passes across the first-order phase transition line, especially when it is close to the CEP. When the evolution trajectory approaches the CEP, the correlation length increases drastically, and the density fluctuation enhances accordingly and reaches its maximum value at the CEP. Studies based on both the hydrodynamic ap-proach [30][31][32] and the microscopic transport model [33] indeed show that the spinodal instability due to the firstorder phase transition between the QGP and hadronic matter can induce large baryon density fluctuations. In the case that such density fluctuations can survive finalstate interactions during the hadronic evolution of heavyion collisions, there should exist strong fluctuations in the nucleon density and thus significant inhomogeneity in the spatial distribution of nucleons at kinetic freeze-out. The baryon density fluctuations is, however, expected to be negligible if the QGP to hadronic matter transition is a crossover. Therefore, the nucleon density fluctuations at kinetic freeze-out in relativistic heavy-ion collisions may provide a unique probe to the critical endpoint in the QCD phase diagram.
In this Letter, we show for the first time that the relative density fluctuation of neutrons (∆n = (δn) 2 / n 2 ) at kinetic freeze-out in relativistic heavy-ion collisions can be encoded in the yield ratio of light nuclei, namely, Our result thus has the advantage of directly measuring the density fluctuation instead of using the number fluctuation to infer the density fluctuation as having been done so far. Our study is based on the coalescence model for light nuclei production [34][35][36][37][38][39][40][41][42][43][44][45][46]. In this model, the probability for the production of a nucleus depends on the nucleon many-body correlations and is thus affected by the fluctuations in the nucleon number or density. From analyzing the very recent data on the proton (p), deuteron (d) and triton (t or 3 H) yields in Pb+Pb collisions at SPS energies measured by the NA49 Collaboration [47], we have observed a possible peak of ∆n in Pb+Pb collisions at √ s N N = 8.8 GeV. This result has further allowed us to estimate that the temperature and baryon chemical potential at which the CEP is located in the QCD phase diagram are T CEP ∼ 144 MeV and µ CEP B ∼ 385 MeV. We start by briefly introducing the newly derived analytical coalescence formula COAL-SH [45] for cluster production in relativistic heavy-ion collisions. In COAL-SH, the yield N c (per unit rapidity) of a cluster at midrapidity and consisting of A constituent particles from the hadronic matter at kinetic freeze-out or emission source of effective temperature T eff (including the effect of transversal flow), volume V , and number N i of the i-th constituent with mass m i reads [45] In the above, M = A i=1 m i is the rest mass of the cluster, l i is the orbital angular momentum associated with the i-th relative coordinate, ω is the oscillator frequency of the cluster's internal wave function and is in-versely proportional to M r 2 rms with r rms being the rootmean-square (RMS) radius of the cluster, and G(l, x) = l k=0 l! k!(l−k)! 1 (2k+1)x 2k with x = (2T eff /ω) 1/2 is the suppression factor due to the orbital angular momentum on the coalescence probability [48]. In addition, g c = (2S + 1)/( A i=1 (2s i + 1)) is the coalescence factor for constituents of spin s i to form a cluster of spin S, g rel is the relativistic correction to the effective volume in momentum space, and g size is the correction due to the finite size of produced cluster.
In central Pb+Pb collisions considered here, V is much larger than the sizes of light nuclei, and we thus set g size = 1. We also set g rel = 1 because the masses of nucleons and light nuclei are much larger than the value of T eff . For light nuclei included in the present study, all the constituent nucleons are in s-state (l = 0), and we thus have G(l, x) = 1. From Eq. (1), the yields of d and 3 H are then simply given by where N p (N n ) is the number of protons (neutrons) in the emission source, the coalescence factor is g d = 3/4 for d and g3 H = 1/4 for 3 H, and we denote  [49]. The effective temperature T eff at the kinetic freeze-out in relativistic heavy-ion collision is typically about 200 MeV and is thus much larger than the oscillator frequencies ω d and ω3 H . Neglecting neutron and proton mass difference (m p = m n = m 0 ) and noting x d , x3 H ≫ 1, we then have Although the coalescence formula COAL-SH is derived by assuming the Bjorken boost invariance [50] for the emission source, Eqs. (4) and (5) turn out to be also valid for an isotropically expanding fireball. This is not surprising as only the effective temperature, volume, and proton and nucleon numbers appear in these equations. Also, the above equations are consistent with the predictions from the thermal (statistical) model [51][52][53][54] if p, n, d and 3 H are assumed to be in thermal and chemical equilibrium and the binding energies of d and 3 H are neglected.
In obtaining Eqs. (4) and (5), we have assumed that nucleons are uniformly distributed in space at kinetic  freeze-out. To take into account density fluctuations of nucleons, we express the neutron and proton density in the emission source as where · denotes the average value over space and δn( r) (δn p ( r)) with δn = 0 ( δn p = 0) denotes the fluctuation of neutron (proton) density from its average value n ( n p ). We can then approximately rewrite Eqs. (4) and (5) as (8) and Assuming δn p ( r) = c( r)δn( r), where the function c( r) can be positive or negative, we can then express the correlation between δn( r) and δn p ( r) as The above equation can also be written as with α being the correlation coefficient and np n accounting for the isospin asymmetry of the emission source. In the case that the neutron and proton density fluctuations are completely correlated, we then have α = 1. By neglecting the term (δn) 2 δn p in Eq. (9), we can rewrite Eqs. (8) and (9) as where ∆n = (δn) 2 / n 2 is a dimensionless quantity that characterizes the relative density fluctuation of neutrons.
Besides depending on ∆n, both d and 3 H yields also depend on T eff , N p and n . The density fluctuation in the emission source can be probed from the following yield ratio: with g = 4/9×(3/4) 3/2 ≈ 0.29. The O p-d-t is constructed in such a way that many effects, such as those due to T eff , N p , n , volume and isospin asymmetry of the emission source, cancel out. Experimentally, one can thus extract ∆n in relativistic heavy-ion collisions by measuring the yield ratio O p-d-t . When α∆n is much smaller than unity, the correction from α in Eq. (14) is second-order, and O p-d-t can be approximated as In this case, O p-d-t has a very simple linear dependence on ∆n. We would like to point out that one may also choose other light nuclei such as 3 He and 4 He to extract the nucleon density fluctuation at kinetic freeze-out. In these cases, however, information on the isospin at freeze-out is needed and also the higher-order density fluctuations may be involved. For example, the yields of 3 He and 4 He are given, respectively, by which further depend on the proton average density n p , its relative density fluctuation ∆n p = (δn p ) 2 / n p 2 and higher-order fluctuations. In Eq. (17), terms like (δn) 2 δn p and (δn p ) 2 δn are neglected. Eqs. (12)- (17) show that large density fluctuations can affect the yields of light nuclei in relativistic heavy-ion collisions and lead to an A dependence different from n A that is expected from the statistical model [55]. Existing experimental data from the Alternating Gradient Synchrotron (AGS) at √ s N N = 4.8 GeV have shown a striking exponential behavior with a penalty factor of about 50 per additional nucleon to the produced nuclear cluster up to A = 7 [55]. Similarly, such a regular exponential behavior is seen at RHIC energies for A ≤ 4 [56]. These results have thus ruled out large nucleon density fluctuations at kinetic freeze-out in heavy-ion collisions at AGS and RHIC top energies.
However, recently published results on light nuclei production in central Pb+Pb collisions at SPS energies [47] show a quite different behavior. This can be seen from the collision energy dependence of O p-d-t and ∆n. Table I summarizes the yields (dN/dy at midrapidity) of p, d, 3 He and 3 H as well as the yield ratio 3 H/ 3 He measured in central Pb+Pb collisions at 20 AGeV (0 − 7% centrality), 30 AGeV (0 − 7% centrality), 40 AGeV (0 − 7% centrality), 80 AGeV (0 − 7% centrality), and 158 AGeV (0 − 12% centrality) by the NA49 Collaboration [47]. In obtaining the yield of 3 H, we have used the relation 3 H= 3 He× 3 H/ 3 He. The derived O p-d-t is also shown in Table I with errors estimated by assuming they are dominated by correlated systematic errors as a result of similar detector acceptance and phase-space extrapolation. It is seen from Table I  Equation (14) shows that for a fixed value of O p-d-t , the extracted value for ∆n depends on the value of α. We note that Eq. (14) has no solution when α is larger than ∼ 0.23 at √ s N N = 8.8 GeV. This feature suggests that a perfect or strong correlation between neutron and proton density fluctuations at kinetic freeze-out (i.e., α = 1 or α > 0.23) cannot appear in collisions at √ s N N = 8.8 GeV. Similar features are also seen at other four collision energies, although the maximum values of α are larger, i.e., 0.32 for 6.3 GeV, 0.28 for 7.6 GeV, 0.32 for 12.3 GeV and 0.29 for 17.3 GeV. Table II shows the extracted values of ∆n for α = −0.2, −0.1, 0, 0.1 and 0.2 at different collisions energies. For all these values of α, a similar non-monotonic behavior is seen in the dependence of ∆n on the collision energy with a peak at √ s N N = 8.8 GeV. Also, the obtained value of ∆n is much larger than that due to the event-by-event statistical fluctuation in the neutron multiplicity, which is expected to be inversely proportional to its mean value and is thus only about a few per cent. To see more clearly the collision energy dependence of ∆n, we plot in Fig. 1 the extracted ∆n as a function of √ s N N for α = −0.2, −0.1, 0, 0.1 and 0.2. The extracted ∆n is seen to increase with increasing value of α, and the increase is faster for larger value of O p-d-t . It is interesting to see that the peak at √ s N N = 8.8 GeV seems to always exist for all values of α considered here. Estimating the statistical significance of the non-monotonic structure of the collision energy dependence of ∆n by the same method as in the analysis of O p-d-t , we find the deviation of the ∆n value at √ s N N = 8.8 GeV from the average value at the other four energies is about 2.3σ, 2.5σ, 2.4σ, 2.4σ and 2.1σ for α = −0.2, −0.1, 0, 0.1 and 0.2, respectively. Given that the present statistical evidence is still weak, it is extremely important to confirm or rule out this possible non-monotonic behavior of the collision energy dependence of ∆n in future measurements with higher precision.
The possible non-monotonic behavior of the collision energy dependence of ∆n can be understood as follows. For central Pb+Pb collisions at higher incident energies (e.g., √ s N N = 17.3 GeV and 12.3 GeV), the reaction system may undergo a crossover rather than a first-order phase transition between the QGP and the hadronic matter, and the density fluctuation in the produced matter is thus insignificant. With decreasing incident energy (e.g., around √ s N N = 8.8 GeV), the reaction system may pass by or approach closely to the CEP and thus develop the largest density fluctuation. With further decrease in the incident energy (e.g., at √ s N N = 7.6 GeV and 6.3 GeV), the reaction system may move away from the CEP and barely cross the first-order transition line, and the density fluctuation decreases as a result of the smaller size and shorter lifetime of the QGP at lower energies. When the incident energy is further lowered, the reaction system may miss the first-order transition line and no QGP to hadronic matter transition occurs in the collisions, thus resulting in negligible density fluctuation at the kinetic freeze-out. Therefore, the possible non-monotonic behavior shown in Fig. 1 is consistent with the scenario that the CEP may be reached or closely approached in the produced QGP during its time evolution in central Pb+Pb collisions around √ s N N = 8.8 GeV. In the above, we have assumed that there is no energy dependence of α in collisions at SPS energies. In general, the correlation between the neutron and proton density fluctuations, characterized by the value of α, near the critical region (e.g., around √ s N N = 8.8 GeV) is likely larger than those at other collision energies, and thus the extracted ∆n from Eq. (14) could be larger and the peak structure would become more pronounced. From the parametrization in Ref. [57] for the chemical freeze-out conditions based on the statistical model fit to available experimental data, the temperature and baryon chemical potential at √ s N N = 8.8 GeV are estimated to be T ∼ 144 MeV and µ B ∼ 385 MeV. It is interesting to note that the estimated µ B ∼ 385 MeV for CEP is close to those predicted from the LQCD [8] and Dyson-Schwinger equation (DSE) [58] as well as that based on the hadronic bootstrap approach [59]. Also, the collision energy √ s N N = 8.8 GeV corresponds to that at which a peak is seen in the measured K + /π + ratio by the NA49 Collaboration [60], which has been interpreted as a signature for the onset of QGP formation [61] or the restoration of chiral symmetry [62] in these collisions.
Although the present study is based on the simple formulas given by Eq. (4) and Eq. (5), the possible nonmonotonic behavior in the relative neutron density fluctuation extracted from the measured yield ratio O p-d-t will still be present if the more accurate formula (with g rel and g size [45]) in Eq. (1) is used. This is because the variation in the value of g in Eq. (14) after taking into account the effects due to g rel and g size is less than 10% for the SPS energies considered here. In addition, although the correlation between neutron and proton density fluctuations influences the value of the extracted ∆n, it does not change the non-monotonic behavior of ∆n as a function of the collision energy. Our study is, however, based on one set of experimental data with large uncertainties and a simplified model. Further experimental and theoretical investigations are needed to verify the present results and eventually establish the yield ratio O p-d-t = N3 H N p /N 2 d as a robust probe to the QCD critical endpoint. These include the experimental BES program at RHIC in the energy range considered here with high luminosity beams as well as detectors of excellent particle identification and large acceptance, and theoretical modeling of light nuclei production and its connection to baryon density fluctuations.
In summary, with a newly derived analytical coalescence formula for cluster production in heavy-ion collisions, we have demonstrated that information on the relative density fluctuation of neutrons (∆n = (δn) 2 / n 2 ) at kinetic freeze-out can be determined directly from the yield ratio O p-d-t = N3 H N p /N 2 d . From measured yields of light nuclei at SPS energies by the NA49 Collaboration, we have extracted the collision energy dependence of ∆n and found a possible non-monotonic behavior with a peak at √ s N N = 8.8 GeV, suggesting that the CEP in the QCD phase diagram may have been reached or closely approached in these collisions with its temperature and baryon chemical potential estimated to be T CEP ∼ 144 MeV and µ CEP B ∼ 385 MeV, respectively. Given that the present statistical evidence for the peak structure of the collision energy dependence of ∆n is still weak, future measurements of light nuclei production in the BES program at RHIC are extremely useful to confirm the present observations and to more precisely determine the location of the CEP in the QCD phase diagram. The