Thermodynamic approach to proton number fluctuations in baryon-rich heavy-ion matter created at moderate collision energies

We develop a framework to relate proton number cumulants measured in heavy-ion collisions within a momentum space acceptance to the susceptibilities of baryon number, assuming that particles are emitted from a fireball with uniform distribution of temperature and baryochemical potential, superimposed on a hydrodynamic flow velocity profile. The rapidity acceptance dependence of proton cumulants measured by the HADES Collaboration in $\sqrt{s_{\rm NN}} = 2.4$ GeV Au-Au appears to be consistent with thermal emission of nucleons from a grand-canonical heat bath, with the extracted baryon number susceptibilities exhibiting an hierarchy $\chi_4^B \gg -\chi_3^B \gg \chi_2^B \gg \chi_1^B$. Naively, this could indicate large non-Gaussian fluctuations that might point to the presence of the QCD critical point close to the chemical freeze-out at $T \sim 70$ MeV and $\mu_B \sim 850-900$ MeV. However, the description of the experimental data at large rapidity acceptances becomes challenging once the effect of exact baryon number conservation is incorporated, suggesting that more theoretical and experimental studies are needed to reach a firm conclusion.

Introduction. The search for the QCD critical point in heavy-ion collisions is centered on the measurements of event-by-event fluctuations of (net) proton number [1,2]. The corresponding measurements have been performed by a number of experiments, including the ALICE experiment at the LHC [3], the NA61/SHINE experiment at SPS [4], and the STAR experiment at RHIC [5,6]. At √ s NN 20 GeV these measurements are consistent with being driven by noncritical physics such as baryon number conservation and repulsive interactions [7,8]. At lower energies, the data indicate an excess of two-proton correlations over a non-critical reference, with deviations increasing as the collision energy is decreased [8]. Recent data at √ s NN = 3 GeV from the STAR fixed target program [9] are consistent with this trend. One potential source of these excess correlations among protons could be the presence of the critical point.
Recently, the HADES Collaboration has presented measurements of the four leading proton number cumulants and integrated correlation functions [10] measured in Au-Au collisions at √ s NN = 2.4 GeV. The data, which are corrected for volume fluctuations, show strong dependence on the phase-space acceptance, approaching the Poisson limit for narrow rapidity bins and indicating the presence of long-range multi-particle correlations in rapidity as the acceptance is expanded. In this work we analyze the HADES data in the context of thermal emission of particles from an interacting fireball. We show that, in this scenario, the measurements in momentum space can be unfolded to extract the grand-canonical baryon number susceptibilities of QCD matter probed by the experiment and provide information on the phase structure of baryon-rich QCD matter. In particular, we discuss the possible presence of the QCD critical point in the baryon-rich region of QCD matter as a driving mechanism behind the behavior of experimental data as well as effect of baryon number conservation.
Cumulants and correlation functions of protons emitted from an expanding fireball. The basis for describing the emission of particles from an expanding fireball is the Cooper-Frye formula [11], which determines their invariant momentum spectrum Here f ∝ exp [−u ν (x)p ν /T (x)] is the phase-space distribution of particles in thermodynamic equilibrium, p µ is their four-momentum 1 , and dσ µ (x) is the invariant element of the freeze-out hypersurface from which the particles are emitted. In Ref. [8] the Cooper-Frye procedure was generalized to describe not only the first moments of particle distributions in a given momentum acceptance, but also the cumulants of their multiplicity distribution that are connected to the grand-canonical (baryon number) susceptibilities. The generalized formalism was used to calculate the proton number cumulants and correlation functions in Au-Au collisions at √ s NN = 7.7 − 200 GeV within numerical relativistic hydrodynamics simulations, as appropriate to the beam energy scan program at RHIC. Here we extend these considerations to lower collision energies, with specific appli-cations to √ s NN = 2.4 GeV Au-Au collisions studied by the HADES experiment at GSI-SIS. Compared to higher collision energies, several modifications and simplifications can be made. First, at √ s NN = 2.4 GeV one can neglect the production of antiparticles, as well as strange particles. Second, instead of using the numerical output from hydro simulations, it is possible to employ a single freeze-out scenario parameterized by an appropriate blast-wave-like hypersurface. Indeed, as recently shown in Ref. [12], the p T spectra of light flavored hadrons in central Au-Au collisions at HADES can be reasonably well described through a spherically-symmetric Siemens-Rasmussen blast-wave model with a Hubble-like radial flow. In this case the freeze-out takes place suddenly at a constant Cartesian time t = const and within a spherical spatial volume of radius R. The hypersurface element reads dσ µ = (1, 0, 0, 0) r 2 sin θ dθ dφ dr, where 0 < θ < π, 0 < φ < 2π, and 0 < r < R. The hydrodynamic flow is spherically symmetric and Hubble-like, u µ = γ(r) (1, v(r)e r ) with γ(r) = (1 − v 2 (r)) −1/2 and v(r) = tanh(Hr).
An important simplification of the single freeze-out scenario is that all elements of the freeze-out hypersurface are characterized by constant values of the temperature and baryochemical potential. This fact is used here to show how the measured fluctuations of baryon (or proton) number are related to the values of baryon number susceptibilities χ B n at the freeze-out point in the QCD phase diagram. Indeed, let us first start with the grand canonical limit, where the emission of baryons from each hypersurface element dσ(x) proceeds independently. The cumulants dκ B n (x) of baryon distribution emitted from dσ(x) is determined by the susceptibilities χ B n , namely Here dV eff (x) = dσ µ (x)u µ (x), thus dκ B n (x) = T 3 χ B n γ(r) r 2 sin θ dθ dφ dr .
The cumulants of the baryon number distribution in the full phase-space in the grand-canonical limit are obtained by integrating Eq. (3) over the freeze-out hypersurface, The experiment measures particles in a restricted part of the full phasespace, however. Let us denote the momentum acceptance by ∆p acc . Every nucleon emitted from dσ(x) ends up in ∆p acc with a certain probability p acc (x; ∆p acc ) which can be evaluated through the Cooper-Frye formula: The momentum acceptance can thus be modeled by convoluting dκ B n (x) with the binomial distribution [13], yielding (see Ref. [8] for details): Here , and B n,l are partial Bell polynomials. As a result, the total cumulants of accepted baryons κ B n (∆p acc ) = x dκ B n (x; ∆p acc ) can be expressed as linear combinations of χ B n : Note that these calculations automatically incorporate the effect of thermal broadening [14], through the integration over the Maxwell-Boltzmann momentum distribution f in Eq. (4). The experimentally measured protons constitute a fraction of all baryons. Based on the isospin content of the colliding nuclei, one can estimate that roughly q iso ≈ 0.4 of all baryons are protons. Furthermore, at low collision energies, a significant fraction q nucl of protons becomes bound into light nuclei, with q nucl ≈ 0.375 based on HADES measurements of proton and light nuclei production [12]. Assuming that this process happens sometime after the chemical freeze-out, for instance via coalescence, this implies that a fraction q = q iso (1 − q nucl ) ≈ 0.25 of all accepted baryons is ultimately measured in the experiment. To leading order, this effect can be modeled by convoluting κ B n (∆p acc ) with a binomial distribution with acceptance q. As a result, the cumulants κ p n (∆p acc ) of measured protons can be expressed as linear combinations of the grand-canonical baryon susceptibilities Here γ(t, q) = ln(1 − q + e t q) and γ . Equation (8) shows that cumulants κ p n (∆p acc ) measured in a particular acceptance ∆p acc are expressed as linear combinations of baryon number susceptibilities χ B l .
As the effect of exact baryon conservation has not been taken into account, this approach should work only when this effect can be neglected, i.e. for a sufficiently small acceptance ∆y acc ∆y beam [15]. Below, we use this fact to extract χ B l from the experimental data of the HADES Collaboration.
Analysis of the HADES data. We analyze the four leading proton cumulants in 0-5% Au-Au collisions at √ s NN = 2.4 GeV measured by the HADES Collaboration in Ref. [10]. To fix the parameters of the freeze-out hypersurface we utilize the results of Refs. [12,16,17], where the chemical and kinetic freeze-out conditions were studied. As shown in [16], a consistent description of hadron yields is obtained for chemical freeze-out temperature of T ≈ 70 MeV. In [12,17] it was shown that the momentum spectra can be reasonably well described within a single freeze-out scenario. We take the following parameters based on Ref. [17]: T = 70 MeV, R = 6.1 fm, The proton cumulants are evaluated in transverse momentum acceptance 0.4 < p T < 1.6 GeV/c around midrapidity, |y| < y cut , for various values of y cut consistent with the experimental acceptance [10]. To fix the values of χ B n we adopt the following strategy: rather than using model assumptions regarding the QCD equation of state at the point of freeze-out, we instead extract the values of χ B n that match the experimental data via Eq. (8). The values ofα n,l (∆p acc ) in Eq. (8) are calculated via Eqs. (9), (7), and (4) where ∆p acc = {0.4 < p T < 1.6 GeV/c, |y| < y cut }. We perform this procedure for values of y cut in a range 0.1 − 0.2, using central values of the experimentally measured cumulants. We do not use the data with y cut > 0.2 because at large acceptance the effect of exact baryon conservation becomes dominant (see the next section below). One obtains the following grand-canonical baryon susceptibility ratios: The extracted values are qualitatively consistent between the two data sets corresponding to different y cut and can be combined as The theoretical uncertainty here comes from the differences in the values of the extracted susceptibilities for y cut = 0.1 and 0.2. Then we use the extracted values of χ B n to evaluate the proton cumulants as a function of y cut . Figure 1 depicts the behavior of the corresponding proton cumulant ratios κ p 2 /κ p 1 , κ p 3 /κ p 2 , and κ p 4 /κ p 2 in the model using the central Note that these data were not used in the procedure to extract the susceptibilities. Remarkably, the model also reproduces semi-quantitatively the experimental data for y cut > 0.2, even though it does not incorporate the exact baryon conservation which is relevant at large acceptances.
Due to the limitations and simplifications of our approach, the extracted baryon number susceptibilities should be taken with a grain of salt from a point of view quantitative estimation of the equilibrium grandcanonical susceptibilities. Qualitatively, however, it is clear that the results indicate the presence of large non-Gaussian fluctuations of the baryon number, characterized by the increasing magnitude of the high-order cumulants. In fact, the extracted baryon susceptibilities exhibit the following hierarchy One known physical mechanism that can generate large non-Gaussian fluctuations of the baryon number in the grand-canonical ensemble is the QCD critical point [18,19]. If this is the correct mechanism, it would indicate that the freeze-out point of √ s NN = 2.4 GeV Au-Au collisions is located in the critical scaling region in close proximity to the critical point, where the large non-Gaussian fluctuations are predicted to emerge. 2 Assuming that the freeze-out indeed takes place close to the critical point, one can use this to estimate its location in the QCD phase diagram. The values of T and µ B at the chemical freeze-out have been estimated by fitting the measured hadron yields in the framework of statistical hadronization model (SHM) in Refs. [12,16]. Although these fits show two distinct minima that describe the data similarly well, it was shown in Ref. [16] using transport model simulations that only one of the two minima corresponds to a reasonable chemical freeze-out picture, and it is located at As discussed in [16], the extracted values are fairly robust with respect to the theoretical uncertainties in the SHM, in particular, to the treatment of nucleon interactions. Note also, that the extract chemical freeze-out temperature of 70 MeV is the same as the T kin = 70 MeV value describing the momentum distributions of protons and pions in the fireball model that we used, consistent with the notion of approximately simultaneous chemical and kinetic freeze-out. If the interpretation of the experimental data presented here is correct, the thermal freeze-out at HADES is located close to the QCD critical point and thus Eq. (14) serves as an estimate for the possible location of the QCD critical point. The estimate, which corresponds to µ B /T ∼ 11, puts the critical point beyond the reach of current lattice QCD methods corresponding to µ B /T 2 − 3.5 [20,21] but in the ballpark of predictions of some effective QCD models like those based on holography [22,23].
The extracted susceptibilities can also be used to improve the analysis presented in Ref. [24], where it was shown that the cumulants can be used to estimate the (isothermal) speed of sound, as well as its logarithmic derivative with respect to baryon density. To do that one needs to determine the true grand-canonical cumulants of baryon number whereas in Ref. [24] direct measurements of the momentum space proton number cumulants were used as a proxy instead. Using Eqs. (10) and (14), we obtain the following for the speed of sound: which is a small value consistent with proximity to the critical point. The logarithmic derivative reads Our results agree with the simplified estimates of Ref. [24] qualitatively, while quantitatively they indicate a smaller c 2 T and a larger d(ln c 2 T )/d(ln n B ) at the HADES collision energy than originally obtained in [24].

Effect of exact baryon number conservation.
Our analysis has so far been based on the picture of thermal particle emission from a grand-canonical heat bath. This does neglect one important effect, namely the exact conservation of the total baryon number, which does not change throughout the heavy-ion collision [25]. In addition, other exact conservation laws like that of electric charge or total energy-momentum can influence the results as well. Correcting the results for conservation laws is challenging, especially at lower collision energies, although theoretical advances have been made recently [15,26]. Here, we use two different methods to incorporate baryon conservation in order to estimate the systematic error due to this effect. First, we consider a simplified correction based on Refs. [15,25] that originally has been derived within the ideal gas model in the canonical ensemble and discussed in the original HADES paper [10]. The proton cumulant ratios subject to baryon conservation are computed in this simplified framework asκ p 2 We denote by tilde the cumulants that are subject to the effect of exact baryon conservation. Here β = 1 − α and is the fraction of measured protons from all baryons. This fraction varies from α ≈ 0.03 for y cut = 0.1 to α ≈ 0.13 for y cut = 0.5. The results for proton cumulant ratios corrected for baryon conservation in this way are shown by the dash-dotted lines in Fig. 1. These results show the same qualitative y cut dependence as the uncorrected results, although they do worsen the quantitative agreement with the experimental data, especially for the scaled variance. It is possible that the agreement could be recovered by readjusting the values of the susceptibilities.
However, it should be noted that the presented correction of proton cumulants for baryon conservation cannot be expected to be accurate. In particular, using the fraction of accepted protons, i.e. excluding other baryons such as neutrons and light nuclei, for calculating α is only justified in the limit of free gas, where the grand-canonical susceptibilities χ B n correspond to Poisson statistics, i.e. if the susceptibilities of all orders are equal. The values of the extracted susceptibilities in Eqs. (10)(11)(12) clearly violate this assumption. In this case, the corrections in Eqs. (17)(18)(19) should only be applied to cumulants of the conserved baryon number rather than that of non-conserved proton number, as discussed in Ref. [15] in the framework of the subensemble acceptance method (SAM-1.0). Therefore, one should apply the correction to the cumulants of accepted baryons in Eq. (6) (rather than protons) using the following expressions from the SAM-1.0: The cumulantsκ p n of accepted protons are then obtained fromκ B n by convoluting the latter with the binomial distribution with acceptance q, as before. The results of this procedure are shown in Fig. 1 by the dotted lines. In addition to SAM-1.0, we also apply an updated framework called SAM-2.0, which has recently been developed in Ref. [27] and used to describe experimental measurements of proton cumulants in Au-Au collisions for √ s NN 20 GeV. The advantage of SAM-2.0 over SAM-1.0 is that its formulation is not restricted to a uniform coordinate space but allows one to perform the corrections in momentum space, which need not be homogeneous. The method involves calculating the joint cumulants proton/baryon distribution both inside and outside the acceptance, and then using a mapping function to evaluate the cumulants of accepted protons subject to baryon conservation. The details of the mapping function and the evaluation procedure for an emission from a Cooper-Frye hypersurface can be found in Refs. [27] and [8], respectively. The results obtained within SAM-2.0 are shown in Fig. 1 by the dashed lines. The results obtained using SAM-1.0 and SAM-2.0 are similar. It is evident that the effect of baryon conservation computed by either of the two methods is strong for y cut > 0.2, where the model disagrees with the experimental data. We checked within SAM-2.0 that the agreement cannot be recovered by refitting the values of χ B n : one can fit χ B n to describe the data for a particular value of y cut > 0.2, however, the model then fails to describe the data for other rapidity cuts. The SAM-2.0 fits for y cut = 0.1 and y cut = 0.2, however, reproduce the hierarchy (13) given by Eqs. (10)- (12), indicating that corrections due to baryon number conservation are indeed subleading for y cut 0.2.
The results indicate that a more involved modeling is needed to describe the cumulants at y cut > 0.2. On one hand, as discussed in Ref. [27], the applicability conditions of the SAM at low collision energies may not be satisfied well, and thus a more involved method for estimating the effect baryon conservation may be needed. Furthermore, the additional effect of exact conservation of electric charge can be relevant as well [26]. Thus, more involved modeling may be warranted than the fireball model considered here. One possibility is to use molecular dynamics with a critical point as recently explored in Ref. [28]. On the other hand, it is also possible that the effects of spectator-participant interactions become sizable at larger y cut , in particular volume fluctuations, thus invalidating the thermodynamic approach for proton cumulants at large rapidities. It is also evident from the results that the baryon conservation effect is moderate and can likely be neglected for y cut < 0.2. The data analysis in that regime can thus be performed in the grand-canonical ensemble, without the theoretically challenging modeling of exact baryon conservation.
Summary. We have shown how experimentally measurable proton number cumulants in heavy-ion collisions at high µ B can be related to baryon number susceptibilities in the picture of thermal emission of particles from approximately uniform fireball superimposed on hydrodynamic flow velocity profile. We applied the formalism to recent experimental data of the HADES Collaboration and showed that the acceptance dependence of proton cumulants at |y cut | 0.2 in Au-Au collisions at