Bayesian analysis of (3+1)D relativistic nuclear dynamics with the RHIC beam energy scan data

This work presents a Bayesian inference study for relativistic heavy-ion collisions in the Beam Energy Scan program at the Relativistic Heavy-Ion Collider. The theoretical model simulates event-by-event (3+1)D collision dynamics using hydrodynamics and hadronic transport theory. We analyze the model’s 20-dimensional posterior distributions obtained using three model emulators with different accuracy and demonstrate the essential role of training an accurate model emulator in the Bayesian analysis. Our analysis provides robust constraints on the Quark-Gluon Plasma’s transport properties and various aspects of (3+1)D relativistic nuclear dynamics. By running full model simulations with 100 parameter sets sampled from the posterior distribution, we make predictions for $p_{\rm T}$-differential observables and estimate their systematic theory uncertainty. A sensitivity analysis is performed to elucidate how individual experimental observables respond to different model parameters, providing useful physics insights into the phenomenological model for heavy-ion collisions.

💡 Research Summary

This paper presents a comprehensive Bayesian inference study of relativistic heavy‑ion collisions in the RHIC Beam Energy Scan (BES) program, focusing on (3+1)‑dimensional dynamical modeling. The authors employ the iEBE‑MUSIC framework coupled with the UrQMD afterburner to simulate Au+Au collisions at √s_NN = 7.7, 19.6, and 200 GeV. Initial conditions are generated with a three‑dimensional Glauber model that includes four “hot‑spots” per nucleon (three quarks and one soft gluon). A total of 20 model parameters are varied, covering rapidity loss (y_loss_2, y_loss_4, y_loss_6), its fluctuation width (σ_y_loss), shadowing (α_shadowing), remnant energy loss (α_rem), string‑tilt (α_string_tilt), pre‑flow (α_preFlow), baryon‑junction parameter (λ_B), string geometry (σ_string_x, σ_string_η), the switching energy density (e_sw), and the transport coefficients of the quark‑gluon plasma (QGP). The shear viscosity η/s is parametrized as a piecewise constant in baryon chemical potential μ_B (η₀, η₂, η₄), while the bulk viscosity ζ is modeled by a Gaussian peak whose height (ζ_max), location (T_ζ₀) and width (σ_ζ) depend on μ_B.

To explore the 20‑dimensional parameter space efficiently, the authors construct Gaussian Process (GP) emulators. They generate 1,000 design points using a Maximum‑Projection Latin Hypercube Design (LHD) and supplement them with 95 high‑probability posterior (HPP) points obtained from an initial Bayesian run. Two GP implementations are used: the PCSK emulator (from the surmise package) and a standard GP from scikit‑learn. Three training configurations are examined: PCSK trained on LHD + HPP, scikit‑learn trained on LHD + HPP, and PCSK trained on LHD only. The emulator predictions replace the expensive full‑model calculations during likelihood evaluation.

Uniform priors are assigned to all parameters. The likelihood assumes a multivariate Gaussian form with a covariance matrix that combines experimental statistical errors (treated as diagonal) and the emulator’s predictive covariance. Posterior sampling is performed with the pocoMC package, which combines a pre‑conditioned normalizing flow with adaptive Sequential Monte Carlo, allowing robust handling of potentially multimodal posteriors.

The posterior distributions are visualized in Figure 1. The authors quantify the information gain relative to the uniform prior using the Kullback–Leibler (KL) divergence. The PCSK emulator trained on LHD + HPP yields the largest KL value (24.6), indicating the strongest constraint, followed by PCSK on LHD only (22.4) and scikit‑learn (20.9). Inclusion of the HPP points modestly sharpens the posteriors.

Key physics results emerge from the posterior analysis. The shear viscosity at zero baryon density η₀ is constrained to ≈ 0.08–0.12, while η₂ and η₄ (values at μ_B ≈ 0.2 GeV and 0.4 GeV) rise to ≈ 0.15–0.25 and ≈ 0.20–0.30 respectively, confirming a μ_B‑dependent increase of η/s at lower beam energies. The bulk viscosity peak height ζ_max is limited to ≈ 0.05–0.12, with the peak temperature T_ζ₀ around 0.18–0.22 GeV and width σ_ζ ≈ 0.03–0.07 GeV. Initial‑state parameters such as the string‑tilt α_string_tilt (≈ 0.4–0.6) and the transverse string width σ_string_x (≈ 0.2–0.5 fm) are also tightly constrained. Rapidity‑loss parameters y_loss_2, y_loss_4, y_loss_6 and their fluctuation σ_y_loss receive non‑trivial posterior peaks, reflecting the importance of longitudinal dynamics at BES energies.

A sensitivity analysis reveals that charged‑hadron multiplicities dN/dy and elliptic flow v₂{2} are most sensitive to the shear viscosity and to geometric parameters (α_string_tilt, σ_string_x), whereas mean transverse momentum ⟨p_T⟩ is driven primarily by bulk viscosity and rapidity‑loss parameters. This mapping clarifies which experimental observables provide the strongest leverage on specific model ingredients.

To demonstrate predictive power, the authors draw 100 parameter sets from the posterior and run full (3+1)D simulations for each. They compute p_T‑differential spectra and flow coefficients v_n{2} for identified hadrons, presenting the spread as a systematic theoretical uncertainty. The predictions agree with existing data within these uncertainties, validating the calibrated model.

In summary, the work showcases how high‑fidelity (3+1)D simulations combined with accurate GP emulators enable robust Bayesian calibration of a complex heavy‑ion collision model. The analysis yields precise constraints on QGP transport properties (shear and bulk viscosities) and on initial‑state dynamics across a wide range of baryon chemical potentials, providing essential input for future studies of the QCD phase diagram, the location of a possible critical point, and the nature of the first‑order phase boundary.