Contents lists available at ScienceDirect

Physics Letters B

www.elsevier.com/locate/physletb

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

Kai-Jia Suna, Lie-Wen Chen3'*, Che Ming Kob, Zhangbu Xuc d

a School of Physics and Astronomy and Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, China b Cyclotron Institute and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA c Brookhaven National Laboratory, Upton, NY 11973, USA

d School of Physics & Key Laboratory of Particle Physics and Particle Irradiation (MOE), Shandong University, Jinan, Shandong 250100, China

CrossMark

A R T I C L E I N F 0

A B S T R A C T

Article history:

Received 15 March 2017

Received in revised form 18 September

Accepted 19 September 2017 Available online 22 September 2017 Editor: W. Haxton

Based on the coalescence model for light nuclei production, we show that the yield ratio 0p-d-t = N3H Np/n2 of p, d, and 3H in heavy-ion collisions is sensitive to the neutron relative density fluctuation An = {(Sn)2)/(n)2 at kinetic freeze-out. From recent experimental data in central Pb + Pb collisions at Vsnn = 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 An as a function of the collision energy with a peak at -,/sNN = 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 ~ 144 MeV and baryon chemical potential of ~ 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.

© 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3.

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-4]. Results from lattice quantum chro-modynamics (LQCD) calculations [5-10] and effective model studies [11-16] have indicated that the QGP to hadronic matter transition is likely a 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 thermody-

* Corresponding author.

E-mail addresses: sunkaijia@sjtu.edu.cn (K.-J. Sun), lwchen@sjtu.edu.cn (L.-W. Chen), ko@comp.tamu.edu (C.M. Ko), xzb@bnl.gov (Z. Xu).

namic observables [12]. The resulting event-by-event fluctuations of conserved quantities in relativistic heavy-ion collisions have thus been extensively studied both theoretically and experimentally. For example, the energy dependence of the fourth-order fluctuation (io2) 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 */sNN = 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 heavy-ion 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 (R20Ut — R2ide) extracted from two-pion interferometry measurements [23-25] in Au + Au (^sNN = 7.7-200 GeV) and Pb + Pb (VsNN = 2.76 TeV) collisions has recently been observed with a maximum value located around „JsNN = 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

https://doi.org/10.1016/j.physletb.2017.09.056

0370-2693/© 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3.

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 hydrody-namic approach [30-32] and the microscopic transport model [33] indeed show that the spinodal instability due to the first-order phase transition between QGP and hadronic matter can induce large baryon density fluctuations. In the case that such density fluctuations can survive final-state interactions during the hadronic evolution of heavy-ion 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 relativis-tic 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 (An = ((Sn)2)/(n)2) at kinetic freeze-out in relativistic heavy-ion collisions can be encoded in the yield ratio of light nuclei, namely, 0p-d-t = N3HNp/N^. 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-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 3H) yields in Pb + Pb collisions at SPS energies measured by the NA49 Collaboration [47], we have observed a possible peak of An in Pb + Pb collisions at „JsNN = 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 tCEP - 144 MeV and /fp - 385 MeV.

We start by briefly introducing the newly derived analytical coalescence formula COAL-SH [45] for cluster production in rel-ativistic heavy-ion collisions. In COAL-SH, the yield Nc (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 Teff (including the effect of transversal flow), volume V, and number Nj of the i-th constituent with mass mj reads [45]

Nc = grel gsize gcM3'2

(4n/M)3/2 ( x

Vx(1 + x2) V1 + x2

G (li, x).

In the above, M = j=i mj is the rest mass of the cluster, lj is the orbital angular momentum associated with the i-th relative coordinate, w is the oscillator frequency of the cluster's internal wave function and is inversely proportional to Mr2ms with rrms being the root-mean-square (RMS) radius of the cluster, and

G(l, x) = J2'k=o why- with x = (2Teff/w)1/2 is the sup-

pression factor due to the orbital angular momentum on the coalescence probability [48]. In addition, gc = (2S + 1)/(f[j=1(2sj + 1)) is the coalescence factor for constituents of spin si to form a cluster of spin S, grel is the relativistic correction to the effective volume in momentum space, and gsize 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 gsize = 1. We also set grel = 1 because the masses of nucleons and light nuclei are much larger than the value of Teff. 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 3H are then simply given by

Nd = gd-

N3h = g3H

(mn + mp)3/2 NpNn (4n/oM)3/2

3/2 3/2

V xd(1 + xd) '

(2mn + mp)3/2 NpNn Vn/rn^3 V2 x2H(1 + x2H)2'

(2) (3)

where Np (Nn) is the number of protons (neutrons) in the emission source, the coalescence factor is gd = 3/4 for d and g3H = 1/4 for 3H, and we denote Xd = (2Teff/«d)1/2 and X3H = (2Teff/«3H)1/2 with the oscillator frequency «d = 8.1 MeV for d and W3H = 13.4 MeV for 3H obtained from their respective RMS radii rrms,d = 1.96 fm and rrms,3H = 1.76 fm [49]. The effective temperature Teff at the kinetic freeze-out in relativistic heavy-ion collision is typically about 200 MeV and is thus much larger than the oscillator frequencies Wd and W3H. Neglecting neutron and proton mass difference (mp = mn = m0) and noting Xd,X3H ^ 1, we then have

_ 3 / 2n \

Nd=2VH moref )

3/2 Np Nn

= 33! ( J^]3 NpN2

' 4 \mo Teff/

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-54] if p, n, d and 3H are assumed to be in thermal and chemical equilibrium and the binding energies of d and 3H 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 densities in the emission source as

n(r) = 1 j n(r)dr + Sn(r) = (n) + Sn(r),

np(r) = 1 j np(r)d? + snp (r) = (np) + Snp(r),

(6) (7)

where (•) denotes the average value over space and 8n(r) (Snp(r)) with (Sn) = 0 ((Snp ) = 0) denotes the fluctuation of neutron (proton) density from its average value (n) ((np)). We can then approximately rewrite Eqs. (4) and (5) as

3 / 2n Y/2 f _ _ ^

Nd=2/2 (moneffj ¡ârn(J)np^

3 ( 2n j3/2

= 232 (mrf (Np(n)+V (SnSnp

33/2 / 2n

33/2 / 2n Vf. . 2

N3H = -(mojf) fdrn(r)2"

2" p (r)

Table 1

Yields (dN/dy at midrapidity) of p, d, 3He and 3H as well as the yield ratio 3H/3He measured in Pb + Pb collisions at SPS energies [47] together with the derived yield ratio Op-d-t. The units for E and -jsNN are AGeV and GeV, respectively.

Centrality

Op-d-t

20 30 40 80 158

0-7% 0-7% 0-7% 0-7% 0-12%

46.1 ± 2.1 42.1 ± 2.0 41.3 ± 1.1 30.1 ± 1.0 23.9 ± 1.0

2.094 ± 0.168 1.379 ± 0.111 1.065 ± 0.086 0.543 ± 0.044 0.279 ± 0.023

3.58(±0.43) ; 1.89 (±0.23) ; 1.28 (±0.15) ; 3.90(±0.50) 1.50 (±0.20) x 10—

10— 10— 10— 10—

1.22 ± 0.10 1.18 ± 0.11 1.16 ± 0.15 1.15 ± 0.19 1.05 ± 0.15

4.37 (±0.64) x 10 2.23 (±0.34) x 10

1.48 (±0.26) ;

4.49 (±0.94) ;

10— 10—

1.58 (±0.31) x 10—

0.459 ± 0.014 0.494 ± 0.020 0.541 ± 0.022 0.458 ± 0.038 0.484 ± 0.037

33/2 / 2n \3 r 2 2

= ~4~ Uf [(<n>2 + <(5n)2>)NP

+ 2V (n)(SnSnp > + V ((Sn)2 Snp >].

Assuming Snp(r) = c(r)Sn(r), where the function c(r) can be positive or negative, we can then express the correlation between Sn(r) and Snp (r) as

(SnSnp> = V j drSn(r)Snp(r)

= V f dr c(r)(Sn(r))2. The above equation can also be written as

(SnSnp> = a {-jPl((Sn)2), (n>

with a being the correlation coefficient and 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 a = 1. By neglecting the term ((Sn)2Snp> in Eq. (9), we can rewrite Eqs. (8) and (9) as

3 f 2n \3/2

Nd = mrf Np (n>(1 + ^

33/2 / 2n \3 2 N3H = [m;^) Np(n>2[1 + (1 + 2a)An], (13)

where An = -(Sn)0>/-n>0 is a dimensionless quantity that characterizes the relative density fluctuation of neutrons.

Besides depending on An, both d and 3H yields also depend on Teff, Np and (n>. The density fluctuation in the emission source can be probed from the following yield ratio:

Op-d-t =

N3h Np N2

1 + (1 + 2a)An (1 + a An)2 ,

with g = 4/9 x (3/4)3/° ^ 0.29. The Op-d-t is constructed in such a way that many effects, such as those due to Teff, Np, (n>, volume and isospin asymmetry of the emission source, cancel out. Experimentally, one can thus extract An in relativistic heavy-ion collisions by measuring the yield ratio Op-d-t. When a An is much smaller than unity, the correction from a in Eq. (14) is second-order, and Op-d-t can be approximated as

Op-d-t ~ g(1 + An).

In this case, Op-d-t has a very simple linear dependence on An. We would like to point out that one may also choose other light nuclei such as 3He and 4He 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 3He and 4He are given, respectively, by

33/2 / 2n \3 2 . N

N3He = — (m^) Nn(np>2 0 + Anp + 2a An),

W 2n \ 2

N4He = mrf Np (np >(n>°

1 + (1 + 4a)An + Anp +

((SnSnp )2> (n>2(np>2

which further depend on the proton average density (np >, its relative density fluctuation Anp = ((Snp)2>/(np>2 and higher-order fluctuations. In Eq. (17), terms like ((Sn)2Snp> and ((Snp)2Sn> 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 „JsNN = 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 Op-d-t and An. Table 1 summarizes the yields (dN/dy at midrapidity) of p, d, 3He and 3H as well as the yield ratio 3H/3He 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 3H, we have used the relation 3H = 3He x 3H/3He. The derived Op-d-t is also shown in Table 1 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 1 that the energy dependence of Op-d-t shows a possible non-monotonic behavior with its largest value at „JsNN = 8.8 GeV. However, it should be pointed out that the evidence for the nonmonotonic behavior may not be statistically significant due to the sufficiently large uncertainty. Indeed, the value of Op-d-t at *ysNN = 8.8 GeV deviates by only about 2.5o from a x2 fit of Op-d-t at VsNv = 6.3 GeV, 7.6 GeV, 12.3 GeV and 17.3 GeV by the constant 0.471 ±0.018.

Equation (14) shows that for a fixed value of Op-d-t, the extracted value for An depends on the value of a. We note that Eq. (14) has no solution when a is larger than ~ 0.23 at ./sNN = 8.8 GeV. This feature suggests that a perfect or strong correlation between neutron and proton density fluctuations at kinetic freeze-out (i.e., a = 1 or a > 0.23) cannot appear in collisions at VsNN = 8.8 GeV. Similar features are also seen at other four collision energies, although the maximum values of a 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 2 shows the extracted values of An for

Table 2

Collision energy dependence of neutron relative density fluctuation An for a =

-0.2, -0.1, 0, 0.1 and 0.2. The units for E and -/s^N are AGeV and GeV, respectively.

E -JsNN Centrality An (a = -0.2) An (a = -0.1) An (a = 0) An (a = 0.1) An (a = 0.2)

20 6.3 0-7% 0.485 ± 0.037 0.526 ± 0.039 0.583 ± 0.048 0.669 ± 0.064 0.816 ± 0.099

30 7.6 0-7% 0.566 ± 0.044 0.623 ± 0.053 0.704 ± 0.068 0.833 ± 0.096 1.093 ± 0.177

40 8.8 0-7% 0.667 ± 0.046 0.746 ± 0.057 0.864 ± 0.076 1.071 ± 0.118 1.620 ± 0.322

80 12.3 0-7% 0.482 ± 0.090 0.523 ± 0.106 0.579 ± 0.130 0.662 ± 0.171 0.807 ± 0.262

158 17.3 0-12% 0.542 ± 0.084 0.594 ± 0.101 0.668 ± 0.127 0.782 ± 0.175 1.002 ± 0.345

-□- a = 0.2 -A- a = 0.1 -■- a = 0 a = -0.1 a = -0.2

i-1-1-1—I—I—I—I—I—I—r

Neutron relative density fluctuation in central Pb+Pb collisions (y=0)

I ^ i - -

- ~ I - H

,0/ 0-7%

0-7% (Centralityj

0-7% 0-12% _i_i_i_i_i_i_i_i

(sNN)la(GeV)

Fig. 1. Collision energy dependence of the neutron relative density fluctuation An in central Pb + Pb collisions at SPS energies based on data from Ref. [47]. Results for a = -0.2, -0.1, 0, 0.1 and 0.2 are shown by various dotted lines.

a = -0.2, -0.1, 0, 0.1 and 0.2 at different collisions energies. For all these values of a, a similar non-monotonic behavior is seen in the dependence of An on the collision energy with a peak at „JsNN = 8.8 GeV. Also, the obtained value of An 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 An, we plot in Fig. 1 the extracted An as a function of „JsNN for a = -0.2, -0.1, 0, 0.1 and 0.2. The extracted An is seen to increase with increasing value of a, and the increase is faster for larger value of Op-d-t. It is interesting to see that the peak at „JsNN = 8.8 GeV seems to always exist for all values of a considered here. Estimating the statistical significance of the nonmonotonic structure of the collision energy dependence of An by the same method as in the analysis of Op-d-t, we find the deviation of the An value at „JsNN = 8.8 GeV from the value obtained by a x2 fit of An at the other four energies is about 2.3a, 2.5a, 2.4a, 2.4a and 2.1a for a = -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 nonmonotonic behavior of the collision energy dependence of An in future measurements with higher precision.

The possible non-monotonic behavior of the collision energy dependence of An can be understood as follows. For central Pb + Pb collisions at higher incident energies (e.g., *,JsNN = 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 „JsNN = 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 „JsNN = 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 nonmonotonic 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 „JsNN = 8.8 GeV.

In the above, we have assumed that there is no energy dependence of a in collisions at SPS energies. In general, the correlation between the neutron and proton density fluctuations, characterized by the value of a, near the critical region (e.g., around „JsNN = 8.8 GeV) is likely larger than those at other collision energies, and thus the extracted An 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 „JsNN = 8.8 GeV are estimated to be T — 144 MeV and /B — 385 MeV. It is interesting to note that the estimated /iB — 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 „JsNN = 8.8 GeV corresponds to that at which a peak is seen in the measured K+/n+ 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 non-monotonic behavior in the relative neutron density fluctuation extracted from the measured yield ratio Op-d-t will still be present if the more accurate formula (with grel and gsjze [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 grel and gsjze 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 An, it does not change the non-monotonic behavior of An 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 Op-d-t = N3HNp/N% 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 (An = ((Sn)2)/(n)2) at kinetic freeze-out can be determined directly from the yield ratio Op-d-t = N3HNp/Njj. From measured yields of light nuclei at SPS energies by the NA49 Collaboration, we have extracted the collision energy dependence of An and found a possible non-monotonic behavior with a peak

at = 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 TCEP ~ 144 MeV and ¡if9 ~ 385 MeV, respectively. Given that the present statistical evidence for the peak structure of the collision energy dependence of An 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 authors thank Vadim Kolesnikov and Peter Seyboth for providing the experimental data. This work was supported in part by the Major State Basic Research Development Program (973 Program) in China under Contract Nos. 2015CB856904 and 2013CB834405, the National Natural Science Foundation of China under Grant Nos. 11625521, 11275125 and 11135011, the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education, China, the Science and Technology Commission of Shanghai Municipality (11DZ2260700), the U.S. Department of Energy under Contract No. DE-SC0015266 and No. DE-SC0012704, as well as the Welch Foundation under Grant No. A-1358 and Shandong University.

References

[1] E. Shuryak, Prog. Part. Nucl. Phys. 62 (2009) 48.

[2] K. Fukushima, T. Hatsuda, Rep. Prog. Phys. 74 (2011) 014001.

[3] P. Braun-Munzinger, V. Koch, T. Schäfer, J. Stachel, Phys. Rep. 621 (2016) 76.

[4] R. Pasechnik, M. Sumbera, Universe 3 (2017) 7.

[5] Y. Aoki, G. Endrodi, Z. Fodor, S.D. Katz, K.K. Szabo, Nature 443 (2006) 675.

[6] T. Bhattacharya, et al., Phys. Rev. Lett. 113 (2014) 082001.

[7] P. de Forcrand, PoS LAT2009 (2009) 010.

[8] Z. Fodor, S.D. Katz, J. High Energy Phys. 0404 (2004) 050.

[9] A. Li, A. Alexandru, K.F. Liu, Phys. Rev. D 84 (2011) 071503.

[10] P. de Forcrand, J. Langelage, O. Philipsen, W. Unger, Phys. Rev. Lett. 113 (2014) 152002.

[11] M. Asakawa, K. Yazaki, Nucl. Phys. A 504 (1989) 668.

[12] M.A. Stephanov, K. Rajagopal, E.V. Shuryak, Phys. Rev. Lett. 81 (1998) 4816.

[13] J. Berges, K. Rajagopal, Nucl. Phys. B 538 (1999) 215.

[14] Y. Hatta, T. Ikeda, Phys. Rev. D 67 (2003) 014028.

[15] M.A. Stephanov, Prog. Theor. Phys. Suppl. 153 (2004) 139.

[16] M. Asakawa, C. Nonaka, Nucl. Phys. A 774 ( 2006 ) 753.

[17] L.P. Csernai, J. Kapusta, L.D. McLerran, Phys. Rev. Lett. 97 (2006) 152303.

[18] R.A. Lacey, et al., Phys. Rev. Lett. 98 (2007) 092301.

[19] C.M. Hung, E.V. Shuryak, Phys. Rev. Lett. 75 (1995) 4003.

[20] D.H. Rischke, M. Gyulassy, Nucl. Phys. A 608 (1996) 479.

[21] P. Batyuk, et al., Phys. Rev. C 94 (2016) 044917.

[22] N.-U. Bastian, et al., Eur. Phys. J. A 52 (2016) 244.

[23] S. Pratt, Phys. Rev. Lett. 53 (1984) 1219.

[24] S. Chapman, P. Scotto, U.W. Heinz, Phys. Rev. Lett. 74 (1995) 4400.

[25] U.A. Wiedemann, P. Scotto, U.W. Heinz, Phys. Rev. C 53 (1996) 918.

[26] R.A. Lacey, Phys. Rev. Lett. 114 (2015) 142301.

[27] X.F. Luo, et al., STAR Collaboration, PoS CP0D2014 (2015) 019.

[28] T. Andrews, Bakerian lecture: on the continuity of the gaseous and liquid states of matter, Proc. R. Soc. Lond. 18 (1869) 42.

[29] B. Berche, M. Henkel, R. Kenna, Rev. Bras. Ensino Fis. 31 (2009) 2602.

[30] J. Steinheimer, J. Randrup, Phys. Rev. Lett. 109 (2012) 212301.

[31] J. Steinheimer, J. Randrup, V. Koch, Phys. Rev. C 89 (2014) 034901.

[32] C. Herold, M. Nahrgang, I. Mishustin, M. Bleicher, Nucl. Phys. A 925 (2014) 14.

[33] F. Li, C.M. Ko, Phys. Rev. C 93 (2016) 035205; F. Li, C.M. Ko, Nucl. Sci. Tech. 27 (2016) 140;

F. Li, C.M. Ko, Int. J. Mod. Phys. E 13 (2017) 1740012; F. Li, C.M. Ko, Phys. Rev. C 95 (2017) 055203.

[34] S.T. Butler, C.A. Pearson, Phys. Rev. Lett. 7 (1961) 69.

[35] H. Sato, K. Yazaki, Phys. Lett. B 98 (1981) 153.

[36] L.P. Csernai, J.I. Kapusta, Phys. Rep. 131 (1986) 223.

[37] C.B. Dover, U.W. Heinz, E. Schnedermann, J. Zimanyi, Phys. Rev. C 44 (1991) 1636.

[38] L.W. Chen, C.M. Ko, B.A. Li, Phys. Rev. C 68 (2003) 017601; L.W. Chen, C.M. Ko, B.A. Li, Nucl. Phys. A 729 (2003) 809.

[39] V. Greco, C.M. Ko, P. Levai, Phys. Rev. Lett. 90 (2003) 202302; V. Greco, C.M. Ko, P. Levai, Phys. Rev. C 68 (2003) 034904.

[40] R.J. Fries, et al., Phys. Rev. Lett. 90 (2003) 202303; R.J. Fries, et al., Phys. Rev. C 68 (2003) 044902.

[41] L.W. Chen, C.M. Ko, Phys. Rev. C 73 (2006) 044903.

[42] R. Fries, V. Greco, P. Sorensen, Annu. Rev. Nucl. Part. Sci. 58 (2008) 177.

[43] L. Zhu, C.M. Ko, X. Yin, Phys. Rev. C 92 (2015) 064911.

[44] K.J. Sun, L.W. Chen, Phys. Lett. B 751 (2015) 272; K.J. Sun, L.W. Chen, Phys. Rev. C 93 (2016) 064909; K.J. Sun, L.W. Chen, Phys. Rev. C 94 (2016) 064908.

[45] K.J. Sun, L.W. Chen, Phys. Rev. C 95 (2017) 044905.

[46] X. Yin, C.M. Ko, Y. Sun, L. Zhu, Phys. Rev. C 95 (2017) 054913.

[47] T. Anticic, et al., NA49 Collaboration, Phys. Rev. C 94 (2016) 044906.

[48] S. Cho, et al., ExHIC Collaboration, Phys. Rev. Lett. 106 (2011) 212001; S. Cho, et al., ExHIC Collaboration, Phys. Rev. C 84 (2011) 064910.

[49] G. Röpke, Phys. Rev. C 79 (2009) 014002.

[50] J.D. Bjorken, Phys. Rev. D 27 (1983) 140.

[51] P. Braun-Munzinger, J. Stachel, J. Phys. G 21 (1995) L17; P. Braun-Munzinger, J. Stachel, Nature 448 (2007) 302.

[52] A. Andronic, P. Braun-Munzinger, J. Stachel, H. Stöcker, Phys. Lett. B 697 (2011) 203.

[53] J. Cleymans, et al., Phys. Rev. C 84 (2011) 054916.

[54] J. Steinheimer, et al., Phys. Lett. B 714 (2012) 85.

[55] T.A. Armstrong, et al., E864 Collaboration, Phys. Rev. Lett. 83 (1999) 5431.

[56] The STAR Collaboration, Nature 473 (2011) 353.

[57] J. Cleymans, H. Oeschler, K. Redlich, S. Wheaton, Phys. Rev. C 73 (2006) 034905.

[58] X.Y. Xin, S.X. Qin, Y.X. Liu, Phys. Rev. D 90 (2014) 076006.

[59] N.G. Antoniou, A.S. Kapoyannis, Phys. Lett. B 563 (2003) 165.

[60] C. Alt, et al., NA49 Collaboration, Phys. Rev. C 77 (2008) 024903.

[61] M. Gazdzicki, M.I. Gorenstein, Acta Phys. Pol. B 30 (1999) 2705.

[62] W. Cassing, A. Palmese, P. Moreau, E. Bratkovskaya, Phys. Rev. C 93 (2016) 014902;

W. Cassing, A. Palmese, P. Moreau, E. Bratkovskaya, Phys. Rev. C 94 (2016) 044912.