Constructing probability density function of net-proton multiplicity distributions using Pearson curve method

TL;DR

This paper constructs probability density functions of proton, anti-proton, and net-proton multiplicity distributions using the Pearson curve method, revealing potential signatures of O(4) criticality in heavy-ion collision data.

Contribution

It introduces a novel application of the Pearson curve method to model multiplicity distributions and compares them with theoretical baselines to identify critical phenomena.

Findings

Some signatures of O(4) criticality are observed.

Constructed PDFs provide a baseline for future cumulant measurements.

Comparison with Skellam and Binomial distributions highlights deviations near criticality.

Abstract

The probability density functions of proton, anti-proton, and net-proton multiplicity distributions are constructed from the Beam Energy Scan results of the STAR experiment using the Pearson curve method. The constructed distributions of proton and anti-proton are compared with Poisson and Binomial distributions. The net-proton probability distributions are compared with Skellam distributions to study the O(4) criticality near the chiral crossover transition. The $C_{6}/C_{2}$ results estimated from the obtained PDFs are compared with Skellam and Binomial baselines for the Beam Energy Scan data. The current study shows some signatures of O(4) criticality, which can be further investigated by precision measurements of the cumulants and understanding the contribution of non-critical fluctuations to them. This study also provides a baseline for the higher-order cumulant measurement in the upcoming RHIC BES II program and future LHC run.

Tables

Table 1. Table 1: Values of the different constant parameters of the shifted and rescaled Beta functions constructed for proton multiplicity distributions for different collision energies.

$\sqrt{s_{NN}}$	$a$	$b$	$\alpha$	$\beta$
11.5	-20.179	236.093	73.615	293.446
19.6	-10.128	121.403	39.351	110.805
27	-10.504	120.919	39.328	122.724
39	-10.944	129.012	39.763	145.223
54	-16.033	235.078	59.736	372.315
62.4	-16.746	374.089	65.315	667.012
200	-11.02	140.058	39.69	161.193

Table 2. Table 2: Numerical values of the different constant parameter of the shifted and rescaled Beta functions constructed for anti-proton multiplicity distributions for different collision energies.

$\sqrt{s_{NN}}$	$a$	$b$	$\alpha$	$\beta$
7.7	-0.13	3.595	0.43	3.353
11.5	-0.531	9.465	1.854	10.552
14.5	-0.816	13.248	2.864	14.181
19.6	-1.443	21.389	5.137	24.078
27	-2.138	29.546	7.646	32.338
39	-3.263	44.796	11.714	51.288
54.4	-4.597	61.003	16.481	72.318
62.4	-5.191	73.611	18.826	93.556
200	-7.926	100.635	28.494	110.911

Table 3. Table 3: Numerical values of the different constant parameter of the shifted and rescaled Beta functions constructed for net-multiplicity distributions for Au-Au collisions at s N ​ N subscript 𝑠 𝑁 𝑁 \sqrt{s_{NN}} = 11.5, 19.6 and 27 GeV.

$\sqrt{s_{NN}}$	$a$	$b$	$\alpha$	$\beta$
11.5	-31.71	398.519	112.311	667.495
19.6	-27.913	168.932	76.134	225.515
27	-40.7822	219.304	115.847	401.39

Table 4. Table 4: Values of N 2 , N 4 subscript 𝑁 2 subscript 𝑁 4 N_{2},N_{4} and N 6 subscript 𝑁 6 N_{6} for most central Au-Au collisions at s N ​ N subscript 𝑠 𝑁 𝑁 \sqrt{s_{NN}} = 14.5. 19.6, 27, 39, 54.4, 62.4 and 200 GeV.

$\sqrt{s_{NN}}$	$N_2$	$N_4$	$N_6$
14.5	48	58	67
19.6	43	53	61
27	38	47	55
39	33	42	50
54.4	31	40	48
62.4	30	39	47
200	26	36	44

Equations

$$\frac{1}{f(x)}\frac{df(x)}{dx}=-\frac{a+x}{b_0+b_{1}x+b_2{x}^2},$$

$$a$$

$$b_0$$

$$b_1$$

$$b_2$$

$$\frac{1}{f(x)}\frac{df(x)}{dx}=-\frac{a_0+a_{1}x}{b_0+b_{1}x+b_2{x}^2},$$

$$a_0$$

$$a_1$$

$$b_0$$

$$b_1$$

$$b_2$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{1197.11+195.271tan{h}^{-1}(0.217+0.012x)}}{(7363.63+36.3573x+x^2{)}^{142.433}}.$$

$$f(x)=\frac{(a-b{)}^{\alpha -1}(x-b{)}^{\zeta }}{B(\alpha -1,1-\zeta )(x-a{)}^{\zeta }},$$

$$f(x)=\frac{x^{\alpha -1}(1-x{)}^{\beta -1}}{B(\alpha ,\beta )};x\in (0,1).$$

$$B(\alpha ,\beta )=\frac{\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )},$$

$$f(x)=\frac{(x-a{)}^{\alpha -1}(b-x{)}^{\beta -1}}{B(\alpha ,\beta )(b-a{)}^{\alpha +\beta -1}};x\in [a,b].$$

$$f(x)=e^{-\mu }\frac{{\mu }^x}{x!}$$

$$f(x)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^x(1-p{)}^{n-x}$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{1186.02+180.332tan{h}^{-1}(0.184+0.012x)}}{(7445.21+31.301x+x^2{)}^{139.926}}.$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{1128.99+246.242tan{h}^{-1}(0.462+0.015x)}}{(5342.11+61.339x+x^2{)}^{146.017}}.$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{1304.06+135.835tan{h}^{-1}(0.266+0.0132x)}}{(6160.13+40.3223x+x^2{)}^{154.198}}.$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{1393.37+105.995tan{h}^{-1}(0.198+0.012x)}}{(7044.98+32.6299x+x^2{)}^{160.123}}.$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{971.675+58.9903tan{h}^{-1}(0.116+0.014x)}}{(5125.62+16.58x+x^2{)}^{115.036}}.$$

$$f(x)=\frac{1}{\sqrt{\pi }}\frac{e^{920.661+16.782tan{h}^{-1}(0.022+0.012x)}}{(6251.8+3.544x+x^2{)}^{105.648}}.$$

$$C_n(\Delta N_p)=c_n(p)+(-1{)}^n{c}_n(\bar{p}),$$

Full text

††journal: Eur. Phys. J. A

1 Introduction

One of the most intriguing subjects in the theory of Quantum Chromodynamics (QCD) is mapping its phase diagram. In this context, we have only a schematic picture of the QCD phase diagram. It was first predicted that at the limit of vanishing quark masses and baryochemical potential ($\mu_{B}=0$), the transition is likely to be second order, belonging to the ${\it O}$(4) universality class of 3-dimensional symmetric spin model Pisarski:1983ms . Current lattice QCD calculation has shown that at vanishing $\mu_{B}$, the chiral and deconfinement transition is a smooth crossover on the temperature ($T$) axis Aoki:2006we . But at non-vanishing $\mu_{B}$, the phase transition will be of first order Ejiri:2008xt ; Bowman:2008kc ; Stephanov:2007fk . Thus, various effective field theory models suggest the presence of a critical point (CP) where the second order phase transition line ends and first order phase transition line starts Asakawa:1989bq ; Hatta:2002sj ; Scavenius:2000qd ; Halasz:1998qr . Hence, locating this QCD critical point in the $T-\mu_{B}$ plane has drawn much attention in recent times. Due to the fermion sign problem at finite $\mu_{B}$, lattice QCD faces tremendous challenges to establishing the fact of the presence or absence of such CP.

Meanwhile, various theoretical works have proposed that if indeed there is a CP in the $T-\mu_{B}$ plane, it can be observed experimentally by varying the collision energy Stephanov:1999zu . Event-by-event measurement of fluctuations of conserved charges multiplicities is an excellent tool for such study Stephanov:1998dy ; Stephanov:2008qz . Then any nonmonotonic behavior of higher-order cumulants and their ratios as a function of $\sqrt{s_{NN}}$ is proposed as the signature of the presence of CP Stephanov:1999zu ; Stephanov:2008qz . Moreover, measurement of higher-order cumulants are potential probes for pseudo-critical behavior near the vicinity of chiral crossover phase transition and determining the freeze-out parameters Karsch:2010ck ; Friman:2011pf ; Bazavov:2012vg ; Borsanyi:2013hza ; Borsanyi:2014ewa . Various experiments, like Beam Energy Scan (BES) program at Relativistic Heavy-ion collision (RHIC) at BNL, the Super Proton Synchrotron (SPS) at CERN and Compressed Baryonic Matter (CBM) experiment at FAIR, GSI aim for such study.

Net-proton multiplicity fluctuations are regarded as a proxy to the net-baryon number fluctuations. The coupling strength of proton with the sigma field is more than the charged particle (pions) Stephanov:2008qz ; Hatta:2003wn . Hence, measurement of higher-order cumulants of net-proton multiplicity fluctuations is an excellent tool to search for critical phenomena. Beam Energy Scan (BES) results of cumulants ($C_{n}$) of net-proton multiplicity distributions up to fourth orders are reported by STAR experiments at RHIC Adamczyk:2013dal ; Luo:2015ewa ; STAR:2020tga ; STAR:2021iop . Here $C_{n}$ represents the $n^{th}$ order cumulant. The most recent energy dependence results of $C_{4}/C_{2}$ shows a nonmonotonic behaviour. In the energy range $\sqrt{s_{NN}}=7.7-62.4$ GeV, the most central collision data show a nonmonotonic variation of a significance of 3.1$\sigma$ with the collision energy STAR:2020tga ; STAR:2021iop . In contrast, a hadronic-transport model (UrQMDBleicher:1999xi ) and different variants of hadron resonance gas (HRG) models, which do not include CP show monotonic energy dependence. Additionally, large values of four-particle correlation function of protons are obtained for $\sqrt{s_{NN}}<$ 19.6 GeV, which could be a signature of a CP or first order phase transition Bzdak:2016jxo . However, data for $\sqrt{s_{NN}}\leq$ 27 GeV are associated with relatively large uncertainties, which will be reinvestigated by the upcoming BES-II program.

Meanwhile, it is demonstrated that the criticality belonging to the O(4) and Z(2) universality class in 3-dimensional symmetric spin model at zero and non-zero $\mu_{B}$ can be studied by knowing the probability distribution of net-baryon (net-proton) number Morita:2013tu ; Morita:2012kt ; Morita:2014fda . It is shown that the characteristic shape of the probability distribution is influenced by the nature of phase transition, and it further reveals the universality class of the transition. For example, the probability distribution appears narrower than the Skellam near the chiral crossover transition belonging to O(4) universality class Morita:2014fda . The narrowing of the probability distributions also leads to a negative value of $C_{6}$. So the O(4) criticality will be reflected both in the probability distributions and in the sixth order cumulants. In Ref. Morita:2014fda , the O(4) criticality is investigated using the efficiency uncorrected net-proton distributions measured by the STAR experiment, which shows some qualitatively similar features as the model. This work further stresses to carry out a similar study with the efficiency corrected data. However, the efficiency corrected net-proton distribution is not obtained by the experiment.

To address this demand, for the first time, the Pearson curve method (PCM) is employed for constructing the probability density function (PDF) of proton, anti-proton and net-proton multiplicity distributions from the efficiency corrected first four cumulants (moments) data. The methodology of the approximation of PDF from the first four cumulants is discussed in Section II. In Section III, this technique is used to obtain the PDF of proton, anti-proton, and net-proton multiplicity distributions for the BES data. Then from the obtained probability density functions, the ratio with respect to Skellam distribution is obtained to test the O(4) criticality. Then from the constructed PDFs the $C_{6}$ and its ratios to $C_{2}$ are estimated as a function of $\sqrt{s_{NN}}$ for most central collision data. The results are also compared with Skellam and Negative Binomial distribution expectations. The summary and outlook are discussed in Section IV.

2 Pearson curve method

A statistical distribution or frequency data is characterized by mean ($\mu$), sigma ($\sigma$), skewness ($S$), kurtosis ($\kappa$) or alternatively by first four moments (cumulants). In 1895, English mathematician Karl Pearson devised a method for deriving the functional form of a PDF from the first four moments of a given frequency data Pearson343 ; Pearson443 . He proposed a family of distributions, also known as a system of curves, which includes famous distributions, like Normal distribution, Beta, Gamma and $\chi^{2}$ distributions. According to Pearson, the PDF $f(x)$, of any of this system of curves can be derived, if it satisfies the following differential equation Pearson343 ; Pearson443 ; statbook .

[TABLE]

where $a,b_{0},b_{1},b_{2}$ are constant parameters of the distribution. These constant parameters are specific functions of the first four moments of a distribution. Generally, the parameters are expressed by central moments as follows statbook ; pearsontypes .

[TABLE]

where $\beta_{1}=m_{3}^{2}/m_{2}^{3}$ = $S^{2}$, $\beta_{2}=m_{4}/m_{2}^{2}$ = $\kappa+3$, and $m_{2},m_{3},m_{4}$ are second, third and fourth central moments, respectively. Eq. (1) yields the PDF with mean value zero. For the PDF with non$-$zero mean, Eq. (1) is modified as follows mathematica .

[TABLE]

where the modified constant parameters are given as follows.

[TABLE]

The PDF obtained using Eq.(6) has the same form as obtained from Eq.(1) except the distribution is shifted by the non$-$zero value of $\mu$. Due to the shift-invariance property of central moments and cumulants, PDF obtained from Eq.(1) and Eq.(6) will yield the same values of central moments and cumulants from second order onwards. Therefore, the PDFs obtained from any of these two can be used conveniently. However, for consistency, PDFs obtained from Eq.(6) will be used for this study.

To determine the PDF of an experimental frequency data, the first step is to calculate, $\mu$, $m_{2},m_{3},m_{4}$ or first four cumulants. Once the parameters are estimated using Eq(7)-(11), the solution of Eq.(6) will yield the PDF. To be noted that the PCM is applicable only when $\beta_{2}-\beta_{1}-1>0$. There are basically seven types of Pearson family of curves exist. Deciding the type of curves depends on the types of roots of the quadratic in the denominator of the differential equations given in Eq.(1) and Eq.(6). Some rigorous examples of deriving the PDF from frequency data using PCM can be found in Ref.Pearson343 ; Pearson443 ; statbook . These families of curves are related to several other known distributions. For example, Beta, Wigner semicircle, Normal, Student’s t, Levy, Fisher-Snedecor, Cauchy distributions are special cases of type 1, 2, 3, 4, 5, 6 and 7 Pearson family of curves, respectively. Nevertheless, the derivation of PDFs using PCM can be done in a straightforward way using modern mathematical tools, like $Mathematica$ mathematica , which is used for this work.

PCM is a very powerful tool to derive the PDF from a set of frequency data. It also eliminates the arbitrariness of using different PDF to the same frequency data. PCM is widely used in Economics, Bio-Science, Engineering for fitting frequency data and to approximate the PDF. A probability distribution can’t describe uniquely by the first four cumulants (moments). PCM can be used to model other probability distributions other than its family of curves. It can exactly reproduce the first four cumulants. However, it does not guarantee to reproduce the higher-order cumulants for $n>4$. In reality, the knowledge about the origin and the exact probabilistic nature of the data is not known as ’a priori’. So far the detector acceptance and efficiency corrected net-proton multiplicity distribution is not obtained. Meanwhile, none of the statistical models (Skellam and Negative Binomial distributions) are successful in describing the experimental results of the individual particle and net-proton multiplicity distributions. Hence, PCM provides an alternative scope for modeling the multiplicity distributions and estimating the higher order cumulants for the baseline study.

3 Results and discussion

The analysis is carried out by using the beam energy dependence results of first four cumulants of proton (p), anti-proton($\bar{\textrm{p}}$) and net-proton ($\Delta\textrm{p}$) multiplicity fluctuations measured at midrapidity $|y|<0.5$, for the transverse momentum ($p_{T}$) range $0.4<p_{T}<2.0$ GeV/$c$ by STAR experiment at RHIC STAR:2021iop . In this article, the analysis is focused only on most central collision events (0-5$\%$ centrality).

First, $\mu$, $m_{2},m_{3}$, and $m_{4}$ are estimated from the detector efficiency corrected results of the first four cumulants of p, $\bar{\textrm{p}}$ and $\Delta\textrm{p}$ multiplicity distributions. The constant parameters of Eq.(6) are estimated by invoking Eq.(7)-(11). Then the functional form of the PDFs are derived using PCM.

3.1 PDF of proton multiplicity distributions

For Au-Au collisions at $\sqrt{s_{NN}}$ = 7.7, the PDF $f(x)$, of proton multiplicity distributions is found to be type 4 Pearson family of curve, which is expressed in Eq. (12).

[TABLE]

Here $x$ represents the proton multiplicity. Type 4 Pearson family of curves is not a standard distribution. For $\sqrt{s_{NN}}$ = 14.5 GeV, the proton multiplicity distribution is constructed as type 6 Pearson family of curves. The PDF is a shifted and rescaled Fisher-Snedecor distribution. It can be expressed as follows.

[TABLE]

for all $x>a$. Here $\zeta=(\beta+\alpha a)/(a-b)$, with $a=-42.6642,b=-302.706,\alpha=920.866$ and $\beta=-24651.3$.

For $\sqrt{s_{NN}}$ = 11.5, 19.6, 27, 39, 54.4, 62.4 and 200 GeV, the proton multiplicity distribution are constructed as type 1 Pearson family of curve, which are shifted and rescaled (SR) Beta distribution. The Beta distribution is expressed as follows.

[TABLE]

Here the normalisation constant $B(\alpha,\beta)$ is Beta function, which is given as,

[TABLE]

where $\alpha$ and $\beta$ are shape parameters. The Beta distribution is the conjugate prior probability distribution for the Bernoulli, (Negative)Binomial Distributions (NBD/BD) and Geometric distributions. The SR Beta distribution is obtained by performing a transformation, $x\rightarrow a+(b-a)x$, which can be written as,

[TABLE]

In Eq.(16), $a$ and $b$ are constant parameters. The values of these constant parameters obtained for $\sqrt{s_{NN}}$ = 11.5, 19.6, 27, 39, 54.4, 62.4 and 200 GeV are given in Table 1.

For illustration purpose, the probability distribution of proton multiplicity is plotted only for most central Au-Au collisions at $\sqrt{s_{NN}}$ = 200 GeV, which is shown in Figure 1. The distribution is also compared with Poisson and BD. A Poisson distribution is a single parameter function of $\mu$, which is expressed as follows.

[TABLE]

For Poisson distribution, the variance, $C_{2}$ is equal to $\mu$. The probability function of a BD distribution of a random variable $x$ is expressed in terms of the number of trials ($n$) and probability of success ($p$) as follows.

[TABLE]

For BD, $\mu=np$ and $C_{2}=np(1-p)$. The ratios of SR Beta distributions with respect to Poisson and BD are shown in the lower panels of Figure 1. From the ratio, it can be seen that the SR Beta distribution obtained for proton multiplicity is wider than the Poisson distribution at very small multiplicity ($N_{p}<$ 10). In the mid-range of multiplicity all the distribution seem indistinguishable. However, for large multiplicity $N_{p}>$ 27, narrowing of the distribution with respect to Poisson distribution occurs. From the ratio it is clear that the SR Beta distribution for proton is closer to BD than the Poisson distributions.

3.2 PDF of anti-proton multiplicity distributions

The PDF of anti-proton multiplicity distributions are constructed using PCM and found to be shifted and rescaled Beta distribution for $\sqrt{s_{NN}}$ = 7.7 - 200 GeV as given in Eq.(16). The numerical values of the constant parameters of the SR Beta distribution for anti-proton multiplicity distributions for different collision energy are given in Table 2.

The probability distribution (SR Beta distribution) of anti-proton obtained from PCM for Au-Au collisions at $\sqrt{s_{NN}}$ = 200 GeV is compared with Poisson and BD distribution, which are illustrated in the upper panel of Figure 2. The difference of SR Beta distribution with respect to Poisson and BD is estimated by taking the ratios and shown in the lower panel of Figure 2. From the ratio, it is observed that the constructed anti-proton distribution has a wider tail than the Poisson and BD at lower multiplicity range ($N_{\bar{p}}<$ 7). With increasing the multiplicity all the distributions almost merge. However, at large value of anti-proton multiplicity ($N_{\bar{p}}>$ 23), it gets narrower with respect to the Poisson distribution. At $\sqrt{s_{NN}}$ = 200 GeV, like the proton multiplicity distribution, the SR Beta distribution for anti-proton is closer to BD than the Poisson distributions.

3.3 PDF of net-proton proton multiplicity distributions

The PDF of net-proton multiplicity distributions in Au-Au collisions at $\sqrt{s_{NN}}$ = 7.7, 14.5, 39, 54.4, 62.4 and 200 GeV are constructed as type 4 Pearson family of curve, which are given in Eq.(19), (20), (21), (22), (23) and (24), respectively.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In these equations, $x$ represents net-proton multiplicity. The PDF of net-proton multiplicity distributions for Au-Au collisions at $\sqrt{s_{NN}}$ = 11.5, 19.6 and 27 GeV are obtained as SR Beta distributions. The constant parameters of the SR Beta functions for the two energies are listed in Table 3 .

The probability distributions constructed for the net-proton multiplicity distributions in the most central Au-Au collisions at $\sqrt{s_{NN}}$ = 7.7 - 200 GeV are illustrated in Figure 3. It can be seen that the net-proton distribution becomes symmetric around zero while going from 7.7 GeV to 200 GeV.

3.4 Test for the criticality of the net-proton multiplicity distributions

The probability distribution of a conserved charge is the fundamental of its cumulants. So the signature of criticality near the QCD phase boundary will be reflected both in its higher order cumulants and in the probability distributions. To demonstrate the Z(2) and O(4) criticality, first the net-quark probability distributions and later the net-baryon probability distributions at zero and non-zero $\mu_{B}$ are investigated Morita:2013tu ; Morita:2012kt ; Morita:2014fda . By comparing the probability distributions with non-singular distribution (Skellam function), it is shown that the deviation with respect to the Skellam function can be used to study the criticality. Using Landau theory of phase transition it is shown that for O(4) universality the probability distribution is characterised by narrower for intermediate values of $|\Delta N_{p}|$ and a wider tail at large value of $|\Delta N_{p}|$ with respect to the Skellam distribution Morita:2012kt . The study done in Ref. Morita:2013tu using the Functional Renormalization Group approach to the Quak-Meson (QM) model shows that the narrowing of the probability distribution relative to the Skellam function can be regarded as a signature of O(4) criticality near the chiral crossover transition. As a consequence, the sixth order cumulant of net-baryon will be negative. Furthermore, it shows that the narrowing is due to the effect Fermi statistics 111The Skellam distribution used in Ref. Morita:2012kt is different than the one used in Ref. Morita:2014fda .

Higher order cumulants are sensitive to the tail of the probability distribution. Hence, the knowledge of the distributions for a sufficiently large value of $|\Delta N_{p}|$ is required for the exact determination of $C_{6}$ and its characteristics properties. The value of $|\Delta N_{p}|$ increases with increasing the order of cumulant. In Ref. Morita:2014fda , the probability distributions ($P(N)^{FRG}$) of net-baryon is obtained from the FRG approach within the QM model at zero and non-zero $\mu_{B}$. Then the ratio of the $P(N)^{FRG}$ and the Skellam distribution ($P(N)^{S}$) is computed as a function of shifted mean ($\delta N_{p}=\Delta N_{p}-\mu$) normalised to $N_{6}$. Here, $N_{6}$ is the maximum value of $|\Delta N_{p}|$ where the $C_{6}$ is exactly reproduced from the Skellam distribution to its exact value Morita:2013tu . In this case, the Skellam distribution has the same mean and variance as the distribution it is compared with. Then the efficiency uncorrected net-proton probability distributions reported by the STAR experiment in the kinematic range $0.4<p_{T}<0.8$ GeV/$c$ and $|y|<0.5$, are used to investigate any possible O(4) criticality. The study shows some qualitatively similar structure as predicted in the model Adamczyk:2013dal ; Morita:2014fda .

Here, an attempt is made to carry out similar exercise to investigate any possible O(4) criticality associated with the constructed probability distributions ($f(\Delta N_{p})^{PCM})$ for the efficiency corrected cumulants results using PCM in most central Au-Au collisions at $\sqrt{s_{NN}}$ = 14.5, 19.6, 27, 39, 54.4, 62.4 and 200 GeV in the kinematic range $0.4<p_{T}<2.0$ GeV/$c$ and $|y|<0.5$. The functional form of the constructed probability distributions are already given in the previous section. First, the maximum value ($N_{max}$) for $2^{nd},4^{th}$ and $6^{th}$ order of cumulants for the Skellam distributions for different collision energies are estimated. Here the Skellam distributions are constructed with same $C_{1}$ and $C_{2}$ value as the net-proton data. For a given order of cumulant, the $N_{max}$ is determined from the value of $|\Delta N_{p}|$ at which the cumulant value converges with its exact values. The convergence is checked by taking the ratio of cumulant estimated from the distribution to its expected value. For illustration purpose the ratio at different value of $|\Delta N_{p}|$ for Au-Au collisions at $\sqrt{s_{NN}}$ = 200 GeV is shown in Figure 4. The open circle, open diamond and open square represent the ratios of the $2^{nd},4^{th}$ and $6^{th}$ order cumulants, respectively, for different value of $|\Delta N_{p}|$. It is clear from Figure 4 that the $N_{max}$ value increases with increasing the order of cumulants. Moreover, it is observed that at relatively smaller range of $|\Delta N_{p}|$ the $C_{4}$ and $C_{6}$ can have negative sign 222In fact with increasing the $|\Delta N_{p}|$, $C_{6}$ oscillates from very high positive value to negative value before converging with its exact value, which is not visible in the current scale of Figure 4. A recent MC study with Skellam distribution reported negative value of $C_{5}$ and $C_{6}$ Pandav:2018bdx . From this study it is evident that their observed negative value is purely due to insufficient event statistics, which fails to explore the required range of $\Delta N_{p}$.

In a similar way, the $N_{max}$ values corresponding to $C_{2},C_{4}$ and $C_{6}$ for other collision energies are also estimated. The $N_{max}$ as a function of collision energy are shown in Figure 5. In Figure 5, the $N_{2}$, $N_{4}$ and $N_{6}$ values correspond to the maximum value of $|\Delta N_{p}|$ at which the $C_{2}$, $C_{4}$ and $C_{6}$, respectively, reproduced from the Skellam distribution as the expected ones. The numerical values of $N_{2}$, $N_{4}$ and $N_{6}$ for different energies are also given in Table 4. From Figure 5 and Table 4 it is clear that for a given order of cumulant, with increasing collision energy the $N_{max}$ decreases. The decrease is faster at the lower energy range. Furthermore, it is observed that for given collision energy the $N_{max}$ shows an increase by almost a constant factor with increasing the order of cumulants. A similar observation is made in Ref. Morita:2013tu . This implies that more event statistics are required to reveal the tail and in particular the higher order cumulants of the distribution.

In the second step, the constructed probability distributions are compared with the Skellam distribution to find out the signature of O(4) criticality. The ratio $f(\Delta N_{p})^{PCM}/P(N)^{S}$ is estimated for different collision energies and illustrated as a function of $\delta N_{p}/N_{6}$ in Figure 6. For Au-Au collisions at $\sqrt{s_{NN}}$ = 14.5 GeV, the ratio $f(\Delta N_{p})^{PCM}/P(N)^{S}$ shows an oscillatory behaviour at very small value of $|\delta N_{p}/N_{6}|$ value, and then increases very rapidly with increasing $|\delta N_{p}/N_{6}|<$ 0.5. For Au-Au collisions at $\sqrt{s_{NN}}$ = 19.6, 27, 39, 54.4 and 62.4 GeV, the ratio $f(\Delta N_{p})^{PCM}/P(N)^{S}$ shows a small increase for very small positive value of $\delta N_{p}/N_{6}$ value. Then it shows a sharp drop for $\delta N_{p}/N_{6}<$ 0.5. This trend is opposite for the negative value of $\delta N_{p}/N_{6}$, which shows the asymmetry of the probability distributions. The drop in the ratio occurs relatively a smaller value of $\delta N_{p}/N_{6}$ with decreasing the energy. This observation is inline with the model study done for the probability distribution constructed based on the FRG approach and with the efficiency uncorrected data Morita:2014fda . This early drop in the ratio and asymmetry of the distribution implies that the criticality will appear at a smaller value of $\delta N_{p}$, and hence in the lower order of cumulants.

For Au-Au collisions at $\sqrt{s_{NN}}$ = 200 GeV, the distributions show very small oscillatory behaviour for small $\delta N_{p}/N_{6}$ value. But with increasing the $\delta N_{p}/N_{6}$ the distribution gets narrower up to 0.5. After that, it again broadens for both values of $\delta N_{p}$, which is reflected in the sharp increase of the ratio. The observed feature for this energy is qualitatively similar to the feature seen in the model study based on the Landau theory of phase transition Morita:2012kt . However, this feature does not agree with the model study based on the FRG approach Morita:2014fda . Interestingly, some of signatures of the O(4) criticality, such as narrowing of the probability distribution for $\delta N_{p}/N_{6}<$ 0.5 is observed. Hence, the possibility of remnants of the chiral phase transition in the constructed probability distributions can not be ruled out. To be noted here that the cumulant results used in the PDF constructions are not corrected for volume fluctuations and the effect of global baryon number conservation. Therefore, more insight can be gained after understanding these effect.

It is also suggested that the narrowing of the distribution due to O(4) criticality results in a negative value of sixth-order cumulant. So for a further test of O(4) criticality, the $C_{6}$ and its ratio $C_{6}/C_{2}$ for different collision energies are estimated from the PDFs, which is discussed in the following section.

3.5 Estimating higher-order cumulants

From the obtained PDFs of net-proton multiplicity distributions the ratio, $C_{6}/C_{2}$, are calculated analytically for $\sqrt{s_{NN}}$ = 7.7 - 200 GeV. As the experimental results are associated with statistical uncertainties, their effects on the estimated $C_{6}/C_{2}$ values need to be evaluated. It should be noted here that the propagation of uncertainties from the first four cumulants to the analytically estimated $C_{6}/C_{2}$ value is rather nontrivial due to their complex relationship with the constant parameters of Eq.(7)-(11), and the constructed PDF in PCM. Therefore, we consider the following sampling method. The value of $C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$ are sampled (5000 times) randomly by varying their values within the 3$\sigma_{c}$ of their mean values, assuming independent Gaussian distributions. Here $\sigma_{c}$ refers to the statistical uncertainty of a given cumulant at a given energy. Then the constant parameters of Eq.(7)-(11) and hence the analytic value of $C_{6}/C_{2}$ are estimated for each set of $C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$ values. Finally, $C_{6}/C_{2}$ and the uncertainty on it is estimated by calculating the mean and standard deviations of all the samples, respectively. The $C_{6}/C_{2}$ result for $\sqrt{s_{NN}}$ = 7.7 - 200 GeV is shown in Figure 7. It can be observed from Figure 7 that the $C_{6}/C_{2}$ result for $\sqrt{s_{NN}}$ = 7.7 GeV has the largest uncertainty, which is due to the largest statistical uncertainties associated with the input cumulants ($C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$), among all the energies. However, it is difficult to draw any direct proportionality relation between the estimated uncertainties and the uncertainties of the lower order cumulants due to the complexity mentioned above.

Furthermore, the baseline estimation is done for the BES data from the Skellam and NBD/BD distributions. When the proton and anti-proton probability density distributions are assumed to be two uncorrelated Poissonian distributions, the net-proton distribution becomes Skellam distribution. Skellam distribution is widely used in HRG model to estimate the cumulants of conserved charges Karsch:2010ck ; Garg:2013ata ; Braun-Munzinger:2014lba . The $n^{th}$ cumulant of Skellam distribution is defined as $C_{n}=c_{1}(p)+(-1)^{n}c_{1}(\bar{p})$. Here $c_{1}(p)$ and $c_{1}(\bar{p})$ are the first cumulant of proton and antiproton distributions, respectively. So the all even order cumulants are same, and hence their ratios, like $C_{4}/C_{2}$ and $C_{6}/C_{2}$ will always yield unity. It can be seen here that values of all orders of cumulants of net-proton distribution depend entirely on the first cumulants of proton and anti-proton distributions.

Recently, NBD/BD has been proposed as one of the baselines for such study Tarnowsky:2012vu . In this case, the proton and anti-proton distributions are assumed two uncorrelated NBD/BD. For a given NBD/BD, one can calculate the cumulants, like $C_{2},C_{4}$ and $C_{6}$ analytically by knowing the value of the first and second cumulants Tarnowsky:2012vu . For the BES data, proton mean value is greater than its variance ($\mu>\sigma^{2}$). This is also true for anti-proton. Hence, BD has to be used to approximate the proton and anti-proton multiplicity distributions. Under this assumption, the $C_{n}$ of net-proton multiplicity distributions are calculated from the cumulants of proton and anti-proton distributions as follows.

[TABLE]

where $c_{n}(p)$ and $c_{n}(\bar{p})$ are the $n^{th}$ order cumulants of the assumed BD for proton and anti-proton multiplicity distributions, respectively.

The $C_{6}/C_{2}$ results of PCM compared with the baselines results estimated from Skellam and BD are illustrated in Figure 7. The baseline result estimated from BD does not show any energy dependence and remains negative for all the energies. Interestingly, the $C_{6}/C_{2}$ result obtained from the constructed PDFs using PCM show clear energy dependence. It can be observed from Figure 7 that the energy dependence of $C_{6}/C_{2}$ show a qualitatively similar trend to the energy dependence behavior of $C_{4}/C_{2}$ reported in Ref.STAR:2020tga ; STAR:2021iop . The BES result of $C_{6}/C_{2}$ deviates from Skellam expectations with oscillating in nature. For the energy range $\sqrt{s_{NN}}\geq$ 39 GeV, the value of $C_{6}/C_{2}$ does not show any significant change. In this energy range, the $C_{6}/C_{2}$ deviates from Skellam similar to the $C_{4}/C_{2}$ data. But in this case, the values are always greater than the Skellam expectations. For Au-Au collisions at $\sqrt{s_{NN}}$ = 19.6 and 27 GeV, the $C_{6}/C_{2}$ becomes negative, and particularly $\sqrt{s_{NN}}$ = 27 GeV result is almost same as the BD expectations. To be reminded here that for these two energies the probability distributions are constructed as SR Beta distributions, which is a conjugate prior probability distribution for the BD. So their cumulants results are close to each other.

Going from 11.5 GeV to 7.7 GeV, the $C_{6}/C_{2}$ values shows a rapid rise and deviate strongly from the Skellam and NBD expectations. At 7.7 GeV, a very large value of $C_{6}/C_{2}$ is observed. A similar high value of $C_{6}/C_{2}$ for 7.7 GeV is estimated from completely different approach Bzdak:2018uhv .

Before making more conclusive remarks on the $C_{6}/C_{2}$ result obtained from PCM, it is imperative to establish the critical behavior from the first four cumulants as they are the sole input parameters of this method. Furthermore, some recent studies show that sources having non-critical fluctuations and their contribution, like baryon stopping, participant fluctuations, and the effect of baryon number conservation, have significant contributions to the measured cumulant results Thakur:2016znw ; Braun-Munzinger:2016yjz . To quantify the effects of such contribution on the estimated $C_{6}/C_{2}$ value from PCM, the following consideration is made. If there will be 5$\%$, 10$\%$ and 20$\%$ uncertainties on $C_{2}$, $C_{3}$, and $C_{4}$ data, respectively, due to such effects for Au-Au collisions at $\sqrt{s_{NN}}$ = 200 GeV, then there will be around 35$-$40$\%$ uncertainty on the $C_{6}/C_{2}$ value. Therefore, it is crucial to have a quantitative idea about the contributions from the sources of non-critical fluctuations before making any conclusive remark on the present result.

Recent measurement of $C_{6}/C_{2}$ from BES-I data by the STAR experiment shows some negative value for 0$-$10$\%$ centrality in Au-Au collisions at $\sqrt{s_{NN}}$ = 200 GeV, which indicates a smooth crossover at top RHIC energy STAR:2021rls . The upcoming RHIC BES-II run aims for studying the $C_{6}/C_{2}$ in larger acceptance with larger event statistics, which will be helpful to validate the baseline predictions of $C_{6}/C_{2}$ obtained from this model Luo:2017faz .

4 Summary

In summary, for the first time Pearson curve method is used to derive the PDFs of proton, anti-proton and net-proton multiplicity distributions using the efficiency corrected results of Au-Au collisions at $\sqrt{s_{NN}}$ = 7.7 - 200 GeV. It is found that for Au-Au collisions at 200 GeV, the constructed PDFs of proton and anti-proton multiplicity distributions are more close to the Binomial distributions than Poisson distribution. The constructed PDFs of net-proton distributions have potential importance in the context of the quark-meson (QM) model within the functional renormalization group approach to study the $O(4)$ criticality. Before testing the O(4) criticality of net-proton distribution, the maximum values of $|\Delta N_{p}|$ for $C_{2}$, $C_{4}$ and $C_{6}$ of Skellam distributions are estimated for the BES data. It is found that the $N_{max}$ value increases with increasing the order of cumulants. Moreover, it is shown that the $C_{4}$ and $C_{6}$ values obtained from a Skellam distribution can have a negative sign at a smaller range of $|\Delta N_{p}|$. The obtained $N_{max}$ value will provide a guideline for MC simulation involving Skellam distribution. The test for O(4) criticality is done by computing the ratio of the constructed PDF to the Skellam distributions for different values of $\delta N_{p}/N_{6}$. For 200 GeV the distribution is found narrower than the Skellam distribution for $\delta N_{p}/N_{6}<0.5$, and gets broader after that. At lower energies (19.6, 27 and 39, 54.4, 62.4 GeV) the distribution is found asymmetric and gets narrower at a very small range of $\delta N_{p}/N_{6}$. From this ratio, some of the signatures of O(4) criticality, such as narrowing of the distribution is observed. For further investigation, the $C_{6}/C_{2}$ is estimated from the constructed PDF of the BES data and compared with Skellam and NBD expectations. The beam energy dependence of $C_{6}/C_{2}$ value has a qualitatively similar trend as the $C_{4}/C_{2}$ data. The ratio deviates from the Skellam and NBD expectation. Only at 19.6 and 27 GeV, the ratio has a negative value. A very large positive value is observed for 7.7 GeV. The lower energy results ($\sqrt{s_{NN}}\leq$ 19 GeV) have large uncertainties. At present, it still remains unclear about the signature of O(4) criticality from the constructed PDFs. It needs further investigation by considering contributions from non-critical fluctuations, like volume fluctuations, baryon stopping and global baryon number conservations. A quantitative comparison between the upcoming RHIC BES-II run data and PCM results will helpful to validate the baseline predictions. In addition to that, the derived PDFs can be used for model study and qualitative comparison with other theoretical predictions.

Acknowledgements.

The authors thanks the STAR Collaboration for providing the preliminary data during the initial phase of this study. The author acknowledges R. Bellweid, S. Dash, K. Morita and B. K. Nandi for their valuable suggestions during the preparation of this manuscript. This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2018R1A5A1025563).