Effects of QCD critical point on light nuclei production

Using the nucleon coalescence model, which can naturally take into account the correlations in the nucleon density distribution, we study the effects of QCD critical point on light nuclei production in relativistic heavy-ion collisions. We find that the yield ratio $N_tN_p/N_d^2$ of proton ($p$), deuteron ($d$) and triton ($t$) increases monotonically with the nucleon density correlation length, which is expected to increase significantly near the critical point in the QCD phase diagram. Our study thus demonstrates that the yield ratio $N_tN_p/N_d^2$ can be used as a sensitive probe of the QCD critical phenomenon. We further discuss the relation between the QCD phase transitions in heavy-ion collisions and the possible non-monotonic behavior of \pdt~in its collision energy dependence.

Introduction. According to the quantum chromodynamics (QCD), for a hadronic matter at sufficiently high temperatures and/or densities, the quarks and gluons inside the hadrons can be liberated to form a new phase of matter called the quark-gluon plasma (QGP) [1]. From calculations based on the lattice quantum chromodynamics (LQCD), the transition between the hadronic matter and the QGP is a smooth crossover if the matter has a low baryon chemical potential (µ B ). This smooth phase transition has been confirmed in experiments on ultrarelativistic heavy-ion collisions [2,3] in which the initially produced matter has a small µ B and high temperature (T ). Due to the fermion sign problem in LQCD [4], it remains, however, an open question whether this smooth crossover changes to a first-order phase transition at high µ B and low T [5], with a possible critical endpoint (CEP), corresponding to a second-order phase transition, on the first-order phase transition line in the µ B − T plane of the QCD phase diagram. Besides its intrinsic interest, knowledge on the properties of QCD phase diagram at finite µ B is also useful for understanding the structure of the inner core of a neutron star [6] and the gravitational wave from neutron-star mergers [7].
Locating the possible CEP and the phase boundary in the QCD phase diagram is one of the main goals for the heavy-ion collision experiments being carried out at the Relativistic Heavy Ion Collider (RHIC), the Facility for Antiproton and Ion Research (FAIR), the Nuclotronbased Ion Collider Facility (NICA), the High-Intensity Heavy Ion Accelerator Facility (HIAF), and the Japan Proton Accelerator Research Complex (J-PARC). For a recent review, see e.g. Ref. [8]. By changing the beam energy ( √ s) in heavy-ion collisions, different regions in the µ B − T plane of the QCD phase diagram can be explored. Although at very high √ s, the evolution trajectory of produced matter in the QCD phase diagram only passes across the crossover line, it could move close to the CEP or pass across the first-order phase transition line as the √ s becomes lower.
Since the hadronic matter, which is formed from the phase transition of the short-lived QGP created in heavyion collisions, undergoes a relatively long expansion, it is a great challenge to find observables that are sensitive to the phase transitions occurred during the earlier stage of the collisions. Based on the generic feature of divergent correlation length (ξ) at the critical point of phase transitions and the resulting properties of self similarities and universality classes [9], observables sensitive to the correlation length, such as the correlations and fluctuations of conserved charges [10][11][12][13], have been proposed. In particular, the fourth-order or kurtosis of net-proton multiplicity fluctuations has been suggested as an observable for the critical point because of its dependence on higher orders in ξ and its non-monotonic behavior as a function of √ s [11,13]. However, the event-by-event fluctuation in the number of net-protons in these theoretical studies [11,13] is for protons and antiprotons in certain spatial volume, which is very different from those within certain momentum cut that are measured in experiments [14,15]. It is unclear how the measured fluctuation in momentum space is related to the fluctuation in coordinate space due to the CEP, especially for heavy-ion collisions at low beam energies because of the lack of boost invariance [16] and the effects of thermal smearing [17].
On the other hand, light nuclei, such as the deuteron (d), triton ( 3 H or t), helium-3 ( 3 He), etc. that are arXiv:2008.02325v1 [nucl-th] 5 Aug 2020 produced in relativistic heavy-ion collisions, provide a promising tool to probe the spatial density fluctuation and correlation in the produced matter as they are formed from nucleons that are located in a very restricted volume of ∆x ∼ 2 fm and ∆p ∼ 100 MeV in phase space [18][19][20][21]. It has been shown that the fluctuation of nucleon density distributions can lead to an enhancement in the yield ratio N t N p /N 2 d through the where ∆ρ n is the average neutron density fluctuation in the coordinate space [20,21]. Also, it has been proposed that this ratio could be enhanced because of the modification of the nucleon-nucleon potential near the CEP [22,23]. Possible non-monotonic behavior of this ratio has indeed been seen in recent experimental data from heavy-ion collisions at both SPS energies [20] and RHIC BES energies [24]. However, no microscopic models have so far been able to explain the data.
In this Letter, we point out that the yield ratio N t N p /N 2 d in relativistic heavy-ion collisions is not only enhanced by a first-order phase transition in the produced matter, it also receives an additional enhancement from the long-range spatial correlation if the produced matter is near the CEP. To quantify the dependence of this ratio on the correlation length ξ in the produced matter in heavy-ion collisions if its evolution trajectory passes close to the CEP of the QCD phase diagram, we derive an expression to relate these two quantities. Because of the intrinsic resolution scale of around 2 fm given by the sizes of deuterons and tritons is comparable to the expected correlation length generated near the CEP [25], the production of these nuclei is thus sensitive to the QCD critical fluctuations. Our finding also suggests that the information on the CEP can be obtained from experiments on heavy-ion collisions by studying the collision energy dependence of the yield ratio N t N p /N 2 d . Effects of critical point on light nuclei production in heavy-ion collisions. Although rarely produced in high energy heavy-ion collisions, light nuclei, such as d, t, 3 He, helium-4 ( 4 He), hypertriton ( 3 Λ H) and their antiparticles, have been observed in experiments at RHIC [26] and the LHC [27]. As shown in Refs. [20,21,28], the effects of nucleon density correlations and fluctuations on light nuclei production can be naturally studied by using the coalescence model [29][30][31][32][33][34]. In this model, the number of deuterons produced from a hadronic matter can be calculated from the overlap of the proton and neutron joint distribution function f np ( 2 are, respectively, the relative coordinates and momentum between its proton and neutron [28], i.e., In the above, g d = 3/4 is the statistical factor for spin 1/2 proton and neutron to form a spin 1 deuteron. As usually used in the coalescence model, we approximate the deuteron Wigner function W d by Gaussian functions in both x and p, i.e., [20,21,28]. For the width parameter σ d in the deuteron Wigner function, it is related to the root-meansquare radius r d of deuteron by σ d = 4/3 r d ≈ 2.26 fm [33,35], which is much smaller than the size of the hot dense hadronic matter created in relativistic heavyion collisions. For the proton and neutron joint distribution function f np (x 1 , p 1 ; x 2 , p 2 ) in phase space, we take it to have the form by assuming that protons and neutrons are emitted from a thermalized source of temperature T and neutron and proton densities ρ n (x) and ρ p (x), respectively. In the above equation, m is the nucleon mass and ρ np (x 1 , x 2 ) is the joint density distribution function of protons and neutrons in the coordinate space, which can be written as in terms of the neutron and proton density correlation function C 2 (x 1 , x 2 ). Substituting f np from Eq. (2) and Eq. (3) in Eq. (1) and integrating over the nucleon momenta, we obtain the deuteron number as In the above, N the proton number in the emission source, denotes the deuteron number in the usual coalescence model studies without density fluctuations and correlations in the emission source, and C np = δρ n (x)δρ p (x) /( ρ n ρ p ) is the correlation between the neutron and proton density fluctuations, with · · · denoting the average over the coordinate space [20,21]. In obtaining Eq. (4), we have used the fact that the width parameter σ d in the deuteron Wigner function is much larger than the thermal wavelength of the nucleons in the emission source. We note that the C np term in Eq. (4) for deuteron production has been carefully studied in Refs. [20,21,28].
The nucleon density correlation becomes important near the CEP, where it is dominated by its singular part given by [9,36] In the above, ξ is the correlation length, η is the critical exponent of anomalous dimension, and λ is a parameter that varies smoothly with the temperature and the baryon chemical potential of the emission source. The nucleon correlation length is similar to the correlation length of the (net-)baryon density because nucleons carry most of the baryon charges in heavy-ion collisions. In the non-linear sigma model [10,37], the correlation length is given by ξ = 1/m σ , where m σ denotes the in-medium mass of the sigma meson and decreases as the system approaches the CEP. Since the value of the anomalous exponent η ≈ 0.04 is small [38], it is neglected in the present study. Using the fact that the nucleon number fluctuation δN 2 ∝ dxC 2 (x) ∝ λξ 2 near the critical point is always positive and enhanced, the λ parameter is positive as well. With Eq. (5), the deuteron number in Eq. (4) then becomes where the function G denotes the contribution from the long-range correlation between neutrons and protons, and it is given by with erfc(z) being the complementary error function. The behavior of G(z) is depicted in Fig. 1. As the correlation length ξ increases, the function G(z) is seen to increase monotonically and saturate to the value 2/π for z 1 or ξ σ d . This means that the divergence of ξ does not lead to a divergence of the deuteron yield, which is different from observables like the kurtosis of net-proton number distribution function. For small ξ, the function G increases as ξ 2 . At ξ ∼ σ d , it increases linearly with ξ, and the increase becomes much slower for ξ > 3σ d .
Similarly, the number of tritons from the coalescence of two neutrons and one proton is given by if the triton Wigner function is also approximated by Gaussian functions in the relative coordinates [28,31,32].
In the above, σ t is related to the root-mean-square radius r t of triton by σ t = r t = 1.59 fm [28,31,32,35]. The three-nucleon joint density distribution function ρ nnp can be expressed as in terms of three two-nucleon correlation functions C 2 (x i , x j ) with i = j, if one neglects the isospin dependence of two-nucleon correlation functions [10], and the three-nucleon correlation function C 3 (x 1 , x 2 , x 3 ). The contribution from the two-nucleon correlation functions to triton production can be similarly evaluated as in Eq. (6) for deuteron production. Keeping only the leading-order term in the function G, the triton number is then given by where N (0) t = 3 3/2 4 2π mT 3 N p ρ n 2 is the triton number in the absence of nucleon density fluctuations and correlations in the emission source, and ∆ρ n = δρ n (x) 2 / ρ n 2 ≥ 0 is the relative neutron density fluctuation [20,21,28]. The term O(G 2 ) in Eq. (10) comes from the contribution of the three-nucleon correlation. We note that the ∆ρ n is closely related to the second-order scaled density moment y 2 by [28], which has been frequently used to describe the density fluctuation or inhomogeneity in coordinate space [18,28,39].
By considering the yield ratio N t N p /N 2 d , which can be considered as the double ratio of N t /N d and N d /N p , one can eliminate the pre-factors of temperature and the proton number in Eqs. (6) and (10), which depend on the beam energy and the collision system. Neglecting the difference between σ d and σ t for simplicity by denoting σ ≈ σ d ≈ σ t and keeping only the leading-order term in G, the yield ratio N t N p /N 2 d can be simplified to Eq. (11) shows that besides its enhancement by the neutron density fluctuation ∆ρ n [20,21], the yield ratio N t N p /N 2 d is also enhanced by the nucleon density correlations characterized by the correlation length ξ. Although the density fluctuations ∆ρ n in a homogeneous system vanishes, the correlation length ξ in the system becomes divergent if the system is close to the CEP. On the other hand, in the presence of a first-order phase transition with the coexistence of two phases, the system could have large density inhomogeneity [18,28,39,40] and thus non-vanishing density fluctuations, even though the correlation length in each phase is significantly smaller than in the case that the system is near the critical point. As a result, the existence of a first-order phase transition and the CEP in the system can both lead to enhancements of the yield ratio N t N p /N 2 d if their effects can survive the hadronic evolution in a heavy-ion collision.
The above result can be generalized to the yield ratios involving the heavier 4 He (α) [20,22,23], e.g. the ratios Compared to the yield ratio N t N p /N 2 d in Eq. (11), these two yield ratios show a larger sensitivity to the correlation length ξ. However, the yield of α particle in heavyion collisions is much more difficult to measure precisely at high collision energies because of its small value due to the large penalty factor e −A(m−µ B )/T for the yield of a nucleus with A nucleons [41].
As can be seen from Eq. (11), the yield ratio N t N p /N 2 d encodes directly the spatial density fluctuations and correlations of nucleons due to, respectively, the first-order QGP to hadronic matter phase transition and the critical point in the produced matter from relativistic heavy-ion collisions. This is in contrast to the higher-order net-proton multiplicity fluctuations [11,13] based on the measurement of the event-by-event fluctuation of the net-proton number distribution in momentum space with its relation to the spatial correlations due to the CEP still lacking [16,17]. Also, the N t N p /N 2 d ratio has a natural resolution scale of around 2 fm, which is comparable to the correlation length ξ 2-3 fm [25] that could be developed in heavy-ion collisions. Such an intrinsic scale is absent in the event-by-event fluctuation of the net-proton number distribution as it only depends on the selected rapidity range in the experiments [8,42].
Collision energy dependence of the yield ratio N t N p /N 2 d . Since both the density fluctuations and longrange correlations of nucleons in the emission source can lead to an enhanced yield ratio N t N p /N 2 d as shown in the above, the effects of QCD phase transitions can thus be studied in experiments from the collision energy dependence of this ratio. It has been demonstrated in Ref. [43] that the correlation length ξ along the chemical freezeout line peaks near the CEP. This indicates that the collision energy dependence of ξ is likely to exhibit a peak structure with its maximum value at certain collision energy √ s H . Due to the critical slowing down [25] in the growth of the correlation length, the value of ξ is, however, limited to be around 2-3 fm at the time the matter produced in realistic heavy-ion collisions is near the CEP. The neutron density fluctuation ∆ρ n is mostly related to the first-order phase transition during which large density inhomogeneity could be developed due to the spinodal instability [18,28,39,40]. It was estimated in Ref. [44] that the largest effect of a first-order phase transition could be developed at around half the critical temperature [44] when the phase trajectory spends the longest time in the spinodal unstable region of the QCD phase diagram. Therefore, the collision energy dependence of ∆ρ n is also expected to have a peak structure with its maximum value at a lower collision energy √ s L .
As a result, the yield ratio N t N p /N 2 d as a function of the collision energy √ s shows two possible nonmonotonic behaviors. The first one has a double-peak structure, with one peak at √ s H due to the critical point and the other peak at a lower collision energy √ s L due to the first-order phase transition. This double-peak structure was conjectured in Ref. [21] without the explicit relation between N t N p /N 2 d and ξ given in Eq. (11). Due to the flattening of the function G given by Eq. (7) for large ξ, the signal from the critical point is broadened. It is thus also possible that √ s L and √ s H are so close that the signals from the CEP and the first-order phase transition overlap, resulting in only one broad peak in the collision energy dependence of the yield ratio N t N p /N 2 d . Possible non-monotonic behavior of the yield ratio N t N p /N 2 d has been seen in recent experimental data from heavy-ion collisions at both SPS energies [20,45] and RHIC BES energies [24]. In particular, the data for N t N p /N 2 d at collision energies from √ s NN = 6.3 GeV to √ s NN = 200 GeV shows a possible double-peak structure. However, due to the large error bars in the data, one cannot exclude the possibility that the double peaks are actually a single broad peak. To extract the collision energy dependence of ∆ρ n and ξ from the experimental data by using Eq. (11), one needs to know the value of the λ parameter in the critical region and its evolution during the hadronic evolution of heavy ion collisions, which requires studies that are beyond the scope of present paper.
In the statistical hadronization model, which assumes that both yields of deuterons and trions remain constant during the hadronic evolution of heavy ion collisions, this ratio would increase with increasing collision energy after including the strong decay contribution to protons [46]. Calculations based on transport models [47,48] without the CEP in the QCD phase diagram all give an essentially energy-independent constant value for N t N p /N 2 d and thus fail to describe the data. However, a recent multiphase transport model study [28], which includes a firstorder QCD phase transition, shows that the density fluctuation or inhomogeneity induced during the first-order phase transition can largely survive the hadronic evolution, because of the fast expansion of the produced matter, and eventually leads to an enhanced yield ratio N t N p /N 2 d at the kinetic freeze out of nucleons when they undergo their last scatterings. One thus expects the long-range correlation to similarly persist until kinetic freeze out and also lead to an enhancement of the yield ratio N t N p /N 2 d . Summary and outlook. Based on the nucleon coalescence model, we have obtained an explicit expression that relates the yield ratio N t N p /N 2 d to the nucleon density correlation length ξ in the hadronic matter produced in heavy-ion collisions, which could be appreciable if the produced matter is initially close to the CEP in the QCD phase diagram. This ratio is found to increase monotonically with the dimensionless quantity ξ/σ where σ ≈ 2 fm denotes the sizes of deuteron and triton. This enhancement is in addition to that due to the large neutron density fluctuation ∆ρ n that could be developed during a first-order QGP to hadronic matter phase transition previously studied in Refs. [20,21,28]. Consequently, the collision energy dependence of this ratio is expected to have a double-peak or a broad one-peak structure depending on the closeness in √ s between the signal of the CEP and that of the first-order phase transition. Such a non-monotonic behavior in the collision energy dependence of the yield ratio N t N p /N 2 d has indeed been seen in the preliminary data from the STAR Collaboration. Our study has thus led to the possibility of extracting the information of the CEP and the phase boundary of QCD phase diagram from comparing the precisely measured data on the yields of light nuclei in heavy-ion collisions with those from theoretical models based on the transport approach [28] and the various hydrodynamic approaches [49][50][51][52][53].
This work was supported in part by the US Department of Energy under Contract No.de-sc0015266 and the Welch Foundation under Grant No. A-1358.