Bayesian inference of nuclear incompressibility from collective flow in mid-central Au+Au collisions at 400–1500 MeV/nucleon

The equation of state (EoS) of cold symmetric nuclear matter (SNM) is a fundamental relationship between the energy per nucleon $ E/A $ and nucleon density ρ. At the saturation density $ \rho_0 = 0.16 $ $ \mathrm{fm}^{-3} $ of SNM, where $ E/A = -16.0 $ MeV and pressure vanishes, the stiffness of the SNM EoS is measured by its incompressibility $ K = 9\rho ^{2}\dfrac{\partial ^{2}E/A }{\partial \rho ^{2} }\mid_{\rho _{0} } $. Precisely determining the value of K has been a longstanding and shared goal of both nuclear physics and astrophysics for the last few decades, due to its strong impact on many aspects of nuclear structure and reactions, as well as the properties of neutron stars, mechanisms of supernovae explosions, and emissions of gravitational waves from mergers of neutron stars (see, e.g., Refs. [1−20]). In fact, thanks to the hard work of many researchers, especially in the last 40 years, much progress has been made in constraining the value of K. In particular, since the pioneering work of Blaizot, who determined K = (210 $ \pm $ 30) MeV by analyzing the giant monopole resonance (GMR) energies in ⁴⁰Ca, ⁹⁰Zr, and ²⁰⁸Pb [1], extensive theoretical studies and systematic measurements of GMR energies [1, 21−26] have led to the community consensus that the K has values in the range of 220–260 MeV [22, 23, 26, 27] or around $ 235\pm 30 $ MeV [28, 29], as indicated by the light blue band in Fig. 1.

Figure 1.  (color online) The light blue band covering $ K = 220-260 $ MeV is the fiducial range from studying giant resonances [1, 22, 23]. The other ranges, as indicated by the legends, represent Refs. [41−49] from left to right. The red points with error bars are the results of K at $ E = 400 $ MeV/nucleon from this work.

Another well-known tool for studying the EoS of nuclear matter is heavy-ion collisions. In particular, various components of nuclear collective flow (e.g., directed and elliptic flow), as well as yields and spectra of multiple particles (e.g., pions, kaons, photons, and dileptons), are useful probes of the EoS, albeit with different sensitivities [30−32]. Information about the EoS is normally extracted by comparing transport model simulations with experimental data in the forward-modeling approach, while in more recent years, data-driven Bayesian inference has become more revealing with quantified uncertainties. In either approach, there are still many interesting issues in narrowing down the remaining uncertainties of various features and parameters of the nuclear EoS using heavy-ion reactions; see, e.g., Refs. [33−35] for reviews. In recent works, heavy-ion collisions have also been widely used to study the EoS of nuclear matter at higher baryon density [36−40].

The purpose of this work is twofold. First, uncertainties remain in the values of K derived both from analyses of GMR energies and certain transport model studies of heavy-ion reactions. For example, as shown by the first point in Fig. 1, analyses of the elliptic flow in Au+Au reactions at beam energies from 0.4 to 1.5 GeV/nucleon from the FOPI collaboration within the Isospin-dependent Quantum Molecular Dynamics (IQMD) model extracted a value of $ K = 190^{+30}_{-30} $ MeV [41]. Furthermore, one analysis of the ratio of kaons in Au+Au over C+C systems within hadronic transport simulations favored a soft EoS with $ K = 200 $ MeV over the hard EoS with $ K = 380 $ MeV [42]. Both studies considered the momentum dependence of isoscalar single-nucleon potentials. For other results on the nuclear incompressibility K, see Fig. 1. It should be noted that the incompressibility K in Refs. [41−45] was obtained from transport models of heavy-ion reactions in which $ \rho > \rho _{0} $. Moreover, in Refs. [46−49], it was derived from studies of the GMR, which is $ \rho \approx \rho _{0} $. The uncertainty of the approximately 40 MeV error band of K from studying GMR and the K from analyzing heavy-ion reactions are still too large for many investigations in both nuclear physics and astrophysics. For example, it has been shown recently that a variation of K between 220 and 260 MeV leads to significant changes in the crust–core transition density and pressure in neutron stars [50, 51]. This subsequently affects the radii, tidal polarizabilities, and crustal moments of inertia of canonical neutron stars, hindering the investigations of the very mysterious high-density cores of neutron stars, as their properties are often strongly correlated with the value of K. Thus, a re-extraction of K with a quantified uncertainty from heavy-ion collisions through a Bayesian approach will be invaluable. Considering that the momentum dependence of single-nucleon potentials is very important for accurately simulating heavy-ion reactions, we have inferred the value of K with or without considering the momentum dependence of single-nucleon potentials.

Second, we investigate the inference of K from different observables under different beam energies. By studying the combined data of rapidity and transverse velocity dependence of proton elliptic flow in mid-central Au+Au reactions at 400 MeV/nucleon, we found an incompressibility of $ K = 188.9^{+2.9}_{-4.5} $ or $ 256.1^{+8.2}_{-8.7} $ MeV at 68% confidence level with or without considering the momentum dependence of single-nucleon potentials, respectively. As shown in Fig. 1, these results are all consistent with the earlier findings of Refs. [41−43]. In contrast, by incorporating the transverse momentum dependence of proton-like directed flow, the inferred K becomes $ 222.3^{+9.0}_{-9.9} $ MeV or $ 285.5^{+6.7}_{-7.3} $ MeV at 68% confidence level with or without considering the momentum dependence of single-nucleon potentials, respectively. Obviously, these results are significantly different from those inferred earlier, indicating that the difference of the extracted EoS is due to the different density regions probed by different observables. While we have not extensively studied all observables in a broad beam energy range yet, the findings reported here are useful for unraveling the physics underlying the incompressibility of nuclear matter and the constraining powers of different observables useful for eventually pinning down the value of K. Furthermore, through the analysis of proton elliptic flow in mid-central Au+Au reactions at beam energies of $ 400-1500 $ MeV/nucleon, with consideration of the momentum dependence of single-nucleon potentials, we found that the value of K increases with beam energy, reflecting the hardening of nuclear matter EoS at higher densities and temperatures in reactions with higher beam energies.

The remainder of this article is organized as follows. In Sec. II, we outline the key ingredients of the IQMD model most relevant for this study. In Sec. III, we recall the key aspects of the Bayesian approach used in this work. In Sec. IV, we present the Gaussian Process (GP) emulator of IQMD and demonstrate the accuracy of its training and testing. The posterior probability distribution functions (PDFs) of the incompressibility K, as well as the corresponding nuclear interaction parameters from our Bayesian inference in different situations, are presented and discussed in Sec. V. In Sec. VI, we summarize the results of this study.

For ease of understanding and completeness, we recall here some key physical ingredients of the IQMD model that are most relevant to this work. This model has been widely and successfully applied to the study of heavy ion collisions at low to intermediate beam energies (see, e.g., Refs. [52−59]). The nuclear effective interactions used in our QMD model include not only a Skyrme term, but also the Yukawa, isospin asymmetric, and Coulomb terms. Optionally, a momentum-dependent interaction (MDI) term can be switched on or off. More specifically, the EoS, pressure, and incompressibility K for cold SNM can be written as follows [52−54, 60−62]:

where α, β, γ, $ t_{4} $, and $ t_{5} $ are parameters, $ P_{F} = \left(3\pi^{2} \hbar^{3} \rho/2\right)^{1/3} $ is the nucleon Fermi momentum at density ρ, $ f\left(\boldsymbol{p}\right ) = \dfrac{3}{4\pi P_{F}^{3} }\Theta \left ( p _{F}-p \right ) $, and Θ is the step function. The values of $ t_{4} = 1.57 $ MeV and $ t_{5} = 5\times 10^{-4} $ $ \mathrm{c}^{2}/\mathrm{MeV}^{2} $ are determined by the experimental nucleon optical potential and isoscalar nucleon effective mass at $ \rho_0 $ [52]. The EoS and its characteristics depend only on α, β, and γ parameters. In this work, default free-space nucleon-nucleon cross-section $ \sigma_{NN}^{\mathrm{free}} $ obtained from experimental data is used for the simulations [63, 64].

For ease of understanding and completeness, we recall here some key physical ingredients of the IQMD model that are most relevant to this work. This model has been widely and successfully applied to the study of heavy ion collisions at low to intermediate beam energies (see, e.g., Refs. [52−59]). The nuclear effective interactions used in our QMD model include not only a Skyrme term, but also the Yukawa, isospin asymmetric, and Coulomb terms. Optionally, a momentum-dependent interaction (MDI) term can be switched on or off. More specifically, the EoS, pressure, and incompressibility K for cold SNM can be written as follows [52−54, 60−62]:

where α, β, γ, $ t_{4} $, and $ t_{5} $ are parameters, $ P_{F} = \left(3\pi^{2} \hbar^{3} \rho/2\right)^{1/3} $ is the nucleon Fermi momentum at density ρ, $ f\left(\boldsymbol{p}\right ) = \dfrac{3}{4\pi P_{F}^{3} }\Theta \left ( p _{F}-p \right ) $, and Θ is the step function. The values of $ t_{4} = 1.57 $ MeV and $ t_{5} = 5\times 10^{-4} $ $ \mathrm{c}^{2}/\mathrm{MeV}^{2} $ are determined by the experimental nucleon optical potential and isoscalar nucleon effective mass at $ \rho_0 $ [52]. The EoS and its characteristics depend only on α, β, and γ parameters. In this work, default free-space nucleon-nucleon cross-section $ \sigma_{NN}^{\mathrm{free}} $ obtained from experimental data is used for the simulations [63, 64].

Bayesian inference is a probabilistic approach employed in Machine Learning (ML) that has been widely utilized in many fields in recent years. Compared to the traditional $ \chi^{2} $ fitting in the forward-modeling approach, it has several advantages in determining model parameters with quantified uncertainties (see, e.g., Refs. [39, 45, 49, 65−81]. For completeness, we recall that Bayes' theorem states

where $ P\left(\theta \mid D\right) $ denotes the posterior probability density function (PDF) of model parameter set θ given the dataset D, while $ P\left(D \mid \theta\right) $ is the likelihood for the model with the parameter set θ to reproduce the dataset D. $ P\left(\theta\right) $ is the prior PDF of the parameter set θ, and $ P\left(D\right) $ serves as the normalization constant. In this study, θ consists of incompressibility K, α, β, and γ constrained by the three SNM saturation conditions at $ \rho_0 $ mentioned earlier. Thus, only one parameter is free, and we choose it to be K. In the Markov Chain Monte Carlo (MCMC) process of our Bayesian analysis, the trial K value is randomly generated uniformly in the prior range of $ K = 170- 420 $ MeV. Once the incompressibility K is selected, the corresponding values of α, β, and γ are then obtained from the saturation conditions. Additionally, two constraints of $ -1000<\alpha <0 $ MeV and $ \beta >0 $ MeV are applied to ensure nuclear matter remains stable and casual at all densities, which in turn raises the lower bound of the prior distribution of K to above 170 MeV. This can be seen in the prior distribution of K shown in Sec. V. Thus, there is a one-to-one correspondence between the incompressibility K and α, β, γ parameters considering the saturation conditions of SNM at $ \rho_0 $ through Eqs. (1)–(3). By inputting a given set of parameters α, β, and γ into the IQMD model, theoretical predictions for the reaction observables corresponding to the selected EoS parameter set can be obtained. These predictions will then be used to evaluate the likelihood function, as we discuss next.

Taking the logarithm of the posterior distribution of parameters enables a transformation of the Bayesian formula [39]:

The logarithm of the likelihood function is evaluated using

In the above, $ \sigma _{i} = \sqrt{\sigma_{\rm exp}^{2} +\sigma_{\rm mod}^{2}} $, where $ \sigma_{\rm exp} $ and $ \sigma_{\rm mod} $ (evaluated from the emulator) represent the experimental and model errors, respectively. $ y _{i}^{\theta } $ and $ y_{i,\rm exp} $ are the model prediction for the observable $ y_i $ and its experimental value, respectively.

Bayesian inference is a probabilistic approach employed in Machine Learning (ML) that has been widely utilized in many fields in recent years. Compared to the traditional $ \chi^{2} $ fitting in the forward-modeling approach, it has several advantages in determining model parameters with quantified uncertainties (see, e.g., Refs. [39, 45, 49, 65−81]. For completeness, we recall that Bayes' theorem states

where $ P\left(\theta \mid D\right) $ denotes the posterior probability density function (PDF) of model parameter set θ given the dataset D, while $ P\left(D \mid \theta\right) $ is the likelihood for the model with the parameter set θ to reproduce the dataset D. $ P\left(\theta\right) $ is the prior PDF of the parameter set θ, and $ P\left(D\right) $ serves as the normalization constant. In this study, θ consists of incompressibility K, α, β, and γ constrained by the three SNM saturation conditions at $ \rho_0 $ mentioned earlier. Thus, only one parameter is free, and we choose it to be K. In the Markov Chain Monte Carlo (MCMC) process of our Bayesian analysis, the trial K value is randomly generated uniformly in the prior range of $ K = 170- 420 $ MeV. Once the incompressibility K is selected, the corresponding values of α, β, and γ are then obtained from the saturation conditions. Additionally, two constraints of $ -1000<\alpha <0 $ MeV and $ \beta >0 $ MeV are applied to ensure nuclear matter remains stable and casual at all densities, which in turn raises the lower bound of the prior distribution of K to above 170 MeV. This can be seen in the prior distribution of K shown in Sec. V. Thus, there is a one-to-one correspondence between the incompressibility K and α, β, γ parameters considering the saturation conditions of SNM at $ \rho_0 $ through Eqs. (1)–(3). By inputting a given set of parameters α, β, and γ into the IQMD model, theoretical predictions for the reaction observables corresponding to the selected EoS parameter set can be obtained. These predictions will then be used to evaluate the likelihood function, as we discuss next.

Taking the logarithm of the posterior distribution of parameters enables a transformation of the Bayesian formula [39]:

The logarithm of the likelihood function is evaluated using

In the above, $ \sigma _{i} = \sqrt{\sigma_{\rm exp}^{2} +\sigma_{\rm mod}^{2}} $, where $ \sigma_{\rm exp} $ and $ \sigma_{\rm mod} $ (evaluated from the emulator) represent the experimental and model errors, respectively. $ y _{i}^{\theta } $ and $ y_{i,\rm exp} $ are the model prediction for the observable $ y_i $ and its experimental value, respectively.

Because it is impractical to invoke the computationally expensive IQMD simulator in Bayesian analyses, its emulator must be used instead. Here, we use the popular and well-tested GP emulator. Below, we provide some details on the training and testing of the emulator.

In this study, we use four observables in Au + Au collisions at $ E = 400 $ MeV/nucleon and two observables at $ E = 600-1500 $ MeV/nucleon, namely, the rapidity dependence of $ -v_{2}\left(y_{0}\right) $ of proton elliptic flow for $ u_{t0}>0.4 $, transverse velocity dependence of $ -v_{2}\left(u_{t0}\right) $ of proton elliptic flow for $ \left | y_{0} \right | <0.4 $ [82], and transverse momentum dependence of $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ of proton-like (all charged particles with $ Z = 1 $ and $ Z = 2 $ are weighted by the charge Z) directed flow in two rapidity ranges of $ y_{0} = 0.5-0.7 $ and $ y_{0} = 0.7-0.9 $ [83]. For proton elliptic flow, the centrality $ 0.25 < b_{0} < 0.45 $ is considered, where $ b_{0} $ denotes the scaled impact parameter defined as $ b_{0} = b/b_{\text{max}} $, with $ b_{\text{max}} = 1.15(A_{P}^{1/3} + A_{T}^{1/3}) $ fm. The range of impact parameter $ 1.9-6.1 $ fm is considered for proton-like directed flow. Here, $ v_{2} = \langle\cos 2\phi\rangle = \langle (p_{x}^{2} - p_{y}^{2}) / (p_{x}^{2} + p_{y}^{2}) \rangle $. The reduced rapidity $ y_{0} $ is defined as $ y/y_{\text{pro}} $, where $ y = \dfrac{1}{2}\ln\left(\dfrac{E+p_{z}}{E-p_{z}}\right) $, and ‘pro’ denotes the incident projectile in the center of mass frame. The transverse (spatial) component $ u_{t} $ of the 4-velocity u is given by $ u_{t} = \beta_{t} \gamma _{c} $. Here, the 3-vector $ \boldsymbol{\beta} $ represents the velocity in units of the speed of light, and $ \gamma _{c} = 1/\sqrt{1-\beta^{2}} $. The unitless transverse velocity $ u_{t0} = u_{t} / u_{\text{pro}} $, where $ u_{\text{pro}} = \beta_{\text{pro}} \gamma_{\text{pro}} $ is used. For the directed flow, $ v_{1} = \langle\cos \phi\rangle = \langle p_{x} / \sqrt{p_{x}^{2} + p_{y}^{2}} \rangle $, we study its dependence on the normalized center-of-mass (c.m.) transverse momentum (per nucleon), defined as $ p_{t}^{\left ( 0 \right ) } = \left ( p_{t}/A \right ) / \left ( p_{P}^{\mathrm{c.m.}}/A_{P} \right ) $. To understand the role of observables in the constrained result of K, we perform comparative studies using two datasets: one group utilizes only the rapidity dependence of proton elliptic flow $ -v_{2}\left(y_{0}\right) $ and its transverse velocity dependence $ -v_{2}\left(u_{t0}\right) $, while the other incorporates the transverse momentum dependence of $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ of proton-like directed flow in two rapidity ranges. Additionally, we utilize only the observables of proton elliptic flow $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ to investigate the beam energy dependence of the constrained results for the incompressibility K.

The initial value of incompressibility K is selected using Latin hypercube sampling (LHS) [39, 45] to cover its entire prior range uniformly and efficiently, enabling the emulator to efficiently learn and reliably predict reaction observables. In generating the training and testing sets with the IQMD simulator, 220 different K values form the training set, while an additional 50 different K values are used for testing purposes. For each incompressibility K, the IQMD model produces 100,000 events to calculate the respective observables $ -v_{2}\left(y_{0}\right) $, $ -v_{2}\left(u_{t0}\right) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $, with each event being associated with a randomly selected impact parameter in the range mentioned earlier.

To facilitate effective learning and reliable predictions of observables by the emulator, the commonly employed Radial Basis Function (RBF) is utilized [84]. To assess the emulator's performance on the testing set, Fig. 2 shows a visual comparison between model calculations (x axis) by the IQMD simulator and the GP emulator predictions with error bars (y axis), where the blue line is the line of perfect match between the GP predictions and IQMD simulations. As evident from Fig. 2, the error bars given by the GP emulator are near or over the blue line, so the GP emulator can basically accurately predict the observables. Clearly, for all the observables in the cases considered, the emulator's performance is satisfactory.

Figure 2.  (color online) Comparison between the predictions for $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ by the IQMD simulator and GP emulator for two randomly selected sets of incompressibility K in the testing dataset. Top: without considering the MDI. Bottom: considering the MDI.

Because it is impractical to invoke the computationally expensive IQMD simulator in Bayesian analyses, its emulator must be used instead. Here, we use the popular and well-tested GP emulator. Below, we provide some details on the training and testing of the emulator.

In this study, we use four observables in Au + Au collisions at $ E = 400 $ MeV/nucleon and two observables at $ E = 600-1500 $ MeV/nucleon, namely, the rapidity dependence of $ -v_{2}\left(y_{0}\right) $ of proton elliptic flow for $ u_{t0}>0.4 $, transverse velocity dependence of $ -v_{2}\left(u_{t0}\right) $ of proton elliptic flow for $ \left | y_{0} \right | <0.4 $ [82], and transverse momentum dependence of $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ of proton-like (all charged particles with $ Z = 1 $ and $ Z = 2 $ are weighted by the charge Z) directed flow in two rapidity ranges of $ y_{0} = 0.5-0.7 $ and $ y_{0} = 0.7-0.9 $ [83]. For proton elliptic flow, the centrality $ 0.25 < b_{0} < 0.45 $ is considered, where $ b_{0} $ denotes the scaled impact parameter defined as $ b_{0} = b/b_{\text{max}} $, with $ b_{\text{max}} = 1.15(A_{P}^{1/3} + A_{T}^{1/3}) $ fm. The range of impact parameter $ 1.9-6.1 $ fm is considered for proton-like directed flow. Here, $ v_{2} = \langle\cos 2\phi\rangle = \langle (p_{x}^{2} - p_{y}^{2}) / (p_{x}^{2} + p_{y}^{2}) \rangle $. The reduced rapidity $ y_{0} $ is defined as $ y/y_{\text{pro}} $, where $ y = \dfrac{1}{2}\ln\left(\dfrac{E+p_{z}}{E-p_{z}}\right) $, and ‘pro’ denotes the incident projectile in the center of mass frame. The transverse (spatial) component $ u_{t} $ of the 4-velocity u is given by $ u_{t} = \beta_{t} \gamma _{c} $. Here, the 3-vector $ \boldsymbol{\beta} $ represents the velocity in units of the speed of light, and $ \gamma _{c} = 1/\sqrt{1-\beta^{2}} $. The unitless transverse velocity $ u_{t0} = u_{t} / u_{\text{pro}} $, where $ u_{\text{pro}} = \beta_{\text{pro}} \gamma_{\text{pro}} $ is used. For the directed flow, $ v_{1} = \langle\cos \phi\rangle = \langle p_{x} / \sqrt{p_{x}^{2} + p_{y}^{2}} \rangle $, we study its dependence on the normalized center-of-mass (c.m.) transverse momentum (per nucleon), defined as $ p_{t}^{\left ( 0 \right ) } = \left ( p_{t}/A \right ) / \left ( p_{P}^{\mathrm{c.m.}}/A_{P} \right ) $. To understand the role of observables in the constrained result of K, we perform comparative studies using two datasets: one group utilizes only the rapidity dependence of proton elliptic flow $ -v_{2}\left(y_{0}\right) $ and its transverse velocity dependence $ -v_{2}\left(u_{t0}\right) $, while the other incorporates the transverse momentum dependence of $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ of proton-like directed flow in two rapidity ranges. Additionally, we utilize only the observables of proton elliptic flow $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ to investigate the beam energy dependence of the constrained results for the incompressibility K.

The initial value of incompressibility K is selected using Latin hypercube sampling (LHS) [39, 45] to cover its entire prior range uniformly and efficiently, enabling the emulator to efficiently learn and reliably predict reaction observables. In generating the training and testing sets with the IQMD simulator, 220 different K values form the training set, while an additional 50 different K values are used for testing purposes. For each incompressibility K, the IQMD model produces 100,000 events to calculate the respective observables $ -v_{2}\left(y_{0}\right) $, $ -v_{2}\left(u_{t0}\right) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $, with each event being associated with a randomly selected impact parameter in the range mentioned earlier.

To facilitate effective learning and reliable predictions of observables by the emulator, the commonly employed Radial Basis Function (RBF) is utilized [84]. To assess the emulator's performance on the testing set, Fig. 2 shows a visual comparison between model calculations (x axis) by the IQMD simulator and the GP emulator predictions with error bars (y axis), where the blue line is the line of perfect match between the GP predictions and IQMD simulations. As evident from Fig. 2, the error bars given by the GP emulator are near or over the blue line, so the GP emulator can basically accurately predict the observables. Clearly, for all the observables in the cases considered, the emulator's performance is satisfactory.

Figure 2.  (color online) Comparison between the predictions for $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ by the IQMD simulator and GP emulator for two randomly selected sets of incompressibility K in the testing dataset. Top: without considering the MDI. Bottom: considering the MDI.

The Affine Invariant MCMC Ensemble sampler algorithm is employed to sample the posterior distribution of incompressibility K. The cumulative mean (running average) plot from the MCMC sampling is typically employed to assess the convergence status of the sampling process. In this work, each of the five MCMC chains independently runs 3 million steps, discarding the first 1.2 million as burn-in steps.

As mentioned earlier, we compare two calculations with and without considering the momentum-dependent interaction (MDI) part of single-nucleon potentials at $ E = 400 $ MeV/nucleon. In both cases, we found that 1 million burn-in steps are sufficient. As an example, Fig. 3(a) and 3(b) show the cumulative mean diagrams of our MCMC sampling without considering the MDI. It can be seen that regardless of the initial states, all chains surely converge to the same equilibrium state after approximately 1 million steps. Moreover, this conclusion is independent of the observables we used. All results presented in the following are obtained by using 3 million steps, discarding the first 1.2 million as burn-in steps.

Figure 3.  (color online) Without considering the MDI: running averages of K as functions of the MCMC steps after 200,000 initial steps in 5 independent chains. Left: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $ Right: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

We first examine the posterior PDFs without considering the MDI. The results obtained from the convergent MCMC chains are shown in Fig. 4(a) (using only $ -v_{2}\left (y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $) and 4(b) (using $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right )} \right ) $). We note that the white dashed line corresponds to the most probable value (MPV) of the posterior distribution. The black dashed line represents the prior distribution of K. In both cases, the posterior PDF is significantly narrower than the prior, demonstrating that the chosen observables provide strong constraints on K.

Figure 4.  (color online) Without considering the MDI: posterior PDFs of K. Left: using the observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Right: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

It is important to note that here we adopt the highest posterior density (HPD) method to calculate the 68% confidence interval [85] because the posterior distribution of K sometimes deviates from a Gaussian distribution. More specifically, we find from the PDF shown in Fig. 4(a) that the mean value of K is $ K = 256.1^{+8.2}_{-8.7} $ MeV at 68% confidence interval when only using the observables $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left ( u_{t0} \right ) $. In contrast, from Fig. 4(b), we find $ K = 285.5^{+6.7}_{-7.3} $ MeV when $ -v_{2}\left(y_{0}\right) $, $ -v_{2}\left(u_{t0}\right) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ are used. In this case, the mean or MPV of K is appreciably larger than that when only the elliptic proton flow data $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ are used. This demonstrates that different combinations of observables, including the proton elliptic flow and proton-like directed flow, impose distinct constraints on the incompressibility K, as discussed in the following.

As previously discussed, given an incompressibility K, the corresponding parameters α, β, and γ can be derived using the equations for $ E/A $, P, and K at saturation density $ \rho_0 $. Utilizing the convergent MCMC chains, we can readily obtain the posterior PDFs of α, β, and γ. As shown in Fig. 5(a)–(f), these posterior PDFs are significantly narrower than their prior distributions, indicating tighter constraints. Such constrained parameters are valuable for testing fundamental interactions, such as the Skyrme nuclear effective interactions.

Figure 5.  (color online) Without considering the MDI: posterior PDFs of α, β, and γ. Top: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Bottom: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

We now turn to the results obtained when incorporating the MDI. We maintain the same prior distributions for incompressibility K and apply identical constraining conditions as in the previous case. The PDFs for K are depicted in Fig. 6(a) and 6(b). The corresponding results for α, β, and γ are illustrated in Fig. 7(a)–7(c) and 7(d)–7(f), respectively. Notably, the mean values of K are now $ 188.9^{+2.9}_{-4.5} $ MeV at 68% confidence level when only the observables $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ are used, and they are $ 222.3^{+9.0}_{-9.9} $ MeV when the observables $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ are incorporated. These results reveal a clear dependence of the inferred K on the choice of observables, similar to the earlier findings. Furthermore, the values of K obtained when including the MDI are significantly smaller than those inferred without considering the MDI. To reproduce the same observables, an equivalent gradient of mean-field potential is required, which can originate from either the density-dependent or momentum-dependent components of the mean-field potential. The MDI potential is repulsive, which provides a larger gradient of mean-field potential for particles without making the density-dependent part of the EoS stiffer [86]. However, the absence of a momentum-dependent potential in the simulation requires compensation from a stiffer density-dependent EoS to reproduce the same observables. Thus, in the case of including the momentum-dependent potential in simulating heavy-ion reactions, regardless of using only $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ or including $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $, the required K values are smaller than without considering the MDI. The results here are consistent with the earlier findings that the MDI reduces the necessary K to reproduce flow data (see, e.g., Refs. [41, 42, 87−89]). However, we acknowledge that the extracted values of K are still subject to uncertainties arising from model dependence. As discussed extensively in the Transport Model Evaluation Project (TMEP) review [34] and references therein, there are systematic differences in various transport models.

Figure 6.  (color online) Considering the MDI: posterior PDFs of K. Left: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Right: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

Figure 7.  (color online) Considering the MDI: posterior PDFs of α, β, and γ. Top: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Bottom: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

The results presented earlier were obtained using proton elliptic flow and proton-like directed flow in mid-central Au+Au collisions at 400 MeV/nucleon. For a systematic study, we now infer the incompressibility solely based on the rapidity and transverse velocity dependence of proton elliptic flow in mid-central Au+Au collisions at beam energies from 600 to 1500 MeV/nucleon, with the momentum-dependent interaction taken into account, maintaining the same prior distribution of incompressibility K and constraining conditions as before. The beam energy dependences of the PDFs of K are illustrated in Fig. 8(a)–8(e). To be more quantitative, the mean values of K as functions of beam energy are shown in Fig. 8 (f), and the values of K are given in Table 1. The incompressibility K increases with beam energy, regardless of whether the observable $ -v_{2}\left(y_{0}\right) $ or $ -v_{2}\left(u_{t0}\right) $ is used. From Eq. (3), it is clear that the incompressibility K is defined at the saturation density $ \rho_0 $, and in principle, K as a property of nuclear matter at $ \rho_0 $ should not depend on the beam energy or specific observable. However, in transport model simulations where the EoS is characterized by a single parameter K, the observable-sensitive density regime probed by heavy-ion collisions increases with beam energy. Therefore, the K values extracted from data at different beam energies effectively reflect the stiffness of the EoS at different densities, rather than the incompressibility at $ \rho_0 $. The observed increasing trend of extracted K values with beam energy should thus be interpreted as an indication that the EoS becomes stiffer at higher densities and temperatures in the energy region from 400 to 1500 MeV/nucleon, which is consistent with the conclusion from the SMASH model in Ref. [90]. The general trend also aligns with Li's recent work [91], despite differences in the models and observables used. Meanwhile, it is opposite to the conclusion that gives a soft EoS with elliptic flow in the same energy region in Le Févre's work [41, 92]. The observable used in this work differs from that in Refs. [41, 92]. Moreover, it should be noted that the description of the $ u_{t0} $ cut condition adopted from Ref. [82] in Ref. [92] is inaccurate, as the descriptions of the $ u_{t0} $ cut in Refs. [82] and [93] are inconsistent. This inconsistency may affect the rigor of the conclusions drawn in Ref. [92]. Notably, the values of K derived from the observables $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ exhibit discrepancies, which arise from the differences in the rapidity and transverse velocity dependence of proton elliptic flow. The reason here is similar to that for the inferred values of K at $ E = 400 $ MeV/nucleon: when different observables probe different density regions, the extracted EoS will differ accordingly. Specifically, elliptic flow may reflect the behavior of particles at low densities, while directed flow might be more sensitive to high-density behavior. As a result, the EoS extracted from proton elliptic flow tends to be softer, whereas including directed flow leads to a stiffer EoS, as illustrated in Figs. 4 and 6.

Figure 8.  (color online) Considering the MDI: posterior PDFs of K at $ E = 600-1500 $ MeV/nucleon and its dependence on beam energy. Observables: $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $.

	noMDI ($ E = 400 $ MeV/nucleon)		MDI ($ E = 400 $ MeV/nucleon)
Observables:	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ) $	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ), $	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ) $	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ), $
		$ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $		$ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $
K Mean	$ 256.1^{+8.2}_{-8.7} $	$ 285.5^{+6.7}_{-7.3} $	$ 188.9^{+2.9}_{-4.5} $	$ 222.3^{+9.0}_{-9.9} $
K MPV	255.5	284.6	187.9	222.3
MDI (observables:$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ) $)
600 MeV/nucleon	800 MeV/nucleon	1000 MeV/nucleon	1200 MeV/nucleon	1500 MeV/nucleon
$ 208.8^{+6.7}_{-7.3} $	$ 271.8^{+7.8}_{-8.7} $	$ 291.4^{+7.8}_{-6.3} $	$ 312.3^{+7.6}_{-6.6} $	$ 314.3^{+6.3}_{-6.5} $
208.5	270.9	291.9	313.1	314.1

Table 1.  Mean (at 68% confidence level) and most probable value of incompressibility K (MeV) inferred with and without considering the momentum-dependent interaction.

The Affine Invariant MCMC Ensemble sampler algorithm is employed to sample the posterior distribution of incompressibility K. The cumulative mean (running average) plot from the MCMC sampling is typically employed to assess the convergence status of the sampling process. In this work, each of the five MCMC chains independently runs 3 million steps, discarding the first 1.2 million as burn-in steps.

As mentioned earlier, we compare two calculations with and without considering the momentum-dependent interaction (MDI) part of single-nucleon potentials at $ E = 400 $ MeV/nucleon. In both cases, we found that 1 million burn-in steps are sufficient. As an example, Fig. 3(a) and 3(b) show the cumulative mean diagrams of our MCMC sampling without considering the MDI. It can be seen that regardless of the initial states, all chains surely converge to the same equilibrium state after approximately 1 million steps. Moreover, this conclusion is independent of the observables we used. All results presented in the following are obtained by using 3 million steps, discarding the first 1.2 million as burn-in steps.

Figure 3.  (color online) Without considering the MDI: running averages of K as functions of the MCMC steps after 200,000 initial steps in 5 independent chains. Left: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $ Right: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

We first examine the posterior PDFs without considering the MDI. The results obtained from the convergent MCMC chains are shown in Fig. 4(a) (using only $ -v_{2}\left (y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $) and 4(b) (using $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right )} \right ) $). We note that the white dashed line corresponds to the most probable value (MPV) of the posterior distribution. The black dashed line represents the prior distribution of K. In both cases, the posterior PDF is significantly narrower than the prior, demonstrating that the chosen observables provide strong constraints on K.

Figure 4.  (color online) Without considering the MDI: posterior PDFs of K. Left: using the observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Right: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

It is important to note that here we adopt the highest posterior density (HPD) method to calculate the 68% confidence interval [85] because the posterior distribution of K sometimes deviates from a Gaussian distribution. More specifically, we find from the PDF shown in Fig. 4(a) that the mean value of K is $ K = 256.1^{+8.2}_{-8.7} $ MeV at 68% confidence interval when only using the observables $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left ( u_{t0} \right ) $. In contrast, from Fig. 4(b), we find $ K = 285.5^{+6.7}_{-7.3} $ MeV when $ -v_{2}\left(y_{0}\right) $, $ -v_{2}\left(u_{t0}\right) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ are used. In this case, the mean or MPV of K is appreciably larger than that when only the elliptic proton flow data $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ are used. This demonstrates that different combinations of observables, including the proton elliptic flow and proton-like directed flow, impose distinct constraints on the incompressibility K, as discussed in the following.

As previously discussed, given an incompressibility K, the corresponding parameters α, β, and γ can be derived using the equations for $ E/A $, P, and K at saturation density $ \rho_0 $. Utilizing the convergent MCMC chains, we can readily obtain the posterior PDFs of α, β, and γ. As shown in Fig. 5(a)–(f), these posterior PDFs are significantly narrower than their prior distributions, indicating tighter constraints. Such constrained parameters are valuable for testing fundamental interactions, such as the Skyrme nuclear effective interactions.

Figure 5.  (color online) Without considering the MDI: posterior PDFs of α, β, and γ. Top: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Bottom: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

We now turn to the results obtained when incorporating the MDI. We maintain the same prior distributions for incompressibility K and apply identical constraining conditions as in the previous case. The PDFs for K are depicted in Fig. 6(a) and 6(b). The corresponding results for α, β, and γ are illustrated in Fig. 7(a)–7(c) and 7(d)–7(f), respectively. Notably, the mean values of K are now $ 188.9^{+2.9}_{-4.5} $ MeV at 68% confidence level when only the observables $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ are used, and they are $ 222.3^{+9.0}_{-9.9} $ MeV when the observables $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $ are incorporated. These results reveal a clear dependence of the inferred K on the choice of observables, similar to the earlier findings. Furthermore, the values of K obtained when including the MDI are significantly smaller than those inferred without considering the MDI. To reproduce the same observables, an equivalent gradient of mean-field potential is required, which can originate from either the density-dependent or momentum-dependent components of the mean-field potential. The MDI potential is repulsive, which provides a larger gradient of mean-field potential for particles without making the density-dependent part of the EoS stiffer [86]. However, the absence of a momentum-dependent potential in the simulation requires compensation from a stiffer density-dependent EoS to reproduce the same observables. Thus, in the case of including the momentum-dependent potential in simulating heavy-ion reactions, regardless of using only $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ or including $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $, the required K values are smaller than without considering the MDI. The results here are consistent with the earlier findings that the MDI reduces the necessary K to reproduce flow data (see, e.g., Refs. [41, 42, 87−89]). However, we acknowledge that the extracted values of K are still subject to uncertainties arising from model dependence. As discussed extensively in the Transport Model Evaluation Project (TMEP) review [34] and references therein, there are systematic differences in various transport models.

Figure 6.  (color online) Considering the MDI: posterior PDFs of K. Left: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Right: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

Figure 7.  (color online) Considering the MDI: posterior PDFs of α, β, and γ. Top: observables $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $. Bottom: observables $ -v_{2}\left ( y_{0} \right ) $, $ -v_{2}\left ( u_{t0} \right ) $, and $ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $.

The results presented earlier were obtained using proton elliptic flow and proton-like directed flow in mid-central Au+Au collisions at 400 MeV/nucleon. For a systematic study, we now infer the incompressibility solely based on the rapidity and transverse velocity dependence of proton elliptic flow in mid-central Au+Au collisions at beam energies from 600 to 1500 MeV/nucleon, with the momentum-dependent interaction taken into account, maintaining the same prior distribution of incompressibility K and constraining conditions as before. The beam energy dependences of the PDFs of K are illustrated in Fig. 8(a)–8(e). To be more quantitative, the mean values of K as functions of beam energy are shown in Fig. 8 (f), and the values of K are given in Table 1. The incompressibility K increases with beam energy, regardless of whether the observable $ -v_{2}\left(y_{0}\right) $ or $ -v_{2}\left(u_{t0}\right) $ is used. From Eq. (3), it is clear that the incompressibility K is defined at the saturation density $ \rho_0 $, and in principle, K as a property of nuclear matter at $ \rho_0 $ should not depend on the beam energy or specific observable. However, in transport model simulations where the EoS is characterized by a single parameter K, the observable-sensitive density regime probed by heavy-ion collisions increases with beam energy. Therefore, the K values extracted from data at different beam energies effectively reflect the stiffness of the EoS at different densities, rather than the incompressibility at $ \rho_0 $. The observed increasing trend of extracted K values with beam energy should thus be interpreted as an indication that the EoS becomes stiffer at higher densities and temperatures in the energy region from 400 to 1500 MeV/nucleon, which is consistent with the conclusion from the SMASH model in Ref. [90]. The general trend also aligns with Li's recent work [91], despite differences in the models and observables used. Meanwhile, it is opposite to the conclusion that gives a soft EoS with elliptic flow in the same energy region in Le Févre's work [41, 92]. The observable used in this work differs from that in Refs. [41, 92]. Moreover, it should be noted that the description of the $ u_{t0} $ cut condition adopted from Ref. [82] in Ref. [92] is inaccurate, as the descriptions of the $ u_{t0} $ cut in Refs. [82] and [93] are inconsistent. This inconsistency may affect the rigor of the conclusions drawn in Ref. [92]. Notably, the values of K derived from the observables $ -v_{2}\left(y_{0}\right) $ and $ -v_{2}\left(u_{t0}\right) $ exhibit discrepancies, which arise from the differences in the rapidity and transverse velocity dependence of proton elliptic flow. The reason here is similar to that for the inferred values of K at $ E = 400 $ MeV/nucleon: when different observables probe different density regions, the extracted EoS will differ accordingly. Specifically, elliptic flow may reflect the behavior of particles at low densities, while directed flow might be more sensitive to high-density behavior. As a result, the EoS extracted from proton elliptic flow tends to be softer, whereas including directed flow leads to a stiffer EoS, as illustrated in Figs. 4 and 6.

Figure 8.  (color online) Considering the MDI: posterior PDFs of K at $ E = 600-1500 $ MeV/nucleon and its dependence on beam energy. Observables: $ -v_{2}\left ( y_{0} \right ) $ and $ -v_{2}\left ( u_{t0} \right ) $.

	noMDI ($ E = 400 $ MeV/nucleon)		MDI ($ E = 400 $ MeV/nucleon)
Observables:	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ) $	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ), $	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ) $	$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ), $
		$ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $		$ v_{1} \left ( p_{t}^{\left ( 0 \right ) } \right ) $
K Mean	$ 256.1^{+8.2}_{-8.7} $	$ 285.5^{+6.7}_{-7.3} $	$ 188.9^{+2.9}_{-4.5} $	$ 222.3^{+9.0}_{-9.9} $
K MPV	255.5	284.6	187.9	222.3
MDI (observables:$ -v_{2}\left ( y_{0} \right ),-v_{2}\left ( u_{t0} \right ) $)
600 MeV/nucleon	800 MeV/nucleon	1000 MeV/nucleon	1200 MeV/nucleon	1500 MeV/nucleon
$ 208.8^{+6.7}_{-7.3} $	$ 271.8^{+7.8}_{-8.7} $	$ 291.4^{+7.8}_{-6.3} $	$ 312.3^{+7.6}_{-6.6} $	$ 314.3^{+6.3}_{-6.5} $
208.5	270.9	291.9	313.1	314.1

Table 1.  Mean (at 68% confidence level) and most probable value of incompressibility K (MeV) inferred with and without considering the momentum-dependent interaction.

In summary, within a Bayesian statistical framework using a Gaussian process emulator for the IQMD simulator, we have inferred the incompressibility K of the SNM from the proton elliptic flow and proton-like directed flow data in mid-central Au+Au reactions at 400−1500 MeV/nucleon from the FOPI collaboration. Compared to previous work, which was mostly based on transport model simulations in the forward modelling approach, Bayesian analyses allow us to infer the underlying transport model parameters with quantified uncertainties. In particular, we infer an incompressibility of $ K = 188.9^{+2.9}_{-4.5} $ MeV or $ K = 256.1^{+8.2}_{-8.7} $ MeV with 68% confidence from the combined proton elliptic flow data at $ E = 400 $ MeV/nucleon in analyses with and without considering the momentum dependence of the single-nucleon potentials. This is consistent with results from previous analyses of the same data using forward modelling but with significantly reduced uncertainties. However, we acknowledge that the derived K values are still subject to model dependence, an issue that warrants further investigation.

Our study also showed that the incompressibility K derived from heavy-ion reactions is related to the observable-sensitive density regime. Specifically, the derived values of K are $ 222.3^{+9.0}_{-9.9} $ MeV or $ 285.5^{+6.7}_{-7.3} $ MeV at $ E = 400 $ MeV/nucleon, depending on whether both proton elliptic flow and proton-like directed flow are considered, with or without accounting for the momentum dependence of the single-nucleon potentials. This finding may suggest that the elliptic flow reflects the behavior of particles at low densities, while the directed flow is more sensitive to high-density behavior. Thus, the EoS extracted from proton elliptic flow tends to be softer, whereas including directed flow leads to a stiffer EoS. Furthermore, the inferred values of K derived from the observables of proton elliptic flow in Au+Au collisions at beam energies ranging from 400 to 1500 MeV/nucleon reveal a general trend of progressive hardening of nuclear matter EoS at higher densities and temperatures.

In summary, within a Bayesian statistical framework using a Gaussian process emulator for the IQMD simulator, we have inferred the incompressibility K of the SNM from the proton elliptic flow and proton-like directed flow data in mid-central Au+Au reactions at 400−1500 MeV/nucleon from the FOPI collaboration. Compared to previous work, which was mostly based on transport model simulations in the forward modelling approach, Bayesian analyses allow us to infer the underlying transport model parameters with quantified uncertainties. In particular, we infer an incompressibility of $ K = 188.9^{+2.9}_{-4.5} $ MeV or $ K = 256.1^{+8.2}_{-8.7} $ MeV with 68% confidence from the combined proton elliptic flow data at $ E = 400 $ MeV/nucleon in analyses with and without considering the momentum dependence of the single-nucleon potentials. This is consistent with results from previous analyses of the same data using forward modelling but with significantly reduced uncertainties. However, we acknowledge that the derived K values are still subject to model dependence, an issue that warrants further investigation.

Our study also showed that the incompressibility K derived from heavy-ion reactions is related to the observable-sensitive density regime. Specifically, the derived values of K are $ 222.3^{+9.0}_{-9.9} $ MeV or $ 285.5^{+6.7}_{-7.3} $ MeV at $ E = 400 $ MeV/nucleon, depending on whether both proton elliptic flow and proton-like directed flow are considered, with or without accounting for the momentum dependence of the single-nucleon potentials. This finding may suggest that the elliptic flow reflects the behavior of particles at low densities, while the directed flow is more sensitive to high-density behavior. Thus, the EoS extracted from proton elliptic flow tends to be softer, whereas including directed flow leads to a stiffer EoS. Furthermore, the inferred values of K derived from the observables of proton elliptic flow in Au+Au collisions at beam energies ranging from 400 to 1500 MeV/nucleon reveal a general trend of progressive hardening of nuclear matter EoS at higher densities and temperatures.

We would like to thank Prof. Kai Zhou for useful communications.

We would like to thank Prof. Kai Zhou for useful communications.