The major goal of the Beam Energy Scan (BES) program at the Relativistic Heavy Ion Collider (RHIC) is to explore the phase diagram of the quantum chromodynamic (QCD) matter. By tuning the collision energies, the QCD phase diagram can be mapped and displayed into a two dimensional plane of temperature (T ) versus baryon chemical potential (μ B ) [1][2][3][4]. Lattice QCD calculations predicted a crossover transition from hadronic matter to a plasma of deconfined quarks and gluons (QGP) at vanishing μ B [5], while QCD-based model calculations suggested that the phase transition is of first-order at large μ B [6]. The critical end point (CEP) is a key feature of the QCD phase diagram, representing the point where the first-order phase transition boundary terminates [6][7][8]. Many efforts have been made to search for the possible CEP in heavy-ion collisions [1,2,[9][10][11], with several measurements from the BES program at RHIC exhibiting a non-monotonic variation with √ s NN , such as the net-proton kurtosis [3,12,13], the Hanbury- Brown-Twiss (HBT) radii [14,15] and the yield ratio of light nuclei production [16].
The aim of this work is to look for critical intermittency induced by the CEP [17,18] in heavy-ion collisions. Upon approaching a critical point, the correlation length of the system diverges and the system becomes scale invariant, or self-similar [19][20][21]. Based on the 3D-Ising universality class arguments [17,[22][23][24], the density-density correlation function for small momentum transfer has a power-law structure, leading to large density fluctuations in heavy-ion collisions [17,[22][23][24]. Such fluctuations can be probed in transverse momentum phase space within the framework of an intermittency analysis by utilizing the scaled factorial moments (SFMs) [17,[22][23][24][25][26]. To achieve this, the D-dimensional phase space is partitioned into M D equal-sized cells and the observable, qthorder SFM or F q (M), is defined as follows [17,18,[27][28][29][30]:
where M D is the number of cells in D-dimensional phase space and n i is the measured multiplicity of a given event in the ith cell. The angle bracket denotes an average over the events.
The intermittency appears as a power-law (scaling) behavior of SFMs [17,18,28,29,31]. If the system features density fluctuations, SFMs will obey a power-law behavior of
where φ q is called the intermittency index quantifying the strength of intermittency [17,22,23,25,27]. In this paper, another expected type of power-law behavior will be used: F q (M) ∝ F 2 (M) β q , M 1, where β q is the scaling index [18,[31][32][33][34][35]. According to the Ginzburg-Landau (GL) theory [18,31], the β q is independent of the details of the critical parameters, allowing for experimental measurement of F q (M) ∝ F 2 (M) β q behavior without the signal being washed out during hadronic evolution. To describe the general consequences of the phase transition, a scaling exponent (ν) is given by β q ∝ (q -1) ν [18,31,32,35,36]. Here, ν also quantifies the strength of intermittency. Near the CEP, the value of ν is predicted to be equal to 1.304 in the entire phase space based on the GL theory [18] and equal to 1.0 from the calculations of the 2D Ising model [31,37]. Over the last decade, the NA49 and the NA61/SHINE experiments have been performing intermittency analyses in heavy-ion systems of various sizes and collision energies [10,25,27,38,39]. The NA49 experiment observed strong intermittency with φ 2 = 0.96 ± 0.16 for protons in Si+Si collisions at 158 A GeV [27]. Two studies, using a Critical Monte Carlo (CMC) model [26] and a cascade ultra-relativistic quantum molecular dy-namics (UrQMD) model with hadronic potentials [40], respectively, suggested that large intermittency could be observed in Au+Au collisions at RHIC energies.
This letter reports the collision energy and centrality dependence of SFMs and intermittency exponents for identified charged hadrons in Au+Au collisions at RHIC/STAR. The data presented here were obtained from Au+Au collisions at √ s NN = 7.7, 11.5, 14.5, 19.6, 27, 39, 54.4, 62.4, and 200 GeV, recorded by the Solenoidal Tracker at RHIC (STAR) experiment from 2010 to 2017 [41]. These energies correspond to μ B values ranging from 20 to 420 MeV at chemical freeze-out [41]. The 7.7,11.5,39,62.4,and 200 GeV data were collected in 2010. The 19.6 GeV and 27 GeV data were collected in 2011, and the 14.5 GeV and 54.4 GeV data were collected in 2014 and 2017. All data were obtained using the Time Projection Chamber (TPC) and the Time-of-Flight (TOF) detectors at STAR [42,43]. Events are selected within a certain Z -position range (|V Z | < 30 cm) from the center of the TPC along the beam line (|V Z | < 50 cm for 7.7 GeV) to optimize for the uniformity in the response of the detectors [13]. Background events, which include interactions with the beam pipe, are rejected by requiring a vertex radius V r less than 2 cm from the center of STAR (V r < 1 cm for 14.5 GeV). To avoid self-correlation [13,44], the centrality is determined from the uncorrected charged particle multiplicity within a pseudorapidity window of 0.5 < |η| < 1, chosen to be outside the analysis window of |η| < 0.5. The centrality is represented by the average number of participating nucleons N part obtained by fitting the reference multiplicity distribution with a Monte Carlo Glauber model [44,45]. The number of events for √ s NN = 7.7, 11.5, 14.5, 19.6, 27, 39, 54.4, 62.4, and 200 GeV, are 3.3, 6.8 13.1, 16.2, 32.2, 89.3, 441.7, 46.7, and 236.0 million, respectively.
Charged hadrons, including protons (p), antiprotons ( p), kaons (K ± ), and pions (π ± ), are identified using the TPC and TOF detectors. TPC particle identification is performed using the measured energy loss (dE/dx), with K ± and π ± requiring a momentum range of 0.2 < p T < 0.4 GeV/c, and p and p requiring a momentum range of 0.4 < p T < 0.8 GeV/c. In addition, the mass squared from the TOF detector is used for particle identification, with K ± and π ± requiring a momentum range of 0.4 < p T < 1.6 GeV/c, and p and p requiring a momentum range of 0.8 < p T < 2.0 GeV/c. A maximum distance of closest approach (DCA) to the collision vertex of 1 cm is required for each candidate track, which helps to suppress contamination due to weak decays and tracks from secondary vertices [12,41]. Tracks must have at least 20 points used in track fitting out of the maximum of 45 hits possible in the TPC. To avoid multiple counting of split tracks, more than 52% of the total possible fit points are required.
When measuring scaled factorial moments, a large number of background effects such as conservation laws, Coulomb repulsion, resonance decays and experimental acceptance, will significantly influence the results [25,27,46,47]. These background contributions must be taken into account in the calculation of the SFMs. We implement the mixed event method to eliminate background contributions in our analysis, following its successful application in the NA49 and NA61/SHINE experiments [25,27,48]. Both the CMC model [28,49] and UrQMD model [49] calculations, have shown that the mixed event method effectively removes background contributions. For this purpose, an additional observable is defined as [25,27,28,48], where the moments from mixed events representing the background contributions are subtracted from the data. Mixed events are constructed by randomly selecting particles from different original events, while ensuring that the mixed events have the same multiplicity and mo-mentum distributions as the original events. We will exclusively use F q (M) instead of F q (M) in the following analysis.
Experimentally, the values of SFMs are influenced by the efficiency of the detector, since they are calculated from the measured multiplicity distribution of particles. To recover the true SFM from the experimentally measured one, the efficiency correction is calculated via the cell-by-cell method [46], which assumes a binomial response of the TPC and TOF detectors [46,[50][51][52]. According to the simulation of the STAR detectors, the detector response is close enough to the binomial distributions within statistical significance up to the 6th-order cumulants [3,13,53]. This also indicates that the impacts of the finite two-track resolution, such as track splitting and merging, are not significant for SFMs as well. The cell-by-cell method corrects the measured q-th moment, n i (n i -1) • • • (n iq + 1) , of SFM for each cell in transverse momentum space, one by one. The efficiency for each cell is calculated by averaging the p T -dependent efficiency of particles located in that cell. This cell-by-cell method has been validated by encoding the tracking efficiency of the STAR detector into the UrQMD event sample [46]. Statistical uncertainty is estimated using the Bootstrap method [54]. Systematic uncertainties are estimated by varying the fit range of M 2 and experimental requirements to reconstruct charged hadrons in the TPC and TOF. These requirements include the distance of closet approach, the track quality reflected by the number of fit points used in track reconstruction, and the dE/dx selection criteria for particle identifications.
We measure the SFMs of identified charged hadrons (h ± ) combining p, p, K ± , and π ± together. Particle identification is required to apply the efficiency correction on the SFMs. The analysis is performed in the measured p T range of 0.2 < p T < 2.0 GeV/c. The domain [-p x,max , p x,max ]⊗ [-p y,max , p y,max ] of the transverse momentum plane with p x,max = p y,max = 2.0 GeV/c is partitioned into M 2 cells to calculate the SFMs according to Eq. (1). Fig. 1 (a)-(d) shows F q (M) data and F q (M) mix corrected for reconstruction efficiency, from the second-order to the sixth-order, in the most central (0-5%) Au+Au collisions at √ s NN = 7.7-200 GeV. The event statistics of BES-I data allow calculating F q (M) up to the sixth order (q = 6) and in the range of M 2 from 1 to 100 2 . It is observed that F q (M) data (q = 2-6) are significantly larger than F q (M) mix in the large M 2 region (M 2 > 1000) at all √ s NN . Therefore, F q (M)
) is significantly larger than zero in the large M 2 region. F q (M) data was observed to overlap with F q (M) mix and F q (M) ≈ 0 from the UrQMD calculations [49], which cannot describe the presented data, since it does not incorporate density fluctuations induced by the QCD phase transition.
In Fig. 1 (e)-(h), the F q (M) (q = 2-6) are shown as a function of M 2 in the most central (0-5%) collisions at four example energies in the √ s NN = 7.7-200 GeV range. We find that F q (M) (q = 2-6) increases with increasing M 2 and reaches saturation when M 2 is large (M 2 > 4000). Therefore, F q (M) (q = 2-6) does not obey a power-law behavior of F q (M) ∝ (M 2 ) φ q over the entire M 2 range. Equivalently, the power-law scaling of F q (M) ∝ (M 2 ) φ q is not valid for the entire M 2 range, and φ q cannot be extracted in a reliable manner (independently of M 2 range). As a result, we will focus on the power-law behavior of F q (M) ∝ F 2 (M) β q and the scaling exponent. Fig. 2 shows
√ s NN = 7.7-200 GeV. We observe that F q (M) (q = 3-6) obey a strict power-law behavior versus F 2 (M) in the most central Au+Au collisions. This power-law scaling of F q (M) ∝ F 2 (M) β q is observed at all collision energies. It is worthwhile to note that one should perform the intermittency analysis only using SFMs in a sufficiently large M 2 region, The scaled factorial moments, F q (M)(q = 2-6), of identified charged hadrons (h ± ) multiplicity in the most central (0-5%) Au+Au collisions at four example energies in the √ s NN = 7.7-200 GeV range. Solid (open) markers represent F q (M) of data (mixed events) as a function of M 2 . (e)-(h) F q (M) (q = 2-6) as a function of M 2 in the most central (0-5%) Au+Au collisions at four example energies in double-logarithmic scale. Statistical uncertainties are obtained from the Bootstrap method.
F q (M) (q = 3-6) as a function of F 2 (M) in the most central Au+Au collisions at √ s NN = 7.7-200 GeV. The solid black lines represent the best power-law fit of F q (M) ∝ F 2 (M) βq with a fitting range of F 2 (M) from M ∈ [30,100]. The value of β q is the slope of the fitting line.
since intermittency is expected to occur at small momentum scale (i.e., small size of cell in momentum space) [17,22,55]. This way one can also avoid the influence of trivial fluctuations at large momentum scale on the determination of the intermittency exponent [25,27,49]. The value of β q is obtained through the best fit as the slope of the straight black line in Fig. 2. We note that β q did not significantly change even when the fitting range was varied. The analysis outlined above was also performed in other central centrality classes (5-10%, 10-20%, 20-30%, 30-40%). Nonetheless, significant statistical uncertainties of higher-order F q (M) (refer to Fig. 8 in Supplemental Material) prevented us to perform the entire chain of analysis in these narrow centrality classes below √ s NN = 19.6 GeV.
Fig. 3 (a) shows β q as a function of q -1 in the most central
√ s NN = 7.7-200 GeV. Note that the extraction of β q for the other narrow central centrality classes (from 5-10% to 30-40%) was prevented by the observed statistical uncertainties, hence we only present this result for the 0-5% centrality class in Fig. 3 (a). In agreement with theoretical expectation, β q also obeys a power-law behavior with q -1. The scaling exponent, ν, can be obtained through a best power-law fit of β q ∝ (q -1) ν . √ s NN = 7.7-14.5 GeV), since the higher-orders F q (M) exhibit statistical uncertainties to such amount that extraction of ν becomes impossible in numerous central centrality classes. As a result, the rest of the results will be presented with a merged centrality class (10-40%) as a baseline for comparison with the most central (0-5%) collisions.
Fig. 4 shows the energy dependence of ν for identified charged hadrons in Au+Au collisions for two different collision centralities (0-5% and 10-40%). In the most central collisions, ν exhibits a nonmonotonic behavior as a function of collision energy and reaches a minimum around √ s NN = 27 GeV. In contrast, ν is consistent with a constant with increasing √ s NN in 10-40% central collisions. The scaling index, β q (q = 3-6), as a function of q -1 in the most central (0-5%) Au+Au collisions at √ s NN = 7.7-200 GeV. The solid lines represent the best power- law fit of β q ∝ (q -1) ν . The statistical uncertainties of β q are smaller than the marker size. (b) The scaling exponent (ν), as a function of average number of participant nucleons ( N part ), in Au+Au collisions at √ s NN = 19.6-200 GeV. The data with the largest number of N part correspond to the most central collisions (0-5%), and the rest of the points are for 5-10%, 10-20%, 20-30% and 30-40% centrality, respectively. The systematic uncertainties of ν are shown in bars and the statistical uncertainties are smaller than the marker size. Both β q and ν at all energies are scaled by different factors. The observed non-monotonic energy dependence of ν in the most central collisions may be due to the signal of density fluctuations induced by the QCD critical point. At √ s NN ≤ 11.5 GeV, there are large systematic and statistical uncertainties for ν. Higher statistics data from the BES-II program [9] will help confirm the energy dependence of ν.
The measured value of ν in Fig. 4 is considerably smaller than the theoretical prediction of the critical ν= 1.30 from GL theory [18] and 1.0 from the 2D Ising model [31,37]. However, these calculations naturally use the entire phase space without any constraint on acceptance, whereas the measurements utilize only the experimentally available region of transverse momentum space within η and p T acceptance. The measured ν is expected to increase in case of a measurement performed in the entire phase space, in particular including higher p T regions, as indicated by the AMPT result which shows a rapid increase in intermittency or fluctuations with the increasing of p T [35]. Moreover, theoretical calculations that consider a reduced transverse momentum phase space and equivalent experimental acceptance, are required to understand the measured scaling exponent. The value of ν = 1.94 ± 0.10 at √ s NN = 200 GeV/c obtained from the AMPT calcula- tion [35], which does not take background subtraction into account and does not incorporate the physics of QCD phase transition, is significantly larger than the measured values. The transport-based UrQMD model is unable to calculate ν due to the absence of the power-law scaling of
new model that exhibits such power-law scaling is required to produce a model baseline for comparison with experimental data.
In summary, we have presented the first measurement of intermittency in heavy-ion collisions at RHIC. The transverse momentum phase space (p x , p y ) scaled factorial moments of identified charged hadrons combining p, p, K ± , and π ± within |η| < 0.5 have been calculated up to the sixth order in Au+Au collisions at √ s NN = 7.7-200 GeV. A distinct power-law scaling of F q (M) ∝ F 2 (M) β q , is observed in Au+Au collisions at all energies after background subtraction. Based on the scaling behavior, the scaling exponent (ν) is extracted and found to decrease monotonically from the peripheral to the central Au+Au collisions. The ν is consistent with a constant for different collisions energies in the midcentral (10-40%) collisions. A non-monotonic energy dependence is observed in the 0-5% most central collisions with ν reaching a minimum around √ s NN = 27 GeV. Whether the observed non- monotonic behavior is related to the CEP or not, further calculations from dynamical modelling of heavy-ion collisions with a realistic equation of state are required. Fig. 5. The scaled factorial moments, F q (M) (q = 2-6), of identified charged hadrons (h ± ) multiplicity in the most central (0-5%) Au+Au collisions at √ s NN = 7.7-200 GeV. Solid (open) markers represent F q (M) of data (mixed events) as a function of M 2 . Statistical uncertainties are obtained from the Bootstrap method.
Fig. 6. The F q (M) (q = 2-6) for identified charged hadrons as a function of M 2 in the most central (0-5%) Au+Au collisions at √ s NN = 7.
7-200 GeV. Statistical uncertainties are obtained from the Bootstrap method. Ministry of Science and Technology of China and the Chinese Ministry of Education, the Higher Education Sprout Project by Ministry of Education at NCKU, the National Research Foundation of Korea, Czech Science Foundation and Ministry of Education, Youth and Sports of the Czech Republic, Hungarian National Research, Development and Innovation Office, New National Excellency Programme of the Hungarian Ministry of Human Capacities, Department of Atomic Energy and Department of Science and Technology of the Government of India, the National Science Centre of Poland, the Ministry of Science, Education and Sports of the Republic of Croatia, German Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie (BMBF), Helmholtz Association, Ministry of Education, Culture, Sports, Science and Technology (MEXT) and Japan Society for the Promotion of Science (JSPS).
Appendix A. Supplemental material Fig. 5 shows F q (M) data and F q (M) mix as a function of M 2 in the most central (0-5%) Au+Au collisions at √ s NN = 7.7, 11.5, 14.5, 19.6, 27, 39, 54.4, 62.4 and 200 GeV. It is observed that F q (M) data (q = 2-6) are larger than F q (M) mix in the large M 2 region (M 2 > 1000) at all √ s NN . Fig. 6 shows F q (M) as a function of M 2 in the most central (0-5%) collisions at √ s NN = 7.7, 11.5, 14.5, 19.6, 27, 39, 54.4, 62.4 and 200 GeV. For all √ s NN , F q (M) (q = 2-6) increase with increasing M 2 , however, it reaches saturation when M 2 is large (M 2 > 4000). The F q (M) does not obey a power-law (scaling) behavior of F q (M) ∝ (M 2 ) φ q over the entire M 2 range. Equivalently, the relationship between F q (M) and M 2 is not linear over the whole M 2 range in double-logarithmic scale. Fig. 7 shows the comparison between the efficiency corrected F q (M) and uncorrected F q (M) for different orders in the most central (0-5%) Au+Au collisions at √ s NN = 27 GeV. The magnitude of the efficiency correction is found to be 19% for the second-order F 2 (M) (M 2 > 1000), and the correction becomes large as the order increases, reaching 51% for the sixth-order F 6 (M).
Fig. 8 (a)-(e) shows F q (M) data and F q (M) mix as a function of M 2 in numerous central centrality classes (0-5%, 5-10%, 10-20%, 20-30%, 30-40%) at √ s NN = 7.7 GeV. In addition, Fig. 8 (f)-(j) shows
F q (M) data as a function of M 2 in these centrality classes at the same √ s NN . The higher-order F 6 (M) has large statistical un-Fig. 7. Efficiency corrected and uncorrected F q (M) for various order (q = 2-6) as a function of M 2 in the most central (0-5%) Au+Au collisions at √ s NN = 27 GeV. The scaled factorial moments, F q (M) (q = 2-6), of identified charged hadrons (h ± ) multiplicity in 0-5%, 5-10%, 10-20%, 20-30%, 30-40% centrality classes at √ s NN = 7.7 GeV. Solid (open) markers represent F q (M) of data (mixed events) as a function of M 2 . (f)-(j) F q (M) (q = 2-6) as a function of M 2 in 0-5%, 5-10%, 10-20%, 20-30%, 30-40% centrality classes at √ s NN = 7.7 GeV.
certainties at larger M 2 regions, beginning from 5-10% centrality.
Therefore, the ν can not be calculated in these centrality classes at lower √ s NN = 7.7-14.5 GeV. Moreover, we can calculate ν in all central centrality classes (from 0-5% to 30-40%) starting from √ s NN = 19.6 GeV. As a result, we only show the centrality depen- dence of ν at higher √ s NN = 19.6-200 GeV in Fig. 3 (b). However, with higher statistics data from the BES-II program, we will be able to present the centrality dependence of ν at all √ s NN .