Energy Dependence of Intermittency for Charged Hadrons in Au+Au Collisions at RHIC

Density fluctuations near the QCD critical point can be probed via an intermittency analysis in relativistic heavy-ion collisions. We report the first measurement of intermittency in Au$+$Au collisions at $\sqrt{s_\mathrm{_{NN}}}$ = 7.7-200 GeV measured by the STAR experiment at the Relativistic Heavy Ion Collider (RHIC). The scaled factorial moments of identified charged hadrons are analyzed at mid-rapidity and within the transverse momentum phase space. We observe a power-law behavior of scaled factorial moments in Au$+$Au collisions and a decrease in the extracted scaling exponent ($\nu$) from peripheral to central collisions. The $\nu$ is consistent with a constant for different collisions energies in the mid-central (10-40\%) collisions. Moreover, the $\nu$ in the 0-5\% most central Au$+$Au collisions exhibits a non-monotonic energy dependence that reaches a possible minimum around $\sqrt{s_\mathrm{_{NN}}}$ = 27 GeV. The physics implications on the QCD phase structure are discussed.

Density fluctuations near the QCD critical point can be probed via an intermittency analysis in relativistic heavy-ion collisions. We report the first measurement of intermittency in Au+Au collisions at √ s NN = 7.7-200 GeV measured by the STAR experiment at the Relativistic Heavy Ion Collider (RHIC). The scaled factorial moments of identified charged hadrons are analyzed at mid-rapidity and within the transverse momentum phase space. We observe a power-law behavior of scaled factorial moments in Au+Au collisions and a decrease in the extracted scaling exponent (ν) from peripheral to central collisions. The ν is consistent with a constant for different collisions energies in the mid-central (10-40%) collisions. Moreover, the ν in the 0-5% most central Au+Au collisions exhibits a non-monotonic energy dependence that reaches a minimum around √ s NN = 27 GeV. The physics implications on the QCD phase structure are discussed.

I. INTRODUCTION
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 firstorder 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 powerlaw 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, qth-order 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 powerlaw behavior of F q (M ) ∝ (M D ) ϕq , M ≫ 1, 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: [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 158A GeV [27]. Two studies, using a Critical Monte Carlo (CMC) model [26] and a cascade ultra-relativistic quantum molecular dynamics (UrQMD) model with hadronic potentials [40], respectively, suggested that large intermittency could be observed in Au+Au collisions at RHIC energies.          lected 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 selfcorrelation [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]. 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 andp 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 andp 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.

II. EXPERIMENT AND DATA ANALYSIS
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 momentum 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 6thorder 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 i −q+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-bycell 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.

III. RESULTS AND DISCUSSIONS
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.
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   It is worthwhile to note that one should perform the intermittency analysis only using SFMs in a sufficiently large M 2 region, 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. 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) ν . Figure 3 (b) shows the extracted ν as a function of ⟨N part ⟩ in various centrality classes in Au+Au collisions at √ s NN = 19.6-200 GeV. It is observed that ν decreases monotonically from mid-central (30-40%) to the most central (0-5%) Au+Au collisions. Figure 3 (b) does not display the centrality dependence of ν at lower collision energies ( √ 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. Figure 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 non-monotonic 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 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 includ- ing 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 calculation [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 ∆F q (M ) ∝ ∆F 2 (M ) βq [49]. Consequently, a new model that exhibits such power-law scaling is required to produce a model baseline for comparison with experimental data.

IV. SUMMARY
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 mid-central (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.  rection 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 ). 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 uncertainties 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 dependence 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 .