Net-Proton Fluctuations at FAIR Energies Using PHQMD Model

1. Introduction

One of the unsolved questions in the study of heavy-ion collisions is understanding how the properties of matter in the Quark-Gluon Plasma (QGP) phase evolves with decreasing beam energy, and determining the energy threshold below which QGP formation no longer occurs. The main goal of relativistic heavy-ion collision experiments is to explore the phase structure of Quantum Chromodynamics (QCD) matter. Lattice QCD calculations predict a smooth crossover from hadronic matter to a deconfined state of quarks and gluons at high temperatures and vanishing net baryon chemical potential (${\mu }_B=0$) [1]. But at higher values of ${\mu }_B$, theoretical QCD-based models suggest the existence of a first-order phase transition boundary [2,3]. This boundary is expected to terminate at a critical point [4,5]. However, due to the sign problem in Lattice QCD at finite baryon density, the location of the critical point remains highly uncertain [4].

Fluctuations of different conserved quantities, such as net-baryon (B), net-charge (Q) and net-strangeness (S) number, are analyzed on an event-by-event basis, where each event corresponds to the collision of two nuclei. Study of these fluctuations allows to probe the system’s correlation length ($\xi$), which characterizes the extent of spatial and temporal correlations among variables. Near the QCD critical point, the correlation length is expected to diverge in an infinite system. As a result, fluctuation measurements provide a sensitive tool for exploring critical phenomena and searching for the QCD critical point [6,7,8,9,10].

The CBM experiment at FAIR, Darmstadt, will offer a unique opportunity to investigate the model-predicted first-order phase transition from hadronic matter to the QGP with high precision, enabled by its exceptionally high collision rate of up to 10 MHz, significantly exceeding that of previous experiments.

In this study, we present the first measurements of the moments of net-proton (a proxy for net-baryon) distributions in Au+Au collisions at $\sqrt{s_{NN}}$ ranging from 3.5 to 19.6 GeV, using the PHQMD model within the Solenoidal Tracker at Relativistic Heavy Ion Collider (STAR) acceptance i.e., $\left|y\right|$<0.5, 0.4$<p_T<$2.0 and CBM detector acceptance i.e., 1.08$<y<$2.08, 0.4$<p_T<$2.0 (y is the rapidity and $p_T$ is the transverse momentum). These results serve as baseline calculations for proton cumulants and correlation functions from the PHQMD model, contributing to the ongoing search for the QCD critical point in heavy-ion collisions. At low collision energies, the net-proton number is primarily dominated by protons, making it a sensitive observable in this regime.

The paper is organized as follow : At first we will introduce the observables of fluctuation study, the relation between the cumulants. Then, we will discuss the PHQMD model. In the results, we will show the energy dependence and centrality dependence of net-proton cumulants in Au+Au collisions at $\sqrt{s_{NN}}$ ranging from 3.5 to 19.6 GeV with the PHQMD model generated data. We will finalize by giving a summary and a conclusion.

2. Cumulants

To quantify the intricacies of the distribution, fluctuations are typically characterized using cumulants. The cumulants of a distribution up to the fourth order are defined as follows,

$$C_1=\langle N\rangle$$

(1)

$$C_2=\langle \delta N^2\rangle$$

(2)

$$C_3=\langle \delta N^3\rangle$$

(3)

$$C_4=\langle \delta N^4\rangle -3{\langle \delta N^2\rangle }^2$$

(4)

where, N is the net-particle number in a single event and $\langle N\rangle$ is the average of N across all the events and $\delta N=N-\langle N\rangle$. The first and second-order moments represent the mean and variance whereas the third and fourth-order moments reflect its skewness and kurtosis, respectively.

$$M=\langle N\rangle$$

(5)

$${\sigma }^2=\langle \delta N^2\rangle$$

(6)

$$S=\langle \delta N^3\rangle /{\left(\langle \delta N^2\rangle \right)}^{3/2}$$

(7)

$$\kappa =\langle \delta N^4\rangle /{\left(\langle \delta N^2\rangle \right)}^2-3$$

(8)

The products of moments $\kappa {\sigma }^2$ and $S\sigma$ can be expressed in terms of cumulant ratios:

$$\kappa {\sigma }^2=C_4/C_2,S\sigma =C_3/C_2$$

(9)

The moment ${\sigma }^2$, skewness S, and kurtosis $\kappa$ have been found to be proportional to the powers of the correlation length $\xi$, specifically ${\xi }^2$, ${\xi }^{4.5}$, and ${\xi }^7$ [7], respectively. The n^th order susceptibilities ${\chi }^{\left(n\right)}$ are related to the cumulant $C_n$ as ${\chi }^{\left(n\right)}=C_n/V{T}^3$ [11], where V and T represent the volume and temperature of the system. The event-by-event fluctuation of the net-proton number (the difference between the number of protons and antiprotons observed) has been proposed as a key observable for signaling the QCD critical end point. In the energy range ($\sqrt{s_{NN}}$=3.5-4.9 GeV) of the CBM experiment, there is a hopeful expectation that the QCD critical point can be discovered.

3. PHQMD

The Parton-Hadron-Quantum-Molecular-Dynamics (PHQMD) model is developed to explore the complex dynamics of heavy-ion collisions and the formation of clusters and hypernuclei over a broad spectrum of beam energies, from several hundred AMeV up to the ultra-relativistic scales. It builds upon the Quantum Molecular Dynamics (QMD) model and the Parton-Hadron-String Dynamics (PHSD). The PHQMD model describes the particle production from SIS up to the LHC energies [12].

• Baryon Propagation (QMD Approach): In PHQMD, baryons are treated as Gaussian wave packets and propagated using the QMD model, which includes density-dependent two-body potential interactions. This approach enables the preservation of full n-body phase-space correlations among the baryons, contrasting with traditional mean-field models that average out such correlations [13,14].
• Collision Dynamics (via PHSD): PHQMD considers the collision integral from PHSD to simulate the full evolution of the system, including hadronic collisions, QGP formation, partonic scatterings, hadronization, propagation of mesons and final-state hadronic interactions. The partonic phase is treated using the Kadanoff-Baym equations [15,16]. The model also incorporates in-medium effects such as collisional broadening and modifications of spectral functions for vector mesons (e.g., $\rho$, $\omega$, $\varphi$) and strange mesons (K, $\overline{K}$, $K^{*}$, ${\overline{K}}^{*}$), enabling a realistic description of hadronic dynamics in dense and hot nuclear matter.
• Dynamic Cluster Formation: PHQMD takes a dynamic approach to cluster formation with QMD. Instead of applying a static coalescence criterion at a specific time, it allows clusters to emerge through ongoing potential interactions throughout the system’s evolution. Cluster recognition is carried out using the Simulated Annealing Clusterization Algorithm (SACA) [17] or the Minimum Spanning Tree (MST) method [12]. The baryon resonances which do not participate in the cluster formation continue to propagate within the PHQMD framework until they undergo their natural decay.
• QGP Phase Identification: In PHQMD, the region where the local energy density exceeds 0.5 GeV/fm³, that regions are considered to be in the quark-gluon plasma (QGP) phase, as hadrons are expected to dissolve into their constituent quarks and gluons. In the QGP phase, the partons-quarks, antiquarks, and gluons are scattered and dynamically propagated within a self-generated scalar mean-field potential. As the system expands and the local energy density decreases to the critical value, the partons undergo hadronization into color-neutral off-shell hadrons—mesons and baryons. This process is described by covariant transition rates that conserve energy, momentum, and quantum numbers on an event-by-event basis.

In the PHQMD model, cluster identification is not inherently performed during the dynamical evolution. Instead, all clusters are initially counted without distinguishing their physical relevance. Not all clusters with mass number A<9 are considered physical; for instance, unbound systems such as a proton-proton pair are excluded from the final list of physical clusters. In our analysis, we have used MST in the PHQMD model for clusterization. In this approach, cluster formation is determined purely based on coordinate space information. Two nucleons are considered part of the same cluster if their spatial separation is less than a threshold distance $r_0=4\mathrm{fm}$. Nucleons connected by this spatial criterion collectively form a cluster.

In our analysis, we have used the PHQMD physical cluster QGP ON and Physical cluster QGP OFF data sets. The QGP ON state refers to the system being in the QGP phase, characterized by a local energy density exceeding 0.5 GeV/fm³. If this condition is not met, the system is considered to be in the QGP OFF state.

4. Results

Two million minimum-bias Au+Au collision events generated using the PHQMD model, are used to conduct this study. We analyzed the net-proton cumulants in Au+Au collisions at $\sqrt{s_{NN}}$ =3.5, 4.0, 4.9, 7.7 and 19.6 GeV for STAR acceptance and $\sqrt{s_{NN}}$ =3.5, 4.0, 4.9 for CBM acceptance.

Although PHQMD offers the flexibility to explore by changing the input parameters, resource and storage constraints prevented the generation of 2 million minimum-bias events for every QGP-ON and QGP-OFF scenario involving physical clustering. To study how the centralitywise $S\sigma$, $\kappa {\sigma }^2$ differs for the QGP ON and QGP OFF case, we have compared 2 million minimum bias events for the $\sqrt{s_{NN}}$ = 7.7 GeV. Our analysis, based on varying model parameters, focuses on results obtained using a reduced maximum impact parameter of $b_{\max }=3.2$ fm ($b_{\max }$ value is obtained as $3.2$ fm by using the Glauber model [18]), corresponding to the 0-5% most central collision events. This data set allows us to investigate and compare the effects of physical clustering and QGP dynamics under top central collision condition.

In the model calculations, we have not applied the volume fluctuations whereas this analysis technique has been used in the real data analysis of STAR. Instead of applying volume fluctuation, we performed the analysis considering smaller bin width. We determined centrality using the multiplicity distribution of all particles, whereas the STAR experiment employed a different method for centrality determination in BES-II [19] for centrality determination. In addition, the Delta theorem method [20,21] has been used to estimate the statistical error in cumulants. In this paper, we focus on the results from the 0-5% most central Au+Au collisions events and also centrality dependent results.

Figure 1 shows the net-proton multiplicity distributions for 0-5% collision centrality for Au+Au collisions at different collision energies from 3.5 GeV to 19.6 GeV considering the STAR acceptance. With increasing energy, proton number as well as net-proton number decrease because the net baryon number tend to be zero at higher energy range.

Figure 2 and Figure 3 shows the centrality dependence of net-proton cumulants ($C_n$) up to the fourth order in Au+Au collisions from PHQMD model for STAR and CBM acceptance, respectively. In Figure 2, the cumulant values tend to deviate towards the central collisions, while for the peripheral collisions, the values are close to zero across all energies. The trend observed in the net-proton $C_1$ (mean) reflects an increase in particle production, favoring matter over antimatter, as collisions go from peripheral to central. Similarly, the $C_2$ (variance) rises due to the same effect. As both proton and antiproton yields increase with the increasing centrality, the variance of their distributions also becomes larger. It is important to note that higher-order cumulants are associated with greater statistical uncertainties.

Figure 4 and Figure 5 shows the same observables as in Figure 2 and Figure 3, but for 0-5% central events and here we replace the centrality with $\sqrt{s_{NN}}$ in the x-axis in the corresponding y and $p_T$ cut for STAR acceptance and CBM acceptance. For $C_1$, $C_2$ and $C_3$, the cumulant values decrease with increase in $\sqrt{s_{NN}}$ but for $C_4$ opposite trend has been seen with more statistical error for lower energies. In Figure 2, Figure 3, Figure 4 and Figure 5 the error bars for $C_1$, $C_2$, and $C_3$ are not visibly prominent, as they lie mostly within the marker size. However, for $C_4$, being a higher-order cumulant, the statistical fluctuations are larger, resulting in clearly visible error bars.

The collision energy dependence of the cumulant ratios, $S\sigma$ (Figure 6a) and $\kappa {\sigma }^2$ (Figure 6b) for net-proton distributions in Au+Au collisions with PHQMD model generated data are presented at the energy range $\sqrt{s_{NN}}$ = 3.5-19.6 GeV considering the STAR acceptance. To eliminate volume dependence, the ratios of cumulants of different orders are utilized. The results are shown for the most central (0-5%) collisions. The STAR data [19] are compared to calculations of PHQMD without a critical point. Figure 7a,b show the $S\sigma$ and $\kappa {\sigma }^2$ of net-proton distribution for top central PHQMD data for the CBM acceptance. $\kappa {\sigma }^2$ values are negative at lower energies but appear to approach zero or slightly positive at higher energies. This might suggest a non-Gaussian distribution, possibly reflecting critical dynamics or non-Poissonian fluctuations.

Figure 8a,b shows the centrality dependence of cumulant ratios of net-proton in PHQMD data for Au+Au collisions at $\sqrt{s_{NN}}$ = 3.5-19.6 GeV for STAR acceptance. In case of CBM acceptance, the centrality dependence of $S\sigma$ and $\kappa {\sigma }^2$ are shown in Figure 9a,b, respectively. Overall $S\sigma$ decrease from central to peripheral collisions and more pronounced negative $\kappa {\sigma }^2$ values at lower energies and in central events. For the energy range $\sqrt{s_{NN}}=3.5$–4.9 GeV, the values of $\kappa {\sigma }^2$ (Figure 8b) initially decrease with centrality and then exhibit a rise in the peripheral region. In contrast, for higher beam energies such as 7.7 and 19.6 GeV, $\kappa {\sigma }^2$ decreases with centrality and remains approximately constant towards the peripheral region.

We present the comparison of centrality dependent PHQMD results for QGP ON and QGP OFF scenarios, showing $S\sigma$ and $\kappa {\sigma }^2$ of the net-proton multiplicity distribution in Au+Au collisions at $\sqrt{s_{NN}}=7.7$ GeV, within the STAR acceptance in the Figure 10. In the Figure 10a, both the blue (QGP ON) and red (QGP OFF) points show a decreasing trend with centrality, indicating that the values of $S\sigma$ are more pronounced in central collisions. The difference between QGP ON and OFF is small but visible at central collisions, where QGP effects are expected to be stronger due to higher energy density.

At the top central collisions (0–5%), there is a noticeable difference between the QGP ON (blue) and QGP OFF (red) scenarios, which we can see in the Figure 10b. QGP ON shows a significantly larger value of possibly reflecting enhanced critical fluctuations or non-Gaussian tails. At mid-central to peripheral collisions, both QGP ON and OFF values remain nearly constant and close to unity. Large error bars at top centrality reflect statistical uncertainty due to a lower number of events.

Figure 11 shows the collision energy dependence of net-proton fluctuations observables $S\sigma$ and $\kappa {\sigma }^2$ obtained from PHQMD simulations for Au+Au collisions under STAR acceptance ($\left|y\right|<0.5$, $0.4<p_T<2.0$ GeV/c) with $b_{\max }=3.2$ fm. The results compare the scenarios with and without a QGP phase (QGP ON and QGP OFF). A noticeable difference is observed in $S\sigma$ at lower energies ($\sqrt{s_{NN}}<10$ GeV), where the QGP ON case shows suppressed values, suggesting a damping of fluctuations due to QGP formation. The $\kappa {\sigma }^2$ exhibits the non-monotonic behavior in both cases, with a slightly stronger suppression in the QGP ON scenario.

5. Discussion

It is important to note that in peripheral collisions, the system formed may not reach the necessary temperature and density to undergo a phase transition or approach the QCD critical point. Therefore, our calculations mainly focus on the most central collisions. The STAR Collaboration [19] observed a non-monotonic energy dependence of net-proton kurtosis ($\kappa {\sigma }^2$) in central Au+Au collisions, suggesting possible signatures of the QCD critical point. In contrast, peripheral collisions showed a monotonic trend consistent with non-critical models. Within uncertainties, the values of $\kappa {\sigma }^2$ of our analysis are consistent with STAR experimental data for central (0-5%) Au+Au collisions energy for 7.7 & 14.5 GeV. The fluctuations decrease with increasing centrality, showing some abrupt behavior in the 5–10% centrality range in Figure 8a and Figure 9a. For $\kappa {\sigma }^2$, similar irregularities are observed in Figure 8b and Figure 9b, particularly within the 0–30% centrality range. Using our generated dataset, we investigated the differences between the QGP ON and QGP OFF scenarios as implemented in the PHQMD model. A clear distinction is observed in the $\kappa {\sigma }^2$ values for central (0–5%) Au+Au collisions at $\sqrt{s_{NN}}=7.7$ GeV within the STAR acceptance in the Figure 10b, highlighting the sensitivity of higher-order fluctuations to the presence of a QGP phase. The results shown in the Figure 11 considering $b_{max}$=3.2 fm indicate that the presence of a QGP phase can significantly influence net-proton fluctuation observables, particularly at lower collision energies, where a suppression in $S\sigma$ and a stronger non-monotonic trend in $\kappa {\sigma }^2$ are observed in the QGP ON scenario.

The exploration of the energy region ( $\sqrt{s_{NN}}$ = 3.0 - 4.9 GeV) with precision measurements is currently underway using STAR’s FXT program. It is important to note that the PHQMD model does not explicitly incorporate the critical point in its current implementation. While a non-monotonic behavior in observables such as $K{\sigma }^2$ is often regarded as a potential signal for the existence of the QCD critical point, its presence alone does not provide conclusive evidence.

Furthermore, conventional transport models like UrQMD do not account for realistic cluster formation, which plays a significant role in the evolution and final-state composition of the system. In contrast, PHQMD includes a dynamic cluster formation mechanism, allowing for the generation of light nuclei and fragments based on nucleon correlations in phase space. This aspect provides a more realistic scenario for studying fluctuations of conserved quantities, particularly in the high baryon density regime relevant to FAIR energies. As a result, PHQMD offers a more comprehensive framework for investigating fluctuation observables in heavy-ion collisions. The upcoming CBM experiment will explore this region and is expected to provide high-statistics data. A significantly higher-statistics dataset is required to investigate the QCD critical point.

6. Conclusions

By maintaining n-body correlations and dynamically modeling interactions from early to late stages, PHQMD enhances the capabilities of earlier models and offers a comprehensive picture of both hadronic and partonic matter, as well as cluster and hypernuclei formation in heavy-ion collisions.

Fluctuations of conserved quantities-such as baryon number, electric charge, and strangeness-are considered sensitive probes in heavy-ion collisions for investigating the QCD phase transition and identifying the critical point. The STAR experiment has reported measurements of the moments of net-proton distributions in Au+Au collisions over a collision energy range of $\sqrt{s_{NN}}$ = 7.7-200 GeV [19]. The PHQMD calculations for $\kappa {\sigma }^2$ are consistent with the trend of STAR data for the most central 0-5% Au+Au collisions within current uncertainties, whereas the same for $S\sigma$ slightly overestimate the data.

QCD-based model calculations suggest that net-proton fluctuations are particularly sensitive observables. In this study, the cumulant ratios of net-proton distributions from PHQMD simulations display a non-monotonic trend as a function of collision energy within the CBM energy range, albeit within the current statistical uncertainties. This analysis offers new insights into the behavior of net-proton fluctuations in the CBM energy domain and contributes to establishing a baseline for interpreting net-proton and net-baryon fluctuations. Such a baseline is essential for future high-statistics measurements with the CBM experiment aimed at probing possible critical phenomena near the QCD critical end point.