Bayesian analysis was employed to constrain the Equation of State (EoS) of nuclear matter with a baryon density of up to six times the nuclear saturation density, using data from heavy-ion collisions at beam energies $\sqrt{s_{\text{NN}}}=2-10$ GeV. The resulting EoS excellently agrees with that constrained by astrophysical observations.

The nuclear EoS describes the relationship between different thermodynamic variables in the bulk nuclear matter. The relationship between the pressure P and energy density ϵ determines the square of the speed of sound in the nuclear matter $c_{\text{s}}^2=dP/dϵ$. The value of $c_{\text{s}}^2$ is approximately 1/3 in the ideal gas limit of the quark-gluon plasma phase (QGP) (or any system with conformal symmetry) and about 1/6 in the hadron resonance gas (HRG) phase of nuclear matter with zero net baryon density. This corresponds to the state of matter existing in the early universe and produced in heavy-ion collisions at the top RHIC and LHC energies, where nearly equal amounts of matter and antimatter are created. Lattice QCD calculations predict that the transition between the QGP and HRG of nuclear matter is a smooth crossover. The nuclear EoS also depends on the net baryon density, such as that produced in low-energy heavy-ion collisions with strong baryon stopping or dense nuclear matter inside neutron stars. Nuclear matter at a finite baryon density can have a different type of phase transition, possibly first order, during which the speed of sound approaches $c_{\text{s}}^2\approx 0$.

Nuclear EoS is essential for describing the dynamic evolution of dense and hot nuclear matter created in heavy-ion collisions because the primary driving force for fluid expansion is the pressure gradient. It is also a vital input for determining the mass-radius relation of neutron stars, as a soft EoS cannot support immense gravitational forces. Nuclear EoS encodes QCD phase transitions, which can be used to identify neutrons and quark stars. Furthermore, the nuclear EoS reveals the emergent properties of numerous nucleons and provides constraints on the nucleon-nucleon interaction, which is crucial for few-nucleon systems such as atomic nuclei.

In a recent publication in Physical Review Letters [1], Manjunath, Steinheimer, Zhou and Stöcker from Goethe University Frankfurt and GSI reported their study on Bayesian analysis to constrain the nuclear EoS in dense regions using data from intermediate-energy heavy-ion collisions. They constructed a density-dependent mean-field potential in UrQMD to simulate heavy-ion collisions at 2-10 GeV. Using published data from heavy-ion collisions, Bayesian analysis was then used to constrain the interaction potential and nuclear EoS. The extracted EoS for nucleon densities ranging from 2n₀ to 4n₀ agrees well with that given by astrophysical observations; the EoS beyond 3n₀ requires more accurate data with higher statistics from the beam energy scan program of STAR-FXT at RHIC, the upcoming CBM experiment at FAIR, and future experiments at HIAF and NICA.

No first-principles calculations are available for the nuclear Equation of State (EoS) at high baryon densities. Although lattice QCD calculations provide the QCD EoS at the vanishing baryon chemical potentials, they cannot accurately predict the EoS at high net baryon densities because of the notorious fermionic sign problem [2]. To investigate the nuclear EoS at high baryon densities and possible critical endpoint between the two phases [3, 4], physicists have turned to a data-driven approach.

The beam energy scan (BES) program at RHIC was designed to accumulate data on heavy-ion collisions at various beam energies with different baryon stopping powers and net baryon chemical potentials in the central rapidity region of the final state. Theoretical approaches, such as the QCD chiral effective field theory (EFT), transport model simulations, and relativistic hydrodynamics, have been employed to extract the nuclear EoS from experimental data on heavy-ion collisions, neutron stars, gravitational waves, and nuclear incompressibility. Fast-moving objects traveling through a medium at speeds greater than the speed of sound create Mach cones. Recent studies also explored the potential use of jet-induced 3D Mach cones to determine the nuclear equation of state [5]. For comparisons between data and theoretical calculations, advanced machine learning methods such as Bayesian analyses [6-10] and deep learning [4, 11-17] are used. For further details, please refer to the following reviews [18-23].

To describe the EoS of dense nuclear matter above two times of the nuclear saturation density 2n₀, the authors of Ref. [1] use a 7-th order polynomial function to parameterize the density dependence of baryon potential energy V(n_B(r)). The seven coefficients in the polynomial function form the parameter space $\theta =\left[{\theta }_1\mathrm{,}{\theta }_2\mathrm{,}\cdots \mathrm{,}{\theta }_7\right]$ to be extracted using Bayesian analysis. The dataset on elliptic flow v₂ and the mean transverse kinetic energy $\langle m_T\rangle -m_0$ of protons at the mid-rapidity of Au+Au collisions at $\sqrt{s}=2-10$ GeV, referred to as $D=\left[v_2\mathrm{,}\langle m_T\rangle -m_0\right]$, were used to constrain the EoS. The goal of their Bayesian analysis was to determine the posterior distributions of the parameters in the polynomial function, given the dataset $P(\theta |D)\propto P(\theta )P(D|\theta )$ The posterior distribution was not normalized. The Markov Chain Monte Carlo (MCMC) method is typically employed to obtain samples from this unnormalized distribution by performing random walks in the parameter space. For instance, the Metropolis-Hastings (MH) algorithm begins at a randomly chosen location ${\theta }_{\text{step}0}$ in the parameter space and proposes the next step ${\theta }_{\text{step}1}={\theta }_{\text{step}0}+\delta \theta$. The algorithm then determines whether to accept this new position based on posterior distribution P(θ | D). If the posterior distribution exhibits a higher probability density at the proposed position than at the current position, ${\theta }_{\text{step}1}$ is accepted as a new sample. Otherwise, it is accepted with a probability ratio between the new and old positions. If ${\theta }_{\text{step}1}$ is rejected, the old position ${\theta }_{\text{step}0}$ is duplicated as a new sample. This process was implemented in PyMC and EMCEE. This study uses PyMC4.0 to sample θ from the posterior distribution P(θ | D). Density estimation provides the distribution of θ, from which one can identify the mean values of θ or parameter combinations that maximize P(θ | D), known as the Maximum A Posteriori (MAP).

In this study, the key equation that relates the potential V(n_B(r)) to the equation of state (pressure) is given by [24] $P\left(n_B\right)=P_{\text{id}}\left(n_B\right)+{\displaystyle {\int }_0^{n_B}}{n}^{\text{'}}\frac{\partial U\left(n^{\text{'}}\right)}{\partial n^{\text{'}}}\text{d}{n}^{\text{'}}$ where $P_{\text{id}}$ is the pressure of the ideal fermionic gas and $U\left(n_B\right)=\frac{\partial \left(n_B\cdot V\left(n_B\right)\right)}{\partial n_B}$. Once the pressure is known, the square of the speed of sound $c_s^2$ in dense nuclear matter can be computed using its derivatives.

As shown in Fig. 1, the extracted $c_s^2$ values using 13 and 15 data points from heavy-ion collisions are compared with those obtained from constraints by Binary Neutron Star Mergers (grey band) [25] and the Chiral Mean Field (CMF) model calculations (purple curve) [24]. The $c_s^2$ values extracted using the Bayesian analysis agreed with the Binary Neutron Star Merger constraints when 15 data points were used. However, reducing the number of data points to 13 leads to a rapid change in $c_s^2$, indicating that the EoS is sensitive to the observables selected in the Bayesian analysis. Therefore, higher statistical data from experimental measurements are required to constrain nuclear EoS in dense regions better.

Fig. 1Full-screen View

The square of the speed of sound $c_s^2$ at the density of 2n₀ - 6n₀, where n₀ is the nuclear saturation density, is extracted from comparisons between experimental data and UrQMD simulations of heavy-ion collisions at the beam energies of 2-10 GeV in Ref. [1]. The Bayesian analysis was employed to constrain the EoS, and the results agree with neutron star merger observations. Notably, the results are sensitive to the number of data points used in the constraining process.

The performance of the Bayesian analysis was verified using closure tests. In these tests, a random potential function was incorporated into UrQMD to compute the physical observables and Bayesian analysis was employed to obtain the posterior distribution of θ in the potential function. The ground-truth potential function is centered around the resulting V(n_B(r)). Closure tests are crucial in Bayesian analysis because many scientific problems do not have well-defined inverse solutions. For example, if two or more potential functions yield the same v₂ and $\langle m_{\text{T}}\rangle -m_0$, then it is impossible to retrieve the used V(n_B(r)). This does not imply that there is no way to distinguish between different EoS but suggests that additional physical observables may be necessary. In this study, the closure tests indicated that the physical observables selected by the authors were sensitive to changes in the EoS.

Constraining the EoS at the high net baryon density is essential for understanding the dynamic evolution of dense nuclear matter during heavy-ion collisions. This knowledge is also crucial for providing key inputs to obtain the mass-radius relationship of neutron stars and interpreting neutron star merger observations. Heavy-ion collisions at various intermediate and low beam energies offer a unique opportunity to study the nuclear EoS under extreme conditions in a laboratory setting [26, 27]. Furthermore, this research can help constrain the fundamental nucleon-nucleon interactions in many-body quantum systems. The incorporation of advanced statistical tools and machine learning methods into data theory comparisons in nuclear physics will accelerate the advancement of this frontier.