Impact of resonance decays on critical point signals in net-proton fluctuations

The non-monotonic beam energy dependence of the higher cumulants of net-proton fluctuations is a widely studied signature of the conjectured presence of a critical point in the QCD phase diagram. In this work we study the effect of resonance decays on critical fluctuations. We show that resonance effects reduce the signatures of critical fluctuations, but that for reasonable parameter choices critical effects in the net-proton cumulants survive. The relative role of resonance decays has a weak dependence on the order of the cumulants studied with a slightly stronger suppression of critical effects for higher-order cumulants.


Introduction
Hot and dense strongly interacting matter is studied experimentally in heavy-ion collisions at different beam energies √ s. One of the major scientific goals of these experiments is to explore the phases of QCD matter as a function of temperature T and baryon chemical potential μ B . At vanishing μ B the transition between the high-T partonic and low-T hadronic phases is an analytic crossover [1], but many models predict that the transition is first order at large μ B [2][3][4]. Experimental searches for the first-order phase transition and the associated critical endpoint are based on the observation that the baryon chemical potential at chemical freeze-out increases with decreasing √ s in the regime probed at the Relativistic Heavy Ion Collider (RHIC) and the SPS fixed target program at CERN. Promising observable signatures for the presence of a critical point include event-by-event fluctuations [5][6][7], in particular fluctuations of conserved quantum numbers such as baryon number or electric charge [8,9]. The influence of a phase transition, or a change of its order, should be reflected in a non-monotonic √ s-dependence of these a e-mail: mbluhm@ncsu.edu observables [10][11][12][13][14][15]. Moreover, non-uniform structures in e.g. the net-baryon number distributions due to domain formation in the spinodal region of the phase diagram may be expected [16,17]. Motivated by these ideas a substantial experimental effort has been made in studying event-by-event multiplicity fluctuations of identified particle species [18][19][20]. Considered as a proxy for net-baryon number fluctuations, net-proton number fluctuations are particularly interesting. For these, intriguing non-monotonic patterns were observed within phase I of the RHIC beam energy scan (BES) program as √ s was decreased [21][22][23]. To deduce from the observed behavior unambiguously the presence of a first-order phase transition or a critical point is, nonetheless, a non-trivial task which requires additional work. On the experimental side the main issue is to reduce systematic and statistical uncertainties, which is one of the central goals of the upcoming phase II of the BES. Moreover, the strong sensitivity of the current data on reconstruction efficiency, kinematic acceptance or collision-centrality needs to be better understood. On the theoretical side a list of challenges has to be resolved before the data can be interpreted meaningfully. For example, net-proton number fluctuations serve, at best, as a proxy for net-baryon number fluctuations. In particular at lower √ s, aspects of exact charge conservation [24,25], repulsive interactions [26] or late stage hadronic phase processes such as resonance regeneration [27,28] or decay have to be considered carefully. Further details and more fluctuation sources were discussed in the recent review [29].
The present work addresses one of the theoretical challenges to data interpretation, namely the role of resonance decays. Resonance decays contribute significantly to the final state particle multiplicities and, consequently, to event-byevent multiplicity distributions. Their impact on net-proton number fluctuations, ignoring the possible presence of a critical point, was previously studied in [30,31]. It was found that the decay processes do not significantly modify the primordial (without resonance decays) ratios of net-proton number cumulants if the probabilistic decay contributions are properly taken into account [31]. This is because a thermal distribution of particles is being folded with a binomial decay distribution yielding again a nearly thermal distribution. However, critical fluctuations lead to sizable deviations from the thermal baseline. In [28] it was argued that the isospin randomization through resonance regeneration processes suppresses strongly the critical fluctuation signals in the netproton number distribution, in particular, in the higher-order cumulants. It is natural to ask if resonance decays similarly obscure critical point signals in fluctuation observables.
Critical fluctuations can be described in terms of an order parameter field σ . In the thermodynamic limit the equilibrium correlation length ξ of the order parameter diverges at the critical point. Under the assumption that QCD belongs to the same universality class as the three-dimensional Ising model [3,4] one finds that the cumulants of the order parameter fluctuations scale as (δσ ) 2 ∝ ξ 2 , (δσ ) 3 ∝ ξ 9/2 and (δσ ) 4 c ∝ ξ 7 , where small anomalous dimension corrections are neglected in these scaling relations [10]. Thus, higher-order cumulants depend more sensitively on the value of the equilibrium correlation length. This is important for heavy-ion collisions, because the actual value of ξ is limited due to finite size and time effects. Most significantly, the temporal growth of the correlation length is controlled by the relaxation time τ ξ , which increases as a power of the equilibrium correlation length. Near the critical point τ ξ ∼ ξ z with z ≈ 3 [32], and critical slowing down places significant limits on the actual growth of the correlation length. In [33], the non-equilibrium evolution of the correlation length was studied and a maximal value of about 2-3 fm was found. Thus, critical fluctuation signals in observables will be weaker than expected based on equilibrium considerations. Within the model of non-equilibrium chiral fluid dynamics [34,35] fluctuations of the order parameter are explicitly propagated in a fluid dynamical evolution of heavy-ion collisions. The effect of such a dynamics on the time evolution of net-proton fluctuations has been investigated recently [36]. It was found that critical fluctuations are affected by the dynamics, and that the location on the hydrodynamic trajectory where fluctuations freeze out determines the magnitude and shape of the fluctuation signals. However, it was also found that a critical point signal remains visible on top of thermal and initial state fluctuations. An analytic model for the time evolution of order parameter fluctuations was recently studied in [37]. The authors found that relaxation time effects introduce a significant lag in the cumulants relative to equilibrium expectations, and that non-equilibrium effects modify the simple scaling relations discussed in [10].
In this work we assume that there is a critical point in the QCD phase diagram, that this point causes critical fluctua-tions in the primordial net-proton number at chemical freezeout, and that these fluctuations are determined by the growth of the correlation length. We study quantitatively the impact of resonance decays on the higher-order net-proton number cumulant ratios. We do not take into account dynamical effects such as the role of the relaxation time and critical slowing down. We assume that (a) the fluctuations reflect thermal equilibrium at chemical freeze-out, and that (b) freeze-out takes place in proximity to the QCD phase transition. These assumptions are motivated by the success of fluid dynamics, based on the assumption of rapid local thermalization, in describing bulk observables [38,39], and by the successful extraction of freeze-out parameters from fluctuation observables [40]. Certainly, results obtained in the limit of rapid local thermalization provide an important benchmark for more sophisticated studies based on non-equilibrium thermodynamics.

Net-proton fluctuations and the impact of a critical point
In this section, we discuss in detail the framework used in our work. We start by describing the thermal, i.e. noncritical, baseline which we assume to be given by a grandcanonical hadron resonance gas model. It was shown that such a model is quite successful in quantitatively describing both lattice QCD results of the thermodynamics in the hadronic phase [41,42] and measured particle yields [43]. Critical fluctuations are then included by making use of universality class arguments and coupling (anti-)particles with the order parameter field, the critical mode, via a phenomenological approach. We discuss how the critical fluctuations can be quantified at the chemical freeze-out and, furthermore, propose different phenomenological ansätze to study, in addition, the impact resonance decays have on the netproton number fluctuations in the presence of a critical point.

Non-critical baseline model
As a baseline theory that exhibits no critical behavior we consider a hadron resonance gas model containing baryons and mesons up to masses of approximately 2 GeV in line with [44]. Starting point for our study is the particle density of species i given by the momentum-integral with degeneracy factor d i and distribution function In Eq. (2), E i = k 2 + m 2 i is the energy of particles of species i with mass m i , and μ i = B i μ B + S i μ S + Q i μ Q is their chemical potential, where μ X are conjugate to the conserved charges N X and B i , S i , Q i represent the corresponding quantum numbers of baryon number, strangeness and electric charge, respectively.
For the grand-canonical ensemble of a non-interacting hadron resonance gas, one defines the mean particle number N i of species i in a fixed volume V as N i = V n i . The cumulants with n > 1, measuring thermal ensemble fluctuations in the particle number, can be related to derivatives of n i /T 3 with respect to μ i /T for fixed T , In this work, we are interested in the first four cumulants C n of the net-proton number N p−p = N p − Np, which read (note that cumulant and central moment are the same for order n ≤ 3) The second equalities in Eqs. (5)-(7) are specific for our baseline model, for which correlations among different particle species vanish.
It is useful to consider particular ratios of these cumulants. Assuming that the fluctuations originate from a source with fixed thermal conditions, we may define and compare model calculations for the cumulant ratios with combinations of mean M = C 1 , variance σ 2 = C 2 , skewness S = C 3 /C 3/2 2 and kurtosis κ = C 4 /C 2 2 of measured event-by-event multiplicity distributions. An advantage of considering ratios is that to first approximation volume fluctuations, which are inherent in the individual cumulants C n , cancel out.

Critical mode fluctuations and coupling to physical observables
In the scaling region near a critical point different physical systems, which belong to the same universality class, exhibit a universal critical behavior [45]. Assuming that QCD belongs to the same universality class as the threedimensional Ising model allows us to relate the order parameter σ of the chiral phase transition in QCD with the order parameter of the spin model, the magnetization M. In the Ising model, the magnetization is a function of reduced temperature r = (T − T c )/T c and reduced external magnetic field h = H/H 0 , where T c is the critical temperature in the spin model and H 0 is a normalization constant. In these variables, the critical point is located at r = h = 0. For r < 0 and h = 0 the transition is of first order, while r > 0 marks the crossover regime. The equilibrium cumulants of the critical mode are then obtained via from appropriate derivatives of M with respect to the classical field keeping r fixed [14,15,37]. A parametric representation of M based on [46] can be found in Appendix A. This gives rise to parametric expressions for the cumulants of the critical mode, which are also given in Appendix A. These expressions are, strictly speaking, only valid for the scaling region close to the critical point. Nevertheless, we employ them for all considered √ s because (a) the size of the scaling region is not known, and (b) the influence of critical fluctuations on observables is suppressed further away from the critical point as we discuss in Sect. 3.
Fluctuations of the critical mode are not directly measurable. Instead, we focus on observables that couple to the critical mode, such as the number of pions or protons [6,7]. In general, the coupling of these observables to the critical mode has to be modeled. In the following we will take the order parameter to be the chiral field σ , and model the coupling of (anti-)particles to the critical mode in a way that is suggested by models in which the chiral field has a linear relation to the dynamically generated mass, δm i = g i δσ [6,12,14]. As a consequence, fluctuations associated with the variations in the particle masses arise in addition to the thermal ensemble fluctuations discussed in Sect. 2.1. In Eq. (10), The coupling of protons and anti-protons to the critical mode with coupling strength g p influences the net-proton number fluctuations near the critical point. According to Eq. (10), critical fluctuations do not affect the mean value C 1 but modify the variance with In these expressions we assumed that non-critical and critical fluctuations are uncorrelated. The coupling to the critical mode induces correlations between protons and anti-protons in C 2 , which reduce the variance compared to that in the case of an independent production. Simultaneously, critical fluctuations enhance the individual proton and anti-proton variances according to Eq. (11). With the help of Eq. (13) we can write C 2 in compact form as and find as net-effect an enhancement of the net-proton variance due to critical fluctuations. Similarly, C 3 and C 4 are modified. For where we neglected subdominant critical fluctuation contributions of order lower than O(ξ 9/2 ). Accordingly, we can write For for m + n = 4, where we also neglected subdominant critical fluctuation contributions, here of order lower than O(ξ 7 ). In compact form one may write We note that if non-critical and critical fluctuations are independent, then Eqs. (17) and (20) will represent the full results for C 3 and C 4 , respectively, in line with the model assumption Eq. (10).

Mapping to QCD thermodynamics
In order to quantify the effect of the critical point on netproton number fluctuations, we have to specify the cumulants of the critical mode (δσ ) n c at chemical freeze-out. In this work we use the chemical freeze-out conditions reported in [40], which are shown in Fig. 1 (squares). These were determined by analyzing the STAR data [18,19] on netproton number and net-electric charge fluctuations. A discussion and practical parametrization of the results (exhibited by the solid curve in Fig. 1) is relegated to Appendix B.
For any set of values of Ising model variables r and h, we can determine the cumulants (δσ ) n c from the parametric representation discussed in detail in Appendix A by using the transformations in Eqs. (A.2) and (A.3). Therefore, we only have to map the thermal variables in QCD, T and μ B , to the variables of the spin model. This mapping is non-universal: The relation between r and h on the one hand, and T and μ B on the other, depends sensitively on model assumptions. The easiest approach, used frequently in the literature [37,47], is a linear map r (T, μ B ) and h(T, μ B ) subject to two conditions: (a) the conjectured QCD critical point (μ The spin model coordinate system (r, h) has its origin in the critical point. Details of the mapping are described in the text. We also show the chemical freeze-out conditions determined in [40] (squares) and a practical parametrization of these results (solid curve), which is discussed in Appendix B The exact location of the critical point and the slope of the adjacent first-order phase transition line are not known. In this work we use lattice information on the location of the pseudo-critical line in the (μ B , T )-plane, indicated by the filled band between the two dashed curves in Fig. 1. The lattice data can be parameterized as where T c,0 = (0.145 · · · 0.163) GeV is the pseudo-critical temperature at vanishing μ B [48,49] and κ c 0.007 · · · 0.059 is the chiral crossover curvature [50][51][52][53][54].
The mapping to h is not well constrained. A typical choice made in the literature is to define the h-axis perpendicular to the r -axis [37,47]. However, since the first-order phase transition line is expected to have some curvature in the (μ B , T )plane it is clear that, in general, the map must be more complicated. In this work, we follow [37] and define T cp (22) such that the h-axis is simply parallel to the T -axis (see Fig. 1). We also define the auxiliary variablẽ The parameters T cp and μ cp B determine the size of the critical region. The Ising variable r is then obtained by rotatingr , where the slope of the r -axis is determined by using the slope of T c (μ B ) at μ cp B .

Resonance decays
The cumulants of a multiplicity distribution influenced by resonance decays were derived in [26,55]. These results take the probabilistic nature of the decay process into account. As critical fluctuations do not alter mean values, cf. Eq. (10), the mean of the final net-proton number readŝ where the notationĈ n refers to a cumulant that includes resonance decays. We also define n i R = r b R r n R i,r as the decay-average number of particles of species i created from resonance R, where n R i,r are created in decay branch r with branching ratio b R r . Note that the sum in Eq. (24) includes all resonances and anti-resonances of the baseline hadron resonance gas model. Analogous formulas for the higher cumulants, which contain critical fluctuation contributions, are given in Appendix C.
To study the impact of resonance decays on the higherorder cumulants of the net-proton number in the presence of a critical point we must quantify the coupling strengths g R between the critical mode and the (anti-)resonances. In line with the chiral model study [56] we may assume that g R is proportional to g p for non-strange baryonic resonances like ++ with a proportionality factor that reflects the mass difference between the resonance and the proton. For strangeness carrying resonances like (1520) we must, in addition, account for the fact that their mass generation is less dominated by the coupling to the chiral field [56]. To accommodate these features we advocate as a phenomenological ansatz for the σ -(anti-)resonance couplings. To analyze the impact of other possible scenarios, we also consider the two cases that either all resonances couple like protons, g R = g p , or that resonances do not couple to the critical mode at all, g R = 0. This allows us to study how sensitively the final results depend on the choice for the couplings g R .

Numerical results
In this section we present numerical results of our study. We start by discussing how to implement kinematic cuts in order to take into account the experimental acceptance. Then, we quantify the effect of resonance decays on critical point signatures in net-proton number fluctuations. We compare our results to preliminary STAR data [23]. We should note that there are sizable uncertainties in both the data and the theoretical analysis. Our goal in this work is not to establish or constrain the existence of a critical point, but to quantify the role of resonance decays.

Kinematic cuts
In the analyses [18,[21][22][23] of net-proton number fluctuations the detector acceptance was limited. Coverage in kinematic rapidity was restricted to the mid-rapidity region |y| ≤ 0.5, with a full coverage in azimuthal angle φ, and in transverse momentum cuts of 0.4 GeV ≤ k T ≤ k max T with k max T = 0.8 GeV [18] or k max T = 2 GeV [21][22][23] were applied. For a smaller experimentally analyzed phase space, one expects critical fluctuation signals in observables to be diminished because the momentum-space correlations, which are caused by the potentially large spatial correlations through thermal smearing, may exceed the considered kinematic acceptance window [57,58]. Indeed, the pronounced structures seen in the preliminary data of the higher-order net-proton number cumulant ratios [21][22][23] disappear to a large extent in the published data [18] for smaller k max T . In this work, we include the kinematic cuts on the level of the momentum-integrals, cf. Eqs. (1) and (13). For this purpose, we rewrite following [30] limit the integrations to −0.5 ≤ y ≤ 0.5, 0 ≤ φ ≤ 2π and 0.4 GeV ≤ k T ≤ k max T and replace E i by cosh(y) k 2 T + m 2 i . This implies that we apply the kinematic cuts to evaluate quantities at the chemical freeze-out and neglect the subsequent modifications of particle multiplicity distributions in a given acceptance window due to elastic scatterings until kinetic freeze-out.

Cumulant ratios
Before we can show our results we have to fix the parameters that enter into our calculation. These include the pseudocritical line parameters T c,0 and κ c , the coordinates of the critical point (μ  Fig. 5 (see Appendix A). We observe that the maximum of the correlation length is about 2 fm, comparable to estimates based on dynamical models.
To quantify the influence of critical fluctuations on the netproton number cumulant ratios, we also need to specify the value of the coupling strength g p between the (anti-)protons and the critical mode, which enters the cumulants via the factors I p and Ip, cf. Eq. (13). In the context of the linear sigma model the coupling constant g p in the ground state can be related to the pion decay constant [7,10], g p m p / f π 10. In chiral models, parameters are typically fixed to reproduce the vacuum masses of baryons and mesons as well as other known properties of nuclear matter at saturation density. In the non-linear chiral model considered in [59] to describe QCD matter in neutron stars, for example, |g p | 10 was used. Quark model calculations [60] of nucleon-meson couplings find g p in the range of 3-7, instead. To illustrate the   14), (17) and (20) including critical fluctuations for the exemplary σ -(anti-)proton couplings g = g p = 3, 5, 7. These results are obtained along the chemical freeze-out curve shown in Fig. 1 and include the kinematic cuts applied in [23]. For guidance, the preliminary data [23] are also shown. See text for more details importance of the actual value of g p , we use three different values g p = 3, 5 and 7 in the following.
In Fig. 2, we show the net-proton number cumulant ratios Sσ and κσ 2 defined in Eq. (8) without including the contributions from resonance decays. The non-critical baseline (thin dashed curves) is determined using Eqs.  The filled bands contain the results of the three different ansätze for g R discussed in this work: the strongest reduction of baseline deviations is seen for g R = 0 when resonances do not couple to the critical mode, while critical fluctuations are least suppressed by resonance decays for g R following the phenomenological ansatz Eq. (25). The case g R = g p lies in between, but closer to the phenomenological ansatz itative features are driven by the sign changes in (δσ ) 3 and (δσ ) 4 c as discussed in Appendix A. The effects are found to be stronger for larger g p and in the region of beam energies, where the correlation length ξ at chemical freeze-out is largest, cf. Fig. 5 in Appendix A. Further away from this region, at small and large √ s, the contributions from critical fluctuations are suppressed and the net-proton number cumulant ratios fall back to the baseline results. At small √ s, this is a consequence of the smallness of ξ further away from the critical point. But more importantly, in particular for large √ s, this is caused by the thermal factor (I p − Ip)/g p , which is positive definite and a monotonically decreasing function with increasing √ s along the chemical freeze-out curve. In Fig. 3 we show the impact of resonance decays, focusing on the case g p = 5. We determine the final net-proton number cumulants including critical fluctuations via the expressions summarized in Appendix C and consider the three different scenarios for the couplings g R between the critical mode and (anti-)resonances discussed in Sect. 2.4. We find that resonance decays influence the net-proton number cumulant ratios visibly, reducing the non-monotonicity induced by the critical fluctuations. Critical fluctuation signals are suppressed in Sσ by about 40% and in κσ 2 by about 50%, i.e. the reduction effect is slightly stronger for higher-order cumulants. This behavior is qualitatively in agreement with the predicted impact of isospin randomization processes on critical fluctuation signals [28].
The results shown in Fig. 3 indicate that resonance decays have significant effects on the cumulant ratios, but that for a sufficiently strong coupling g p between critical mode and (anti-)protons critical fluctuations survive final state processes such as resonance decays. Details in the couplings g R play a rather minor role as evident from Fig. 3. The strongest suppression of critical fluctuation signals is found for g R = 0, while the least is observed for g R in line with the phenomenological ansatz Eq. (25).

Discussion
The results discussed in Sect. 3.2 are obtained for fixed values of g p . In general, it is reasonable to assume that the coupling strengths between the critical mode and (anti-)particles depend on √ s. For example, in [6] it was argued that at the critical point the σ -pion coupling could be significantly smaller than in vacuum. Moreover, quark-meson model [61] and Nambu-Jona-Lasinio (NJL) model [62] studies find that meson-nucleon couplings decrease both with increasing T and/or μ B . Thus, by phenomenologically adjusting g p ( √ s), as done in [63], a quantitative model-to-data comparison would certainly be feasible, which is, however, beyond the scope of the present work.
The specific behavior of the skewness, which shows a reduction of Sσ compared to the baseline result, depends on our choices for the direction of the h-axis relative to the T -axis, the relative sign of δσ and h, and the sign of g p . For simplicity we have taken all these signs to be positive. It is important to note that the sign of the critical contribution changes, leading to an enhancement of Sσ relative to the noncritical baseline, if either of these signs is changed. Indeed, in the case of the liquid-gas phase transition in ordinary liquids, the order parameter is δn = n −n cr , which is mapped to the positive h-direction in the Ising model. This supports the idea of aligning h with the positive T -axis in the QCD phase diagram, but it also suggests, based on Eq. (10) and g > 0, that δσ is aligned with the negative h-axis. The assumption g > 0 is consistent with the simple sigma model picture that g controls the dynamically generated mass, m ∼ gσ . An NJL model study [64] of the σ -mode finds that in the critical region positive fluctuations in σ are associated with STAR prel. HRG baseline crit net−p + reso reduced cuts Fig. 4 Cumulant ratio σ 2 /M as a function of √ s including resonance decays and critical fluctuation contributions for g p = 5 (both filled bands). As in Fig. 3, for g R = 0 the results are closest to the noncritical baseline (thin dashed curve). For guidance, we also show the preliminary data [23]. Reducing the kinematic acceptance, for example by decreasing the value for k max T from 2 GeV (upper band) to 0.8 GeV (lower band), reduces the non-monotonic behavior caused by critical fluctuations slightly positive fluctuations in the net-baryon density. In our simple phenomenological approach this corresponds to g < 0, cf. Eq. (10). A different NJL model study [65] shows that the baryon skewness is negative on the high (μ B , T ) side of the transition line. Again, this result corresponds to changing one of the signs in our model.
It is also important to quantify the influence of a critical point on the lowest-order cumulant ratio σ 2 /M. According to Eq. (14), σ 2 /M is, compared to the non-critical baseline, increased through the critical fluctuation contributions in C 2 . In contrast, neither the published [18] nor the preliminary STAR data [23] show within errors any sign of non-monotonic behavior in σ 2 /M. This puts significant constraints on the values of g p , μ cp B , T cp , M 0 and H 0 , and the location of the critical point relative to the chemical freezeout in the QCD phase diagram.
In Fig. 4, we show σ 2 /M including the contributions from resonance decays for the case g p = 5 along the parametrized freeze-out curve. The filled bands comprise again the results of the three different ansätze for g R , cf. Fig. 3. For comparison also the preliminary STAR data [23] and the non-critical baseline (thin dashed curve) are shown. The non-monotonic structures seen in Fig. 4 are caused by the critical fluctuations. By decreasing k max T from 2 GeV (upper band) to 0.8 GeV (lower band), the non-monotonic behavior is only slightly reduced. The observed reduction is a consequence of the fact that the thermal factors I i in the cumulants correlate (anti-)particles of different momenta in the distribution with each other. By decreasing the considered kinematic acceptance, these correlations are reduced, as is also discussed e.g. in [57,58]. The effect is found to be stronger for higher-order cumulant ratios. Moreover, increasing the considered rapidity range from |y| ≤ 0.5 to |y| ≤ 0.75 instead of changing k max T from 0.8 to 2 GeV has a quantitatively comparable influence, cf. also [57].

Conclusions
In this work, we have quantified the impact of resonance decays on higher-order net-proton number cumulant ratios for the case that the primordial net-proton number distribution contains critical fluctuations induced by the presence of a QCD critical point. The critical fluctuations were included on top of a non-critical background described by a thermal hadron resonance gas model. The coupling of particle and anti-particle resonances to the critical mode was modeled using a sigma model coupling. Fluctuations of the critical mode were determined, applying universality class arguments, from derivatives of the order parameter in the threedimensional Ising model. We employed a linear map to relate QCD chemical freeze-out parameters (μ fo B , T fo ) to the Ising variables (r, h).
We find that resonance decays reduce, but do not completely obscure, the non-monotonic beam energy dependence of cumulant ratios induced by a critical point. The suppression is found to be slightly stronger for higher-order cumulants, in qualitative agreement with the influence of isospin randomization processes in the hadronic phase [28]. This is fairly independent of the strength of the coupling between (anti-)resonances and the critical mode. Shrinking the kinematic acceptance reduces the critical point signals as expected. The coupling to the critical mode results in correlations between protons and anti-protons both from the direct coupling and from the indirect coupling via resonances. The correlations turn out to be fairly small, modifying the cumulant ratios by a few percent at most.
Clearly, there are many issues that need to be addressed in order to make more quantitative model-to-data comparisons. This includes, in particular, the role of dynamical effects, both in the evolution of the order parameter and in the coupling to resonances.
to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix A: Critical mode fluctuations based on universality class arguments
A suitable parametric representation of the magnetization M(r, h) in the three-dimensional Ising model, which is valid in the scaling region around the critical point, is given in terms of the auxiliary variables R and θ as [46] withh(θ ) = cθ 1 + aθ 2 + bθ 4 being an odd function of θ which is regular near θ = 0 and θ = ±1. The parametric representation [46] uses the universal critical exponents β = 0.3250 and δ = 4.8169, while the coefficients inh(θ ) read a = −0.76145, b = 0.00773 and c = 1 such that the relevant non-trivial roots ofh(θ ) are located at θ 0 = ± √ 1.33128. Furthermore, the representation Eqs. (A.1)-(A.3) is uniquely defined for R ≥ 0 and −|θ 0 | ≤ θ ≤ |θ 0 |. For R > 0 and θ = 0, we have h = 0 and positive r , which corresponds to the crossover transition, while θ = θ 0 corresponds to the first-order phase transition line. The normalization constants M 0 and H 0 [66,67] are of mass dimension one and three, respectively, and the magnetization M is positive for (r < 0, h → 0 + ) and for (r → 0, h > 0).
We obtain the equilibrium cumulants of the critical mode within this parametrization from Eq. (9). They read as functions of R and θ 6 5 . (A.6) The variance (δσ ) 2 is, by definition, a positive function. The coefficients in (δσ ) 3 are (A.11) and in (δσ ) 4 c they read (A.14) (A. 20) In the limit of the so called linear parametric representation with critical exponents β = 1/3 and δ = 5 and coefficients a = −2/3, b = 0 and c = 3 these expressions simplify to the cumulants reported in [14,37].
One can study the dependence of the cumulants of the critical mode on the spin model variables by analyzing, for example, the dimensionless skewnessS = (δσ ) 3 / (δσ ) 2  Beam energy dependence of ξ along the considered chemical freeze-out curve (see Fig. 1). The correlation length falls off slower in the crossover regime, i.e. at larger √ s, as is expected from the underlying universality class. The depicted result will be described qualitatively with the phenomenological ansatz for ξ introduced in [12], if a comparable size of the critical region is considered θ . The skewness is an odd function of h, cf. also [15,37], which is positive for h < 0 and changes sign continuously (discontinuously) for positive (negative) r when h becomes positive. The kurtosis, instead, is an even function of h, which for a given r > 0 is negative within an interval |h| < h * (r ) which depends on r . For r < 0, one finds positiveκ as function of h with a cusp at h = 0. As a function of r ,S exhibits a maximum (minimum) for negative (positive) h at an r > 0 which itself depends on h, whileκ exhibits with decreasing positive r first a minimum followed by a maximum whose locations depend on the non-zero value of |h|.
Let us compare the parametric expressions for the cumulants of the critical mode in Eqs. (A.4)-(A.6) with those obtained from the corresponding Ginzburg-Landau effective potential. Including cubic and quartic interactions one finds [10,37] (δσ ) 2 where the correlation length dependence of the cubic and quartic coupling strengths as known in the scaling region has already been taken into account, such thatλ 3 andλ 4 are dimensionless parameters. Comparing (δσ ) 2 in Eq. (A.4) with Eq. (A.21) allows us to determine the √ s-dependence of ξ along the chemical freeze-out curve considered in this work. This is shown in Fig. 5. One finds that the correlation length falls off slower in the crossover regime of the QCD phase diagram in line with the underlying universality class. Details of this behavior depend very sensitively on the chosen values for the normalization constants M 0 and H 0 , but also on the details of the mapping such as the size of the critical region and the proximity of the chemical freeze-out to the QCD phase transition. For the parameter choices specified in Sect. 3.2 we find that ξ becomes as large as 2.3 fm for √ s 19.2 GeV, which corresponds to a freeze-out baryon chemical potential μ fo B significantly smaller than μ cp B . A similar comparison of the higher-order cumulants can be made to determine the behavior of the dimensionless parametersλ 3 andλ 4 . These are found to be independent of R and, thus, remain finite in the domain −|θ 0 | ≤ θ ≤ |θ 0 | of the parametric representation. In particular,λ 3 turns out to be negative for h < 0. This implies that the sign of the skewnessS depends sensitively on the probed values of h.

Appendix B: Parametrization of the chemical freeze-out curve
The chemical freeze-out conditions reported in [40] were determined by analyzing the published STAR data [18,19] of σ 2 /M for event-by-event net-proton number and netelectric charge distributions. The qualitative behavior of these observables is dominated by charged pions and (anti-) protons. The published results for net-proton number fluctuations [18] were obtained by considering a significantly smaller kinematic acceptance in transverse momentum k T compared to [21][22][23]. Within errors, no pronounced deviation from non-critical baseline expectations was found in the data [18,19] for the lowest-order cumulant ratios. This justified the application of the non-critical baseline model in the analysis [40].
In particular for larger √ s, the chemical freeze-out temperatures found in [40] are significantly lower than those obtained in traditional approaches in which yields or yield ratios are analyzed, cf. e.g. [68]. Moreover, for small μ B a positive slope of the chemical freeze-out curve is found, cf. Fig. 1. This is opposite to the curvature of the crossover line obtained in lattice QCD, as can be seen from Eq. (21) and the cited values for κ c . Nonetheless, both small values for the chemical freeze-out temperature [69,70] and a positive slope of the freeze-out curve for small μ B [70] have also been obtained in lattice QCD analyses of the same fluctuation observables [18,19].
A practical parametrization of the freeze-out conditions [40] can be given in form of functions of the beam energy √ s. We define analytic functions for the chemical potentials at chemical freeze-out μ fo X similar to [68] reading   [40] and hold for beam energies √ s ≥ 7 GeV. The non-vanishing values for μ Q and μ S for a given √ s are a consequence of physical side conditions for first cumulants, namely overall strangeness neutrality and an initial proton to baryon ratio of about 0.4 as realized in Au + Au or Pb + Pb heavy-ion collisions. Due to the smallness of μ fo B at high beam energies, the latter condition is still quantitatively appropriate at high √ s even though the created system is almost isospin symmetric at mid-rapidity.

Appendix C: Resonance decay contributions to netproton cumulants including critical fluctuations
The variance of the final net-proton number is given bŷ C 2 = ( N p ) 2 + ( Np ) 2 −2 N p Np , where the twoparticle correlators including resonance decays read [26,55] ( with ( n i ) 2 R = r b R r (n R i,r ) 2 − n i 2 R . Note again that the sums over R include all resonances and anti-resonances of the non-critical baseline model such that if R denotes an anti-resonance like¯ (1520), the correspondingR is the resonance (1520).
According to Eq. (C.2), resonance decays can induce correlations between different particle species. Nonetheless, without critical fluctuations protons and anti-protons would remain uncorrelated in our grand-canonical hadron resonance gas model approach because of baryon charge conservation in each decay. The individual contributions ( N R ) 2 and N R NR in Eqs. (C.1) and (C.2), which include critical fluctuations, follow from Eqs. (11) and (12) analogous to ( N i ) 2 and N i N j for (anti-)protons. Critical point induced correlations between a resonance and its anti-resonance are, therefore, inherited by their decay products. This induces additional correlations between protons and anti-protons on top of the primordial correlations discussed in Sect. 2.2.
From Eqs. (C.1) and (C.2) we find thatĈ 2 contains noncritical fluctuation contributions, which are independent of the correlation length, and critical fluctuation contributions of order O(ξ 2 ). This implies that the critical contribution to the resonance terms scales with the same power of ξ as in the primordial result Eq. (14). The situation is more complicated in the case of the cumulants of third and fourth order, see [26]. Consider the third-order cumulant. The critical contribution to the primordial term, Eq. (15), scales as ξ 9/2 . The resonance terms contain a contribution, induced by fluctuations ( N R ) 3 , that scales with the same power of ξ , but they also contain contributions, induced by the probabilistic nature of branching, that scale with lower powers of ξ . Since we do not keep subdominant terms, we drop these contributions in the resonance terms. The same statement applies to the fourth-order cumulant.
N p ( Np ) 2 the three-particle correlators including resonance decays were derived in [26]. Taking only non-critical fluctuation contributions and critical fluctuation contributions of order O(ξ 9/2 ) into account, we find for these correlators where ( n i ) 3 R = r b R r (n R i,r ) 3 − 3 n i R r b R r (n R i,r ) 2 + 2 n i order cumulants follow from the four-particle correlators including resonance decays, cf. [26]. Taking forĈ 4 only non-critical fluctuation contributions and critical fluctuation contributions of order O(ξ 7 ) into account, we find The individual contributions in Eq. (C.5) follow from Eqs. (3) and (18) and in Eq. (C.6) from Eq. (19).