Higher-Order Cumulants and Correlation Functions of Proton Multiplicity Distributions in $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV Au+Au Collisions at the RHIC STAR Experiment

The dataset analyzed in this paper was collected in 2018 by the Solenoidal Tracker at RHIC (STAR) using a fixed-target configuration. The gold target of the thickness of 1.93 g/cm² ($0.25$ mm) corresponding to a 1% interaction probability was located 200.7 cm from the center of the Time Projection Chamber (TPC) Ackermann et al. (2003). A beam, consisting of 12 bunches of $7\times 10^{9}$ gold ions is circulated in the RHIC ring at a frequency of 1 MHz with an energy of 3.85 GeV per nucleon.

Proton multiplicities were recorded in the TPC and Time-of-Flight detectors (TOF) Llope (2012), which are located inside STAR’s solenoidal magnet. The magnet provides a uniform 0.5 T field along the beam axis. A total of $1.4\times 10^{8}$ Au+Au events at $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV were used in this analysis. The minimum bias events required a hit in either the Beam-Beam Counter (BBC) Bieser et al. (2003) or the Event Plane Detector (EPD) Adams et al. (2020) and at least three hits in the TOF. To remove collisions between the beam and the beam pipe, event vertices are required to be less than 1.3 cm from the Au target along the beamline and less than 1.5 cm from the target radially from the mean collision vertex from the TPC center along the beam line. Events are also checked on the average of variables for different run periods: charged particle multiplicity, vertex position, and track’s pseudo-rapidity $\eta$ ($\eta=0.5*\ln[\frac{p+p_{\rm z}}{p-p_{\rm z}}]$, where $p$ and $p_{\rm z}$ are total momentum and its fraction in beam direction), the distance of closest approach (DCA), and transverse momentum ($p_{\rm T}$). The outlier runs which deviate more than $\pm$3 $\sigma$ are excluded in the analysis where $\sigma$ is the standard deviation of run-by-run distributions of variables listed above.

The Au+Au collisions are characterized by their centrality. Here, the centrality is a measure of the geometric overlap of the two colliding nuclei and can be determined by measuring charged particle multiplicity in the TPC. To maximize the centrality resolution and minimize the self-correlation effect Chatterjee et al. (2020); Adam et al. (2021a); Abdallah et al. (2021); Chatterjee et al. (2021), the reference multiplicity includes charged particles except protons in the full TPC acceptance (TPC covers the pseudo-rapidity $\eta$ of $-2<\eta<0$ in lab frame). Protons 3 $\sigma$ away from theoretical expectation are excluded by a TPC particle identification cut. The anti-proton production is negligible at $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV and does not affect the centrality determination ($\overline{p}/p\sim\exp({-2\mu_{\rm B}/T_{\rm ch}})<10^{-6}$) Andronic et al. (2018). The reference multiplicity distribution in Fig. 1 is fitted with a Monte Carlo Glauber model (GM) Miller et al. (2007) coupled with a two-component model Kharzeev and Nardi (2001). The two-component model assumes multiplicity $n$ in nuclear collisions has two components which are respectively proportional to the number of participants $N_{\rm part}$ and the number of binary collisions $N_{\rm coll}$:

where $x$ and $n_{\rm pp}$ denote the fraction of multiplicity from $N_{\rm coll}$ and the mean multiplicity measured in $pp$ collisions per unit of pseudo-rapidity due to $N_{\rm coll}$, respectively. Then the simulated multiplicity per event is obtained to sample $n_{\rm pp}$ times of the negative binomial distribution (shown in Eq. 2, where $\mu$ is the mean, and is set to $n_{\rm pp}$).

The fit is performed by minimizing a $\chi^{2}$ between the measured multiplicity and the GM from the reference multiplicity from 20 to 80. The parameters for the best fit are: $n_{\rm pp}$ = 0.62, $x$ = 0.06, and $k$ = 5.56. At reference multiplicities below $10$, the data and the GM disagree due to peripheral event trigger inefficiency. At multiplicities above $80$, double collision (pileup) events dominate the multiplicity distribution. The collision centrality is determined by fitting the Glauber calculation of charged particle multiplicity distribution to that of data. According to the normalized distribution from the Glauber model, one can extract the collision parameters such as $\langle N_{\rm part}\rangle$ and the fraction of the collision centrality, 0-5%, 5-10%, …, 50-60%. In addition to a pileup correction discussed in Sec. II.6, events above the reference multiplicity of 80 are removed from the 0–5% centrality class. The selection cuts for each centrality class, $\langle N_{\rm part}\rangle$ as well as the pileup fraction are shown in Tab. 1.

The TPC measures both the trajectory and the energy loss (${\rm d}E/{\rm d}x$) of a particle. TPC spatial hits are fitted with helices to determine the charge and momentum of each charged particle. To ensure track quality, tracks are required to meet selection criteria which are at least 10 hits and more than five ${\rm d}E/{\rm d}x$ measurements. Additionally, to prevent double-counting reconstructed tracks from a single particle, a selected track is required to have more than 52% of the maximum-possible fit points, which peaks at 45 possible hits. To suppress the contamination from spallation in the beam pipe and secondary protons from hyperon decays, DCA $<$ 3 cm criterion is placed at the distance of the closest approach (DCA) in 3-dimensions of the reconstructed track’s trajectory to the primary vertex position. The results presented here are within the kinematics $-0.9<y<0$ and $0.4<p_{\rm T}<2.0$ GeV/$c$.

Particle identification (PID) is performed by measuring the ${\rm d}E/{\rm d}x$ and the time of flight in the TPC and TOF, respectively. Figure 2 (a) shows the $\langle{\rm d}E/{\rm d}x\rangle$ as a function of rigidity (the ratio of total momentum over electric charge, $|p|/q$ GeV/$c$) for the positively charged tracks. To select proton candidates, the measured values of ${\rm d}E/{\rm d}x$ are compared to a theoretical prediction Bichsel (2006) (red line). The quantity $N_{\sigma,p}$ for charged tracks in the TPC is defined as

where $\left\langle{\rm d}E/{\rm d}x\right\rangle$ is the truncated mean value of the track energy loss measured in the TPC, $\left\langle{\rm d}E/{\rm d}x\right\rangle^{\rm th}$ is corresponding theoretical prediction, and $\sigma_{\rm R}$ is the track length dependent ${\rm d}E/{\rm d}x$ resolution. The $N_{\sigma,p}$ distribution appears as a standard Gaussian distribution with a mean close to zero. The offset from zero is measured as a function of momentum in 0.1 GeV/$c$ bins and the $N_{\sigma,p}$ distribution is re-centered. The proton tracks are selected within three standard deviations of the re-centered $N_{\sigma,p}$ distribution ($|N_{\sigma,p}|<3.0$).

Figure 2 (b) shows the mass-squared ($m^{2}$) versus rigidity of charged particles in the TPC and TOF. The $m^{2}$ is given by

where $p$, $t$, and $L$ are the momentum, time of flight, and path length of the particle, respectively. The speed of light in vacuum is denoted by $c$. The protons are identified by selecting charged tracks with mass-squared values between $0.6<m^{2}<1.2$ GeV²/$c^{4}$. While the mass-squared cut provides high proton purity, it introduces a 60% matching efficiency. The proton purity is required to be higher than 95% at all rapidities and momenta for the subsequent cumulant analysis.

Figure 2 (c) shows the transverse-momentum versus rapidity for protons selected in the TPC within $|N_{\sigma,p}|<3.0$. The tracks above a momentum of 2.0 GeV/$c$ in the lab frame are required to have a mass-squared cut. The kinematic acceptance of the analysis ($-0.9<y<0$ and $0.4<p_{\rm T}<2.0$ GeV/$c$) is indicated by a black box.

Here, the definition of cumulants is provided. Let $N$ be the number of particles measured in each event. Then mean value of $N$ is given by $\left\langle N\right\rangle$ and $\delta N=N-\left\langle N\right\rangle$ is the deviation from the mean value where the $\langle.\rangle$ symbol indicates the average over events. The $r^{\rm th}$-order central moment of a distribution is described by

In terms of central moments, the cumulants are defined as

The cumulants can also be expressed in terms of raw moments (Eq. 25 in Appendix A). Some commonly used moments and ratios are given as

where $M$, $\sigma^{2}$, $S$, and $\kappa$ are mean, variance, skewness, and kurtosis, respectively. The products $S\sigma$ and $\kappa\sigma^{2}$ can be expressed in terms of the cumulant ratios as

In case there are no intrinsic correlations among the measured particles, all ratios of the cumulants are unity, thus Poisson statistics is a non-trivial baseline for experimentally measured cumulant ratios.

The probability distribution of Poisson statistics is $P(N)=\lambda^{N}e^{-\lambda}/N!$, where $\lambda$ is an average of the number of measured particles per event. The cumulants are equal to the mean value: $C_{1}=C_{2}=...=C_{n}=\lambda$, where $C_{n}$ is the $n^{\rm th}$-order cumulant. Furthermore, all cumulant ratios equal one. More discussion can be found in Ref. Asakawa and Kitazawa (2016) and Appendix B.

As discussed in Refs. Bzdak et al. (2017); Abdallah et al. (2021), the cumulants $C_{r}$ can be algebraically converted to the integrals of the corresponding multi-particle correlation functions. These integrated correlation functions also known as factorial cumulants will be simply called correlation functions (denoted by $\kappa_{i}$) henceforth. In terms of cumulants, the correlation functions up to $6^{\rm th}$-order are

A compact form of the above equations can be seen in Eq. 28 of Appendix A.

We define $C_{i}/{C_{1}}-1,\quad{i}=1,2,\cdots$ as reduced cumulant ratio. The reduced cumulant ratio of different orders can be displayed in terms of correlation function ratios as

It is clear that the $n^{\rm th}$-order reduced cumulant ratio is a combination of all multi-particle correlation functions up to the $n^{\rm th}$-order.

In this analysis, the proton tracks are corrected for detector inefficiency in the TPC and TOF. The TPC efficiency is calculated by placing Monte Carlo tracks into a GEANT FINE and NEVSKI (2000) detector simulation. The GEANT detector response is mixed with the detector response from real data and a track reconstruction process is performed. The TPC efficiency is calculated by counting the fraction of successfully reconstructed Monto Carlo tracks. Figure 3 (a) shows the TPC detector efficiency for proton tracks with respect to kinematic acceptance ($p_{\mathrm{T}}$ versus $\eta$). The TOF matching efficiency is estimated directly from data by measuring the fraction of TPC tracks that satisfy the TOF-matching criteria. Recall, a TOF-matched proton candidate requires $p>2$ GeV/$c$, $|N_{\sigma,p}|<3.0$, and $0.6<m^{2}<1.2$ GeV²/$c^{4}$. Figure 3 (b) shows TOF matching efficiency for proton tracks in a window of $p_{\mathrm{T}}$ versus $\eta$. The detector efficiency corrections are performed on a “track-by-track” basis Nonaka et al. (2017); Luo and Nonaka (2019), where the proton reconstruction efficiency as a function of $p_{\rm T}$ and rapidity is applied as a weight to each track.

The proton cumulants are evaluated in an event-by-event manner. To extract an averaged value of the cumulants from a range of the measured reference multiplicities, or in other words, from a centrality bin, a proper procedure called centrality bin width correction (CBWC) Luo et al. (2013) method is applied. The number of events from each multiplicity bin is used as a weight in the averaging procedure. Figure 4 shows centrality dependence of proton cumulants up to $6^{\rm th}$-order in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV. The black circles show the multiplicity dependence of the cumulants, while red circles and blue squares are centrality binned cumulants with and without CBWC, respectively. The CBWC is necessary to extract properly averaged cumulants in a given centrality bin.

Double collisions (“pileup”), seen in Fig. 1, are the largest source of background in the STAR fixed-target experiment. Here, we discuss the correction technique Nonaka et al. (2020) used to remove the pileup statistically. The correction technique requires an estimate of the pileup contribution as a function of reference multiplicity. Thus, an unfolding method Zhang et al. (2022) used to estimate the event-averaged pileup fraction is discussed.

The correction method assumes a pileup event is the superposition of two independent single-collision events. Let $P_{m}(N)$ be the probability distribution function to find an event with $N$ particles at multiplicity $m$. If the probability of a pileup event at the $m^{\rm th}$ multiplicity bin is $\alpha_{m}$, then $P_{m}(N)$ is

where $P^{t}_{m}(N)$ and $P^{pu}_{m}(N)$ are the single-collision and pileup probability distribution functions, respectively.

The pileup events can be decomposed into the sub-pileup event probability distribution function $P^{sub}_{i,j}(N)$ as

where $w_{i,j}$ is the probability to observe the sub-events among all pileup events at multiplicity $m$, where $m=i+j$. Additionally, the sum over $i$ and $j$ runs over non-negative integers and $\displaystyle\sum_{i,j}\delta_{m,i+j}w_{i,j}=1$, where $w_{i,j}=w_{j,i}$.

Following the procedure outlined in Ref. Nonaka et al. (2020), the single collision moments can be recursively expressed in terms of the measured moments of lower multiplicity bins as

where $\beta^{({r})}_{m}$ is defined as

and

The correction requires both $\alpha_{m}$ and $w_{i,j}$ to be determined with a high level of precision. Both parameters can be expressed in terms of the multiplicity of the single collision events $T(m)$ as

where $\alpha$ is the total pileup fraction overall reference multiplicities. Therefore, the accuracy of $\alpha_{m}$ and $w_{i,j}$ is determined by one’s ability to extract the single collision distribution from the measured reference multiplicity.

For this analysis, an unfolding technique Esumi et al. (2021) is used to estimate $T({m})$. An overview of the unfolding procedure and a closure test of simulated events can be found in Ref. Zhang et al. (2022). The unfolding is performed by generating both a pileup distribution and single collision distribution from Monte-Carlo (toy-MC) events. The difference between the toy-MC (single + pileup) distribution and the data multiplicity distribution is measured and propagated back to the toy-MC single collisions. The process is repeated until the toy-MC and data agree. The bottom panel of Fig. 1 shows the ratio of the data and toy-MC after 100 iterations. In the top panel of Fig. 1, the single collision and pileup distributions are represented by blue and green dashed lines, respectively. The procedure has one free parameter, which is the total pileup probability $\alpha$ in Eq. 15. The procedure is run for various $\alpha$ parameters and a $\chi^{2}$ test is performed. The pileup probability $\alpha$ is determined to be ($0.46\pm 0.09$)% for all events and ($2.10\pm 0.40$)% in the 0–5% centrality class. With the unfolded single collision distribution and the $\alpha$ parameter, the response matrix $w_{i,j}$ can be simulated as shown in Fig. 5. As stated, $w_{i,j}$ is the probability to observe a sub-pileup event at multiplicity $m$ with $m=i+j$. The pileup corrected cumulants are shown in Fig. 6. Additionally, the event-averaged pileup corrected (red) and uncorrected (blue) cumulants are displayed. For all cumulants, only results from the top centrality class (0-5%) are affected. Figure 7 are the pileup corrected and uncorrected cumulant ratios. Similar to the cumulants, the cumulant ratios are only affected in the most central collisions. Pileup correction will increase uncertainties in the high multiplicity region, especially for reference multiplicity larger than 60. After the pileup correction, higher-order cumulant ratios, $C_{4}/C_{2}$, $C_{5}/C_{1}$ and $C_{6}/C_{2}$, are consistent with zero within uncertainty for the most central multiplicity bins.

Physics results will be discussed for a given event centrality class. Since the physics of higher-order cumulants and their ratios are supposed to be sensitive to collision dynamics including the centrality, it is important to understand the correlation between the experimentally measured reference multiplicity distribution and the extracted class of collision centrality. It is well-known that quantum fluctuations in particle production and fluctuation of the participating nucleon pairs will affect the final centrality determination, especially at low-energy collisions. The microscopic hadronic transport model UrQMD (v3.4) Bass et al. (1998); Bleicher et al. (1999), which does not contain critical phenomena physics, has been used to show the volume fluctuation effect. As an illustration, the UrQMD model results on the correlation of the reference multiplicity and participating nucleons $N_{\rm part}$ is shown in the left panel of Fig. 8. The right panel shows the root-mean-square (RMS) values of the $N_{\rm part}$ distribution at a given fixed reference multiplicity.

As one can see, the correlation is broad and the dispersion (RMS) of $N_{\rm part}$ is as large as 30 in the mid-central collisions for 3 GeV Au+Au collisions. Even in the most central 5% collisions, the dispersion is in the range of 15. Primarily, the large dispersion is due to the fluctuation of the $N_{\rm part}$ in addition to the variation in the charged particle production for a given pair of nucleons. The variation of the initial number of participants for fixed reference multiplicity is also called initial volume fluctuation (IVF) and its implications on the results of the higher moments of proton distributions will be discussed later in the paper. In fact, as indicated by the red dashed lines in the plot(a), the top 5% central collisions are largely different events with a small overlap.

The $N_{\rm part}$ distributions are not measurable experimentally but obtained from the calculations of the Glauber and UrQMD models. Figure 9 shows the distributions from both Glauber (black solid line and red shaded area) and UrQMD (blue solid line and blue dashed lines) models. The hatched areas are corresponding to various collision centralities determined from the charged particle reference multiplicity. It is obvious that the overall distributions are quite different. More so, the widths of the top 5% are dramatically different. As discussed in Refs. Skokov et al. (2013); Braun-Munzinger et al. (2017), the initial volume fluctuation can be partly suppressed within the framework of the Wounded Nucleons Model (WNM) Bialas et al. (1976). The WNM model assumes that produced particles in nucleus-nucleus collisions are generated from inelastic scattered wounded nucleons. Each wounded nucleon (or participating nucleon) is treated as an independent source and contributes to the total number of produced particles. However, the difference in the $N_{\rm part}$ distributions of UrQMD and Glauber Model would imply a strong model dependence.

In order to demonstrate the effect of volume fluctuations, a volume fluctuation correction (VFC) method proposed in Ref. Braun-Munzinger et al. (2017) has been applied to proton cumulants from the transport UrQMD model. As seen in Fig. 9, the $N_{\rm part}$ distributions are different using different centrality determination methods. As a result, there are sizable differences in the corrected proton cumulants results shown in Fig. 10 where black dots are the ratios of proton cumulants from the UrQMD model. The results of the corrected ratios, using $N_{\rm part}$ determined from the UrQMD model directly or from Glauber fits to the charged multiplicity distributions, done exactly as in data analysis in the experiment, are shown by blue squares and red triangles, respectively. Although the correction with the UrQMD $N_{\rm part}$ is supposed to be the answer, the results with Glauber are quite different except for the most central collisions. The maximum effect, due to the volume fluctuations, is around the mid-central centrality bins. This is qualitatively consistent with the mid-central peak of the dispersion in $N_{\rm part}$ presented in Fig. 8 (b). The negligible impact on the most central Au+Au events is due to the constraint of the total number of participating nucleons $N_{\rm part}^{\rm max}=394$.

Calculations from the UrQMD model with the fixed impact parameter range $b\leq 3$ fm, corresponding to the top 5% collisions, are performed for the 3 GeV Au+Au collisions. The result is shown as blue open crosses in Fig. 10. As mentioned earlier, although the selection of centrality with impact parameter or $N_{\rm part}$(Tab. 2) eliminated the IVF in the model calculations, it is not experimentally measurable. In addition, the approach collected a different class of events as shown in Figs. 8 and 9.

The statistical uncertainties are obtained using the Bootstrap approach Efron and Tibshirani (1994) in which events are re-sampled with replacement and the analysis is re-run. The Bootstrap procedure is repeated 200 times and the statistical uncertainty is the standard deviation of the Bootstrapped observable values, such as the cumulants and their ratios.

The systematic uncertainty of the cumulant calculation can be subdivided into three categories: pileup correction, centrality determination (Tab. 1), and track selection. The track selection includes the track reconstruction requirements (TPC spatial hits, DCA), the mass-squared cut, and the efficiency in the TPC and TOF. The effect of lowering the ${\rm d}E/{\rm d}x$ cut to $|N_{\sigma,p}|<2$ was tested but did not affect the final result.

To estimate the systematic uncertainty, the analysis was repeated with different analysis requirements which are outlined in Table 3. The final total value is the quadratic sum of uncertainties from centrality, pileup, and the dominant contribution from TPC points, DCA, and PID $m^{2}$ and efficiency $\epsilon$. The difference between the systematic analyses and nominal analysis in $C_{2}/C_{1}$, $C_{3}/C_{2}$, and $C_{4}/C_{2}$ in 0-5% central Au+Au collisions is listed in Table 4.

Experimental data of proton cumulants and their ratios as a function of the reference multiplicity from 3 GeV Au+Au collisions are shown in Fig. 11 and Fig. 12. The reference multiplicity dependence of data without initial volume fluctuation correction is shown as grey open squares while data with volume fluctuation correction using $N_{\rm part}$ distributions from Glauber and UrQMD model are shown as black filled circles and black open triangles, respectively. The centrality binned results with CBWC are shown as blue open squares, red filled circles and orange filled triangles correspondingly. By definition, $C_{1}$ is not affected by the volume fluctuation correction while strong model dependence for higher order cumulants in the initial volume fluctuation corrections is clear in the figure.

For higher-order cumulants, a maximum difference between results with and without VFC is seen around mid-central Au+Au collisions and the difference slightly depends on the order of the cumulants. In the most central bin, the corrected and uncorrected proton cumulants $C_{i}$ ($i>3$) are very similar. One can see that cumulants show strong multiplicity dependence. Rapid decreases are seen from mid-central (5-10%) to most central collisions (0-5%) in $C_{3}$ to $C_{6}$. And in the high reference multiplicity region ($>50$) there is an increase with multiplicity. Indeed, as one can see that in the central collision region, 0-5% and 5-10%, the values of $C_{4}$ to $C_{6}$ (Fig. 11) and corresponding ratios (Fig. 12) change from positive to negative and to positive value again at the multiplicity larger than 60. Again, when multiplicity is larger than 50, the VFC shows no effect on either the proton cumulants ($C_{4},C_{5}$ and $C_{6}$ in Fig. 11) or their ratios ($C_{4}/C_{2},C_{5}/C_{1}$, and $C_{6}/C_{2}$ in Fig. 12).

In later discussions, the experimental data will be compared with model calculations. For the UrQMD model, around 80 million events are produced in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 3 GeV with the cascade mode. The same kinematic cuts and centrality bins used in the data analysis are applied in the model calculations. The collision centrality is determined using the charged particle multiplicity excluding protons in the acceptance of the TPC, the same procedure as used in the data analysis. The CBWC procedure is also applied to get the properly weighted centrality binned model results.

The UrQMD model calculations as a function of reference multiplicity are shown in Fig. 13. Cumulants show a strong dependence on reference multiplicity. The strong multiplicity dependence in the higher order of the proton ratios is very similar to that observed in experimental data, see Fig. 12. Stronger variations are seen in the higher order ratios.

In this section, we discuss the centrality dependence of the proton cumulants and correlation functions along with the corresponding ratios. Assuming that collisions are a superposition of independent sources, one expects the cumulant values to increase with $\left\langle N_{\rm part}\right\rangle$. The centrality classes are related to the average number of participating nucleons $\left\langle N_{\rm part}\right\rangle$ shown in Tab. 1.

Figures 14 and 15 show proton cumulants and correlation functions as a function of $\left\langle N_{\rm part}\right\rangle$ with VFC using two different models. In the cumulant and correlation function ratios, the corrections suppress the large values of cumulant ratios in mid-central and peripheral collisions. In addition, one observes large variations in the VF corrected results between the UrQMD and Glauber Model. However, as can be seen in these figures, in the most central 0-5% Au+Au collisions, the difference for higher order ones ($C_{i},\quad i>2$) between data without the VFC and with the VFC using different model inputs is small. This implies that the results of cumulants as well as the correlation functions in the most central collisions are least affected by the volume fluctuations. Because of the strong model dependence, starting from Fig. 16 the VFC method is not adopted and the results from 3 GeV experimental data are only applied with pileup correction and CBWC. For UrQMD results, only CBWC is applied.

Figure 16 shows the centrality dependence of the proton cumulants and their ratios extracted from the kinematic acceptance $-0.5<y<0$ and $0.4<p_{\rm T}<2.0$ GeV/$c$. The experimental data are compared with the UrQMD calculations (gold bands). As one can see, only $C_{1}$ and $C_{2}$ values increase with $\left\langle N_{\rm part}\right\rangle$. For the higher order cumulants ($C_{i},\quad i>2$), the cumulants increase with $\left\langle N_{\rm part}\right\rangle$ ($\left\langle N_{\rm part}\right\rangle<200$) but change rapidly in the more central centrality region. For all cumulant ratios, the values are above unity in the peripheral and are closer to unity for mid-central collisions. Figure 17 shows the centrality dependence of the proton correlation functions and ratios with the same acceptance as Fig. 16. The correlation functions also deviate from a monotonic increase around $\left\langle N_{\rm part}\right\rangle\sim 200$. The trends of cumulants and correlation functions shown in Figs. 16 and 17 can be qualitatively reproduced by the UrQMD calculations.

In this section, we discuss the rapidity and $p_{\rm T}$ dependence of the cumulants and cumulant ratios.

Figures 18 and 19 depict the rapidity (left panels) and transverse momentum (right panels) dependence of ratios of proton cumulants and correlation functions. Data from most central ($0-5\%$) and peripheral ($50-60\%$) centrality classes are shown in the figures. The measured rapidity window covers $y_{\rm min}<y<0$, where $y_{\rm min}$ changes from -0.2 to -0.9 and $p_{\rm T}$ window is $0.4<p_{\rm T}$ (GeV/$c)<p_{\rm T}^{\rm max}$, where $p_{\rm T}^{\rm max}$ is varied from $0.8$ to $2.0$ GeV/$c$. Corresponding results from the UrQMD calculations are shown as colored bands in the figures.

As one can see, in the most central collisions, the cumulant ratio $C_{2}/C_{1}$ in Fig. 18, remains above unity at all rapidities. The $C_{3}/C_{2}$ ratio is slightly above unity for the smallest rapidity window ($y_{\rm min}=-0.1$) and decreases as the rapidity window increases.

As is expected, the $C_{4}/C_{2}$ ratio is close to unity in the smallest rapidity window and seems to go back to unity with large uncertainty when the rapidity window is larger than $y$ = -0.5. Similarly, the ratios of the correlation function in Fig. 19 (e) are also close to zero (Poisson baseline) at the smallest rapidity window but show deviations from zero when it goes to a larger rapidity window. These rapidity dependencies are reproduced by the UrQMD calculations.

Overall, the results of the hadronic transport model UrQMD calculations qualitatively reproduce both the rapidity and transverse-momentum dependence. As discussed earlier, in addition to the genuine collision dynamics in the model, the effect of the volume fluctuations is also present in the calculation.

In the first order, taking the ratio of cumulants cancels the effect of volume but not the fluctuations in volume. Figure 20 depicts the collision energy dependence of the cumulant ratios from 0-5% central (top panels) and 50-60% peripheral (bottom panels) collisions. The new result of protons from 3 GeV Au+Au collisions data, shown as filled squares, is compared to that of protons (open squares) and net-protons (filled circles) from Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7 – 200 GeV.

UrQMD results with $|y|<0.5$ of proton are shown as gold bands while that of net-proton are shown as green dashed lines or green bands. At 3 GeV, the model result for proton ($-0.5<y<0$) are shown as blue crosses. While the net-proton ratios show a clear energy dependence, the proton $C_{2}/C_{1}$ and $C_{3}/C_{2}$ ratios are relatively flat and around unity as a function of collision energy except for the 3 GeV data. The new proton data from 3 GeV does not follow this trend in the most central collisions.

Notably, both proton (open squares) and net-proton (filled circles) cumulant ratios converge at collision energies below 20 GeV. This implies that at high baryon density region, the anti-proton production becomes negligible. At the center of mass energy of 2.4 GeV, HADES reported the values for proton cumulant ratios in 0-10% central Au+Au collisions: $C_{3}/C_{2}$ = -1.63 $\pm$ 0.09 (stat) $\pm$ 0.34 (sys) and $C_{4}/C_{2}$ = 0.15 $\pm$ 0.9 (stat) $\pm$ 1.4 (sys) from $|y|<0.4$, $0.4<p_{\rm T}<1.6$ GeV/$c$ Adamczewski-Musch et al. (2020). While the value of $C_{4}/C_{2}$ from the HADES experiment is consistent with the new 3 GeV data, the sign of $C_{3}/C_{2}$ is opposite to what we observed here.

Except for the $C_{3}/C_{2}$ ratio from 3 GeV central collisions, the UrQMD results reproduce the energy dependence trend well for both proton and net-proton, see green and gold bands in the figure. For the peripheral 50-60% collisions, the $C_{4}/C_{2}$ ratio from 3 GeV is larger than that from higher energy collisions, by a factor of five. A rapid increase in the energy dependence seems confirmed by the UrQMD model calculations, see both the blue cross and gold band in the figure. In the 3 GeV most central collisions, unlike all higher energy collisions, the value of $C_{4}/C_{2}$ is negative. The UrQMD model calculations, again, reproduce the trend well: due to baryon number conservation, the $C_{4}/C_{2}$ is dramatically suppressed in the high baryon density region.

Hydrodynamic calculations are shown as red dashed lines in Fig. 20 for the 0-5% Au+Au collisions. The hydrodynamic evolution is made with the open-source code MUSIC (v3.0) Denicol et al. (2018). The initial condition is taken from Ref. Shen and Alzhrani (2020) and the particlization is given by the Cooper-Frye formula Cooper and Frye (1974) with non-ideal hadron resonance gas model Vovchenko et al. (2017). At the grand canonical limit, including both effects of excluded volume and global baryon number conservation, the net-proton cumulants are evaluated on the Cooper-Frye hypersurface. One may find more details of the model calculations in Ref. Vovchenko et al. (2022). Unlike the commonly used transport model approach, here all calculations, starting from initial condition to hydro-evolution to hadronization, are all performed using averaged ensembles. Cumulant ratios $C_{2}/C_{1}$, $C_{3}/C_{2}$, and $C_{4}/C_{2}$ in hydrodynamic calculations are all below unity. Interestingly, the UrQMD results with a fixed impact parameter are also suppressed, see open blue crosses. Qualitatively, the results from the fixed impact parameter ($b\leq 3$ fm blue open crosses) UrQMD calculations follow that of the hydrodynamic calculations (red dashed lines) with a canonical ensemble.

The fact that the negative $C_{4}/C_{2}$ in the most central Au+Au collisions at 3 GeV is reproduced by the hadronic transport UrQMD model although $C_{3}/C_{2}$ is over-predicted, implies that the system is dominated by hadronic interactions. This conclusion is also consistent with the measurements of the collectivity of light hadrons Abdallah et al. (2022a) as well as the strange hadron production Abdallah et al. (2022b) at the same collision energy.

Due to baryon number conservation, proton multiplicity distributions are also modified by the formation of light nuclei in the same collision Fecková et al. (2015). The effect is especially strong in the high baryon density region where the production of light nuclei is expected to be relatively high Andronic et al. (2011). The influence of the light nuclei on the proton cumulants and their ratios will be analyzed in the future when the yields of light nuclei become available.