Criticality of the net-baryon number probability distribution at finite density

We compute the probability distribution $P(N)$ of the net-baryon number at finite temperature and quark-chemical potential, $\mu$, at a physical value of the pion mass in the quark-meson model within the functional renormalization group scheme. For $\mu/T<1$, the model exhibits the chiral crossover transition which belongs to the universality class of the $O(4)$ spin system in three dimensions. We explore the influence of the chiral crossover transition on the properties of the net baryon number probability distribution, $P(N)$. By considering ratios of $P(N)$ to the Skellam function, with the same mean and variance, we unravel the characteristic features of the distribution that are related to $O(4)$ criticality at the chiral crossover transition. We explore the corresponding ratios for data obtained at RHIC by the STAR Collaboration and discuss their implications. We also examine $O(4)$ criticality in the context of binomial and negative-binomial distributions for the net proton number.


Introduction
Fluctuations of conserved charges are promising observables for exploring critical phenomena in relativistic heavy ion collisions [1][2][3][4]. A particular role is attributed to higher order cumulants of the net baryon number and electric charge fluctuations, which in a QCD medium can be negative near the chiral transition [5][6][7].
At physical values of quark masses, the phase transition in QCD is expected to change from a crossover transition at small values of the baryon chemical potential to a first-order transition at large net baryon densities. The first-order chiral phase transition, if it exists, then begins in a second-order critical point, the critical end point (CEP) [8]. Owing to the divergent correlation length at the CEP [7], and the spinodal phase separation in a non-equilibrium first-order transition [9], one expects large fluctuations of the netbaryon number in heavy ion collisions, at beam energies where the system passes through the first-order phase boundary or close to the CEP.
The conjectured existence of a CEP in the QCD phase diagram has so far not been confirmed by lattice QCD calculations (LQCD) [10,11]. At small values of the quark chemical potential, μ q /T < 1, LQCD exhibits a chiral crossover transition. There are indications that, in the chiral limit for light quarks, the QCD transition belongs to the universality class of 3-dimensional O (4) spin systems [12,13]. Thus, a promising approach for probing the phase boundary in heavy ion collisions, is to explore the fluctuations of the chiral phase transition, assuming O (4) criticality. Owing to the proximity of the chemical freeze-out to the chiral crossover at small values of the baryonic chemical potential, one may expect that the critical fluctuations are reflected in the data on conserved charges [14]. A baseline for the cumulants of charge fluctuations is provided by the hadron resonance gas (HRG) model, which reproduces the particle yields at chemical freeze-out in heavy ion collisions [15], as well as the LQCD equation of state in the hadronic phase [16,17].
At the CEP, which is expected to belong to 3d Z (2) universality class, the second and higher order cumulants diverge. By con- class, at vanishing baryon chemical potential, low-order cumulants remain finite, while the sixth and higher order cumulants diverge. 1 For non-zero quark masses, the divergences are replaced by a rapid variation of the cumulants near the crossover temperature, including changes of sign [18].
The fluctuations of the net-baryon number, more precisely of the net-proton number, were measured in heavy ion collisions by STAR Collaboration at RHIC [19][20][21][22]. Data on the mean (M), variance (σ ), skewness (S) and kurtosis (κ ) of the net-proton number were obtained in a broad energy range and for different centralities. These observables are linked to the cumulants χ n of the net-baryon number, and are accordingly modified by the critical chiral dynamics.
The deviations in the latter are small at the top RHIC energy, increase with the order of the cumulants at fixed collision energy, and show a non-monotonic dependence on the energy with a maximum at √ s N N 19 GeV. In the O (4) universality class one expects 2 χ 4 /χ 2 < 1, while in Z (2) this ratio is expected to be larger than unity [23]. Thus, the systematics of the ratios of cumulants in central Au-Au collisions, as measured by STAR [22], indicate that the observed deviations of the net proton number fluctuations from the HRG values may be attributed to O (4) criticality at the phase boundary, at least for The cumulants of a conserved charge are given by appropriate combinations of moments of the corresponding probability distribution. Thus, the behavior of cumulants near criticality must be reflected in the properties of the probability distribution. Moreover, it is expected that the critical behavior of the probability distribution depends on the universality class. Indeed, we have recently shown, that the structure of the probability distribution for the net baryon number depends on the properties of the critical chiral fluctuations [23,24]. In particular, we have argued, that at vanishing chemical potential, the residual O (4) critical fluctuation at physical pion mass leads to narrowing of the probability distribution relative to the Skellam function. This corresponds to a negative structure of the sixth order cumulant at the chiral crossover transition [24].
In this paper, we extend our previous studies to non-zero chemical potential and propose a method for identifying the characteristic properties of the net baryon probability distribution, which are responsible for the critical behavior of the cumulants at the chiral transition. We apply this method to the net proton probability distributions obtained by the STAR Collaboration in central Au-Au collisions at √ s N N ≥ 19 GeV. We also critically examine the question whether O (4) criticality can be captured by assuming that the baryon and antibaryon multiplicities are described by binomial or negative binomial distributions.
In this paper, we show that the (suitably rescaled) ratio of the net baryon probability distribution to the corresponding Skellam function reveals the critical narrowing of the probability distribution, which is characteristic for the O (4) scaling.

The net-baryon number probability distribution
In the grand canonical ensemble specified by temperature T , subvolume V and chemical potential μ, the probability distribution for the conserved charge N, is given by 1 For μ = 0, diverging cumulants appear already at third order. 2 Although the results of [23] were obtained in the chiral limit, it is plausible that they remain valid also for physical values of the quark masses. (1) where the canonical partition function Z (T , V , N) is obtained e.g. by a projection of the grand partition function Z(T , V , μ), In the HRG the probability distribution of the net-baryon number is, within the Boltzmann approximation, given by the Skellam function [21,25] where b = N b and b = Nb are the thermal averages of the number of baryons and anti-baryons, respectively.
The HRG model reproduces the particle yields in heavy ion collisions in a broad energy range from SIS to LHC. Furthermore, it describes the equation of state obtained in LQCD, as well as the first and second order cumulants of the net baryon number for temperatures below the chiral crossover temperature. On the other hand, as suggested in [18], the deviation of higher order cumulants and their ratios from the HRG results provides a potential signature for criticality at the phase boundary.
These considerations indicate, that the probability distribution of the HRG, the Skellam function, offers an appropriate baseline for P (N). Indeed, for small N, where the probability distribution is fixed by the non-critical lowest order cumulants the Skellam function provides a good approximation to P (N). On the other hand, for large N the two distributions differ, since the critical fluctuations modify the tail of the distribution, which in turn determines the higher cumulants. Thus, it is natural to consider the Skellam function as a reference for identifying criticality in the probability distribution of the net-baryon number [21]. Specifically, we show that the ratio of P (N) to the Skellam function exposes the effect of critical fluctuations.
We extract the characteristic features of the probability distribution near the chiral crossover transition within the O (4) universality class by applying the Functional Renormalization Group (FRG) approach to the quark-meson (QM) model [26][27][28]. The QM model exhibits the relevant chiral symmetry of QCD, and belongs to the same O (4) universality class [29,30].
The Lagrangian density in the QM model reads where q and q denote the quark and anti-quark fields coupled with the O (4) chiral meson multiplet φ = (σ , π). The last three terms in Eq. (4) constitute the mesonic potential with the symmetry breaking term. The thermodynamic potential is calculated in the QM model, within the FRG approach [26]. Applying the optimized regulator to the exact flow equation for the effective average action in the local potential approximation [26], the flow equation for the scale dependent thermodynamic potential density Ω k reads [5]  where ρ = (σ 2 + π 2 )/2 is the reduced field variable and n F and n B are the Fermi and Bose distribution functions, respectively.
The single particle energies of π , σ and q are given by: where Ω k and Ω k denote the first and the second derivatives of The full thermodynamic potential is given by the minimum of Ω k→0 (ρ). We solve the flow equation (5) numerically by making use of the Taylor expansion method [5,31]. At the ultraviolet cutoff scale k = Λ = 1 GeV, the initial condition Ω Λ (ρ) is fixed so as to reproduce the physical pion mass m π = 135 MeV, and the sigma mass m σ = 640 MeV. The strength of the Yukawa coupling is fixed to be g = 3.2 by the constituent quark mass To avoid the unphysical behavior of thermodynamic quantities at high temperatures, we include the higher momentum contributions, beyond the cutoff scale Λ, by accounting for the μand T -dependent thermodynamic potential obtained through the flow equation for a non-interacting gas of quarks and gluons [5,32].
In Ref. [24], the FRG approach was applied to compute P (N) of the net-baryon number within the QM model at μ = 0. In the present paper, we extend these studies to finite chemical potential and identify the qualitative structures of the net baryon probability distribution which are due to O (4) criticality at the chiral crossover transition. We also evaluate the ratio of the data on the net proton probability distribution [22] to the Skellam function and discuss the results in the light of the theoretical considerations.
In general, the probability distribution P (N) depends on the volume parameter. However, as shown in [24], the volume dependence of the re-scaled distribution reduced. This approximate scaling is valid for both the Skellam function, and the P (N) calculated within the QM model for sufficiently large V T 3 , and is exact for a Gaussian distribution. Thus, in the ratios of P (N) and Skellam, the leading volume dependence is canceled. At finite density where M > 0, the scaling property holds for δN = N − M, after shifting the mean.
In order to compute the cumulants χ n reliably, knowledge of the probability distribution P (N) for sufficiently large |δN| = |N − M| is needed. For a given order n, it is sufficient to know the distribution P (N) for |δN| ≤ N n , where N n grows with n and with the volume of the system [24]. In the O (4) universality class and at μ = 0, χ 6 is the first cumulant which exhibits criticality.
Thus, to be able to identify criticality in the distribution, we need to know P (N) in the range needed to obtain a converged result for χ 6 , i.e. for |δN| < N 6 . For a given volume, N 6 is determined by requiring that the sixth cumulant of the Skellam function is reproduced. As expected, we find that N 6 to a good approximation scales with √ V . Thus, in a plot of the ratio of P (N) to the Skellam function as a function of δN/N 6 , criticality is characterized by deviations from unity for |δN/N 6 | 1.
In Fig. 1(a) we show the ratio of P FRG (N), computed in the quark-meson model at μ = 0 within the FRG approach [24], and the Skellam distribution P S (N), with the same variance, as a function of δN/N 6 . The results are shown for different temperatures T /T pc , where T pc is the chiral crossover or pseudocritical temperature. This ratio exhibits a characteristic dependence on temperature, as T pc is approached from below. For T T pc , the ratio P FRG (N)/P S (N) is less than unity, indicating a narrowing of the probability distribution for larger |δN|, owing to O (4) criticality. Indeed, the decrease of the probability ratio for δN/N 6 1 near T pc is responsible for the negative values of χ 6 , which are characteristic of the chiral crossover transition in the O (4) universality class [24]. We note, that the narrowing of P (FRG) (N) for smaller values of δN/N 6 , seen in Fig. 1(a), can be partly attributed to the non-critical reduction of χ 4 . However, the smoothly decreasing tail of P (N), close to |δN/N 6 | 1, is entirely due to O (4) criticality, resulting in the characteristic shape and the negative structure of χ 6 . Consequently, shrinking probability distribution, relative to the Skellam one, can indeed be considered as a necessary condition for O (4) criticality [24]. At finite chemical potential, the probability distribution P (N) of the net-baryon number is no longer symmetric around the mean. Thus, it is not a priori clear, how the distribution is modified by The asymmetry of P (N) at μ = 0 appears due to the fugacity factor e μN/T in Eq. (1), which suppresses the contribution from N < 0 and enhances that from N > 0. Consequently, at finite chemical potential, the tail of the probability distribution P (N) is enhanced, and criticality is expected to appear at smaller |δN|, and thus also in lower order cumulants.  consistent with the fact that at finite μ, already the third cumulant exhibits O (4) critical behavior [18]. On the other hand, for negative δN the ratio of the distributions exhibits the opposite behavior, reflecting the asymmetry of the probability distribution at non-zero net baryon density.
At finite chemical potential and at large |δN|, the calculations of P (N) are difficult due to the oscillating nature of the integrand in the projection on the canonical partition function (2). The numerical integration yields reliable results only up to δN/N 6 ≤ 0.6.
Consequently, the complete χ 6 cannot be reconstructed due to insufficient information on the tail of the distribution. Nevertheless, the narrowing of the probability ratio shown in Fig. 1(b) clearly exhibits the characteristic features of P (N), which are due to O (4) criticality. Evidently the deviation of P FRG (N)/P S (N) from unity near T pc (μ) will grow with increasing |δN/N 6 | and μ.

O (4) criticality in heavy ion collisions
In heavy ion collisions, particle yields, charge densities and their variance are described consistently by the HRG model on the same chemical freeze-out line in the (T , μ B )-plane [14,22,33]. For a given collision energy one can identify a unique point on the freeze-out line. If the freeze-out takes place sufficiently close to the chiral crossover transition, the critical fluctuations are expected to leave a characteristic imprint on the cumulants and on the corresponding probability distribution.
In Fig. 2 we illustrate the expected structure of the probability distribution at chemical freeze-out by showing the QM model results for P FRG (N)/P S (N) at a set of points in the (T , μ) plane.
They lie on the approximate freeze-out line, defined by requiring the same variance per unit volume of the net baryon number as in the μ = 0 point. The μ dependence of the ratios, with a narrowing of the distribution for positive and a broadening for negative δN with increasing μ is characteristic for the critical region. As shown in Fig. 1, the distribution in a non-critical system exhibits the opposite trend, with a broadening for positive and a narrowing for negative δN.
In general, the measurement of higher order cumulants, which are particularly sensitive to criticality, need high statistics owing to the increasing importance of the tail of the distribution. Furthermore, the experimental conditions, such as acceptance corrections, must be under control in order to make a meaningful comparison of the measured cumulants and their probability distribution with theoretical predictions [34][35][36].
Recently the STAR Collaboration has published extensive results on the probability distribution of the net-proton number N p = Fig. 3. Ratios of the efficiency uncorrected probability distributions of the net-proton number P ( N p ) by STAR Collaboration [22] to the Skellam function P S ( N p ) with the same mean and variance as P ( N p ). The data are for the most central Au-Au collisions, with the number of events N ev > 100.
N p − Np and the corresponding cumulants up to the fourth order, obtained in heavy ion collisions [22]. Also preliminary results on the sixth order cumulant have been presented by STAR for several collision energies and centralities [36]. While the cumulant ratios measured by STAR [22] were efficiency corrected and tested against possible modifications due to volume fluctuations and accepted kinematical windows, the data on the probability distributions of the net proton number are uncorrected. Furthermore, in the model calculations, the net baryon number rather than the net proton number measured by the STAR Collaboration is considered. We assume, that the criticality due to chiral symmetry restoration, which appears in the net baryon number is also reflected in fluctuations of the net proton number. However, as shown in Refs. [37] and [38], the quantitative differences between cumulants of the net baryon and the net proton number are not excluded. Therefore, the significance of a direct comparison of model predictions with the measured probability distribution is a priori not clear-cut. Here we assume the isospin invariance such that the net baryon number fluctuations are equivalent to those of net proton number and invoke the generic properties of the probability distribution of the net baryon number due to O (4) criticality found in the QM model persist in the net proton number probability distribution.
Nevertheless, we have verified that the data on P ( N p ) obtained by STAR [22] are dominated by physics. The contribution of volume fluctuations to the data is small, as demonstrated by the approximate scaling of P ( N p ) with the standard deviation σ in central and semi-central collisions, for GeV. This scaling holds also for the Skellam and P (FRG) (N) distributions for sufficiently large volumes. Moreover, the ratios of cumulants computed directly from the uncorrected P ( N p ) data [22] exhibit similar systematics as the efficiency corrected ratios. These tests indicate that the data may yield at least a qualitative indication whether the measured distribution exhibits criticality or not.
In Fig. 3 we show the probability ratio P ( N p )/P S ( N p ) obtained from the uncorrected data [22] in the highest centrality bin, for √ s N N = 200, 62.4, 39, 27 and 19.6 GeV. The probability ratio is constructed using the same method as in Fig. 1. In order to avoid large uncertainties, we have restricted the data to those with more than 100 events. Consequently, the probability distributions are limited to |δN/N 6 | < 0.5. This implies, that the present statistics does not allow for a reliable estimate of the sixth order cumulant. Nonetheless, the ratios in Fig. 3 clearly exhibit a structure qualitatively similar to that shown in Figs. 1 and 2 which is a reflection of the underlying O (4) criticality. In particular, the characteristic narrowing of the probability distribution relative to the Skellam function for positive δN, and early drop of the ratio below unity for √ s = 19.6 GeV, are characteristic signatures for O (4) criticality at non-zero chemical potential. There are several potential contributions to the cumulants and the probability distribution from sources other than critical fluctuations [37,[39][40][41], as well as experimental issues e.g. regarding efficiency corrections [42]. Thus, a final conclusion on the criticality of P ( N p ) can be drawn only once the role of these effects has been sorted out.
In [22], the cumulant ratios Sσ and κσ 2 are analyzed with efficiency and centrality bin width corrections. By constructing Sσ and κσ 2 from the uncorrected P (N) data discussed above, we have found that the deviations from Skellam distribution are slightly smaller than that seen in the corrected ratios. However the systematics and the energy dependence are almost the same. Therefore, we regard the results shown in Fig. 3 as the lower limit for possible deviations from the Skellam function. We stress, that the present method provides a transparent framework, where such corrections can be included. If the narrowing of P ( N p ) relative to the Skellam, as seen in Fig. 3, is still observed after these corrections are included, then this will provide potential evidence for remnants of the chiral crossover transition in experimental data.

O (4) criticality and binomial distributions
In order to reveal the O (4) criticality in the net baryon number probability distribution and in the corresponding cumulants we have used the Skellam distribution as a reference. The Skellam function is the natural choice, since in heavy ion collisions data on particle yields, as well as the QCD thermodynamics, are well reproduced by the hadron resonance gas partition function. In the HRG, baryons multiplicity is distributed according to Poisson and the net-charge distribution is then given by the Skellam function.
However, it has been shown that the lowest cumulants of the net proton fluctuations and the corresponding probability distributions obtained by the STAR Collaboration are consistent also with negative binomial (NBD) or binomial (BD) distributions [22,43]. Thus, it is of interest to verify to what extent these distributions can describe critical fluctuations at the chiral transition. This study can be done within the QM model, where the O (4) critical structure of the cumulants and the corresponding probability distribution are manifest.
In the case of NBD or BD, the net-baryon probability distribution is constructed assuming independent emission of baryons and antibaryons, with n being the number of baryons or antibaryons. For μ = 0, the property of the net-baryon P (N) is uniquely determined by two parameters (r, p), characterizing the NBD or BD. Using the additive property of the cumulants, one finds We construct the reference NBD/BD, which has the same χ 2 and χ 4 as the O (4) distribution obtained in the QM model within the FRG method. 3 From Eqs. (8)- (13), it is clear, that χ 4 /χ 2 > 1 for NBD and χ 4 /χ 2 < 1 for BD. For χ 4 /χ 2 = 1, both distributions are reduced to the Skellam function for r → ∞ and p → 0. Therefore, we use NBD for temperatures where χ 4 /χ 2 > 1 and BD for χ 4 /χ 2 < 1.
This implies, in particular, that NBD cannot describe fluctuations around T pc , where χ 4 /χ 2 < 1. Fig. 4 shows the ratios χ 4 /χ 2 and χ 6 /χ 2 obtained in the QM model as functions of temperature near T pc and at μ = 0. Also shown are the corresponding ratios obtained from the NBD/BD distributions, where the parameters (r, p) are fixed so as to reproduce χ 2 and χ 4 . Clearly the ratio χ 6 /χ 2 , in particular the negative values near T pc , is not reproduced by the binomial distribution. We therefore conclude that the NBD/BD distributions clearly fail to describe the critical fluctuations near the chiral transition.
At μ = 0, the baryon and antibaryon distributions are different, so that there are two sets of parameters (r, p) in the NBD/BD, one for protons and one for antiprotons. Thus, in principle, one may construct a reference distribution for the net baryon number, P (N), which reproduces the four leading cumulants, 4 i.e. χ n with n = 1, 2, 3, 4, ignoring the fact that at non-zero μ, χ 3 and χ 4 may be affected by criticality. We stress that it is unclear whether these parameters also yield a good description of the baryon and antibaryon distributions. Moreover, also in this case, the NBD/BD distributions that reproduce the leading cumulants, cannot describe higher order ones, χ n with n > 4, nor the tail of distribution, P (N).
Thus, an unambiguous verification of critical fluctuations in heavy ion collisions requires knowledge of the tail of net proton distribution P (N), so that the sixth order cumulant can be determined reliably, although at non-zero μ, an indication of criticality may be exposed in cumulants of lower order.

Concluding remarks
We have discussed the properties of the net-baryon number probability distribution P (N) near the chiral crossover transition at vanishing and at finite baryon chemical potential. The critical properties of P (N) in the quark-meson model were obtained within the functional renormalization group approach. In the chiral limit this model exhibits a second order phase transition belonging to O (4) universality class in three dimensions.
At the physical value of the pion mass, the O (4) criticality is reflected in the tail of the distribution P (N). We have shown that the ratio of P (N) to the Skellam function P S (N), constructed with the same mean M and variance as P (N), clearly exhibits the influence of the O (4) criticality on the probability distribution.
We have shown that at vanishing chemical potential there is a characteristic reduction of this ratio below unity near the phase boundary. At finite chemical potential, the ratio P (N)/P S (N) exhibits a characteristic asymmetry in δN = N − M. For δN < 0, the probability ratio is enhanced near the O (4) pseudocritical point, while for δN > 0 it is suppressed. The asymmetry of the distribution is enhanced with increasing μ along the freeze-out line.
The relevance of our results for heavy ion experiments was discussed. In particular, we have computed the corresponding probability ratios for the efficiency uncorrected net proton number obtained by the STAR Collaboration, and discussed their interpretation. We have also demonstrated that O (4) criticality, in particular its reflection in higher cumulants of the net baryon number, is not consistent with a description of the baryon and antibaryon multiplicities in terms of binomial or negative-binomial distributions.
Finally we stress that an unambiguous identification of O (4) chiral criticality on the phase boundary requires high statistics data on the net proton probability distribution over a range in δN that allows a reliable determination of the sixth order cumulant.