Net-proton number fluctuations in the presence of the QCD critical point

TL;DR

This paper models net-proton fluctuations near the QCD critical point, reproducing some experimental fluctuation patterns and suggesting the critical point's influence on observed heavy-ion collision data.

Contribution

It introduces a phenomenological model coupling the critical mode to protons, providing insights into net-proton fluctuations and their relation to the QCD critical point.

Findings

Reproduces the non-monotonic kurtosis observed in experiments

Shows the variance-to-mean ratio changes smoothly with beam energy

Highlights discrepancies in skewness behavior compared to critical expectations

Abstract

Event-by-event fluctuations of the net-proton number studied in heavy-ion collisions provide an important means in the search for the conjectured critical end point (CP) in the QCD phase diagram. We propose a phenomenological model in which the fluctuations of the chiral critical mode couple to protons and anti-protons. This allows us to study the behavior of the net-proton number fluctuations in the presence of the CP. Calculating the net-proton number cumulants, $C_n$ with n=1,2,3,4, along the phenomenological freeze-out line we show that the ratio of variance and mean $C_2/C_1$, as well as kurtosis $C_4/C_2$ resemble qualitative properties observed in data in heavy-ion collisions as a function of beam energy obtained by the STAR Collaboration at RHIC. In particular, the non-monotonic structure of the kurtosis and smooth change of the $C_2/C_1$ ratio with beam energy could be due to the CP located near the freeze-out line. The skewness, however, exhibits properties that are in contrast to the criticality expected due to the CP. The dependence of our results on the model parameters and the proximity of the chemical freeze-out to the critical point are also discussed.

Figures

Tables

Table 1. Table 1: Considered locations of the QCD critical point in the ( μ B , T subscript 𝜇 𝐵 𝑇 \mu_{B},T )-plane. The parameters T c , 0 subscript 𝑇 𝑐 0 T_{c,0} and κ c subscript 𝜅 𝑐 \kappa_{c} of the crossover (pseudo-critical) line Eq. ( 30 ) needed to determine T c ​ p subscript 𝑇 𝑐 𝑝 T_{cp} and the slope of the first-order phase transition line at the critical point for a given μ c ​ p subscript 𝜇 𝑐 𝑝 \mu_{cp} are also listed. The locations of these critical points in the QCD phase diagram are shown in Fig.  3 .

$C{P}_i$	${\mu }_{cp}$[GeV]	$T_{cp}$[GeV]	$T_{c,0}$[GeV]	${\kappa }_c$
1	0.390	0.149	0.156	0.007
2	0.420	0.141	0.155	0.010
3	0.450	0.134	0.155	0.012

Equations

$$n_i(T,{\mu }_i)=d_i\int \frac{d^{3}k}{(2\pi {)}^3}{f}_i^0(T,{\mu }_i),$$

$$f_i^0=\frac{1}{(-1{)}^{B_i}+e^{(E_i-{\mu }_i)/T}}$$

$$C^{i}_{n}=VT^{3}\frac{\partial^{n-1}(n_{i}/T^{3})}{\partial(\mu_{i}/T)^{n-1}}\bigg{|}_{T}\,,$$

$$C_1=⟨N_{p-\bar{p}}⟩=⟨N_p⟩-⟨N_{\bar{p}}⟩,$$

$$C_2=⟨(\Delta N_{p-\bar{p}}{)}^2⟩=C_2^p+C_2^{\bar{p}},$$

$$C_3=⟨(\Delta N_{p-\bar{p}}{)}^3⟩=C_3^p-C_3^{\bar{p}},$$

$$C_4=⟨(\Delta N_{p-\bar{p}}{)}^4{⟩}_c=C_4^p+C_4^{\bar{p}},$$

$$\Delta N_{p-\bar{p}}=N_{p-\bar{p}}-⟨N_{p-\bar{p}}⟩,$$

$$⟨(\Delta N_{p-\bar{p}}{)}^4{⟩}_c=⟨(\Delta N_{p-\bar{p}}{)}^4⟩-3⟨(\Delta N_{p-\bar{p}}{)}^2{⟩}^2.$$

$$\frac{C_2}{C_1}=\frac{{\sigma }^2}{M},\frac{C_3}{C_2}=S\sigma ,\frac{C_4}{C_2}=\kappa {\sigma }^2,$$

$$\delta f_i=\frac{\partial f_i}{\partial m_i}\delta m_i=-\frac{g_i}{T}\frac{v_i^2}{{\gamma }_i}\delta \sigma ,$$

$$C_2=C_2^p+C_2^{\bar{p}}+⟨(V\delta \sigma {)}^2⟩(I_p-I_{\bar{p}}{)}^2,$$

$$C_3=C_3^p-C_3^{\bar{p}}-⟨(V\delta \sigma {)}^3⟩(I_p-I_{\bar{p}}{)}^3$$

$$C_4=C_4^p+C_4^{\bar{p}}+⟨(V\delta \sigma {)}^4{⟩}_c(I_p-I_{\bar{p}}{)}^4,$$

$$I_i=\frac{g_i{d}_i}{T}\int \frac{d^{3}k}{(2\pi {)}^3}\frac{v_i^2}{{\gamma }_i}$$

$$⟨(V\delta \sigma {)}^n{⟩}_c={\left(\frac{T}{V{H}_0}\right)}^{n-1}{\frac{{\partial }^{n-1}{M}_I}{\partial h^{n-1}}|}_r,$$

$${\chi }_B\simeq {\chi }_B^{\mathrm{reg}}+{\sigma }^2\chi .$$

$$C_2^{\mathrm{sing}}=⟨(V\delta \sigma {)}^2⟩m_p^2(J_p-J_{\bar{p}}{)}^2,$$

$$J_i=\frac{g_i{d}_i}{T}\int \frac{d^{3}k}{(2\pi {)}^3}\frac{v_i^2}{E_i}.$$

$${\gamma }_i\to E_i/(g_i\sigma ),$$

$$C_2=C_2^p+C_2^{\bar{p}}+g_p^2{\sigma }^2⟨(V\delta \sigma {)}^2⟩(J_p-J_{\bar{p}}{)}^2,$$

$$C_3=C_3^p-C_3^{\bar{p}}-g_p^3{\sigma }^3⟨(V\delta \sigma {)}^3⟩(J_p-J_{\bar{p}}{)}^3$$

$$C_4=C_4^p+C_4^{\bar{p}}+g_p^4{\sigma }^4⟨(V\delta \sigma {)}^4{⟩}_c(J_p-J_{\bar{p}}{)}^4.$$

$$M_I=M_0{R}^{\beta }\theta ,$$

$$r$$

$$h$$

$$w(\theta )=c\theta (1+a{\theta }^2+b{\theta }^4)$$

$$\tilde{r}=\frac{{\mu }_B-{\mu }_{cp}}{\Delta {\mu }_{cp}},$$

$$\tilde{h}=\frac{T-T_{cp}}{\Delta T_{cp}},$$

$$T_c({\mu }_B)=T_{c,0}\left[1-{\kappa }_c{\left(\frac{{\mu }_B}{T_c({\mu }_B)}\right)}^2+...\right],$$

$$\int d^{3}k⟶\int k_T\sqrt{k_T^2+m_i^2}\cosh yd{k}_{T}dyd\varphi ,$$

Full text

• May 2019

Delineating the phase diagram of quantum chromodynamics (QCD) at finite temperature $T$ and baryon-chemical potential $\mu_{B}$ is a challenging problem in theoretical and experimental studies [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. The Beam Energy Scan (BES) at the Relativistic Heavy Ion Collider (RHIC) [12, 14] has been dedicated to the search for the conjectured QCD critical end point (CP) through the systematic studies of various observables. In particular, a non-monotonic behavior of the fluctuations of conserved charges with beam energy ($\sqrt{s}$) is considered as a conceivable experimental signature of the chiral critical behavior and of the CP in heavy-ion collisions [11, 15, 16, 17, 18, 19, 20]. Such a behavior is typically associated with the divergence of the correlation length and the fluctuations of the critical mode ($\sigma$) at the CP [15, 16, 21, 22, 23, 24, 25]. Another, indirect way of verifying the existence of the CP is given by searching for non-uniform structures in multiplicity distributions due to domain formations in the region of the first-order phase transition adjacent to the critical point [26, 27, 28].

Measurements of event-by-event fluctuations in the net-proton number [29, 30, 14, 32], as a proxy for net-baryon number, the net-electric charge [33] and the net-kaon number [34], as a proxy for net-strangeness, have been performed in heavy-ion collisions at RHIC and LHC energies. While non-monotonic structures in the still preliminary STAR Collaboration data on higher-order net-proton fluctuations were indeed observed [30, 31, 32], their unambiguous interpretation as the consequence of the presence of the CP has not been achieved yet [19, 20, 35, 36].

Different QCD-like effective models indicate, that the chiral CP belongs to the Z(2) static universality class of the $3$-dimensional Ising model, and that the chiral critical mode can be associated with the order parameter, the chiral condensate [23, 24, 25, 37]. For a static and infinite system in thermodynamic equilibrium, the scaling of the fluctuations of $\sigma$ modes with the diverging correlation length at the critical point is governed by the Z(2) critical exponents [38]. However, a medium created in a heavy-ion collision is neither static nor infinite and due to its expansion dynamics, the non-equilibrium effects can play an important role. The temporal growth of the correlation length is already dynamically limited by the phenomenon of critical slowing down [39, 40]. In addition, non-equilibrium effects can lead to the retardation and damping of the critical signals [41, 42, 43, 44]. Furthermore, even in thermodynamic equilibrium, the exact charge conservation [45, 46, 47], volume fluctuations [48, 49, 46, 50], and further sources of non-critical fluctuations in the data [51], as well as, late hadronic stage processes [52, 53, 54, 55, 56], can modify signals of the critical fluctuations.

Given the above mentioned challenges it becomes clear that the theoretical description and interpretation of data on fluctuation observables require special care and eventually a dynamical framework to match the conditions expected in heavy-ion collisions. Nevertheless, numerical simulations of the critical dynamics require also input from static, equilibrium models to provide an analytic benchmark [44, 50, 57, 58]. One of such phenomenological models, that accounts for the critical fluctuations in heavy-ion collisions has been formulated in Ref. [56]. There, the fluctuations of the critical mode were linked to (anti-)protons and resonances by allowing the particle masses to fluctuate as a consequence of the (partial) mass generation through the coupling to $\sigma$. As was demonstrated in Ref. [56], the non-monotonicity observed in the higher-order fluctuation data of STAR Collaboration [30, 31, 32] can be qualitatively described within such formulation. However, in contrast to the experimental data a pronounced peak structure was also found in the lowest-order cumulant ratios of the net-proton fluctuations even for small values of couplings [56].

An additional complication arises from the fact that the relevant quark masses $m_{q}$ are finite but small. In the limit of vanishing $m_{q}$, there is a second-order phase transition for small and even zero $\mu_{B}$ belonging to the O(4) static universality class. This line of second-order phase transitions terminates in a tricritical point [59]. Recent lattice QCD studies [60, 61] and measurements of higher-order cumulants of net-particle multiplicity distributions [62] have revealed a significant impact of the hidden O(4) criticality at larger $\sqrt{s}$ even in the case of explicit chiral symmetry breaking. Moreover, in [59] it was shown that for reasonable light quark masses the physics near CP is strongly affected by the presence of the tricritical point. In fact, only very close to CP, for distances smaller than $1.5$ MeV in $\mu_{B}$, criticality in line with Z(2) has been observed. Outside, in a wider region around CP, traces of the hidden tricritical point become visible as a difference in the scaling behavior of chiral and net-baryon number susceptibilities. From a phenomenological point of view it is, therefore, reasonable to attempt an implementation of the overlap of the two static universality classes relevant for the QCD phase transition, O(4) and Z(2).

The objective of this work is to re-examine the model assumptions introduced in Ref. [56] to improve the discrepancies between the model predictions and STAR data on the variance of net-proton number fluctuations obtained in heavy-ion collisions. The general idea is based on the universality for the critical scaling behavior of the net-quark (net-baryon) number susceptibility [59] in the vicinity of CP that is influenced by the tricritical point. We will show that the consistent implementation of critical scaling leads to a much weaker singularity than that seen for the net-proton number variance in Ref. [56]. In the mean field approximation of the chiral effective models, the singular part of the net-quark number susceptibility is proportional to the chiral susceptibility multiplied by the order parameter squared [63, 64]. This is correct both for O(4) and the tricritical point. In line with these findings, we propose the modification of the model introduced in Ref. [56] to further consistently identify the influence of criticality due to the existence of the CP on the net-proton number fluctuations in such a way that the impact of the O(4) criticality is taken into account. We will quantify different cumulant ratios of these fluctuations along the chemical freeze-out line obtained in heavy-ion collisions. We will also identify systematics of the net-proton number cumulant ratios for different locations of the CP in the $(T,\mu_{B})$-plane, relative to the freeze-out line.

The theoretical tools used in our study are introduced in Section 1. In particular in Section 1.2 we outline our procedure for capturing the overlap between the O(4) and Z(2) criticalities phenomenologically. In Section 2, we present results on the properties of the net-proton number fluctuations in the presence of the CP and for different assumptions on its location or the coupling strength between particles and the critical mode. We conclude our findings in Section 3.

1 Modeling critical fluctuations near CP

In the following, we describe the theoretical tools to be used in order to quantify the net-proton number fluctuations in the presence of the QCD critical point. We first define the baseline model which does not contain contributions from chiral critical mode fluctuations. Then we explain how the critical fluctuations can be coupled to particles near the chiral CP. Our approach is motivated by the observation that the critical contribution to the variance of the net-proton number should obey a certain scaling behavior [63, 64]. On the mean-field level, this idea can be extended to higher-order fluctuations. The critical mode fluctuations are obtained by using universality class arguments between QCD and the $3$-dimensional Ising spin model. The necessary mapping between the corresponding variables is also discussed.

1.1 Baseline model

As a baseline model for calculating the net-proton number cumulants we employ the hadron resonance gas (HRG). In this model, the pressure of the interacting hadron gas is approximated by the sum of the partial pressures of non-interacting hadrons and their resonances [9, 65]. In the HRG model, the particle density of each particle species is given by the ideal gas expression

[TABLE]

where $d_{i}$ is the degeneracy factor and

[TABLE]

is the thermal equilibrium distribution function. In Eq. (2), $E_{i}=\sqrt{k^{2}+m_{i}^{2}}$ and $\mu_{i}=B_{i}\mu_{B}+S_{i}\mu_{S}+Q_{i}\mu_{Q}$ are the energy and chemical potentials of a particle with mass $m_{i}$, baryon number $B_{i}$, strangeness $S_{i}$ and electric charge $Q_{i}$, and $\mu_{X}$ is the conjugate to the conserved charge number $N_{X}$. For a grand canonical ensemble, the average number in a constant volume $V$ is $\langle N_{i}\rangle=Vn_{i}$.

In the thermal medium the particle number fluctuates around its mean on an event-by-event basis. These fluctuations can be quantified in terms of cumulants, where the n-th order cumulant is defined as

[TABLE]

at constant temperature $T$. In the following, we consider the first four cumulants of the net-proton number $N_{p-\bar{p}}=N_{p}-N_{\bar{p}}$, which are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

The second equalities in Eqs. (5)-(7) hold only when correlations between different particle species vanish which is the case in the HRG baseline.

In the following, we will ignore the contributions stemming from resonance decays and concentrate solely on the primary proton and anti-proton numbers and their fluctuations. Since the cumulants are volume-dependent, it is useful to consider their ratios in which, ignoring volume fluctuations, the volume dependence cancels out. Consequently,

[TABLE]

where $M=C_{1}$ is the mean, $\sigma^{2}=C_{2}$ the variance, $\kappa=C_{4}/C_{2}^{2}$ the kurtosis, and $S=C_{3}/C_{2}^{3/2}$ the skewness. Under the assumption that experimentally measured event-by-event multiplicity fluctuations originate from a thermal source with given $T$ and $\mu_{X}$, one can compare the model results for the cumulant ratios with the experimental data to deduce features of the QCD phase diagram.

1.2 Coupling to critical mode fluctuations

Fluctuations of the chiral critical mode $\sigma$ near the QCD critical point are expected to affect various experimentally measurable quantities, including the net-proton number cumulants [15, 16, 66]. Currently, there is no general prescription of how to model the effect of critical fluctuations on observables. In [56], inspired by the way in which the particle mass is generated near the chiral transition in the sigma models, the critical mode fluctuations were incorporated into the HRG model by allowing the particle mass to fluctuate on an event-by-event basis around its mean value. Consequently, this leads to fluctuations in the distribution function, $f_{i}=f_{i}^{0}+\delta f_{i}$. Here, the change of the distribution function due to critical mode fluctuations reads

[TABLE]

where $g_{i}$ is the coupling strength between $\sigma$ and the particle of type $i$, which in principle can also depend on $T$ and $\mu_{X}$, $v_{i}^{2}=f_{i}^{0}((-1)^{B_{i}}f_{i}^{0}+1)$ and $\gamma_{i}=E_{i}/m_{i}$.

Due to the above modification of the distribution function, proton and anti-proton number fluctuations are no longer independent. Considering only the most singular contributions to the fluctuations (see [67] for a discussion of the impact of less critical contributions), one obtains the following expressions for the net-proton number cumulants influenced by the conjectured critical point,

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

and $\langle(V\delta\sigma)^{n}\rangle_{c}$ are the critical mode cumulants.

To calculate the cumulants of the critical mode, we apply universality class arguments which state, that close to the critical point different physical systems belonging to the same universality class, exhibit the same critical behavior characterized by the corresponding critical exponents [68]. Assuming that the QCD critical end point belongs to the same universality class as the $3$-dimensional Ising model, we identify the order parameter in QCD, $\sigma$, with the magnetization $M_{I}$, i.e. the order parameter of the Ising spin model. This allows us to define the critical mode cumulants as

[TABLE]

where $r=(T-T_{c})/T_{c}$ and $h=H/H_{0}$ are the reduced temperature and magnetic field in the spin model, respectively, and the critical point is located at $r=h=0$.

In our approach, the second cumulant $C_{2}$ of the net-proton number receives critical contributions through the coupling of (anti-)protons to the critical mode via the first derivative of the magnetization $\partial M_{I}/\partial h$. In the spin model, this quantity is related to the magnetic susceptibility, and because of universality, to the chiral susceptibility ($\chi$) in QCD. For Z(2), the chiral susceptibility is known to diverge with a stronger amplitude but the same critical exponent as the net-baryon number susceptibility $\chi_{B}$. This is, however, correct only extremely close to the CP [59]. Further away, in a larger region around CP, traces of the hidden tricritical point lead to different critical exponents, making $\chi$ to diverge stronger than $\chi_{B}$ [59, 63, 64]. Therefore, the model in its form introduced in Ref. [56] requires some modifications. These can be accomplished by using the following chiral relation between the net-baryon number and the chiral susceptibilities, which was found in the effective model calculations [63, 64] for O(4) and the tricritical point,

[TABLE]

Here $\chi_{B}^{\rm{reg}}$ is the regular part of the net-baryon number susceptibility, whereas the second term in the above equation constitutes the singular part, $\chi_{B}^{\rm{sing}}$. Although this relation holds on the mean-field level, we will still use it in the present model to capture the correct dominant scaling behavior near but not exactly at the QCD critical end point.

Following Eq. (12), the singular part of the second-order cumulant is written, as

[TABLE]

where

[TABLE]

In the linear-$\sigma$ models, the proton mass is related to the chiral condensate. By replacing the factor $m_{p}$ in Eq. (18) by $g_{p}\sigma$, the critical contribution to the second-order cumulant in Eq. (12) has the same form as the one in Eq. (17). We note that this replacement is only done at leading-order in the derivation, i.e. factors of $E_{i}$ remain unaffected. This phenomenological procedure results in the following modification

[TABLE]

which is applied also to the higher-order cumulants of the net-proton number. Consequently, we arrive at

[TABLE]

[TABLE]

and

[TABLE]

By evaluating $\sigma$ and its cumulants via Eq. (16) from applying the scaling equation of state of the $3$-dimensional Ising model, see Section 1.3, we aim at capturing the crossover between the O(4) and Z(2) universality classes near the CP. As a result, $C_{2}$ now scales as $R^{\beta(3-\delta)}$ instead of $R^{\beta(1-\delta)}$ with the distance $R$ from the CP in the region that is dominated by the influence of the tricritical point. It should only approach the scaling $\propto R^{\beta(1-\delta)}$ extremely close to the critical end point. Nonetheless, we note that we will never consider a situation in which we are substantially close enough to see the pure Z(2)-criticality.

In QCD with finite quark masses the proper order parameter of Z(2) is a linear combination of scalar density, quark number density and energy density. Along this direction the grand potential exhibits zero curvature at the Z(2) CP [69]. The soft mode emerges as a space-like collective excitation in the scalar spectral function [70]. A massive $\sigma$-mode appears in the time-like sector and decouples from the critical fluctuations near the CP. In our study, the $\sigma$ in Eq. (17) must be interpreted as an effective soft mode for the overlap between the O(4) and Z(2) criticalities even in the presence of finite quark masses. In the following, we use $g$ instead of $g_{p}$ and apply Eqs. (4) and (21)-(23) together with Eq. (19) to calculate the net-proton number cumulant ratios.

1.3 Magnetic equation of state

For the magnetic equation of state we use a parametric representation,  [71]

[TABLE]

which is strictly speaking valid only in the scaling region close to the critical point. Here $R$ and $\theta$ are auxiliary variables which depend on $r$ and $h$, from which $R$ measures the distance from the critical point. They are determined by solving the following equations,

[TABLE]

where $\beta$ and $\delta$ are critical exponents and

[TABLE]

is an odd polynomial in $\theta$. The parameters entering Eqs. (24)-(27) are determined numerically by the Monte-Carlo simulations or by other theoretical tools such as the $\epsilon$-expansion or functional renormalization group methods. In [71], renormalization group and field theoretical methods were used to determine the critical exponents and coefficients of $w(\theta)$, reading, $\beta=0.325$ and $\delta=4.8169$, and $a=-0.76145$, $b=0.00773$ and $c=1$. For the normalization constants $H_{0}$ and $M_{0}$ in Eqs. (16) and (24) we follow [56] and set exemplarily, $M_{0}=5.52\times 10^{-2}\,\rm{GeV}$ and $H_{0}=3.44\times 10^{-4}\,\rm{GeV}^{3}$. Differentiation of Eqs. (24)-(27) with respect to $h$ as defined in Eq. (16) allows us to determine the cumulants of the critical mode. These expressions are summarized in the appendix of Ref. [56], and will not be repeated here.

The parametrization of the magnetic equation of state used in this work provides an accurate description of the order parameter close to the critical point. In a simpler, linear parametric representation the coefficients and critical exponents read, $a=-2/3$, $b=0$ and $c=3$, $\beta=1/3$ and $\delta=5$. In this representation, the cumulants of the critical mode are considerably simpler and can be found in [43, 72]. We note, that the numerical results presented in [56] were obtained using this simpler representation. Qualitatively, features of the results obtained in either of these parametrizations are very similar, however the non-monotonic structures due to the critical point are quantitatively more pronounced in the parametrization which we employ in this study.

1.4 Mapping between spin model and QCD

In order to utilize the universality class argument, one needs a mapping between the reduced temperature $r$ and magnetic field $h$ in the spin model and the QCD temperature $T$ and baryon-chemical potential $\mu_{B}$. Such a mapping is non-universal. Moreover, it is sensitive to the model assumptions. One of the frequently used models is a linear mapping between the spin model and QCD phase diagrams [43, 73]. There, the main assumptions in the mapping are: (i) The conjectured QCD critical point at $(\mu_{cp},T_{cp})$ is located at $r=h=0$ in the Ising model coordinate system, and (ii) the $r$ axis is tangential at the critical point to the first-order phase transition line in QCD, where the positive $r$ direction points towards the QCD crossover region. The orientation of the $h$ axis is not well constrained. In this work we assume that this axis is perpendicular to the $r$ axis and its positive direction points towards the hadronic phase of QCD (see Fig. 1).

To obtain $(r,h)$ corresponding to a given $(\mu_{B},T)$ pair it is convenient to introduce an auxiliary coordinate system $(\tilde{r},\tilde{h})$, originating at the QCD critical point and oriented such that the $\tilde{r}$ axis is parallel to the $\mu_{B}$ axis. Then the mapping is defined as

[TABLE]

and

[TABLE]

where $\Delta T_{cp}$ and $\Delta\mu_{cp}$ are parameters which are connected to the size of the critical region. Following [56] we set $\Delta T_{cp}=0.02\,$GeV and $\Delta\mu_{cp}=0.42\,$GeV. The corresponding point in the spin model coordinate system is obtained by rotation of the auxiliary coordinate system, where the angle is determined by the slope of the first-order phase transition line of QCD at the critical point.

The exact location of the QCD critical point and the slope of the first-order phase transition line are not known. Input provided by lattice QCD calculations may be used to constrain these parameters. The crossover line can be parametrized as

[TABLE]

where $T_{c,0}=(0.145\,\dots\,0.163)\,$GeV is the critical temperature at vanishing chemical potential [3, 4] and $\kappa_{c}\,\simeq\,0.007\,\dots\,0.059$ is the chiral crossover curvature [74, 75, 76, 77]. Then, for given $\mu_{cp}$, $T_{c,0}$ and $\kappa_{c}$ the temperature $T_{cp}$ and the slope of the first-order phase transition line at $(\mu_{cp},T_{cp})$ are obtained from Eq. (30).

Finally, to make contact between our model calculations and the experimental data on net-proton fluctuations we calculate the net-proton number cumulants at chemical freeze-out. In this work we use the chemical freeze-out conditions which were determined by analyzing the measured hadron yields [78, 79, 80, 81, 82, 83]. The solid blue line in Fig. 1 indicates the recent parametrization from Ref. [9].

2 Numerical results on the net-proton number fluctuations

In the following we discuss the numerical results on different fluctuations of the net-proton number obtained in the above phenomenological approach that correctly embeds the expected scaling behavior of the net-proton variance [63, 64]. To link our results to STAR data we first discuss how the experimentally applied kinematic acceptance cuts can be included into the framework. Then, the impact of our phenomenological modification introduced in Eq. (20) on fluctuation observables is studied. This includes also a discussion of the effects of modifying the coupling strength $g$ and the proximity of the thermal conditions at which the cumulant ratios are evaluated near the QCD critical point.

2.1 Kinematic cuts

In the fluctuation measurements, the investigated phase-space coverage is limited by the detector design and specific demands from the experimental analysis, like e.g. to optimize efficiency. Since the observables depend on the implemented kinematic acceptance, (see e.g. refs. [29] and [30, 31, 32, 84]) it is important to incorporate the experimental cuts into our theoretical framework. Following [85], we include restrictions in kinematic rapidity $y$, transverse momentum $k_{T}$ and azimuthal angle $\phi$ by replacing

[TABLE]

and $E_{i}\to\sqrt{k_{T}^{2}+m_{i}^{2}}\cosh y$ in the momentum-integrals. In line with [32], we consider the following phase-space integrations: $-0.5\leq y\leq 0.5$, $0\leq\phi\leq 2\pi$ and $0.4\,\rm{GeV/c}\leq k_{T}\leq 2\,\rm{GeV/c}$. We note, that this procedure cannot account for scattering of particles in and out of the acceptance window during the late stage evolution of a medium created in heavy-ion collisions.

2.2 Net-proton number cumulant ratios

The contributions of critical fluctuations to the net-proton number cumulants are sensitive to the value of the coupling $g$ between (anti-)protons and the critical mode. This value may, in principle, depend on $T$ and $\mu_{X}$, i.e. on the position in the QCD phase diagram where the cumulant ratios are evaluated. In fact, in the quark-meson [86] and NJL [87] models, the meson-nucleon couplings are found to decrease both with increasing $T$ and/or $\mu_{B}$. In the following, we consider fixed values for $g$ along the chemical freeze-out curve depicted in Fig. 1, i.e. independent of the beam energy $\sqrt{s}$. Typical values for $g$ may be inferred from various effective model calculations. In the linear sigma model this parameter can be related in the ground state to the pion decay constant, $g\simeq m_{p}/f_{\pi}\simeq 10$ [38, 66]. Similar values can be found in non-linear chiral models [88] describing QCD matter in neutron stars. On the other hand, based on different quark-meson models the value of $g\simeq 3-7$ for the nucleon-meson couplings is well conceivable [89]. To highlight the features of our model results we will use $g$ in the range between $3$ and $5$.

In the previous work [56], the critical point was exemplarily located at $\mu_{cp}=0.39$ GeV and $T_{cp}=0.149$ GeV (see Fig. 1). There, even for rather small values of $g\simeq 3$, the cumulant ratio $C_{2}/C_{1}$ exhibited a clear peak structure compared to the non-critical baseline and in contrast to the STAR data [29, 30, 31, 32]. This is illustrated in Fig. 2.

The effect of the modified scaling, discussed in Sec. 1.2, is the substantial reduction of the critical contribution to $C_{2}$ in Eq. (21) implying the disappearance of the maximum in $C_{2}/C_{1}$ even for large values of $g\simeq 5$, as seen in Fig. 2. Within error bars, the model results are in agreement with data of the STAR Collaboration [32] on the $C_{2}/C_{1}$ ratio. We note, that in order to see a similar maximum in $C_{2}/C_{1}$ in this model calculations, as found in [56] for a given location of the critical point, a significantly larger value of $g$, outside the expected range discussed above, would be necessary. Furthermore, the differences between the $C_{2}/C_{1}$ ratios shown in Fig. 2 are independent of the particular choice for the phenomenological freeze-out conditions. Indeed, in these studies we have adopted the freeze-out line from Ref. [9], whereas in [56] the $C_{2}/C_{1}$ ratio was calculated along the freeze-out line from Refs. [90, 91].

The higher-order cumulant ratios of the net-proton number fluctuations are also shown in Fig. 2. The model introduced in [56] exhibits clearly pronounced non-monotonic structures of higher-order cumulant ratios for $g=3$ and $5$. We note, that while in the model [56] the ratios $C_{2}/C_{1}$ and $C_{4}/C_{2}$ are independent of the choice made for the orientation of the $h$-axis, the behavior of $C_{3}/C_{2}$ is sensitive to this choice. In [56], the positive direction of the $h$-axis was defined to point upward towards larger values of $T$. In our work, we choose the opposite direction as discussed in Sec. 1.4 such that $g\sigma$ is positive along the chemical freeze-out curve.

As seen in Fig. 2, in the present model calculations, the non-monotonic structures of the higher-order cumulant ratios become strongly suppressed even for $g=5$. In fact, for the considered setup, the model results show rather small deviations from the non-critical baseline. Moreover, in contrast to [56], the behavior of $C_{3}/C_{2}$ in the present model (blue lines in Fig. 2) does not depend on the orientation of $h$. This is because the combined $\theta$-dependence in $C_{3}$ of the critical mode fluctuations $\langle(V\delta\sigma)^{3}\rangle$ and the factor $\sigma^{3}$ is even, see Eq. (22). Thus, although $h$ is an odd function in $\theta$ as seen from Eqs. (26) and (27), a re-orientation of the $h$-axis would have no effect.

The substantial reduction of the critical signal in the net-proton number cumulant ratios seen in our model results is a consequence of (i) the reduced scaling of the critical contributions to $C_{n=2,3,4}$ in Eqs. (21)-(23) and (ii) the magnitude of the factors $(g\sigma)^{n}$. As a result of the phenomenological implementation in Eq. (20), the scaling of $C_{n}$ at the CP is weakened by an additional $n\beta$ factor in the critical exponents. Moreover, the factor $g\sigma$ differs from the vacuum proton mass $m_{p}$ in Eq. (18) as employed in [56]. In the present calculations, $g\sigma$ is of the order of $0.2-0.3$ GeV for most $\sqrt{s}$. Its actual values depend on the parameters in the magnetic equation of state, most notably on the value of $M_{0}$ (see Sec. 1.3), and the mapping between spin model and QCD (see Sec. 1.4). These values receive support from recent works on the origin of the baryon masses in both lattice QCD [92, 93] and effective models based on parity doubling [94, 95]. There, the baryon masses are found to be given to a large extent by $\sigma$-independent contributions.

The properties of different cumulant ratios shown in Fig. 2 can lead to the conclusion, that the monotonic beam energy dependence seen experimentally in $C_{2}/C_{1}$ together with the non-monotonicity in the higher-order cumulant ratios cannot be explained simultaneously by a model that (i) includes critical mode fluctuations through the coupling of $\sigma$ with the particles and (ii) obeys the connection between the chiral and net-baryon number susceptibilities observed in effective models [63, 64]. However, one notes that the model results depend not only on the values of the coupling $g$ but also on the non-universal details of the mapping between QCD and the spin model discussed in Sec. 1.4. One of them is the unknown distance of the QCD critical point from the chemical freeze-out conditions at which the fluctuations are determined. To study this effect we keep the chemical freeze-out conditions fixed but vary the location of the critical point in the QCD phase diagram as summarized in Tab. 1 and depicted in Fig. 3.

In Fig. 4 we show the net-proton number fluctuations along the phenomenological freeze-out line at fixed value of the coupling $g=5$ and assuming different locations of the CP. As evident from this figure, moving the critical point closer to the chemical freeze-out curve leads to an increase of non-monotonic structures in the net-proton number cumulant ratios. While deviations from the non-critical baseline (black dotted lines) remain moderately weak in $C_{2}/C_{1}$, they become more pronounced with increasing order of the fluctuations in $C_{n}/C_{2}$ ratios. We note that even for CP3 we remain always at a distance of about or larger than $3$ MeV in $\mu_{B}$ along the chemical freeze-out line.

In Fig. 5 we show the influence of the critical point located at the closest distance CP3 on the energy dependence of the net-proton number fluctuations for different values of the coupling $g$. As can be seen, all cumulant ratios depend strongly on the actual value of $g$. This behavior is expected from Eqs. (21)-(23) and Eq. (19) since cumulants $C_{n}$ scale as $g^{2n}$. This differs from the model introduced in [56] where, as seen from Eqs. (12)-(15), they scale only as $g^{n}$. When compared to the STAR data, the results for $C_{2}/C_{1}$ and $C_{4}/C_{2}$ shown in Fig. 5 are in qualitative agreement with data, given the uncertainties in the model assumptions. However, the $C_{3}/C_{2}$ ratio also increases beyond the non-critical baseline towards the lower beam energies in contrast to the STAR data.

From the results shown in Figs. 2, 4 and 5 it is clear that for small couplings $g\simeq 3$ deviations from the non-critical baseline are negligible in all cumulant ratios irrespective of the studied location of the critical point. By increasing $g$, non-monotonic structures in the $\sqrt{s}$-dependence of the net-proton cumulant ratios develop and are sensitive to the relative distance between the critical point and the chemical freeze-out curve. This behavior is also stronger in the higher-order cumulant ratios.

It is therefore conceivable that, by an appropriate choice of the location of the CP and the model parameters, it is possible to describe the energy dependence of some ratios of the net-proton number fluctuations as seen in the preliminary STAR data. Consequently, the rather strong increase of the $C_{4}/C_{2}$ ratio beyond the HRG baseline and the smooth dependence of $C_{2}/C_{1}$ observed in heavy-ion collisions at energies $\sqrt{s}<20$ GeV could be due to the contribution from the CP located near the chemical freeze-out line. However, in this case the $C_{3}/C_{2}$ ratio should also exceed the non-critical baseline, which is not seen experimentally.

Thus, based on the presented equilibrium model results one concludes, that the energy dependence of the $C_{4}/C_{2}$, $C_{3}/C_{2}$ and $C_{2}/C_{1}$ ratios observed in heavy ion collisions at $\sqrt{s}<20$ GeV by the STAR Collaboration does not follow the systematics expected from the contribution of the CP to the net-proton number fluctuations alone. This conclusion is consistent with the previous analysis of different fluctuation observables based on lattice QCD and PNJL model results [19, 20]. We note, however, that the above statement requires further theoretical and empirical justifications due to current uncertainties in the model assumptions and the experimental data. It remains to be seen, for example, whether non-equilibrium effects can, in fact, push the $C_{3}/C_{2}$ ratio below the non-critical HRG model baseline. Moreover, the fireball evolution or late hadronic stage processes such as resonance decays or isospin randomization, which have not been included in our study, may influence the theoretical results.

3 Conclusions

We have studied the influence of the QCD critical end point (CP) on the properties of the n-th order cumulants ($C_{n}$) of the net-proton number and their ratios. The results were addressed in the context of the recent data from the STAR Collaboration on the energy dependence of the net-proton number fluctuations in Au-Au collisions obtained within the Beam Energy Scan program at the Relativistic Heavy Ion Collider (RHIC).

To calculate the net-proton number cumulants we have proposed a phenomenological model where non-critical fluctuations are obtained from the hadron resonance gas (HRG) statistical operator, which is known to describe data on particle yields in heavy-ion collisions and the lowest-order fluctuation observables from lattice QCD. For simplicity, the baryonic sector of the HRG was approximated by contributions from primary protons and anti-protons. To quantify the role of the chiral criticality due to the CP, the phenomenological model was introduced to describe the non-analytic part of the statistical operator in which the fluctuations of the chiral critical mode $\sigma$ are coupled to the (anti-)protons. This was achieved by linking their masses to the $\sigma$ mode, as suggested by different chiral models. Consequently, the (anti-)proton mass and its momentum distribution function fluctuate on an event-by-event basis around its mean or equilibrium value, respectively.

The critical mode fluctuations were determined by applying universality class arguments between QCD and the $3$-dimensional Ising spin model. We have extended the model introduced in Ref. [56] by accounting for the critical scaling behavior of the net-baryon variance $\chi_{B}$ suggested by effective chiral models. There $\chi_{B}$ is linked to the product of the chiral susceptibility and the chiral order parameter squared. In this way, we have phenomenologically embedded the overlap between a hidden O(4) and the Z(2) criticality, which is theoretically expected in the vicinity of CP up to a tiny region in $\mu_{B}$ [59], into the model.

We have found a substantial reduction of the critical mode contributions to the net-proton number fluctuations compared to the results in Ref. [56]. This is a consequence of the reduced critical scaling imposed by respecting the proper scaling relation between net-baryon number and chiral susceptibilities and the size of the proton mass modification due to the coupling to the $\sigma$ mode. This brings our results for different n-th order cumulants $(C_{n})$ of the net-proton number, calculated along the phenomenological chemical freeze-out line, closer to the experimental observations made by the STAR Collaboration for the energy dependence of the cumulant ratios in heavy-ion collisions at RHIC. In particular, with an appropriate choice of the model parameters and the location of the CP relative to the chemical freeze-out line, the model can reproduce the smooth energy dependence of $C_{2}/C_{1}$ and the increase and non-monotonic variation of $C_{4}/C_{2}$ towards lower beam energies, as is observed by the STAR Collaboration. However, the decrease of the $C_{3}/C_{2}$ ratio towards lower beam energies seen in the STAR data is inconsistent with the systematics expected in the present model from the contribution of the CP, which would predict an access of this ratio beyond the non-critical baseline. Thus, our conclusion is that it is rather unlikely that the properties observed in the low energy behavior of different ratios of the net-proton number cumulants in heavy-ion collisions are due to the existence of the critical point near the phenomenological chemical freeze-out line alone.

The work of M. Szymański and M. Bluhm is funded in parts by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska Curie grant agreement No 665778 via the Polish National Science Center (NCN) under grant Polonez UMO-2016/21/P/ST2/04035. M. Bluhm acknowledges also the partial support by the program “Etoiles montantes en Pays de la Loire 2017”. This work is partly supported by the Polish National Science Center (NCN) under Maestro grant no. DEC-2013/10/A/ST2/00106 and Opus grant no. 2018/31/B/ST2/01663. K. Redlich also acknowledges partial support of the Polish Ministry of Science and Higher Education. We thank F. Geurts and J. Thäder for providing the preliminary STAR data [32] shown in Figures 2 and 5. We also thank Bengt Friman and Marlene Nahrgang for interesting discussions and comments.