From the theory side, the heavy-ion collisions are usually described by relativistic hydrodynamics [14][15][16][17]. At a socalled particlization stage [18], the QCD fluid is transformed into an expanding gas of hadrons and resonances, based on the picture of grand-canonical hadron resonance gas (HRG). Indeed, this picture works quite well to describe bulk properties of measured hadrons, such as spectra and flow, in a broad range of collision energies [19][20][21][22][23]. However, a quantitative theoretical analysis of fluctuations in this picture is challenging. In contrast to mean hadron yields, the event-by-event fluctuations, especially the high-order cumulants, are influenced by several effects which make direct application of the grand-canonical statistical mechanics questionable. The most important effects are the global conservation laws [24][25][26][27][28] and the smearing of fluctuations due to momentum cuts [29,30]. Other mechanisms include volume fluctuations [31][32][33] or diffusion in the hadronic phase [34,35].

The two main issues mentioned above were recently addressed in Ref. [36] at LHC energies via a generalized Cooper-Frye procedure called the subensemble sampler, utilizing the approximately boost invariant nature of heavyion collisions at the LHC parametrized by the blast-wave model. In this paper we extend this method to the RHIC-BES energies. This is achieved in the following way. First, we relax the assumption of boost invariance. Instead, we employ realistic particlization hypersurfaces obtained from (3+1)-dimensional viscous hydrodynamics simulations using code MUSIC [37][38][39]. Second, we calculate the cumulants of (anti)proton number distribution emitted from the hypersurface subject to global baryon conservation analytically (rather than using Monte Carlo as in Ref. [36]). For this purpose we use a recently developed subensemble acceptance method 2.0 (SAM-2.0) [40], which allows one to perform a correction of cumulants of accepted protons for the effect of exact global baryon conservation analytically. As a result, we are able to calculate cumulants of (anti)proton number distribution emerging from hydrodynamics to high orders without the need to generate large numbers of Monte Carlo events. The results are then confronted with the experimental data of the STAR Collaboration.

The paper is organized as follows. Section II presents the method to calculate cumulants of proton and baryon number distribution at particlization. The calculation results for Au-Au collisions at RHIC-BES energies are presented in Sec. III. Discussion and summary in Sec. IV close the paper.

We employ relativistic viscous hydrodynamics to simulate the evolution of a system created in 0-5% Au-Au collisions at RHIC, using the open-source code MUSIC v3.0 [19,37,41]. We perform hydrodynamic simulations with event averaged initial density profiles at each collision energy, which describes the expansion of quark-gluon plasma created in the earlier stage of the collision and its transition to a hadron gas. Cumulants of proton and baryon number distributions are computed at the end of hydrodynamic evolution at particlization.

The three-dimensional initial conditions are taken from Ref. [39]. They are based on the Glauber collision geometry, employ local energy and momentum conservation, and are calibrated to reproduce the measured proton transverse momentum distributions and midrapidity yields at different collision energies. This makes it suitable for the analysis of second-and higher-order proton cumulants, which we perform here.

The hydrodynamic evolution is calculated in MUSIC v3.0 by numerically solving the equations corresponding to energymomentum and baryon number conservation, as well as Israel-Stewart-type relaxation equations describing the viscous stress tensor. We include shear viscous corrections but neglect bulk viscous corrections and baryon diffusion. We employ a NEOS-BSQfoot_0 equation of state from Ref. [42] which interpolates between the lattice QCD equation at large temperatures [43], described via the Taylor expansion, and hadron resonance gas at low temperatures. This equation of state imposes vanishing net strangeness, n S = 0, and a fixed charge-to-baryon ratio, n Q /n B = 0.4, in every fluid cell. The shear viscosity to entropy baryon ratio is temperature and chemical-potential dependent; the details can be found in Fig. 4 of Ref. [39]. The hydrodynamic equations are solved in Milne coordinates, (τ, x, y, η s ), and evolved in τ until all computational cells reach a threshold energy density ε sw for particlization.

The hydrodynamic evolution ends at a particlization hypersurface σ (x) of constant "switching" energy density ε sw . The value of ε sw = 0.26 GeV/fm 3 has been adjusted in Ref. [39] to fit bulk observables and it is used here throughout. Note that a further improvement of the proton spectra at √ s NN 39 GeV can be achieved by increasing the value of ε sw at those energies [44]. In our calculations we observed a mild effect of a changing ε sw on proton cumulants (this is detailed in Appendix A), mainly in the form of stronger excluded volume effects at larger densities; thus we keep ε sw = 0.26 GeV/fm 3 at all energies for consistency. The QCD fluid is transformed at particlization into a system consisting of hadrons and resonances. The momentum spectrum for hadron species j emerging from hydrodynamics is given by the Cooper-Frye formula [45] foot_1

with

Here λ j (x) describes deviations from the ideal gas distribution function due to interactions. We assume that these deviations are momentum independent. In the ideal HRG limit one has λ j (x) = 1. On the other hand, this factor is different from unity in case of a nonideal HRG. Here we employ the excluded volume HRG (EV-HRG) model with repulsive baryon-baryon interactions [47,48], which has been observed to provide an improved description of the lattice QCD data on baryon number susceptibilities near the pseudocritical temperature T PC ∼ 155 MeV at µ B = 0 compared to the standard ideal HRG. This model has been used in our previous study of proton fluctuations for the LHC energies [36] and we refer to that work for further details on the EV-HRG model. We also perform calculations in the ideal HRG limit to establish the relevance of baryon repulsion in the investigated observables.

We use the open source package THERMAL-FIST v1.3 in all our HRG model calculations [49].

The particlization hypersurfaces, consisting of a list of fluid elements each containing the normal four-vector dσ µ , the fluid four-velocity u µ , as well as energy and baryon densities, are available via Ref. [50]. For each hypersurface element we recalculate the values of the temperature T (x) and the chemical potential µ B (x) such that the HRG model equation of state at these T -µ B values matches the energy and baryon densities corresponding to this hypersurface element. We also enforce n Q /n B = 0.4 and n S = 0 for each fluid element to determine the electric charge and strangeness chemical potentials µ Q and µ S .

With the numerical output from MUSIC the Cooper-Frye integral becomes a sum over all fluid elements:

In what follows we neglect the modification of the shape of (anti)baryon spectra due to resonance decays and evolution in the hadronic phase. All baryons are modeled as thermal particles with nucleon mass m N = 0.938 GeV emitted from a Cooper-Frye hypersurface. These simplifications make it feasible to evaluate proton number cumulants analytically. We relax these assumptions in the Appendix B using a Monte Carlo approach and show that they have only small influence on the resulting proton number cumulants, at least up to the third order. Equation ( 3) is sufficient to calculate the number of (anti)baryons in a given momentum acceptance by integrating over the momenta. To calculate fluctuations, however, we need a generalization beyond the standard Cooper-Frye prescription.

Let us first calculate the cumulants in the grand-canonical limit, i.e., neglecting the exact global conservation of the baryon number. We shall correct the cumulants for the baryon number conservation via the recently developed method of Ref. [40] afterwards. We further assume that the dynamical correlation length ξ that defines the range of correlations is smaller than any other relevant scale, such that one can assume ξ → 0. This is in analogy to the model of critical fluctuations at freeze-out studied in Refs. [29,51]. In our case, where particle number correlations in the grand-canonical ensemble are attributed purely to the excluded volume effect, the emission of particles from all the hypersurface elements proceeds independently. To calculate the cumulants of (anti)baryon number distribution inside a particular momentum acceptance it is thus sufficient to sum up contributions from all the volume elements independently. The number of (anti)baryons emitted from a hypersurface element i fluctuates in accordance with the grand-canonical susceptibilities χ B ± of (anti)baryon number fluctuations. The corresponding cumulants, therefore, read

Here

is the effective volume of a hypersurface element i. The probability p acc (x i ; p acc ) that an (anti)baryon emitted from a fluid element i ends up in a momentum acceptance p acc can be calculated from the Cooper-Frye formula (3):

. (5) The contribution δκ B ± ,GCE n (x i ; p acc ) of element i to the nth-order cumulant of the accepted (anti)baryons is obtained by convoluting the cumulants {δκ B ± ,GCE l (x i )}, l = 1 . . . n with the binomial distribution with probability p acc (x i ; p acc ). One obtains (see, e.g., Ref. [52])

Here φ ≡ φ(t, p acc ) = ln(1 -p acc + e t p acc ) and B n,l are partial Bell polynomials.

In the same way we can also obtain the (anti)baryon cumulants δκ B ± ,GCE n (x i ; p acc ) outside the acceptance p acc , i.e., for p ∈ p acc . This is achieved by substituting p acc (x i ; p acc ) → 1 -p acc (x i ; p acc ) in Eq. (6).

To obtain (anti)proton number cumulants one can apply the arguments of Kitazawa and Asakawa [53,54]: Based on the isospin randomization in the hadronic phase, the (anti)proton cumulants are obtained by the binomial filtering of the (anti)baryon cumulants. Note that this argument does not necessarily require a long-lasting hadronic phase with many pion-nucleon scatterings to randomize the isospin. The binomial filtering is valid also in the case of models where primordial correlations between baryons do not depend on the isospin. This is the case for the EV-HRG model studied here, where the same excluded volume parameter is used for all baryon pairs, regardless of their isospins. The probability that a randomly chosen (anti)baryon is an (anti)proton is simply the ratio between the mean numbers of (anti)protons and (anti)baryons, q(x i ) = N p ± (x i ) N B ± (x i ) . The value of q(x i ) is calculated using the HRG model in each hypersurface element. To be consistent with the experimental conditions realized in the STAR experiment, we include all weak decay contributions [44]. In practice, this yields q(x i ) ≈ 1/2 in most cases.

We shall use the binomial distribution argument to obtain the joint baryon-proton cumulants δκ B ± ,p ± n,m (x i ; p acc ) in terms of baryon cumulants δκ B ± ,GCE n (x i ; p acc ) and the proton-tobaryon ratio q(x i ), for each hypersurface element x i . We calculate the joint cumulants because these quantities will later be needed to apply the correction for baryon number conservation using the method of Ref. [40]. Given the probability P(N B ) to have N B baryons the joint probability P(N B , N p ) to have N B baryons and N p protons is

where

is the binomial distribution. The joint cumulant generating function for N B and N p reads

Here G B is the cumulant generating function of the N B distribution and

To obtain Eq. ( 11) we used an identity

of the baryon-proton distribution correspond to the Taylor expansion coefficients of the generating function G(t B , t p ) around t B = t p = 0. The corresponding derivatives of G(t B , t p ) are evaluated with the help of Faà di Bruno's formula and expressed in terms of the cumulants κ B n of the N B distribution:

The same procedure applies for the joint cumulants of the antiproton-antibaryon distribution. Rewriting Eq. ( 12) for the cumulants corresponding to the accepted particles emitted from volume element x i we get

). (13) The joint cumulants of all accepted (anti)baryons/(anti) protons are obtained by summing over the contributions of all the hypersurface elements:

The joint net-baryon/net-proton cumulants can be obtained straightforwardly in the case when there are no grandcanonical correlations between baryons and antibaryons. This is the case for the EV-HRG model used here. The cumulants read

We also list here, for completeness, the joint net-baryon/ proton and net-baryon/antiproton cumulants

( p acc ), ( 16)

The joint net-baryon/(net-)(anti-)proton cumulants outside the acceptance are obtained in the same fashion, by substituting p acc (x i ; p acc ) → 1 -p acc (x i ; p acc ).

To account for the exact global baryon conservation we apply a generalized version of the subensemble acceptance method (SAM) developed in Ref. [40]. The SAM-2.0 allows one to calculate the effect of global conservation of a conserved quantity, such as net baryon number, on the cumulants of arbitrary nonconserved quantity, such as (net) proton number. The original SAM [55,56] was formulated for uniform thermal systems in the thermodynamic limit and acceptances in the coordinate space. The SAM-2.0 extends the method to nonuniform systems and momentum space acceptances. The method takes joint net-baryon/(net-)(anti-)proton number cumulants calculated inside and outside the acceptance without the account of exact baryon conservation and produces the cumulants that are subject to global baryon conservation, i.e., it provides a mapping:

Here p acc corresponds to particles outside the momentum acceptance p acc . The details of the procedure to calculate the mapping S are presented in Ref. [40]. The grand-canonical cumulants κ B,p,GCE n,m ( p acc ) and κ B,p,GCE n,m ( p acc ) entering the right-hand side of Eq. ( 19) were calculated in the previous subsection.

It is assumed in SAM-2.0 that the system is sufficiently large, such that the means of all the relevant quantities correspond to the maximum of probability distribution (see Ref. [40] for details). This assumption is realized in central Au-Au collisions at RHIC-BES. The second assumption of the method is that the distributions inside and outside the acceptance are independent, i.e., that cumulants κ B,p,GCE n,m ( p acc ) and κ B,p,GCE n,m

( p acc ). This assumption is satisfied exactly in the ideal HRG model and to a high accuracy within the EV-HRG model at RHIC-BES energies. As discussed in Ref. [40], even in an extreme case where the additivity of cumulants does not hold, the SAM-2.0 results exhibit only small deviations from the exact result, thus the possible deviations from the true results for the EV-HRG model applications considered in the present paper are expected to be negligible. The results presented in the next section that incorporate the effect of global baryon conservation have all been obtained using SAM-2.0. A MATHEMATICA notebook which calculates the mapping S is available via Ref.

[57] and used in this paper. on sw . The net proton, proton, and antiproton cumulants are calculated in the relevant momentum acceptances analytically, following the method presented in the previous section. We perform separate calculations employing EV-HRG and ideal HRG models, and study the behavior of cumulants with and without the correction for baryon number conservation. These different configurations allow us to establish the relevance of repulsive interactions and global baryon conservation. In Appendix B we perform a cross-check of the analytic results for the case of the ideal HRG model by means of Monte Carlo sampling of hadrons at particlization.

First we look at the net proton rapidity distributions. To calculate the net proton dN/dy we partition the rapidity axis into bins. The rapidity density in a given bin then corresponds to the first net proton cumulant evaluated for that bin. As dN/dy is determined by the mean numbers of particles, it is unaffected by the correction for global baryon conservation. We observe that net proton rapidity distributions calculated in the EV-HRG and ideal HRG models virtually coincide. This is attributed to the fact that we match the net baryon density at particlization to the MUSIC output in both scenarios, which leads to virtually identical mean numbers of net protons.

The resulting rapidity distributions are depicted in Fig. 1 for various collision energies. The results agree qualitatively with earlier MUSIC calculations in Ref. [39]. The calculations also agree within errors with the midrapidity net proton yields measured by the STAR Collaboration [58][59][60] The net proton rapidity distributions show peaks in the forward-backward rapidity regions at all collision energies except for 7.7 GeV. This reflects the fact that large rapidities probe baryon-rich matter. It is observed, for instance, that larger longitudinal space-time rapidities are characterized by larger values of the baryochemical potential µ B at particlization. This underlines the fact that it is impossible to characterize the whole fireball by a single pair of temperature T and baryon chemical potential µ B but instead one has to integrate over different µ B -T values encompassing the hypersurface.

The leading four cumulants of net proton distribution have been measured and presented by the STAR Collaboration in Ref. [9] as a function of collision energy. The measurements were performed in the momentum acceptance |y| < 0.5 and 0.4 < p T < 2.0 GeV/c. Here we calculate these cumulants in the same momentum acceptance. The results are presented in Fig. 2. The calculations are performed within the EV-HRG model with (solid red lines) and without (dotted black lines) the effect of exact baryon number conservation. We also perform a calculation within the ideal HRG model but including the effect of baryon number conservation (dash-dotted red line). The dashed blue lines correspond to the uncorrelated (anti)proton baseline, which is given by the Skellam distribution, i.e., κ

The first net-proton cumulant is unaffected by the correction for global baryon conservation. It is also virtually unaffected by the excluded volume effects due to the matching of the EV-HRG model equation of state to the net baryon density at particlization. The model provides a reasonable description of the experimental data, with the possible exception of the 19.6-and 27-GeV points. The description of these data points can be improved by fine-tuning the simulation parameters.

The second, third, and fourth cumulants are affected by both the excluded volume and baryon conservation, the latter effect being the stronger of the two. Both effects suppress the cumulants, and the suppression is stronger for higher-order cumulants and lower collision energies. The effect of excluded volume is stronger at lower collision energies because they probe baryon-rich matter with smaller interparticle distances between baryons at particlization. The baryon conservation plays a larger role at smaller energies because a larger fraction of the full baryon number ends up in the midrapidity region, which is where the measurements are performed. Compared to the STAR data, the calculation with simultaneous excluded volume and baryon conservation effects generally yields the best agreement. The agreement is not perfect everywhere; in particular κ 2 [pp] is notably overestimated by the model at √ s NN 19.6 GeV. This is a reflection of an overestimated mean number of protons and antiprotons produced by the model compared to the data. There are different possible reasons for this. For instance, if weak decay contributions are overestimated in the model calculation, this may explain the discrepancy. Although it has been argued that the integrated yields of (anti)protons measured by STAR contain essentially all weak decay contributions [44], the situation might be different in the measurements of fluctuations. We performed calculations of the cumulants neglecting all weak decay contributions, and indeed obtain an improved description of the data at some of the energies, although in this case the data are generally underestimated by the model. We did observe, however, that weak decays affect only mildly the various volume-independent ratios of cumulants, thus we keep all weak decay contributions in all our further results throughout this paper. Other potential reasons that may contribute to the overestimation of (anti)proton yields include neglecting baryon annihilation in the hadronic phase [64], or a reflection of the so-called thermal proton yield anomaly in the HRG model [65,66]. As all these possible mechanisms are not linked directly to the QCD equation of state, we shall not study them in detail here but instead look at observables where the effect of describing the total (anti)proton yield is minimized.

We thus analyze the following ratios of cumulants:

These two ratios have a baseline of unity in the (Skellam) limit of an uncorrelated production of protons and antiprotons, at any collision energy. Deviations from unity can be a signal of new physics; in particular, the QCD critical point has been predicted to have a strong influence on these non-Gaussian measures of net proton number fluctuations [7,67,68]. In this sense, Sσ 3 /M is more convenient than the commonly adopted normalized skewness sσ 2 ≡ κ 3 /κ 2 which shows strong collision energy dependence even in the Skellam limit. Note that the absolute yields of (anti)protons drop out in the ratios of cumulants, thus the ratios are not very sensitive to possible inaccuracies in the description of the overall yields discussed above. The ratios, however, are sensitive to both the excluded volume and baryon number conservation. Figure 3 depicts the collision energy dependence of Sσ 3 /M and κσ 2 calculated in the model along with the experimental data of the STAR Collaboration from Ref. [9]. Both the data and the model calculations show a suppression of Sσ 3 /M with respect to the Skellam baseline of unity. When baryon excluded volume, but not global baryon conservation, is incorporated (dotted black lines), this leads to an improved agreement with the data compared to Skellam, although the suppression of Sσ 3 /M due to the excluded volume is not sufficient to obtain a quantitative agreement. Calculations that incorporate global baryon conservation but not excluded volume (dash-dotted red lines) indicate that the former effect is stronger than the latter one. In this case the data at √ s NN 20 GeV are described but not at higher collision energies. Finally, when both the baryon conservation and excluded volume are incorporated, the experimental data at √ s NN 20 GeV are described on a quantitative level. On the other hand, the data at lower collision energies are underestimated. It should be noted that the magnitude of the excluded volume effects in the EV-HRG model that we use has been constrained using lattice QCD data at µ B = 0 in Ref. [48]. Thus, we expect the model to be most reliable at the highest collision energies that probe the QCD phase diagram close to the vanishing net baryon density. Deviations from the data at √ s NN 20 GeV may be an indication of a breakdown of the EV-HRG model that we use. We explore this issue further in the next subsection by studying proton correlation functions. The behavior of the net proton kurtosis κσ 2 is qualitatively similar to Sσ 3 /M, although the error bars are considerably larger, especially at the lower collision energies. This precludes making strong conclusions from the available data on κσ 2 from RHIC-BES-I; it should however be possible to use the more precise future data from RHIC-BES-II for this purpose. Figure 4 depicts our predictions for the net proton hyperkurtosis, κ 6 /κ 2 . This quantity is strongly suppressed by both the excluded volume and baryon conservation, and it turns negative as the collision energy is decreased to below √ s NN 40-60 GeV. These predictions can be probed by future high-statistics measurements at RHIC.

We also compare our results with a noncritical baseline of Ref. [27], which is based on the ideal HRG model and the empirical rapidity distributions. These results, shown in Figs. 3 and4 by the blue points, agree fairly well with our ideal HRG model results with exact baryon conservation (dash-dotted red lines). Thus, there is a consistency between our ideal HRG model based calculations and the prior literature.

More information can be obtained by analyzing cumulants of proton and antiproton distributions separately. In particular, one can construct the so-called correlation functions (factorial cumulants) Ĉn from the ordinary cumulants κ n [69]. The correlation functions characterize genuine multiparticle correlations. For n > 1 they vanish in the limit of Poisson statistics (uncorrelated particle production). The correlation functions are linear combinations of the cumulants:

Experimental measurements of both the cumulants κ n and the correlation functions Ĉn have recently been presented by the STAR Collaboration in Ref. [10]. Figure 5 depicts the comparison of our model calculations with the experimental data, in terms of normalized quantities, κ n /κ 1 -1 and Ĉn / Ĉ1 . 3Deviations of these quantities from zero signal physics beyond the uncorrelated gas of hadrons.

The normalized second-order cumulants and correlation functions are equivalent and characterize the two-particle correlations. The experimental data clearly establish the existence of negative two-particle correlations, among both the protons and the antiprotons. The data for protons at √ s NN 20 GeV are adequately described when both the baryon conservation and excluded volume effects are taken into account. The baryon conservation exerts a stronger influence on Ĉ2 / Ĉ1 than the excluded volume although both effects are necessary to obtain a fair agreement with the data at √ s NN 20 GeV.

At lower energies the suppression of Ĉ2 / Ĉ1 is overestimated by the model, especially at √ s NN = 7.7 GeV. This can be due to different mechanisms. For instance, we have neglected the effect of volume fluctuations which would increase Ĉ2 / Ĉ1 [70]. The STAR Collaboration has applied the centrality bin width correction [71] to minimize the effects of volume fluctuations in the data, which, however, does not remove volume fluctuations completely [33]. To leading order, the volume fluctuations lead to an additive correction to Ĉ2 / Ĉ1 [32], namely, Ĉvol. fluc.

Ĉvol. fluc.

Here v 2 is a normalized variance of volume fluctuations. The 7.7-GeV STAR data point could be described with volume fluctuations for v 2 ≈ 0.002, however one would require v 2 < 0.001 to not spoil the agreement at higher collision energies. Thus, the volume fluctuations could only explain the deviations from experimental data if their effect is considerably larger at √ s NN = 7.7 GeV than at higher energies. Qualitatively, such a behavior has been indicated before [72] and it remains to be seen if it can provide a quantitative resolution.

A potentially more intriguing explanation for the disagreement at √ s NN = 7.7 GeV is the presence of sizable attractive interactions among protons. This effect is not included in our model and would lead to an increase in Ĉ2 / Ĉ1 . If this is the case, the data would indicate a transition from repulsive to attractive proton interaction regime as the collision energy is decreased below √ s NN 20 GeV. One possible mechanism for this would be approaching the QCD critical point. It should be noted however that approaching the QCD critical point would be also expected to generate multiparticle correlations [29], which has not yet been established by the data. At lower collision energies one can also expect sizable effects due to the nuclear liquid-gas transition [73][74][75].

It should be noted that at √ s NN = 7.7 GeV the proton cumulants are expected to also be affected by exact conservation of electric charge. We estimate this effect in Appendix C and show that this would lead to a further suppression of Ĉ2 / Ĉ1 by about 20%. This would increase the disagreement with the experimental data.

The trends in the data for antiprotons are reproduced by the model, although, in contrast to the protons, the data point to considerably stronger anticorrelation among the antiprotons than predicted by the model. This difference between the protons and antiprotons may be related to their possibly different production mechanisms. While the measured protons consist of both the stopped and the newly produced protons, all the measured antiprotons are the newly produced particles only. If the newly produced particles are affected by a different mechanism compared to the stopped protons, for instance by local rather than global baryon conservation, this may lead to a difference in the behavior of two-particle correlations of protons and antiprotons. A more involved modeling, however, is required to shed light on this possibility. It should also be noted that the agreement of our present model with the data is already much better than that of the UrQMD model calculations [76,77] shown in the STAR paper [10]. The results based on the noncritical baseline of Ref. [27] are shown in Fig. 5 by the blue points. They agree well with our ideal HRG model based results and thus show a similar quantitative disagreement with the STAR data for the antiprotons.

The higher-order correlation functions Ĉ3 / Ĉ1 and Ĉ4 / Ĉ1 show only small deviations from zero in our model. This is consistent with the fact that our model has no critical point and the associated critical fluctuation dynamics which would be expected to generate strong multiparticle correlations among protons in the vicinity of the critical point [29]. The result is also consistent with an earlier observation made in Ref. [70] that baryon conservation, which is included in our model, has a modest effect on three-and four-proton correlation functions. Our results for Ĉ3 / Ĉ1 and Ĉ4 / Ĉ1 are consistent with the STAR data, if the experimental error bars are to be taken seriously. This also means that the statistically significant deviations of the third-and fourth-order cumulants from the Skellam baseline are indeed driven by the two-particle correlations, i.e., by the contributions of Ĉ2 to κ 3 and κ 4 rather than by genuine multiparticle correlations. Thus, these data presently show no indication for the existence of the QCD critical point in the studied collision energy regime. As the experimental uncertainties at √ s NN 14.5 GeV are sizable, however, the present data also do not rule out a possible presence of notable multiparticle correlations among protons in that collision energy regime. The high statistics data coming from RHIC-BES-II will thus be able to shed light on this issue. 5. Collision energy dependence of scaled (anti)proton cumulants and factorial cumulants (correlation functions) in 0-5% Au-Au collisions up to fourth order. The solid lines depict calculations that incorporate both the baryon conservation and excluded volume effects (EV-HRG model) while the dashed lines correspond to baryon conservation only (ideal HRG model). The red squares and gray triangles correspond to the experimental data of the STAR Collaboration [10] for protons and antiprotons, respectively. The blue circles correspond to the canonical ensemble ideal HRG model calculation based on (anti)proton acceptance fractions from Ref. [27].

The cumulants and correlation functions have been measured by the STAR Collaboration as a function of acceptance in rapidity. Here we compare our model predictions for the acceptance dependence of cumulants with the STAR data. As neither the model nor the STAR data show conclusive notable deviations from zero for the higher-order normalized correlation functions Ĉ3 / Ĉ1 and Ĉ4 / Ĉ1 , we focus the analysis of the acceptance dependence on the second normalized correlation function Ĉ2 / Ĉ1 .

The results for proton and antiproton number Ĉ2 / Ĉ1 as a function of the rapidity cut y max (i.e., |y| < y max ) are FIG. 6. Rapidity acceptance dependence of the second normalized factorial cumulant (correlation function) Ĉ2 / Ĉ1 of protons (blue lines) and antiprotons (black lines) calculated from hydrodynamics and compared to the experimental data of the STAR Collaboration [10]. The calculations incorporate both the excluded volume effect and global baryon conservation.

shown in Fig. 6. The magnitude of Ĉ2 / Ĉ1 increases with y max for both the protons and antiprotons, at all collision energies. This is the expected result which reflects that (i) the effect of baryon conservation becomes stronger when a larger fraction of particles is accepted [26] and (ii) the thermal smearing of local correlations diminishes for larger y max [29].

For protons, the STAR data at √ s NN = 19.6 GeV and above are described by the model quite well, including the slope of the y max dependence. For √ s NN = 14.5 GeV the data are described up to y max = 0.3, whereas at y max = 0.4 and 0.5 the model predictions are below the data, i.e., the slope in the data changes faster than in the model. For √ s NN = 7.7 GeV the model is below the data for all y max , the largest deviations being observed at the maximum measured y max = 0.5. Interestingly, the data at √ s NN = 7.7 and 14.5 GeV (as well as at √ s NN = 11.5 GeV [10] not shown here) show indications that the slope of the y max dependence of proton Ĉ2 / Ĉ1 may flip sign at y max > 0.5. Such a qualitative feature is not observed in our model, i.e., it cannot be attributed to baryon conservation or baryon repulsion.

The STAR data for antiprotons are described at the lowest two energies, √ s NN = 7.7 and 14.5 GeV, as well as at the top RHIC energy √ s NN = 200 GeV. At the intermediate energies the magnitude of the negative two-particle correlation among the antiprotons is underestimated by the model, for all values of y max . The data indicate a larger negative slope of antiproton Ĉ2 / Ĉ1 than predicted by the model at these energies.

Our calculations have been focused on 0-5% Au-Au collisions, as in that regime the assumptions of our formalism are most appropriate. Namely, the degree of thermal and chemical equilibration is expected to be the highest in most central collisions, the effect of volume fluctuations is expected to be smaller than in peripheral collisions, and the applicability conditions of SAM-2.0 are expected to be satisfied with the highest accuracy. Nevertheless, for the sake of completeness, we have also performed calculations of proton cumulants at different centralities within our event-averaged hydrodynamics framework. We find that all cumulant ratios stay essentially flat as a function of centrality at a given collision energy, i.e., our framework predicts essentially no centrality dependence for all volume-independent measures of event-by-event fluctuations. This observation agrees well with the measurements of the STAR Collaboration [10] for √ s NN 20 GeV and N part 100, while at lower energies and in peripheral collisions deviations from this picture are observed, indicating that a more involved modeling is warranted in those regimes.

Here we discuss an important issue which affects many theory-to-experiment comparisons, namely, the difference between baryon and proton number fluctuations. The experiment has direct access to the latter but it is notoriously difficult to measure all baryons. On the other hand, proton number is inaccessible in many (effective) QCD theories, for instance lattice QCD. Cumulants of net baryon number are computed instead, and often directly compared to net proton cumulants measured in the experiment.

In our model it is possible to compute both the proton and baryon number cumulants. This then allows us to establish to what extent the two correspond to each other for conditions realized in heavy-ion collisions at RHIC-BES. Figure 7 depicts the beam energy dependence of Sσ 3 /M, κσ 2 , and κ 6 /κ 2 of net protons (red lines) and net baryons (black lines) calculated within our formalism including the excluded volume and baryon conservation effects. It is seen that net proton and net baryon cumulants are quite different, with baryons generally exhibiting larger deviations from the Skellam baseline. In particular, the net baryon Sσ 3 /M disagrees with the experimental data on net protons at all collision energies, whereas the net proton calculation agrees much better. The error bars in the data for κσ 2 are still too large to make a clear distinction between the model predictions for net protons and net baryons. However, from the difference between red and black lines it is clear that the issue will persist once precision measurements of κσ 2 become available.

Qualitatively, Sσ 3 /M and κσ 2 exhibit similar trends as functions of collision energy when compared between protons and baryons. As seen in the bottom panel of Fig. 7, this is no longer the case for κ 6 /κ 2 , where the net proton hyperkurtosis monotonically increases with √ s NN while the net baryon hyperkurtosis exhibits a nonmonotonic dependence. Note that this nonmonotonicity here is unrelated to critical phenomena, which our model does not have, but is caused by an interplay V NV N N FIG. 7. Collision energy dependence of Sσ 3 /M, κσ 2 , and κ 6 /κ 2 of net protons (red lines) and net baryons (black lines) in 0-5% central Au-Au collisions. The calculations incorporate both the baryon conservation and excluded volume effects. The experimental data of the STAR Collaboration [9] are depicted by the red circles. between baryon repulsion and conservation which is sensitive to the collision energy.

Note that a method exists to reconstruct baryon number fluctuations from the measured proton number fluctuations [53,54]. This method assumes that the isospin is randomized at the final stage of the collision, for instance due to resonance formation and regeneration reactions during the hadronic phase. In that case baryon number cumulants can be obtained by performing a binomial unfolding of the measured proton cumulants. This essentially corresponds to an additional binomial efficiency correction with the efficiency probability q ≈ 1/2, as discussed in Refs. [53,54,78]. Experiment can do this procedure, however this requires very precise measurements of high-order cumulants, because the binomial unfolding increases the error in κ n by a factor of FIG. 8. Collision energy dependence of the second scaled factorial cumulant Ĉ2 / Ĉ1 of protons (red line) and baryons (black line) in 0-5% central Au-Au collisions. The calculations incorporate both the baryon conservation and excluded volume effects. The red squares depict the experimental data of the STAR Collaboration for protons [10] while the black circles correspond to baryons reconstructed from the proton data using the binomial distribution. order q -n ≈ 2 n . Thus applying the unfolding may not be useful using the presently available data on the skewness and kurtosis of (net-)proton distributions which have sizable error bars.

On the other hand, the second-order cumulants have already been measured fairly precisely, making the baryon unfolding procedure feasible to do. To illustrate this, we apply the unfolding to the STAR data on the ratio Ĉ2 / Ĉ1 of proton number factorial cumulants. Following Refs. [53,54], the baryon and proton factorial cumulants are related by ĈB n = Ĉ p n /q n . Thus, ĈB

We apply Eq. ( 25) with q = 1/2 to the STAR data to reconstruct two-particle correlations of the baryon number. The result is depicted in Fig. 8. The baryon number Ĉ2 / Ĉ1 constructed from the STAR data agrees well with our model calculation at √ s NN 20 GeV, similarly to the agreement for the proton number Ĉ2 / Ĉ1 . The results presented in this section indicate that there are essential quantitative differences between proton and baryon number cumulant ratios, and that the two may not be compared directly. Either baryon number cumulants should be reconstructed from the data on proton number cumulants via the method of Refs. [53,54], or the proton number cumulants, rather than baryon number cumulants, should be calculated in the theory.

In this paper we calculated the leading six (net-) (anti)proton cumulants in heavy-ion collisions at RHIC beam energy scan energies in the framework of relativistic viscous hydrodynamics. The cumulants have been calculated in the momentum acceptance where experimental measurements have been performed by the STAR Collaboration. For the first time, effects of exact global baryon conservation and repulsive interactions among baryons, modeled by excluded volume, have been incorporated simultaneously and in a dynamical description of heavy-ion collisions. The excluded volume parameter has been chosen such that deviations from the Skellam distribution in net baryon cumulants computed in lattice QCD at µ B = 0 are reproduced by the excluded volume HRG model used in our analysis. The model does not incorporate any critical point effects at finite baryon density; thus, our calculations correspond to the no-critical-point scenario.

We obtain good agreement with the experimental data of the STAR Collaboration on net-proton cumulant ratios Sσ 3 /M and κσ 2 at √ s NN 20 GeV. It is observed that the effect of baryon conservation has a stronger influence on the proton cumulants than excluded volume, although incorporating both effects simultaneously is required in order to reproduce the experimental data on Sσ 3 /M quantitatively. At lower collision energies, √ s NN 20 GeV, the data are underestimated by the full model. The model is in fair agreement with the data on net proton κσ 2 , although the available data have rather large error bars. Our model predicts that the net proton hyperkurtosis κ 6 /κ 2 in the STAR acceptance |y| < 0.5, 0.4 < p T < 2.0 GeV/c turns negative at √ s NN 40 GeV, mainly as a consequence of strong effect of baryon conservation.

We explored the behavior of cumulants and factorial cumulants of proton and antiproton distributions. It is observed that our model produces notable negative two-particle correlations among protons and antiprotons, but only mild three-and four-particle correlations, characterized by small values of the third and fourth scaled factorial cumulants, Ĉ3 / Ĉ1 and Ĉ4 / Ĉ1 . In this case the behavior of the high-order cumulants such as skewness and kurtosis is driven by the two-particle correlations rather than by multiparticle correlations that would have been expected near the critical point. The experimental data are consistent with small, if not vanishing, threeand four-particle proton correlations within errors, thus these data do not contain statistically significant evidence for the existence of the QCD critical point in the studied collision energy regime. Note, however, that the error bars on Ĉ3 / Ĉ1 and especially Ĉ4 / Ĉ1 are quite large, and it is possible that more precise measurements may find evidence for multiparticle correlations. The upcoming data from RHIC-BES-II will be essential to shed light on this possibility.

The experimental data for the second normalized correlation function Ĉ2 / Ĉ1 of protons are described well by the model at energies √ s NN = 19.6 GeV and above. At lower energies the model predictions are notably below the data. We discussed volume fluctuations and/or attractive interactions as a possible explanation for this discrepancy but further studies are required to shed more light on this issue. It has also been observed that negative two-particle correlations among antiprotons seen in STAR data are notably underestimated by the model at √ s NN = 11.5 and 19.6-62. 4 GeV, whereas at √ s NN = 7.7, 14.5, and 200 GeV the data are described

well. An explanation of this observation is presently an open question.

Compared to other quantitative theoretical predictions for the cumulants that are available in the literature, our model demonstrates a much better agreement with the STAR data than the UrQMD transport model calculations shown in Ref. [10]. Our calculations that incorporate baryon conservation but not the excluded volume repulsion (ideal HRG model) are consistent with the data-driven approach of Ref. [27] within the same ideal HRG model framework. We do observe, however, that the quantitative agreement with the STAR data on net proton Sσ 3 /M and the proton number normalized correlation function is obtained at √ s NN 20 GeV only when the baryon repulsion is incorporated in addition to baryon conservation. The presence of baryon repulsion is in line with the behavior of baryon number susceptibilities at µ B = 0 observed in lattice QCD.

Comparisons between (net-)proton and (net-)baryon cumulants revealed essential quantitative differences between the two in the RHIC-BES regime. The higher-order net baryon cumulants generally exhibit larger deviations than the net proton ones from the Skellam distribution baseline of an uncorrelated particle production. This is due to the fact that protons form a subset of all baryons, thus the strength of the measured correlations is diluted compared to the full baryon set. As a consequence, for instance, the calculated net baryon Sσ 3 /M underestimates significantly the measured net proton Sσ 3 /M, whereas the calculated net proton Sσ 3 /M agrees well with the data. It is thus essential that the same quantities are employed for theory-to-experiment comparisons; in particular, we find no justification for direct comparisons between the measured (net-)proton and calculated (net-)baryon cumulants that have often been performed in the literature.

One way to address the issue of the difference between proton and baryons cumulants is to unfold the baryon cumulants from the proton ones using the method of Kitazawa and Asakawa [53,54]. In the present paper we have demonstrated the feasibility of this procedure by unfolding the scaled factorial cumulant ĈB 2 / ĈB 1 of baryons from the corresponding scaled factorial cumulant Ĉ p 2 / Ĉ p 1 of the protons that was measured by the STAR Collaboration. The resulting data on baryon number ĈB 2 / ĈB 1 agree reasonably well with our model calculations of this quantity. The challenge in applying the unfolding to high-order cumulants lies in the fact that this procedure significantly increases the resulting experimental uncertainties. For this reason the method may not be viable for the presently available data from RHIC-BES phase I but should be viable to do using the more precise data coming from phase II. We hope that results presented here will serve as a motivation for this procedure to be done.

We have not incorporated any critical fluctuation dynamics associated with the QCD critical point in our paper. In that sense, our results can be viewed as a baseline that incorporates essential noncritical contributions to (net-)proton number cumulants stemming from baryon number conservation and repulsive baryonic interactions. Unambiguous signals of the QCD critical point in the beam energy scan regime, if there is one to be found, shall manifest themselves as deviations from our model calculations. We view the three-and four-particle correlation functions (factorial cumulants) of proton number to be particularly promising in that regard, perhaps more so than the ordinary cumulants. Our model, which has no critical fluctuations, predicts these scaled factorial cumulants to be small. On the other hand, the multiparticle correlations among protons are expected to be strong in the vicinity of the critical point [29] and may well be reflected in a sizable magnitude of the corresponding scaled factorial cumulants such as Ĉ3 / Ĉ1 and Ĉ4 / Ĉ1 . We note that the development of a quantitative hydrodynamics framework to describe critical fluctuations is in progress [79][80][81] which should eventually be able to provide more robust predictions of the critical point signals in both ordinary and factorial cumulants.

FIG. 10. Collision energy dependence of scaled (anti)proton cumulants κ 2 /κ 1 -1 (top) and κ 3 /κ 1 -1 (bottom) in 0-5% Au-Au collisions, evaluated at Cooper-Frye particlization using the ideal HRG model analytically (dashed lines) and via Monte Carlo sampling (bands).

(3) The event is rejected if the sampled yields violate the exact global baryon number conservation, i.e., if the sampled total net baryon number does not coincide with the expected baryon number computed in the first step. If the sampled yields satisfy the global conservation, we go to the next step.

(4) Momenta and coordinates of each hadron are sampled, independently from all other hadrons. To do that first we determine the hypersurface element from which the given hadron is sampled; this is done via the multinomial distribution where each volume element is weighted by its grand-canonical mean yield for the given hadron species. Then the momenta and coordinates of the hadron emitted from the chosen hypersurface element are sampled via the standard procedure.

(5) The chain of all strong, electromagnetic, and weak decays is performed.

We sample 100 000 events for each collision energy. Figure 10 depicts the collision energy dependence of the subtracted scaled second and third cumulants κ 2 /κ 1 -1 and κ 3 /κ 1 -1 for protons and antiprotons, evaluated in the STAR momentum acceptance. The analytic calculations (dashed lines) are the same that are shown in Fig. 5 of the main text. The Monte Carlo results are depicted by the bands. The Monte Carlo and analytic results are in good agreement with each other. This validates the analytic calculations and also indicates that the simplifying assumptions made in the analytic calculation, like neglecting the difference between the masses of nucleons and baryon resonances, have negligible influence on the cumulants. The Monte Carlo results also validate the accuracy of the generalized subensemble acceptance method of Ref. [40] used to correct the grand-canonical cumulants for baryon number conservation.

Note that although here we explicitly tested only the cumulants up to third order, we expect the analytic calculations to be accurate also for the high-order cumulants. This is due to the fact that both the baryon conservation and excluded volume generate only small multiparticle correlations (see Fig. 5), thus the high-order cumulants are mainly determined in this setup by two-particle correlations that we checked to be calculated accurately. The situation may change if one incorporates effects that generate strong multiparticle correlations like the QCD critical point. In that case the question of accuracy of the analytic calculations for high-order cumulants may have to be revisited.

In our analysis we have incorporated the effect of exact global baryon conservation but not of other conserved charges like electric charge and strangeness. While the effect of baryon conservation is expected to be the dominant one, the (anti)proton cumulants are also affected by other conserved charges [36,56]. This can be particularly relevant at lower collision energies where protons form a considerable fraction of the total electric charge. Thus here we evaluate the effect of multiple exactly conserved charges through Monte Carlo sampling within the ideal HRG model.

The sampling algorithm in Appendix B is adjusted in the following way: all events that do not satisfy exact conservation of all three conserved charges are rejected. The total net strangeness is constrained to be exactly zero in all events while the total electric charge is constrained to reproduce the charge-to-baryon ratio of Q/B = 0.4. 4 To speed up the sampling procedure, we employ the multistep method of Becattini and Ferroni [82]. Sampling with multiple conservation laws is more time consuming than with only baryon number conservation. We generate about 20 000 events for each collision energy and restrict the analysis to the second-order moments only.

The collision energy dependence of the subtracted scaled proton cumulant κ 2 /κ 1 -1 in 0-5% Au-Au collisions evaluated with the simultaneous conservation of baryon number, electric charge, and strangeness is shown in Fig. 11 by the blue band. The result is compared to the calculations with exact conservation of only baryon number (red band). The two calculations agree within statistical errors at √ s NN 20 GeV, suggesting that the restriction of the global conservation laws to only baryon number might be sufficient in that regime. At the lower collision energies the conservation of electric charge and strangeness leads to a notable further suppression of the two-particle correlation function of the protons. Thus, accounting for exact conservation of multiple conserved charges is important for analyses of fluctuations at √ s NN 20 GeV.

Our analytic calculations neglect rescatterings in the hadronic phase. The hadronic phase may affect the cumulants in a couple of different ways: (i) it modifies the p T spectrum of (anti)protons, thus the number of protons in the acceptance may change, and (ii) B B annihilations may decrease the numbers of protons and antiprotons. To evaluate the possible effect of the hadronic phase we use the Monte Carlo sampling of the ideal HRG from the previous subsection and, instead of performing the chain of decays, we run the output through the hadronic afterburner UrQMD [83,84]. is achieved by replacing step 5 the Monte Carlo sampling in the previous subsection by the following: all resonances which are not recognized by UrQMD are decayed until only hadrons and resonances recognized by UrQMD are left and then the hadronic phase is simulated by UrQMD via the hadronic afterburner toolkit from Ref. [85]. We evaluate the effect of the hadronic afterburner at √ s NN = 27 GeV by sampling around 200 000 events. The relevance of the hadronic phase is established by comparing the results with the afterburner (hydro + UrQMD) to the results without applying the afterburner (hydro + decays). Figure 12 depicts the rapidity acceptance dependence of the second scaled factorial cumulant Ĉ2 / Ĉ1 of protons and antiprotons in the STAR transverse momentum range, 0.4 < p T < 2.0 GeV/c. For antiprotons the difference between the two scenarios is within the statistical uncertainty. For protons the difference is also mild, with indications that Ĉ2 / Ĉ1 is slightly more suppressed when the hadronic phase evolution is included. This appears to be due to a larger fraction of protons ending up in the STAR acceptance and hence a larger effect of the baryon number conservation. The hadronic phase leaves the y max dependence of the scaled cumulants essentially unchanged. Thus, incorporating the time-consuming hadronic afterburner would seem to only be necessary for very precise studies of cumulant ratios, at least for √ s NN = 27 GeV.