Net-proton number cumulant ratios as function of beam energy from an expanding nonequilibrium chiral fluid

The beam energy scan program at RHIC provides data on net-proton number fluctuations with the goal to detect the QCD critical end point and first-order phase transition. Interpreting these experimental signals requires a vital understanding of the interplay of critical phenomena and the nonequilibrium dynamics of the rapidly expanding fireball. We study these aspects with a fluid dynamic expansion coupled to the explicit propagation of the chiral order parameter sigma via a Langevin equation. Assuming a sigma-proton coupling through an effective proton mass, we relate cumulants of the order parameter and the net-proton number at freeze-out and obtain observable cumulant ratios as a function of beam energy. We emphasize the role of the nonequilibrium first-order phase transition where a mixed phase with gradual freeze-out can significantly alter the cumulants. We find that the presence of a critical end point is clearly visible in the cumulant ratios for a relatively wide range of center-of-mass energies.


I. INTRODUCTION
Experiments at SPS and STAR have found that nuclear matter at high temperatures undergoes a transition from a hadronic phase to a quark-gluon plasma [1,2], a state of deconfinement and chiral symmetry restoration. While this transition is a continuous crossover for zero and small baryochemical potential µ B [3][4][5], a critical endpoint (CEP) and first-order phase transition (FOPT) are expected at large µ B . Here, a variety of effective models of quantum chromodynamics (QCD) [6][7][8] as well as functional techniques [9,10] yield widely different results regarding existence and location of CEP and FOPT in the space of T and µ B .
Besides these theoretical efforts, considerable workforce is invested in ongoing experimental programs aiming at understanding the QCD phase diagram, e.g. the beam energy scan at STAR [11], NA49/61 [12,13], HADES [14], or the upcoming facilities NICA [15] and FAIR [16]. Important observables in this context are cumulants of conserved quantities, namely baryon number, strangeness, and electric charge. As shown in lattice QCD [17,18], various studies of effective models [19][20][21], and functional techniques [22], these cumulants diverge at the CEP and behave characteristically in the critical region around the CEP. Even though these calculations are based on equilibrium thermodynamics, it is widely believed that a measurement of cumulants or cumulant ratios in heavy-ion collisions should reveal nonmonotonic behavior as function of center-of-mass energy in the presence of a CEP. Here, however, a thorough understanding of the nonequilibrium dynamics and a design of proper experimental methods are crucial. For an overview of CEP physics at STAR, see [23]. Since cumulants of various order are directly proportional to certain powers of the correlation length, they diverge in an infinitely large equilibrated medium close to the CEP. In a heavy-ion collision, their growth is limited not only by the finite system size but also by critical slowing down [24]. The importance of a proper dynamical description of the critical dynamics has been emphasized by a variety of publications in recent years [25][26][27][28][29][30][31][32][33][34]. The direct impact on experimental observables, e.g. the net-proton number fluctuations, has been argued and demonstrated in [35][36][37][38]. Equally important than understanding the complicated dynamics near the CEP is understanding and modeling of evolutions passing through a FOPT where spinodal decomposition enforces density inhomogeneities within single events [39][40][41][42][43][44][45] and leads to an increased production of entropy [46,47] which could be observed, e.g., via an enhancement of the pion-to-proton ratio.
In this paper, we apply the nonequilibrium chiral fluid dynamics model [48] to a longitudinal Bjorken expansion along the beam axis [47]. This model describes the dynamics of the sigma field as the chiral order parameter with a Langevin equation interacting with a locally thermalized expanding quark fluid. Field and fluid exchange energy-momentum via a source term. We extract cumulants of the sigma field on an event-by-event basis along a parametrized freeze-out curve. As shown in [36], we relate these cumulants to net-proton number fluctuations by assuming a superposition of standard Poisson and critical fluctuations. We correct our results by the effect of volume fluctuations as detailed in [49]. In the present work, we aim at revealing possible signatures in the net-proton number cumulants that would confirm or rule out a CEP and FOPT in stronly-interacting matter.
After a description of the model in Section II, we report our results on various net-proton number cumulant ratios as function of beam energy in comparison to data from STAR and HADES in Section III, and finally conclude with a summary in Section IV.

II. MODEL DESCRIPTION
The model is based on the Lagrangian of the widely studied quark-meson model [6,7,50] with a CEP at (T CEP , µ CEP ) = (100, 200) MeV. Although arguably simple, the model provides a description of a chiral phase diagram with the generic features of a crossover, CEP, and FOPT. Currently, there is no agreement within the theoretical physics community on where to find the CEP in the phase diagram [51], yet recent QCD-based calculations from functional renormalization group techniques [10,52] suggest a location similar to the one in the quarkmeson model. We use the fluctuation of the sigma field as the critical mode of the QCD CEP, characterized by a vanishing sigma mass which has also provided one of the fundamental motivations for the beam energy scan program at RHIC [36,53].
Here and in the following, µ denotes the quark chemical potential, thus µ = µ B /3. The Lagrangian for light quarks q = (u, d) and the chiral order parameter σ with potential U reads with standard parameters f π = 93 MeV, m π = 138 MeV and U 0 such that U (σ) = 0 in the ground state. The value of the pion fields has already been set to its vacuum expectation value of zero. The quark-sigma coupling constant g is fixed requiring that gσ equals the nucleon mass of 940 MeV in vacuum. The grand potential Ω in mean-field approximation resembles that of a Fermi gas of quarks with energies E = p 2 + g 2 σ 2 and is evaluated as Here, N f = 2, N c = 3 denote the number of light quark flavors and the number of colors.

A. Equations of motion
We evolve the zero mode or volume-averaged sigma field defined as σ(τ ) = 1 V d 3 xσ(τ, x) using a Langevin equation of motion, neglecting spatial fluctuations. We describe the expanding fluid using a Bjorken model, and consequently use proper time τ rather than coordinate time t, starting from an initial thermalization at τ 0 = 1 fm. Consequently, the dots in Eq. (5) denote derivatives with respect to τ . For our case of purely longitudinal hydrodynamic flow, we set D = 1 in the Hubble term. The full and proper nonequilibrium dynamics of sigma is encoded in the dissipation coefficient η and the stochastic noise ξ which are related by a dissipation-fluctuation relation [48], Here, ξ is assumed Gaussian and white, i.e. it is not correlated over time. The coefficient η includes effects from various processes: • Mesonic interactions, i.e. σσ ↔ σσ (scattering of a condensed sigma with a thermal sigma), and σ ↔ ππ (two-pion decay) [54], described by a phenomenological damping coefficient of η = 2.2/fm [55] wherever kinematically allowed.
• Meson-quark interactions, σ ↔ qq, leading to a Tand µ-dependent coefficient [48], We assume an ideal fluid of quarks and antiquarks described by the energy-momentum tensor T µν q = (e + p)u µ u ν − pg µν . Due to energy-momentum conservation the divergence of the total energy-momentum tensor T µν q +T µν σ vanishes, which leads to the evolution equation for the energy density, while the net-baryon density simply followṡ In Eqs. (5) and (8), the pressure is given by p = −Ω qq . It is an explicit function of σ which during the evolution is not fixed to its equilibrium value. The fireball volume which appears in Eq. (6) and also later in the description of the freeze-out cumulants, is given by V = πR 2 τ , with a gold nucleus radius of R = 7.3 fm, assuming central Au+Au collisions.
B. Freeze-out and mapping to beam energies We investigate higher order cumulant ratios of the netproton number along a freeze-out curve which has been obtained from thermal model fits to experimental data from SIS, AGS, SPS, and RHIC for a wide range of beam energies from √ s = 2.24 to 200 AGeV [56]. The parametrization reads with a = 0.166 GeV, b = 0.139 GeV −1 and c = 0.053 GeV −3 . Since the phase boundary of the quarkmeson model would lie below the thus obtained freezeout line, we scale T and µ in this parametrization with a common factor T crossover (µ = 0)/T freeze (µ = 0) yielding a freeze-out curve which coincides with the crossover line at µ = 0 and consistently lies below the crossover and phase transition for µ > 0. Here, T crossover (µ = 0) = 145 MeV is determined by a maximum of the quark number susceptibility.
We define initial conditions for the evolution by adopting initial values T i and µ i from previous studies [47]. During the evolution of the fluid according to Eqs. (5), (8), (9), the trajectory in T -µ space will eventually hit the freeze-out curve. This initial hit point is then used to map the evolution to a corresponding beam energy via: with parameters d = 1.308 GeV and e = 0.273 GeV −1 determined in accordance with the freeze-out curve above [56]. The baryochemical potential in Eq. (11) is obtained from averaging over events with the same initial condition. Note that the thus obtained beam energies are to be understood as guidelines to put our results into the context of experimentally obtained data from STAR and HADES. Even though in our present model, the CEP is passed for evolutions with an initial √ s ∼ 5 GeV, this does not mean that we necessarily expect the physical CEP to be found there. However, it puts us in a position to estimate a CEP's impact on measurable observables in a fully dynamical nonequilibrium setup if it indeed exists in this low-energy range.

C. Sigma and net-proton number cumulants
To relate the fluctuations in the chiral order parameter σ to fluctuations or cumulants of the net-proton number, we follow the strategy outlined in [36]. Consider an infinitesimal change of the chiral field, δσ, leading to a change of the effective proton mass by δm = gδσ. Assuming a sigma-proton coupling gσpp, we may write fluctuations of the momentum space distribution function for protons, f k , as The first term δf 0 k is the purely statistical fluctuation, and in the second term, n FD denotes the Fermi-Dirac distribution for a particle of a given mass m. The fluctuation of the net-proton multiplicity N = V d d 3 k (2π) 3 f k is then given by where d = 2 is the spin degeneracy factor. The first term δN 0 can be assumed Poisson distributed. Consequently, all of its cumulants are equal to the expectation value N . In leading order and assuming no correlations between δN 0 and δσ, we can express cumulants of order n as In this notation, σ V = d 3 xσ = σV as we neglect spatial fluctuations, and · c is the respective cumulant, which is equal to the expectation value for n = 1 and to the corresponding central moment for n = 2, 3. For n = 4, we have and similarly for δN n c . In the next section III, we will use the shorter notation C n for the net-proton number cumulants δN n c which has also been commonly used in experimental studies of recent years.

III. RESEARCH PROCEDURE AND RESULTS
We initialize the fluid at a set of fixed values (T i ,µ i ) and define initial sigma field, energy density, quark number density, and pressure as the corresponding equilibrium values. The pairs of initial values are hereby adopted form earlier works [47,57] and will be matched to centerof-mass energies using the freeze-out condition, Eq. (11). The coupled system evolves according to the equations of motion until the freeze-out curve is hit. Notably, for expansions at high baryochemical potential, the freeze-out curve can be hit more than once due to the nonequilibrium evolution of the expanding plasma. This effects occurs due to the sudden release in latent heat that drives the system back into the chirally symmetric phase, visible in a back-bending of the trajectories in T and µ [42,47]. In contrast to that, an equilibrated hydrodynamic system with constant S/A would pass along the phase boundary for a finite amount of time [58]. To take into account this effect of a possible mixed phase in a heavy-ion collision and its impact on the observed cumulants, we evolve the system for these cases until a second crossing of the freeze-out curve and subsequently calculate cumulants at both crossing or hit points. In the following figures and text, we consequently refer to the "first hit" as the cumulants evaluated at the first crossing of the freeze-out line and the "second hit" as those evaluated at the second crossing, where applicable. This is relevant for evolutions near the CEP in the phase diagram or crossing the FOPT line [57]. Ultimately, it will provide us with a range of possible cumulants for the respective energies. The characteristic back-bending after passing the phase boundary becomes most pronounced for evolutions passing the FOPT due to significant energy dissipation pushing the system back into the chirally restored phase. Physically, this process could manifest in a gradual freeze-out where the medium undergoes droplet formation [39,59]. Another suggested signal for such a delayed transition process is an enhancement in the dilepton production [60].
We simulate 10 7 events and calculate event-by-event fluctuations in terms of cumulants of σ V . Since these are subject to significant fluctuations of the freeze-out volume, we include the corresponding corrections as derived in [49]. This is most relevant for evolutions at the lowest energies with large variations in the time at which different events from the same initial condition hit the freeze-out curve [57]. Eq. (14) allows us to determine the net-proton number cumulants. The such obtained values are compared with a Poisson baseline. For the net-proton number (p −p), the cumulants assuming Poisson distributed proton and antiproton numbers are calculated by C n,p−p = C n,p + (−1) n C n,p , where C n,p , C n,p are equal to the expectation values of the Poisson distribution for all orders n, cf. [11]. We furthermore provide comparison to the equilibrium netbaryon number susceptibilities which are obtained as Figure 1 (top) shows the ratio C 2 /C 1 of the net-proton number compared to results from STAR for energies √ s NN ≥ 7.7 GeV [11] and HADES for √ s NN = 2.4 GeV [14]. Since the HADES collaboration reported a strong dependence of cumulant ratios on the chosen rapidity window, we depict results for both |y| < 0.4 and |y| < 0.5 for comparison, the latter one also being applied to STAR data. In [14], a reliable rapidity window of |y| < 0.46 is quoted. We see that for high energies, our results to- gether with STAR data lie close to the Poisson baseline. With decreasing energy, our model yields values that are gradually enhanced until they reach a maximum at around 5 GeV for the evolution that passes through or close by the CEP and freezes out close to a spinodal line of the FOPT [57]. Along the spinodal lines, all susceptibilities diverge in nonequilibrium with critical exponents that become larger with order of the susceptibility [43,61]. This is clearly reflected in the peak of the corresponding susceptibility ratio visible on the bottom of figure 1. Further lowering the energy results in a gradual return to the baseline for the first hit of the freeze-out curve, while for the second hit, significantly larger fluctuations are observed, possibly due to an enhancement of spinodal instabilities. Notably, the data from HADES lies within the thus obtained range of cumulant ratios at the lowest beam energy. Obviously, a significant gap in the experimental data will have to be filled by future experiments such as FAIR and NICA that aim at exploring the high-µ B region of the QCD phase diagram. The susceptibility ratio on the right hand side of figure 1 shows a similar trend for high energies. The freeze-out close to the spinodal line, where susceptibilities in the presence of spinodal instabilities diverge and change sign [43,61], results in the strongly negative value of χ 2 /χ 1 which is reflected in the peak of C 2 /C 1 , similar to what we found for the other cumulant ratios as we will discuss below.
The cumulant ratio C 3 /C 2 is shown and compared to experimental data in figure 2 (top). The most notable feature is, again, the strong impact of the CEP around √ s = 5 GeV. Besides that, the obtained points from our model are close to the baseline for high energies and the second hit at the lowest energy is close to the data point from HADES, where a suppression of C 3 /C 2 was found, possibly a result of the dynamics at the FOPT. This is also found in the net-baryon number susceptibility ratio on bottom part of the same figure, approaching zero for the first hit and remaining negative for the second hit of the parametrized freeze-out curve. The most apparent difference is the sign at the CEP evolution which is positive for the susceptibilities, but negative for the cumulants. As mentioned before, the evolution passing close to the CEP freezes out very close to the spinodal line where some of the susceptibilities change sign [43]. Therefore, finite-time effects can here dramatically influence the final values at freeze-out.
Finally, the ratio C 4 /C 2 is depicted in figure 3, on the top we see qualitative similarities between our model results and the experimental data. A slight suppression in the STAR data around 20 GeV is also present in our model calculations where the cumulant ratio lies below the Poisson baseline, although with smaller significance. Then, as the beam energy is lowered, the notable point at 7.7 GeV where C 4 /C 2 is enhanced is reflected in an enhancement, albeit orders of magnitude larger, of the ratio from our calculation. Here, it is necessary to emphasize that the aforementioned freezing out near the spinodal line leads to larger and larger cumulants at higher orders. Lowering the beam energy even further, our calculations approach the HADES results which for this cumulant ratio show the strongest dependence on the applied experimental cut. Within error bars, both points lie within our range defined by the first and second hit of C 4 /C 2 ∼ 1 − 10. The susceptibility ratios, shown on the bottom of the same figure, are close to zero at these low energies which could indicate that an enhancement of the net-proton number cumulants occurs through a prolonged evolution in the mixed-phase region for the FOPT. The positive peak for the CEP evolution is also found in the susceptibilities, however, for larger energies, the values are slightly negative, which is only partly reflected in the obtained cumulants.

IV. SUMMARY AND CONCLUSIONS
We have studied cumulant ratios C 2 /C 1 , C 3 /C 2 , C 4 /C 2 of the net-proton number at STAR and HADES energies within a nonequilibrium chiral Bjorken expansion. Here, a sigma model served as input for a generic chiral phase structure and net-proton cumulants have been calculated event-by-event from cumulants of the sigma field at a parametrized freeze-out curve. Volume fluctuations have been accounted for and were properly corrected. Although admittedly crude and neglecting effects of an inhomogeneous medium, the dynamical description nevertheless shows some qualitative resemblance to the experimental data, in the approach of the Poisson baseline for high energies far away from the critical region, but also the enhancement or suppression of certain cumulant ratios at a speculated CEP or FOPT. We have demonstrated the general impact of a CEP and FOPT on cumulant ratios as key observables for probing the QCD phase structure. If indeed a CEP is present for center-of-mass energies below 7.7 GeV, it should be clearly visible in a relatively wide range energy range and manifest itself and the adjacent FOPT through an enhancement and/or suppression of cumulant ratios. Clearly, the current gap in beam energies from 2.4 for 7.7 GeV requires filling by future experiments.
Possible future improvements of our current model include the consideration of a spatially inhomogeneous fluid and an extension of the study to full (3+1) dimensional hydrodynamics.