High cumulants of conserved charges and their statistical uncertainties^*

One of the main motivations of high energy heavy-ion collisions is to study the structure of the Quantum Chromodynamics (QCD) phase diagram [1]. The cumulants of conserved charges are powerful observables that are sensitive to the location of the QCD critical point [2–8]. Theoretical calculations predict that the third order cumulant (C₃) is proportional to correlation length (ξ) as C₃∼ξ^4.5, and the fourth order cumulant grows more rapidly as C₄∼ξ⁷ [2]. Hence, the cumulants of conserved charges exhibit large fluctuations near the QCD critical point. Besides locating the QCD critical point, the extraction of the freeze-out temperature and baryon-chemical potential directly from QCD first principles becomes achievable, by comparing the experimental and lattice QCD measured cumulants [9–12].

With the data collected from the first beam energy scan (BES I) program at the Relativistic Heavy Ion Collider (RHIC), cumulants of net-proton [13, 14] and net-charge [15, 16] multiplicity distributions have been measured. For κσ² of net-charge multiplicity distributions, where κσ²≡C₄/C₂, no energy dependence has been observed by either the STAR or PHENIX experiment. The κσ² of net-proton multiplicity distributions, especially for the recent preliminary results with larger transverse momentum range, 0.4 < p_T < 2.0 GeV/c, shows a possible signal of non-monotonic variation at RHIC/STAR. When comparing experimental measurements to theoretical calculations, we should keep in mind that the non-critical background contribution can affect the measured results in the experiment. Currently, there are many discussions on techniques for reducing the non-critical contributions for cumulants, such as statistical error estimation [17–19], centrality bin width effect [20–22], and efficiency correction [23–26].

With the current statistics collected at RHIC BES I, errors of the high cumulants are significantly larger than other event-by-event observables, such as dynamical net-charge fluctuations [27], K/π, p/π and K/p fluctuations [28]. It has been found that the estimated errors are reasonable based on delta theorem and bootstrap methods [17, 19]. However, with limited statistics, we did not notice the relationship between the measured values and their corresponding errors. In fact, a study using the method of centrality bin width correction (CBWC) has indicated that the estimated error is influenced by the measured cumulant [18]. When applying the CBWC method with two different schemes, weighting the cumulant by the number of events or error, Ref. [18] shows that values of κσ² based on these two schemes are not consistent with each other. It indicates that the estimated error is not only proportional to , where n is the number of events, but also to the measured value of κσ².

For N independent observations (x₁, x₂, ···, x_N) with the same unknown expectation (μ), they are approximated with a normal distribution as x_i∼(μ,σ²), if they are obtained from independently generated samples with the same statistics and conditions. The error of x_i is the width of this normal distribution (σ). Experimentally, we cannot directly obtain the width σ from the normal distributions, since what we can obtain is just one measurement of the N observations. In this case, the same data sample is utilized to obtain the value of the measurement and estimate the corresponding error. Consequently, the estimated error is influenced by the measured value. In this paper, we will comprehensively study the properties of errors of κσ² in relativistic heavy-ion collisions. We will start with discussions of the correlation between measured values of κσ² and estimated uncertainties. Generally, with a given number of events, the larger κσ², the larger the statistical uncertainties that are estimated. In Section 3, we will discuss the effectiveness of the measured averaged κσ² in two schemes: the event- and error- weighted. The probabilities of κσ² that lie outside three and four standard deviations of the expectation are discussed in Section 4. Finally, the results are summarized in Section 5.

Supposing the particle and anti-particle are produced independently according to a Poisson distribution, then the net-particle follows a Skellam distribution [29–31]. The parameters of the Skellam distribution are the means of those for the particles and anti-particles. Referring to the means of protons and anti-protons at low energies of RHIC/STAR, the parameters of the Skellam distribution are set to m₁=14.5,m₂=0.6 [13].

The upper panel of Fig. 1 shows probability distributions of κσ² with different numbers of events: 0.01 million (0.01M), 0.1M and 1M. In each sub-figure, the probability distribution is obtained from 5M independent randomly generated samples. The means of these three distributions are 0.9895, 0.9989 and 0.9998, respectively. Those values are on the brink of unity, which is the theoretical expectation. The widths of these three distributions are 0.91, 0.29 and 0.092, respectively. According to the definition of the statistical error, the width of the distribution is the statistical error with the given statistics. In our simulations, the same data are used to measure κσ² and estimate their errors. Consequently, in Fig. 1(a), the estimated errors of all κσ² are of course not completely the same as the widths of the distributions. Not only that, the estimated error should be different and fluctuate even for the same values of κσ². To generally see the correlation between the values of κσ² and their associated errors, Figs. 1(d) to (f) show 〈error〉 as a function of κσ². 〈error〉 is the average of the estimated errors of κσ² in the same histogram bin. The binning method of κσ² in x–axis in the lower panel is exactly the same as that in the upper panel. With the given statistics, the larger the cumulant we measured, the larger the 〈error〉 we obtained. 〈error〉 is influenced by the measured κσ².

Zoom In Zoom Out Reset image size

When κσ² is unity, the vertical lines from Figs. 1(d) to (f) show the values of 〈error〉 are 0.91, 0.29 and 0.092, respectively. Those values are in agreement with the widths of the corresponding probability distributions in the upper panel. Consequently, 〈error〉 is under-estimated when κσ² is smaller than its expectation, while 〈error〉 is over-estimated when κσ² is larger than its expectation. As the statistics increase, Figs. 1(d) to (f) demonstrate that the increment of 〈error〉 decreases. Consequently, the larger the statistics we use to measure κσ², the smaller the correlation between the measured κσ² and its estimated statistical error.

Here we only show influence of the measured κσ² on its statistical error. We want to give a reminder that this kind of relationship cannot be eliminated and it is not dependent on the details of the method of error propagation as long as the same collected data are used to extract the measured value and its corresponding error.

Since this influence cannot be ignored in cumulant analysis, we have to be careful when extracting some other results based on cumulants and their errors. Experimentally, the rule of three sigma (or even higher) standard deviations is the most important indication of a new phenomenon, where σ is the error of measurement. Currently, in order to reduce the initial size fluctuation, it is suggested to calculate the cumulant at each of N_ch. The cumulant is averaged over all multiplicities in a given centrality, where the average is always weighted by the number of events in each N_ch. This is called the centrality bin width correction [18, 32]. The difference in the cumulant obtained by the CBWC method based on error-weighted average and event-number average should be studied. In the following sections, we will discuss κσ² obtained by these two different schemes.

For N independent observations (x₁, x₂, ···, x_N) with the same unknown expectation (μ), if the error of each observable is e_i, the maximum likelihood estimations of their expectation and error [33] are

and

In this case, the expectation is evaluated by the error-weighted average. In general, we can choose this method to obtain 〈x〉. However, the same data are always used to obtain x_i and its error e_i in experiment, which will lead to a correlation between x_i and e_i. In theory, if the difference of e_i is only proportional to , their expectation and error can be determined by weighting the average by the number of events,

and

Figure 2 demonstrates the statistics dependence of 〈κσ²〉 obtained by event- and error- weighted averages, respectively. 〈κσ²〉 first increases as the number of events increases based on these two schemes. As we have mentioned in Ref. [34], κσ² is sensitive to the tail of the net-particle multiplicity distributions. When the number of events is insufficient, the statistics are lower, and the information we can detect about the tail of the net-particle multiplicity distributions is poorer. That is why 〈κσ²〉 obtained from the event-weighted average firstly increases as the number of events increases, and then finally saturates to unity.

Zoom In Zoom Out Reset image size

Figure 2 shows that the effect of underestimation for 〈κσ²〉 obtained from the error-weighted average is stronger than that obtained from the event-weighted average. The minimum statistics required for the event-weighted method is around 0.01M events, while it is about 0.1M events for the error-weighted method. Currently, in most central collisions (0-5% centrality), the number of events in each of N_ch, or in each δ1% centrality bin width, is significantly lower than 0.1M. This is why the event-weighted average is in favor in current RHIC/STAR experiments [13, 15, 20]. On the other hand, we still need to know that the event-weighted method is effective only when the number of events in each of N_ch is larger than 0.01M. Currently, the analyzed statistics for cumulants of net-proton multiplicity distributions is 6.6M at in Au + Au collisions. If we simply assume that of events are within 0-5% centrality and there are 100 multiplicity bins in 0-5% centrality, then the number of events in each N_ch is just around 4100. It means that κσ² would be under-estimated if the event-weighted CBWC method is applied in each of N_ch at .

Figure 2(b) shows qualitatively similar statistical dependence to Fig. 2(a). However, it requires more statistics to approach the saturated value of unity. This means more statistics are needed to get a stable expectation as the parameters of the Skellam distribution become larger. So the required statistics depend on the distribution, or the mechanism of particle production.

Here, we only show the statistics dependence of 〈κσ²〉. For other cumulants, such as Sσ and C₆/C₂, the required statistics should be observable dependent. It should be studied carefully case by case.

Since the estimated error is influenced by the measured κσ², we should examine that if three-sigma rule of thumb is still applicable. To see the deviation of measured κσ² from its expectation, the probability distributions of , , is demonstrated in Fig. 3. With three cases of statistics shown in Fig. 3, decreases rapidly as increases. The tail of is dependent on the statistics. The lower the statistics, the longer the tail of the probability distribution.

Zoom In Zoom Out Reset image size

In a normal distribution, 0.3% and 0.007% of the values lie outside 3σ and 4σ standard deviations of the expectation respectively. Table 1 shows that when statistics are lower than 1M, the values of and are all slightly larger than that in a normal distribution. Therefore, three sigma and/or four sigma confidence intervals for κσ² are applicable when the number of the events is above 1M.

Table 1. Probabilities of measured κσ² that lie outside three and four standard deviations of the expectation.

statistics (in millions)
0.01	0.74	0.079
0.1	0.34	0.015
1	0.28	0.0077

In summary, we have studied high cumulants of conserved charges and their statistical uncertainties in relativistic heavy-ion collisions. With a given number of events in Monto Carlo simulations, the measured cumulants randomly fluctuate with an approximately normal distribution. The mean of the distribution is equal to the theoretical expectation and the width of the distribution is proportional to . The estimated uncertainty is not only determined by , but also influenced by the measured values of the cumulants. Generally, with a given number of events, the larger the cumulants measured, the larger the estimated statistical uncertainties. As the number of events increases, the influence of the measured cumulant on its statistical uncertainty is reduced.

Since the estimated statistical uncertainties are influenced by their measured values of cumulant, we should be careful if we re-evaluate other measurements based on the measured cumulants and their uncertainties. It is found that values of 〈κσ²〉 are different depending on whether the event-weighted or error-weighted method is used. Consequently, the values of κσ² are different when using the CBWC method based on the event-weighted and error-weighted methods. With limited available STAR BES I data, the error-weighted average scheme is not suitable. At , in most central collisions, κσ² is also under-estimated if the event-weighted CBWC method is applied in each of N_ch.

We have also studied the probability distributions of with different statistics, and found that the three sigma rule of thumb is applicable when the statistics are above 1M.