Measurements of charged-particle multiplicity dependence of higher-order net-proton cumulants in p + p collisions at √ s = 200 GeV from STAR at RHIC

We report on the charged-particle multiplicity dependence of net-proton cumulant ratios up to sixth order from √ s = 200 GeV p + p collisions at the Relativistic Heavy Ion Collider (RHIC). The measured ratios C 4 / C 2 , C 5 / C 1 , and C 6 / C 2 decrease with increased charged-particle multiplicity and rapidity acceptance. Neither the Skellam base-lines nor PYTHIA8 calculations account for the observed multiplicity dependence. In addition, the ratios C 5 / C 1 and C 6 / C 2 approach negative values in the highest-multiplicity events, which implies that thermalized QCD matter may be formed in p + p collisions.


Introduction
Exploring the quantum chromo-dynamics (QCD) phase structure is one of the ultimate goals in heavy-ion collision experiments.The QCD phase diagram is characterized in terms of baryon chemical potential (µ B ) for the x-axis and temperature (T ) for the y-axis.In the conjectured structure of the QCD phase diagram, there is a quark-gluon plasma (QGP) phase in the high temperature region [1] where the thermalized QCD matter is formed, while a hadron gas phase exists at lower temperature and higher baryon density [2].However, the transition between the two phases is not well-understood.According to lattice QCD calculations, the phase transition at µ B /T < 2 is likely a smooth crossover [3,4].Model calculations predict a first-order phase transition in the finite µ B region [5] as well as a QCD critical point [6].
The C 4 /C 2 value from the fixed-target experiment for 3 GeV collisions is consistent with hadron transport model calculations [17].For precise measurements at 3 < √ s NN < 20 GeV, the Beam Energy Scan program phase II (BES-II) was carried out from 2019 to 2021 both in the collider mode (7.7 ≤ √ s NN ≤ 19.6 GeV) and the fixed-target mode (3.2 ≤ √ s NN ≤ 7.7 GeV).
Higher-order cumulant ratios have been measured by STAR for Au+Au collisions at √ s NN = 200 GeV [18].In particular, the ratios C 6 /C 2 were found to be systematically negative from peripheral to central collisions.Albeit with large uncertainties, the negative values are qualitatively consistent with QCD models and lattice QCD calculations [19,20,21] which seems to favor a smooth crossover transition.This Letter reports the net-proton cumulant ratios,  [22,23,24,25].This scenario can be tested by measuring the multiplicity dependence of higher-order cumulant ratios.Although collectivity [22,23], collective motion of final particles with different masses such as the anisotropy flow v 2 , and strangeness enhancement [25], due to abundant production of gluons, have been observed in p+p collisions, it does not necessarily indicate that the thermalized QCD drops of matter have been created.On the other hand, higher-order cumulants of net-proton distributions offer more sensitive probe for the nature of thermalization in such collisions.This will be done by comparing various orders of cumulants between data and calculations from both thermal model and transport model [26].

Dataset
The data set was taken in 2012 using a minimum-bias trigger, with 220 million events analyzed.All data were measured using the Time Projection Chamber (TPC) and Time of Flight (TOF) detector at the STAR experiment [27].The collision vertex is required to be within 30 cm of the detector center and within 2 cm in the transverse direction relative to the beamline.Pileup events are suppressed by requiring the difference in the vertex position along the beamline measured by the TPC and the Vertex Position Detector to be within ±3 cm, and by requiring that tracks are matched to the TOF hits.Protons and antiprotons are analyzed in the transverse momentum range 0.4 < p T < 2.0 GeV/c and in the midrapidity acceptance |y| < 0.5.Figure 1 (a) shows the ionization energy loss (dE/dx) as a function of the momentum of charged particles, as measured by the TPC. Figure 1 (b) depicts correlations between mass squared (m 2 ) measured by the TOF and the momentum measured by the TPC.Protons and antiprotons are identified by using both dE/dx and m 2 .The distance of closest approach of a reconstructed track to the collision vertex is required to be less than 1 cm to suppress the contribution from secondary protons from weak decays.Cumulants are calculated at each bin of measured charged-particle multiplicity, m TT ch , which is defined as charged tracks having hits both in the TPC and TOF at midrapidity.Protons and antiprotons are excluded from m TT ch in order to suppress self-correlations [28].The m TT ch distribution is shown in Fig. 1 (c).The event-by-event net-proton number distributions for two ranges of m TT ch are shown in Fig. 1 (d).The TOF hit requirement is designed to suppress the contribution from collision pileup events in high-luminosity p+p collisions.Since the efficiency of the TOF is well understood in the experiment [11], m TT ch is converted to the corresponding multiplicity in the TPC, m TPC ch , in subsequent discussions, as the new results will be finally compared with previous measurements from Au+Au collisions [29].Correlations between the mass squared, m 2 , measured by the TOF and the ratio of momentum to the electric charge, p/q.Colored bands represent the identified protons and antiprotons used in the analysis.(c) Charged-particle multiplicity distribution.(d) Event-by-event net-proton multiplicity distributions for |y| < 0.5 and 0.4 < p T < 2.0 GeV/c at two ranges of charged particle multiplicity as indicated in the legend.

Observables
The nth-order cumulant is defined by the nth derivatives of the cumulant-generating function.Explicitly, cumulants are expressed in Eqs. ( 1)- (6).
where ⟨δN⟩ = N − ⟨N⟩, N is the number of particles in one event, and the angle brackets indicate the average over all events.The cumulant ratios, and C 6 /C 2 , are employed to cancel the trivial volume dependence of cumulants [30].Neutrons cannot be measured in the STAR experiment, hence net-proton number is measured as a proxy for net-baryon number [31].A Skellam distribution, which is defined as the difference between two independent Poisson distributions, is employed as a statistical baseline.The odd-and even-order cumulants of the Skellam distribution are expressed by the difference and the sum of the averaged values of protons and antiprotons, respectively.Consequently, the Skellam baselines for C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 are always unity, while the corresponding baselines for C 2 /C 1 and C 3 /C 2 depend on the averaged value of protons and antiprotons.
Statistical uncertainties are estimated by the bootstrap method [41,37].Systematic uncertainties are estimated by varying the selection ranges for track quality, particle identification, luminosity, and by varying the detector efficiencies for the efficiency corrections.The luminosity ranges from 2 kHz to 15 kHz for the coincidence rates for the Zero Degree Calorimeters, which is divided into 10 groups to estimate the variations of the cumulants.A Barlow test was performed in order to remove contribution from statistical uncertainties [42].The systematic uncertainties from each source are 0.29%, 0.32%, 0.38%, and 0.69%, respectively, for multiplicity-averaged results of C 4 /C 2 .The total systematic uncertainties are 0.90%, 7.8%, and 8.4% for C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 , respectively.

Results and discussions
Charged-particle multiplicity dependence for the net-proton cumulant ratios, C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 are shown in Fig. 2 as red dots.Vertical bars represent statistical uncertainties while the red bands represent systematic uncertainties.The main source of the large systematic uncertainties for each charged-particle multiplicity bin is the luminosity.It constitutes 50% out of total systematic uncertainty for C 4 /C 2 at m TPC ch = 4, which increases up to above 90% at m TPC ch > 18.The effect is more prominent for C 5 /C 1 and C 6 /C 2 , greater than 90% in most m TPC ch bins.Calculations from the Skellam baseline and PYTHIA8 model [43] are shown as black dashed lines and purple bands, respectively.Approximately 800 million PYTHIA8 events, with the option of SoftQCD, are generated for the comparison.Figure 2 indicates all ratios are decreasing with increasing multiplicity, except for C 2 /C 1 .Although the PYTHIA8 model shows multiplicity dependence, it does not provide quantitative description of data.For the higher-order ratios C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 , the model (purple squares) also overpredicts the multiplicity-averaged data values (cyan solid circles).As reported in Ref. [44], thermodynamic model calculations from both lattice QCD [3,21] and functional  renormalization group (FRG) [20] predict a special ordering of the higher-order baryon-number susceptibility ratios: The multiplicity-averaged ratios in Fig. 2 show the hierarchy where the significance of the first and second inequalities is 8.7 σ and 0.9 σ, respectively.At high charged-particle multiplicity bins m TPC ch > 12, both values of C 5 /C 1 and C 6 /C 2 are consistent with zero within uncertainties, which hints at a possible sign change in these ratios at even higher multiplicity.These observations are consistent with the thermodynamic model expectation for the formation of thermalized QCD matter in high-multiplicity p+p collisions.We discuss this point later when comparing with results from 200 GeV Au+Au collisions.Note that the PYTHIA8 calculations fail to reproduce the hierarchy in the ratios and no hint of a sign change is observed.At LHC energies, due to multi-parton interactions in p+p collisions, the color-reconnection (CR) mechanism is used to mimic collective excitation for strangeness and heavy-quark production [45,25,46,47,48,23].We confirmed that PYTHIA8 calculations with CR included do not have any significant effect on the net-proton higher-order cumulant ratios from the 200 GeV p+p collisions.
Studying the rapidity dependence of the cumulant ratios may shed light on the time evolution of the collision dynamics [49].The rapidity-window dependence of the net-proton cumulant ratios vs. multiplicity is shown in Fig. 3, with the colored markers representing data from three rapidity bins.In order to avoid overlapping uncertainties, data are plotted up to m TPC ch = 14.The higher-order ratios C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 show that the larger the rapidity acceptance, the more the ratios deviate from the corresponding ratio of the Skellam distribution (dashed lines) and the values are largely decreasing with increasing multiplicity.Both features imply the strongest correlation in the largest rapidity window.A similar rapidity-window dependence has also been observed in Au+Au collisions at RHIC [18].In addition, in the widest rapidity bin, |y| < 0.5, the ratios C 5 /C 1 and C 6 /C 2 approach negative values at the highest multiplicity, indicating that thermalized QCD matter may be created in very high-multiplicity events in p+p collisions at RHIC, as we noted for Fig. 2. None of the observations are reproduced by the QCD-inspired event generator PYTHIA8, which implies a lack of collective excitation dynamics in the model.Please note that the baryon number conservation could also lead to the decrease of higher-order cumulant ratios, especially for C 5 /C 1 and C 6 /C 2 , and the effect depends on the detector acceptance [50].Using the PYTHIA8 model and STAR acceptance, we find that both the ratios decrease as charged multiplicity increases, but the rate of the decrease is slower than that of data.More importantly, these ratios remain to be positive even for the highest multiplicity events.(More details can be found in the supplemental materials.)We also note that the baryon number conservation is present in PYTHIA8.Nevertheless, the ratios from PYTHIA8 does not show negative signs of C 5 /C 1 and C 6 /C 2 and the decreasing trend with respect to the multiplicity is not as prominent as the experimental results.Net-proton ratios C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 are compared with those from Au+Au collisions at √ s NN = 200 GeV [11,44] in Fig. 4, shown as a function of the charged-particle multiplicity.The kinematic acceptance 0.4 < p T < 2.0 GeV/c and |y| < 0.5 is used in both p+p and Au+Au collisions.The five data points for Au+Au collisions correspond to 0-40%, 40-50%, 50-60%, 60-70%, and 70-80% centrality bins.The multiplicity averaged ratios from the p+p collisions seem to follow the trend in the peripheral Au+Au collisions.However, the slope of the multiplicity dependence in the p+p range, 5 < m TPC ch < 20, is steeper than that in the range of 25 < m TPC ch < 460, for Au+Au collisions.The steeper slope indicates that thermalized QCD matter are created more efficiently in the high multiplicity p+p collisions than that in heavy-ion collisions.In particular, both the ratios C 5 /C 1 and C 6 /C 2 from p+p collisions decrease with multiplicity, approach negative values for the highest-multiplicity events.The ratios from Au+Au collisions at similar multiplicity region are positive, which could be attributed as the initial volume fluctuations in heavy-ion collisions [51,52,53].The negative ratios in the context of lattice QCD calculations imply that even in the small p+p system, thermalized QCD matter may be created in the highest-multiplicity collisions, if the trend in multiplicity dependence of the cumulant ratios reflects the thermalization of the bulk matter as that in Au+Au collisions.Similar phenomena have been observed in the multiplicity dependence of strange hadron and J/ψ production [25,46], and long-range collective motion [22,23] in p+p collisions at LHC energies.Overall, the high-order net-proton cumulant ratios from both 200 GeV p+p and Au+Au collisions show a clear decreasing trend from low to high chargedparticle multiplicity, and eventually reach values consistent with lattice QCD calculations assuming µ B = 25 MeV and T = 155 MeV [21], within uncertainties.Measurements of the cumulant ratios at high center-of-mass energies or in larger collision system than p+p, (e.g.p+Au, d+Au, Zr+Zr, and Ru+Ru collision systems) could provide systematic information on the dynamics of the observed multiplicity dependence.[21], where the multiplicity range is chosen arbitrary.

Summary
In summary, we report the first measurements of higher-order cumulant ratios C 2 /C 1 , C 3 /C 2 , C 4 /C 2 , C 5 /C 1 and C 6 /C 2 of net-proton multiplicity distributions in p+p collisions at √ s = 200 GeV, as measured by the STAR detector at RHIC.Both charged-particle multiplicity and rapidity-cut dependencies are reported.It is found that the ratios are all below Skellam expectations.Overall, the net-proton cumulant ratios C 4 /C 2 , C 5 /C 1 and C 6 /C 2 decrease progressively between the 200 GeV p+p collisions and Au+Au central collisions.Calculations from PYTHIA8 [43] fail to reproduce the multiplicity dependence and the hierarchy observed in net-proton ratios C 4 /C 2 , C 5 /C 1 and C 6 /C 2 .The C 5 /C 1 and C 6 /C 2 values are consistent with zero within uncertainties at the highest multiplicity bins with the maximized rapidity acceptance.Negative values of C 5 /C 1 and C 6 /C 2 are predicted, according to lattice QCD calculations, for thermalized nuclear matter [3,4].These new results imply that the thermalized QCD drops of matter could be created in truely high multiplicity events of the p+p collisions.Systematic measurements of the cumulant ratios from higher collision energies, in larger colliding systems, and wider acceptance in rapidity will provide important information on the underlying dynamics of high-order net-proton cumulants and the process of thermalization in high-energy collisions.

Effect of baryon number conservation
As discussed in Ref. [50], it is known that the baryon number conservation artificially decreases the higher-order cumulants of net-particle distribution.To study the effect on our measurements, we performed numerical calculations with the procedures as follows. .This is because the mean values of (anti)protons within the rapidity acceptance is much smaller for m ch = 6 compared to the higher m ch events.All ratios decrease with increasing the acceptance factor.We find that, at the highest multiplicity events in the experiment (m ch = 20), the ratios are C 5 /C 1 ≈ 0.08 and C 6 /C 2 ≈ 0.2.These ratios are not negative even in the further higher multiplicity, m ch = 30.In Fig. 6, we overlaid

Figure 1 :
Figure 1: (a) Correlations between the energy loss, dE/dx, of charged tracks measured by the TPC and momentum divided by electric charge.(b)Correlations between the mass squared, m 2 , measured by the TOF and the ratio of momentum to the electric charge, p/q.Colored bands represent the identified protons and antiprotons used in the analysis.(c) Charged-particle multiplicity distribution.(d) Event-by-event net-proton multiplicity distributions for |y| < 0.5 and 0.4 < p T < 2.0 GeV/c at two ranges of charged particle multiplicity as indicated in the legend.

Figure 2 :
Figure 2: Net-proton cumulant ratios, (a) C 2 /C 1 , (b) C 3 /C 2 , (c) C 4 /C 2 , (d) C 5 /C 1 , and (e) C 6 /C 2 as a function of charged-particle multiplicity from √ s = 200 GeV p+p collisions.Black solid lines and red bands represent the statistical and systematic uncertainties, respectively.Cyan points represent event averages for 3 < m ch < 30, and they are plotted at the corresponding value of m TPC ch .The uncertainties on the cyan points are smaller than the marker size.The Skellam baselines are shown as dashed lines.The results of the PYTHIA8 calculations are shown by hatched-golden bands.The golden bands at m TPC ch ≈ 6 are the results from the PYTHIA8 calculations averaged over multiplicities.

Figure 3 :
Figure 3: Rapidity-acceptance dependence of the net-proton cumulant ratios, (a) C 4 /C 2 , (b) C 5 /C 1 , and (c) C 6 /C 2 shown as a function of the charged-particle multiplicity, m TPC ch , from √ s = 200 GeV p+p collisions.Solid squares, triangles, and circles represent the results with |y| < 0.1, 0.3, and 0.5, respectively.Black solid lines and colored bands show the statistical and systematic uncertainties.The data points and their uncertainties for |y| < 0.1 and 0.3 are slightly shifted for clarity.The Skellam baselines are shown by dashed lines.

Figure 4 :
Figure 4: Net-proton cumulant ratios, (a) C 4 /C 2 , (b) C 5 /C 1 , and (c) C 6 /C 2 as a function of charged-particle multiplicity for p+p collisions and Au+Au collisions at 200 GeV for |y| < 0.5 and 0.4 < p T < 2.0 GeV/c.Cyan markers represent event averages for 3 < m TPC ch < 30 from the p+p collisions.Results from Au+Au collisions are shown as triangles for the 0-40%, 40-50%, 50-60%, 60-70%, and 70-80% centrality bins.Red and grey bands show the systematic uncertainties for p+p collisions and Au+Au collisions, respectively.The m TPC ch values are not corrected for the reconstruction efficiency of the TPC.The golden bands at m TPC ch ≈ 6 represent the results from PYTHIA8 calculations averaged over multiplicity.The Skellam baselines are shown in dotted lines.The purple bands show corresponding susceptibility ratios of baryon number from lattice QCD calculations[21], where the multiplicity range is chosen arbitrary.

1 .
Generate baryons and antibaryons based on the Poisson distribution, where their mean values are taken from PYTHIA8 calculations for the reference multiplicity m ch = 6.The m ch is defined within the same acceptance for the experimental analysis.2. Events are proceeded only if the net-baryon number is 2, to incorporate the baryon number conservation.3. Produced (anti)baryons are randomly sampled based on the acceptance factor with the binomial response to determine the (anti)protons in the experimental acceptance.The acceptance factor is determined for |y| < 0.1, 0.3, and 0.5, by PYTHIA8 calculations.4. The above 1)-3) are repeated for many events, then we calculate net-proton C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 .5. 1)-4) are repeated by varying the reference multiplicity for m ch = 14, 20, and 30.Note that m ch = 20 corresponds to the highest multiplicity that we can access in our experimental data.The results are shown in Fig. 5. Different lines/markers represent the reference multiplicity m ch .Each line has three points, corresponding to |y| < 0.1, 0.3, and 0.5 from left to right.It is seen that the results from different m ch are mostly on the same trend, except for m ch = 6

Figure 5 :
Figure 5: Results from numerical calculations on net-proton C 4 /C 2 , C 5 /C 1 , and C 6 /C 2 as a function of acceptance factor.Different markers represent four charged particle multiplicity bins chosen in PYTHIA8 calculations.Three data points for each line corresponds to the proton rapidity acceptance of |y| < 0.1, 0.3, and 0.5.

Figure 6 :
Figure 6: Net-proton cumulant ratios as a function of reference multiplicity.Blue points represent the results from the numerical calculations incorporating Poisson fluctuations with the baryon number conservation.