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

Syed Afrid Jahan, 1, ∗ Hendrik Roch, 1, † and Chun Shen 1, ‡

1 Department of Physics and Astronomy, Wayne State University, Detroit, Michigan 48201, USA

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 eventby-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 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.

## I. INTRODUCTION

Relativistic heavy-ion collision programs at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) have been a key tool in studying the properties of the Quark Gluon Plasma (QGP), a state of matter exhibiting quarks and gluon degrees of freedom [1, 2]. These experiments allow high-energy nuclear physicists to quantitatively characterize the strongly interacting nuclear matter under extreme temperatures and densities, thereby exploring the phase structure of Quantum Chromodynamics (QCD). Especially with the systems studied in the RHIC Beam Energy Scan (BES) program, where the center-of-mass energy is varied, it is possible to study the QGP over a wide range of temperatures and densities [3-6]. With precision measurements from this experimental program, it is exciting to investigate the transition between the QGP and hadronic matter, search for a critical point and the associated first-order phase boundary in the QCD phase diagram, and study the emergent properties of the produced state of matter (see reviews [7-11]).

Modeling heavy-ion collisions event-by-event is a complex and computationally demanding task [12-17]. Because the relevant length scales change dynamically as the collision systems evolve, it requires matching different types of physical models in a multi-stage approach. Fluid dynamics serves as an effective long-wavelength description of the collective behavior of the strongly coupled QGP produced in heavy-ion collisions. Relativistic viscous hydrodynamics can well describe the time evolution of the produced bulk medium, efficiently transforming the initial-state spatial inhomogeneities to anisotropies in the final-state hadrons momentum distributions. This theoretical description allows us to access information

∗ hm0746@wayne.edu

† Hendrik.Roch@wayne.edu

‡ chunshen@wayne.edu

about the initial conditions of heavy-ion collisions and transport properties of the QGP from the experimental measurements.

It is challenging to compute the QGP transport coefficients from first principles [18]. There have been several attempts using lattice QCD techniques to compute the plasma's shear viscosity for a pure gluon system [19-22]. On the other hand, extensive phenomenological studies were able to show that the hadronic observables measured from heavy-ion collisions are sensitive to the viscosities in the QGP [14, 23-28]. There have been many studies using hydrodynamics in heavy-ion collision simulations to estimate the specific shear viscosity η/s for the QGP [23, 24, 29-31]. Nowadays, there is more effort put into the direction of extracting the QGP-specific shear and bulk viscosities, together with their uncertainties from varying other aspects of the theoretical model, based on large-scale model-to-data comparisons with Bayesian inference analysis in a high dimensional model parameter space [32-44].

The RHIC BES program generates a wealth of data for heavy-ion collisions at different collision energies. These measurements enable phenomenological extraction of the QGPtransport properties at finite baryon density via the global Bayesian inference analysis. For collisions with energies around O (10) GeV, full (3+1)D simulations are required to perform precision comparisons [45, 46].

This work expands upon Ref. [47] and performs a systematic Bayesian inference analysis of RHIC BES measurements using full (3+1)D event-by-event simulations. In addition to the constraints on the QGP-specific viscosity, we will report and analyze the posterior distribution for all model parameters in Sec. IV. We improve the constraints on the model parameters by applying a more accurate Gaussian Process (GP) model emulator [48] than the one used in [47]. Sec. V will further provide a selection of model predictions by performing numerical simulations with the parameter sets drawn from the posterior distribution. The relative variations in the model predictions from different parameter sets provide the theoret-

ical uncertainty from the Bayesian calibration. We will perform a detailed sensitivity analysis between the model parameters and experimental observables in Sec. VII.

## II. THE HYBRID THEORETICAL FRAMEWORK AND MODEL PARAMETERS

The theoretical simulations for Au+Au collisions at √ s NN = 7 7 19 6, and 200 GeV used for this Bayesian . , . analysis were performed using the iEBE-MUSIC framework [47]. To simulate the full (3+1)D dynamical evolution, a 3D-Glauber model is coupled with relativistic viscous hydrodynamics ( MUSIC ) [49-51] and hadronic transport ( UrQMD ) [52, 53]. The 3D-Glauber model considers the finite thickness of nuclei along the beam direction as they pass through each other, which becomes important as the center of mass collision energy decreases to O (10) GeV [54-56]. The transverse collision geometry is controlled by the shape of nuclei and their relative impact parameter. Each nucleon inside the colliding nuclei contains four hotspots related to three constituent quarks and one soft gluon. These hotspots lose energy and longitudinal momentum by being decelerated by strings between them during nucleon-nucleon collisions [17]. The averaged rapidity loss is parametrized as a function of the incoming rapidity y init in the collision pair rest frame [47]:

$$\langle y _ { \text{loss} } \rangle = \begin{cases} y _ { \text{loss}, 2 } \frac { y _ { \text{init} } } { 2 }, \, 0 < y _ { \text{init} } \leq 2 & \text{led} \, v \\ y _ { \text{loss}, 2 } + ( y _ { \text{loss}, 4 } - y _ { \text{loss}, 2 } ) \frac { y _ { \text{init} } - 2 } { 2 }, \, 2 < y _ { \text{init} } < 4 & \text{individual} \\ y _ { \text{loss}, 4 } + ( y _ { \text{loss}, 6 } - y _ { \text{loss}, 4 } ) \frac { y _ { \text{init} } - 4 } { 2 }, \, y _ { \text{init} } \geq 4 & \text{(1)} \\ \text{where } v _ { \text{need} } \, \text{ is the average amount of radidity loss for} \end{cases}$$

̸

where y loss ,n is the average amount of rapidity loss for y init = n . The event-by-event fluctuations of the rapidity loss are modeled by the parameter σ y loss , which controls the variance of the fluctuations [17]. The parameter α shadowing takes into account coherence effects when the incoming nucleons go through multiple collisions. With α shadowing = 0, all binary collisions will produce strings while α shadowing = 0 sets a finite probability to shadow string production from secondary collisions between nucleon pairs. The parameter α rem controls the amount of energy-momentum lost by the collision remnants [17]. These wounded partons are connected via strings and decelerated in the longitudinal direction until τ hydro = 0 5 fm . /c in the collision rest frame. Then, they are fed as source terms for the system's energymomentum tensor and net baryon current,

$$\partial _ { \mu } T ^ { \mu \nu } = J ^ { \nu }, \quad \ \ \ ( 2 ) \quad \text{of}$$

$$\overset { \cdot } { \partial } _ { \mu } J _ { B } ^ { \mu } = \rho _ { B }. \quad \quad \ \ ( 3 ) \quad \overset { \text{sp} } { \text{d} \epsilon }$$

The definitions of the source terms were explained in detail in Ref. [17]. Here, we introduce a longitudinal tilted string profile by parameterizing the transverse profile as,

$$\text{el} \\ f ( \vec { x } _ { \perp }, \eta _ { s } ) & = \frac { 1 } { 2 \pi ( \sigma _ { x } ^ { \text{string} } ) ^ { 2 } } \\ & \exp \left [ - \frac { [ x - x _ { c } ( \eta _ { s } ) ] ^ { 2 } + [ y - y _ { c } ( \eta _ { s } ) ] ^ { 2 } } { 2 ( \sigma _ { x } ^ { \text{string} } ) ^ { 2 } } \right ], \quad ( 4 )$$

where the position of the string center ( x , y c c ) depends on the space-time rapidity,

$$\underset { \underset { \underset { \underset { \underset { \underset { \underset { \underset { \underset { \underset { \text{by} } } { \text{by} } } { \text{by} } } } } } { \text{ns at} } } \quad x _ { c } ( \eta _ { s } ) = \frac { x ^ { P } + x ^ { T } } { 2 } + \alpha _ { \text{string tilt} } \frac { x ^ { P } - x ^ { T } } { \eta _ { s } ^ { P } - \eta _ { s } ^ { T } } \left ( \eta _ { s } - \frac { \eta _ { s } ^ { P } + \eta _ { s } ^ { T } } { 2 } \right ), \\ \text{evolu} \quad \quad \sim \quad \sim \quad \sim \quad \sim.$$

$$\text{evou} \\ \text{ic vis} \\ \text{trans} \\ \text{siders} \\ \text{ection}$$

Here, ( x P/T , y P/T , η P/T s ) are the space-time coordinates for the projectile ( P ) and target ( T ) hotspots after the string deceleration.

Motivated by the baryon junction model [57], we allow for the incoming baryon numbers to be located at the string ends or fluctuate along the string with a parameter λ B [17],

$$\text{The} \quad P ( y _ { P / T } ^ { B } ) = ( 1 - \lambda _ { B } ) y _ { P / T } + \lambda _ { B } \frac { e ^ { y _ { P / T } ^ { B } - ( y _ { P } + y _ { T } ) / 2 } } { 4 \sinh [ ( y _ { P } - y _ { T } ) / 4 ] }. \ ( 7 )$$

The pre-equilibrium evolution of the system is modeled via a blast-wave-like transverse flow profile for the individual strings developed in the time ≤ τ hydro with a transverse flow rapidity [47, 58]

$$\eta _ { \perp } ( { \mathbf x } _ { \perp } ) = \alpha _ { \text{preFlow} } | \tilde { \mathbf x } _ { \perp } |,$$

where ˜ x ⊥ = ( x -x , y c -y c ), with x c and y c being the string's center in the transverse plane in Eqs. (5) and (6).

The hydrodynamic equation of motion, Eqs. (2) and (3), are evolved with the lattice-QCD-based NEOS-BQS equation of state at finite densities, which assumes strangeness neutrality and n Q = 0 4 . n B for Au+Au collisions [59]. The system's viscous stress tensors are evolved with the Denicol-Niemi-Molnar-Rischke (DNMR) theory [60]. The µ B dependence of QGP-specific shear viscosity is parametrized by

$$\text{$\mathfrak { n }$ un} _ { \mathfrak { n } $ un} _ { \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $ \mathfrak { n } $(\mathfrak { n } _ { B } ) = \begin{cases} \eta _ { 0 } + ( \eta _ { 2 } - \eta _ { 0 } ) \frac { \mu _ { B } } { 0. 2 }, \ 0 < \mu _ { B } \leq 0. 2 \, \text{GeV} \\ \eta _ { 2 } + ( \eta _ { 4 } - \eta _ { 2 } ) \frac { ( \mu _ { B } - 0. 2 ) } { 0. 2 }, \ 0. 2 < \mu _ { B } < 0. 4 \, \text{GeV} \\ \eta _ { 4 }, \ \mu _ { B } \geq 0. 4 \, \text{GeV} \end{cases}$$

(9)

where ˜ η ≡ ηT/ e ( + P ) and η 0 , η 2 and η 4 are the values of the specific shear viscosity at µ B = 0 0 2 0 4 GeV re-, . , . spectively. We did not introduce an explicit temperature dependence in the parameterization to keep the number of model parameters manageable. A previous analysis showed that a temperature-independent η/s is compatible with experimental data at high energy collisions [37].

We parameterize the specific bulk viscosity by

$$\tilde { \zeta } ( T, \mu _ { B } ) & = \begin{cases} \zeta _ { \max } \exp \left [ - \frac { [ T - T _ { \zeta } ( \mu _ { B } ) ] ^ { 2 } } { 2 \sigma _ { \zeta, - } ^ { 2 } } \right ], \, T < T _ { \zeta } ( \mu _ { B } ) & \text{For} \\ \zeta _ { \max } \exp \left [ - \frac { [ T - T _ { \zeta } ( \mu _ { B } ) ] ^ { 2 } } { 2 \sigma _ { \zeta, + } ^ { 2 } } \right ], \, T \geq T _ { \zeta } ( \mu _ { B } ) & \text{space} \\ \text{withor} \ \tilde { \zeta } \ = \ \tilde { \ } T / \slash e \, \rightharpoonup \, P \, \text{ with the bill, wrote one book} \quad \end{cases} \quad \begin{array} { c } \text{For} \\ \text{comp} \end{array}$$

where ˜ ζ ≡ ζT/ e ( + P ) with the bulk pressure peak temperature T ζ ( µ B ) ≡ T ζ, 0 -0 15 . 1 GeV µ 2 B . This parameterization follows the constant energy density curve e = e T ( ζ, 0 , µ B = 0) of the NEOS-BQS equation of state closely and ensures that the parameterization peaks close to the phase crossover at µ B &gt; 0 [11, 61, 62].

We choose to normalize the QGP shear and bulk viscosity by T/ e ( + P ) instead of the local entropy density s to better control the variations of the viscous relaxation times in numerical simulations. Take the shear relaxation time, for example,

$$\tau _ { \pi } \equiv C _ { \eta } \frac { \eta } { ( e + P ) } = C _ { \eta } \tilde { \eta } \frac { 1 } { T } & & \text{bias} \\ & & \text{(11)} & & \text{all} y, \\ & & \text{where} ^ { + } }$$

$$\Phi = C _ { \eta } \frac { \eta } { s } \frac { 1 } { T } \frac { 1 } { 1 + \frac { \mu _ { B } n _ { B } } { T s } }. \quad ( 1 2 ) \quad \text{d}$$

One can see that the shear relaxation time τ π in Eq. (11) remains well behaved when we vary ˜ within the prior η range. If we parameterize η/s as in Eq. (9), then τ π in Eq. (12) could potentially reduce to very small values and violate causality conditions [63, 64] when the factor µ B n B / Ts ( ) becomes large during the dynamical evolution.

All fluid cells dropping below an energy density e sw are particlized using a Cooper-Frye freeze-out procedure including δf corrections for multiple conserved charge currents ( B,Q,S ) using the Grad moment method [12, 65]. The hadronic rescatterings and decays in the afterburner phase are modeled using UrQMD [52, 53].

All 20 model parameters, together with their prior ranges, are listed in Table I.

TABLE I. The 20 model parameters and their prior ranges.

| Parameter                                                                                                 | Prior                                                                                                | Parameter                                                                                       | Prior                                                                                                                                                    |
|-----------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------|
| B G [GeV - 2 ] α shadowing y loss , 2 y loss , 4 y loss , 6 σ y loss α rem λ B σ string x [fm] σ string η | [1 , 25] [0 , 1] [0 , 2] [1 , 3] [1 , 4] [0 . 1 , 0 . 8] [0 , 1] [0 , 1] [0 . 1 , 0 . 8] [0 . 1 , 1] | α string tilt α preFlow η 0 η 2 η 4 ζ max T ζ, 0 [GeV] σ ζ, + [GeV] σ ζ, - [GeV] sw [GeV / fm 3 | [0 , 1] [0 , 2] [0 . 001 , 0 . 3] [0 . 001 , 0 . 3] [0 . 001 , 0 . 3] [0 , 0 . 2] [0 . 15 , 0 . 25] [0 . 01 , 0 . 15] [0 . 005 , 0 . 1] [0 . 15 , 0 . 5] |

## III. MODEL EMULATION AND BAYESIAN INFERENCE SETUP

For the exploration of the 20-dimensional parameter space shown in Table I with Bayesian inference, we need efficient Gaussian Process (GP) emulators to replace the computationally expensive high-fidelity model calculations with the iEBE-MUSIC framework. In this work, we will make use of the PCSK emulator implemented in the surmise package [66] and a standard GP emulator from the Scikit-learn Python package [47]. For the details of the implementation of these emulators, we refer to our previous comparison of these GP emulators in Ref. [48].

For an efficient exploration of the parameters in Table I, we generate model simulations at 1,000 different parameter sets in 20 dimensions using the Maximum Projection Latin Hypercube Design (LHD). At every design point, the model simulations consist of 1,000 minimum bias Au+Au collisions at 200 GeV and 2,000 minimum bias events each at 19 6 GeV and 7 7 GeV. . . Additionally, we make use of 95 additional full model simulations, where the parameters were sampled from the posterior distribution. This additional dataset is denoted as High Probability Posterior (HPP) points. In our standard setup, we train the GP emulator with the LHD + HPP datasets combined, which improves the emulator predictions in the parameter region with higher probability.

The fast prediction of model outputs at different points in the parameter space using GP emulators allows us to make use of Bayes' rule:

$$\mathcal { P } ( \theta | _ { \mathbf Y _ { \exp } } ) = \frac { \mathcal { P } ( \mathbf Y _ { \exp } | \theta ) \mathcal { P } ( \theta ) } { \mathcal { P } ( \mathbf Y _ { \exp } ) }. \quad \quad ( 1 3 )$$

P ( θ | y exp ) represents the posterior probability density function, P ( θ ) is the prior probability density function, P ( y exp ) stands for the total evidence, and P ( y exp | θ ) is the likelihood function. In this work, we choose uniform prior distributions for all the model parameters listed in Table I. The likelihood function is defined as,

$$\text{prior} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad \text{***} \quad
\text{iges.} \quad \times \exp \left [ - \frac { 1 } { 2 } ( y _ { \sin } ( \theta ) - y _ { \exp } ) ^ { T } \Sigma ^ { - 1 } ( y _ { \sin } ( \theta ) - y _ { \exp } ) \right ], \\ \quad. \quad - \quad - \quad - \quad. \quad. \quad. \quad. \quad.$$

where Σ ≡ Σ model + Σ exp is the covariance matrix. The GP emulator provides the model covariance matrix Σ model . Lacking the knowledge of the full covariance matrix for the experimental measurements, we assume a diagonal form for Σ exp using the square of the statistical errors for its elements.

To tackle the computation of the posterior distribution numerically, we employ Markov Chain Monte Carlo (MCMC) techniques [67, 68]. In Ref. [48], we have shown that the affine invariant MCMC implemented in the emcee Python package [69] and the Parallel Tempering Langevin Monte Carlo (PTLMC) implemented in

surmise lead to consistent results. The parallel tempering methods generally work better when there are multiple modes in the posterior distributions. This was the case for the closure test performed there, where the uncertainty of the test point was sufficiently small. For real experimental data with a larger uncertainty, we will therefore use a method that can deal with multi-modal posterior distributions. Since the previously tested PTLMC method does not compute the evidence P ( y exp ) (this factor is essential for computing the Bayes factor, which we use to compare different training processes in our analysis), we will employ another efficient MCMC sampler routine implemented in the pocoMC package [70, 71]. This sampler implements Preconditioned Monte Carlo (PMC), where first normalizing flow is applied to decorrelate the parameters and the distributions, and then adaptive Sequential Monte Carlo (SMC) is applied.

The Bayesian inference analysis in this paper is performed using the STAR and PHOBOS measurements for Au+Au collisions at the RHIC BES program. The details are listed in Table II.

TABLE II. The experimental measurements for Au+Au collisions used in our Bayesian inference analysis. We include identified d N/ y d and ⟨ p T ⟩ from central to 60% in collision centrality and the charged hadron v n { } 2 from 0 to 50% centrality.

|   √ s NN [GeV] | STAR                                                                                                        | PHOBOS                              |
|----------------|-------------------------------------------------------------------------------------------------------------|-------------------------------------|
|          200   | d N/ d y ( π + , K + , ¯ p, p ) [72] ⟨ p T ⟩ ( π + , K + , p, ¯ p [72] v ch 2 { 2 } [75], v ch 3 { 2 } [76] | d N ch / d η [73] v ch 2 ( η ) [74] |
|           19.6 | d N/ d y ( π + , K + , p ) [77] ⟨ p T ⟩ ( π + , K + , p, ¯) p [77] v ch 2 { 2 } [75], v ch 3 { 2 } [76]     | d N ch / d η [73]                   |
|            7.7 | d N/ d y ( π + , K + , p ) [77] ⟨ p T ⟩ ( π + , K + , p, ¯) p [77] v ch 2 { 2 } [75], v ch 3 { 2 } [76]     |                                     |

In Ref. [48], we systematically studied the performance of different model emulation setups on these experimental measurements and quantified their impacts on the closure test.

## IV. POSTERIOR ANALYSIS

In this section, we perform Bayesian analyses on the RHIC measurements using three different GP models: PCSK and Scikit GP trained with LHD + HPP simulations, and PCSK trained only with the LHD dataset. We will first evaluate and compare the general aspects of the posterior distributions from the three GP models and then discuss the physics insights from the constraints on model parameters.

In Fig. 1, we show the posterior marginal distributions for the 20 model parameters using the PCSK emulator (blue), the Scikit GP emulator (green), and the PCSK emulator trained only with the 1,000 LHD points (magenta). The former two GP emulators were trained with the 1,095 LHD + HPP training points. We showed that the PCSK emulator performed much better than the Scikit GP in the closure test in Ref. [48]. However, when applying the Bayesian inference analysis to the real experimental measurements, the two emulators result in similar posterior distributions for the model parameters. For most of the parameters, we find that the PCSK emulator still gives slightly more peaked distributions than those from the Scikit GP.

Comparing the blue and magenta posterior distributions, we find that including the additional 95 HPP points in training the GP emulator results in slightly stronger constraints on the model parameters.

The amount of information gained in our Bayesian inference analysis can be quantified using the Kullback-Leibler (KL) divergence, which is defined as,

$$D _ { \text{KL} } ( p | | q ) = \int \text{d} \theta \, p ( \theta ) \ln \left ( \frac { p ( \theta ) } { q ( \theta ) } \right ). \quad ( 1 5 )$$

Here, we compute the KL divergence with p ( θ ) as the posterior distribution P ( θ | y exp ) in Eq. (13) and q ( θ ) as the uniform prior distribution P ( θ ). Using Eq. (13), we can compute the KL divergence from the sampled posterior chains as

$$D _ { \text{KL} } ( p | | q ) = \frac { 1 } { N _ { \text{samples} } } \sum _ { i } \ln \left ( \frac { \mathcal { P } ( \mathbf y _ { \text{exp} } | \boldsymbol \theta _ { i } ) } { \mathcal { P } ( \mathbf y _ { \text{exp} } ) } \right ), \quad ( 1 6 )$$

where the index i runs through the N samples parameter samples in the posterior chains.

TABLE III. The Kullback-Leibler divergence of posterior distributions to the uniform prior in our Bayesian inference analysis using different GP emulators.

| GP Model         |   D KL |
|------------------|--------|
| PCSK (LHD + HPP) |   24.6 |
| Scikit GP        |   20.9 |
| PCSK (LHD only)  |   22.4 |

Table III reports the KL divergence for the three posterior distributions presented in Fig. 1. The ordering of the D KL is consistent with our observation in Fig. 1 discussed above. The more accurate model emulator, PCSK, gives stronger constraints on model parameters than those from the Scikit GP. Including 95 high probability posterior (HPP) points in training the GP emulator further tightens the model constraints.

Another metric to compare the posterior distributions is the Bayes factor. The Bayes factor is defined as the ratio of posterior probabilities from two models, A and

FIG. 1. Comparison of the marginal posterior distributions of the 20 model parameters using the MCMC sampling with the PCSK and Scikit GP emulators. Values above the plots show the parameters' median and 90% confidence intervals. The x -axes represent the ranges of the uniform prior distributions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-6720d3e8-7e3f-4b26-b0e9-d77311ee25ab.jpg)

B , given the experimental data y exp , set across the entire parameter space:

$$\mathcal { B } _ { A / B } \equiv \frac { \mathcal { P } ( A | \mathbf y _ { \exp } ) } { \mathcal { P } ( B | \mathbf y _ { \exp } ) } = \frac { \mathcal { P } ( \mathbf y _ { \exp } | A ) \mathcal { P } ( A ) } { \mathcal { P } ( \mathbf y _ { \exp } | B ) \mathcal { P } ( B ) }. \quad ( 1 7 )$$

Here, we compare the three GP emulator models shown in Fig. 1. Because the three GP models are emulating the same theoretical model, we assume the prior model probability is the same for all the models, namely P ( A ) = P ( B ) in Eq. (17). Thus, the Bayes factor reduces to the ratio of marginalized likelihoods, also known as Bayesian evidence:

$$\mathcal { B } _ { A / B } = \frac { \mathcal { P } ( \mathbf y _ { \exp } | A ) } { \mathcal { P } ( \mathbf y _ { \exp } | B ) }. \quad \quad ( 1 8 ) \quad \text{ac} \quad \text{co} \\ \text{dail}$$

The Bayesian evidence for a model A on the given experimental data y exp , P ( y exp | A ), can be computed by integrating the likelihood function over the model parameter

$$\mathcal { P } ( \mathbf Y _ { \exp } | A ) = \int \mathrm d \theta _ { A } \, \mathcal { P } ( \mathbf Y _ { \exp } | \theta _ { A }, A ) \mathcal { P } ( \theta _ { A } ). \quad ( 1 9 )$$

We compute the Bayesian evidence for the three models using the pocoMC package. The Bayes factors are reported in Table IV.

Following the Jeffreys' Scale of the Bayes factor in Ref. [78], we find strong evidence in favor of the analysis with the PCSK emulator than the one with the Scikit GP. This result emphasizes the importance of getting an accurate model emulator in Bayesian studies. We also compare the two model setups with different training data sets using the PCSK emulator. In this case, we find a Bayes factor with a magnitude between 1 and 2.5, which suggests that the PCSK trained with only the LHD points is favored but only with weak evidence. Naively,

one would expect the GP emulator with more training data to be more favored. The reversed hierarchy here is because the PCSK emulator trained with the LHD + HPP dataset gave a smaller model covariance matrix Σ model , which leads to an overall lower likelihood in Eq. (14) compared to the results from PCSK trained with LHD only.

Considering the KL divergence and Bayes factor results, we will perform further analysis on the posterior distribution from the PCSK emulator (LHD + HPP) from now on.

TABLE IV. Bayes factors for different combinations in the model setup.

| Model A          | Model B         | ln ( B A/B )      |
|------------------|-----------------|-------------------|
| PCSK (LHD + HPP) | Scikit GP       | 6 . 94 ± 0 . 03   |
| PCSK (LHD + HPP) | PCSK (LHD only) | - 1 . 41 ± 0 . 03 |

Now, we discuss the physics implications of the constraints on model parameters from the posterior distribution.

- · The parameter B G is the width of the nucleon to sample its hot spots [17]. The posterior range B G = 15 99 . +7 14 . -9 66 . GeV -2 translates to the nucleon mass width between 0.5 and 0.95 fm at 90% confidence level.
- · The parameter α shadowing controls the probability of producing strings from secondary binary collisions in the 3D-Glauber initial state model [17, 55]. Our Bayesian analysis prefers α shadowing ∼ 0 15, in-. dicating that only 15% of the binary collisions are shadowed for soft string production. The numbers of produced strings scale closer to the number of binary collisions than the number of participants.
- · The rapidity loss fluctuation parameter σ y loss ∼ 0 3 . in our analysis, which is smaller than the value σ y loss = 0 6 used in Ref. [17]. . The latter was calibrated using the multiplicity distribution in p + p collisions, which has a stronger sensitivity to this model parameter than the 10% centrality binned particle yields in Au+Au collisions [17]. The particle multiplicity fluctuations in Au+Au collisions include the number of participant fluctuations, further diluting the sensitivity to the rapidity loss fluctuation parameter.
- · The parameter for the nucleons' remnant rapidity loss α rem = 0 54 . ± 0 15 is . constrained by the shape of the pseudo-rapidity distribution d N ch / d . η Its value is consistent with the previous simulations [17].
- · The parameter λ B controls the amount of initialstate baryon charge transport to the string junction [17]. Its constrained value λ B ∼ 0 15 is consis-. tent with previous analyses [17, 79].
- · Our Bayesian analysis prefers α preflow ∼ 0, sug-

FIG. 2. Posterior distributions for the average rapidity loss ⟨ y loss ⟩ (top) as a function of the initial rapidity, the QGP specific shear viscosity (middle) and bulk viscosity (bottom) as functions of µ B and temperature respectively. Bands indicate the 68% (dashed) and 90% (dotted) confidence levels. The 90% prior is shown as a gray band. For all curves, the PCSK (LHD + HPP) emulator is used in the MCMC.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-0f3ca29e-13af-494e-bbe0-b588501ad358.jpg)

gesting small pre-equilibrium flow in the transverse plane for the string source terms before depositing into hydrodynamics.

- · The posterior distribution for the hotspot transverse size σ string x peaks at the lower bound of the prior, suggesting a small hotspot size is preferred.
- · The parameter σ string η controls the string profile's fall-off along the longitudinal direction. Its marginal distribution has a double-peaked structure, indicating this parameter may have a √ s NN dependence.

To quantify the constraints on the functional dependence for the initial-state rapidity loss and the QGPspecific viscosities, we display the 90% confidence level of the prior distributions (gray) with the 68% and 90% confidence levels of the posterior distributions (green) for these parameters in Fig. 2.

The average rapidity loss in the initial state is, compared to the prior distribution (gray), tightly constrained by the particle yield measurements at multiple collision energies. Despite using different types of GP emulators and MCMC algorithms, our results are consistent with those in Ref. [47]. We also obtain strong constraints on the specific shear and bulk viscosities of QGP as functions of µ B and T , respectively. The values of QGPspecific shear viscosity show a significant increase with µ B within [0, 0.2] GeV. This constraint primarily comes from the pseudo-rapidity dependence of elliptic flow v 2 ( η ) measurements [80]. The values of ˜ at larger η µ B remain around 0.3 closer to the upper bound of the prior for the parameter η 4 . For the QGP's specific bulk viscosity, we find a peaked distribution around a temperature of ≈ 200 MeV with a maximum value of 0.12. The values of ˜ fall off to zero quickly when departing from the ζ maximum value.

## V. MODEL SIMULATIONS WITH POSTERIOR DISTRIBUTIONS

After using the GP emulators to obtain the posterior distributions of model parameters displayed in Fig. 1, we will verify these constraints with full model simulations and further make model predictions. We randomly choose 100 samples { θ post } from the posterior distributions generated with the PCSK emulator (LHD + HPP) and then perform event-by-event simulations with the iEBE-MUSIC framework. At each posterior model parameter set θ post i , we simulate 1,000 minimum bias Au+Au collisions at √ s NN = 200 GeV, 2,000 events at √ s NN = 19 6 GeV, and 3,000 collisions at . √ s NN = 7 7 GeV. With . these simulations, we compute the averages of the observables from the 100 posterior parameter sets and denote their standard deviations as the remaining systematic uncertainties after the Bayesian analysis. Here, we assume the statistical errors of observables at individual posterior parameter sets are much smaller compared to the standard deviations of the 100 posterior parameter sets. In Fig. 3, we show the comparison to the experimental observables (see Tab. II) used in the Bayesian inference.

In general, the full model calculations with parameters sampled from the posterior distribution can provide a good description of the 544 data points in the Bayesian inference. The agreement is generally better for central collisions and near the mid-rapidity region. This behavior is related to the fact that the GP emulator has smaller emulation errors for those observables than the observables in peripheral centralities and forward rapidity regions [48].

Looking into details of the mid-rapidity observables at √ s NN = 200 GeV (right columns of Fig. 3), our model slightly overpredicted the particle yields of K + mesons. The centrality dependencies of proton and anti-proton yields are somewhat flatter than the STAR measurements. All these observables are sensitive to the switching energy density from hydrodynamics to hadronic transport, the current model discrepancy could be reduced by allowing the parameter e sw to be centrality dependent. The identified particles' mean p T and mass hierarchy reflect the amount of radial flow developed in the hydrodynamic phase. Our model provides a good description of the STAR measurements for pions and kaons, while slightly underestimating the proton and anti-proton mean p T but still within the experimental uncertainties. The averaged transverse momenta of protons is more sensitive to the scattering cross sections and baryon-anti-baryon annihilation processes in the hadronic transport model [28].

For the mid-rapidity observables at √ s NN = 19 6 GeV, . the description of kaon's particle yields is slightly improved compared to that at 200 GeV. At this collision energy, the averaged net baryon chemical potential can reach up to 150-200 MeV at mid-rapidity [77]. Our model gives a ∼ 50 MeV larger mean p T for anti-protons than those for protons at finite baryon density. The current STAR measurement can not distinguish this difference. In our model, the mean p T difference between protons and anti-protons primarily comes from the UrQMD hadronic scattering phase as anti-protons emit relatively earlier than protons from the hydrodynamic hypersurface and experience more hadronic interactions [51]. This phenomenon is more pronounced at √ s NN = 7 7 GeV, . at which the anti-proton mean p T is about 100 MeV larger than that of protons in our calculations. At √ s NN = 7 7 GeV, our model underestimates the mean . p T of pions and kaons, suggesting the systems do not develop enough radial flow in the hydrodynamic phase. Introducing an explicit µ B dependence on the QGP-specific bulk viscosity would help to reduce the tension here. Because the hydrodynamic phase at √ s NN = 7 7 GeV is limited, . we observe that the triangular flow coefficients v ch 3 { 2 } start to deviate from the experimental data in peripheral collisions, similar to the finding in Ref. [83].

FIG. 3. Observables computed with the full theoretical model at 100 posterior samples { θ post } . For the three different collision energies √ s NN = 7 7 GeV, 19.6 GeV, and 200 GeV, we simulate 3,000, 2,000, and 1,000 minimum bias Au+Au collisions at each . posterior point θ post i , respectively. The full line indicates the mean value of the observable averaged from the 100 parameter sets, and the band indicates the one standard deviation interval. The experimental data are summarized in Table II.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-94db335c-3ff2-4efd-b0c8-a62f722db30e.jpg)

FIG. 4. Model predictions for p T -differential observables from the 100 posterior parameter samples { θ post } . The three columns from left to right represent observables at √ s NN = 7 7 GeV, 19 6 GeV, and 200 GeV, respectively. The first two rows show . . p T -differential π + and p spectra compared with the STAR measurements [72, 77], while the bottom row shows the p T -differential elliptic flow v 2 { SP } compared with the STAR data [81, 82]. The band indicates the one standard deviation interval computed from the 100 parameter sets.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-25c1ae34-9ec5-486c-ae44-9ab6f7d6640b.jpg)

PHOBOS charged hadron d N ch / d η measurements across collision centrality at 200 and 19.6 GeV. Our model gives a flatter shape for v ch 2 ( η ) compared to the PHOBOS measurements. This is because our parameterization for the QGP shear viscosity is temperature-independent. Introducing a η/s T ( ) will improve the model-to-data comparison [84].

The full model simulations allow us to go beyond the calibrated experimental measurements and make model predictions with controlled systematic uncertainties. In this work, we select p T -differential observables to demonstrate the predictive power of our model. We will report a systematic extrapolation of our model to different collision systems such as ( p , d , 3 He)+Au, O + O , Ru+Ru,

Zr+Zr collisions, and Au+Au collisions at other collision energies in the RHIC BES program in future work.

Since our model is calibrated using only the p T -integrated observables in Table II, the model-to-data comparison for p T -differential observables can effectively test the consistency of our Bayesian analysis.

Figure 4 shows our model predictions for the p T -differential spectra for pions and protons and charged hadron v 2 ( p T ) in Au+Au collisions at 7.7, 19.6, and 200 GeV compared with the STAR measurements. We find a good description of the identified particle p T spectra at 200 and 19.6 GeV, while the p T spectra are slightly steeper than the STAR data at 7.7 GeV. These results are consistent with our description of the identified parti-

cle mean p T in Fig. 3. Our results show that the first two moments of the particle p T spectra, namely the particle yield and mean p T , are sufficient to produce reasonable p T spectra in the Bayesian calibration [32].

For the charged hadron p T -differential elliptic flow, our model prediction can reproduce the experimental data fairly well up to 2 GeV. We notice that the systematic uncertainty in v 2 ( p T ) grows with p T in the calculations. The effects from the out-of-equilibrium δf corrections in the Cooper-Frye particlization become important at high p T . The v 2 ( p T ) at 200 GeV shows a strong suppression in p T &gt; 2 GeV from the δf based on the Grad's moment method.

## VI. HIGH LIKELIHOOD PARAMETERS

Although our approach in Sec. IV is a preferred way to make model predictions with controlled systematic errors, it is computationally demanding to perform eventby-event simulations at O (100) parameter sets sampled from the posterior distribution, especially for observables which require high statistics.

The Bayesian analysis usually reports the parameter set θ MAP at the MAP (maximum a posteriori probability) estimate using the GP emulator. In our case, we find the MAP parameter set in the entire prior space gives a reverse ordering for some rapidity loss parameters, namely y loss 4 , &gt; y loss 6 , . In the Bayesian analysis above, we did not impose any constraints on the prior distribution of the rapidity loss parameters. As shown in Fig. 2, the experimental measurements at 200 GeV give strong constraints on the initial-state rapidity loss near its beam rapidity, y beam ≡ arccosh ( √ s NN / (2 m N ) ) = 5 36. . Therefore, the values of y loss 4 , and y loss 6 , are anti-correlated (see Fig. 6 below), and both orderings are allowed in the posterior distribution. Although the MAP parameter set with y loss 4 , &gt; y loss 6 , is sufficient to perform simulations for the RHIC BES program, we expect it to fail when extrapolating to the higher LHC energies. Therefore, we

TABLE V. The highest likelihood parameter set obtained from the MCMC with the PCSK emulator (LHD + HPP), imposing the monotonic constraint y loss 2 , ≤ y loss 4 , ≤ y loss 6 , .

| Parameter       |   Value | Parameter          |   Value |
|-----------------|---------|--------------------|---------|
| B G [GeV - 2 ]  |  17.095 | α string tilt      |   0.884 |
| α shadowing     |   0.145 | α preFlow          |   0.004 |
| y loss , 2      |   1.467 | η 0                |   0.045 |
| y loss , 4      |   1.759 | η 2                |   0.28  |
| y loss , 6      |   2.26  | η 4                |   0.287 |
| σ y loss        |   0.356 | ζ max              |   0.148 |
| α rem           |   0.611 | T ζ, 0 [GeV]       |   0.214 |
| λ B             |   0.129 | σ ζ, + [GeV]       |   0.018 |
| σ string x [fm] |   0.113 | σ ζ, - [GeV]       |   0.04  |
| σ string η      |   0.156 | e sw [GeV / fm 3 ] |   0.35  |

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-d75948bb-3e4f-4160-9337-fe77d3c08ad0.jpg)

FIG. 5. Posterior distributions for the average rapidity loss ⟨ y loss ⟩ as a function of y init . The band indicates the 90% confidence level imposing the monotonic constraint y loss 2 , ≤ y loss 4 , ≤ y loss 6 , . The black dotted line shows the 90% confidence interval without the constraint from Fig. 2, and the 90% prior is shown as a gray band.

FIG. 6. Single parameter (diagonal) and two-parameter joint (off-diagonal) distributions of the posterior using the PCSK emulator (LHD + HPP) in the MCMC for the parameters related to the initial-state rapidity loss. The upper right part (magenta) shows the parameter correlation without the monotonic constraint. The lower left part (green) represents the results imposing the constraint y loss 2 , ≤ y loss 4 , ≤ y loss 6 , .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-989e00b9-9f78-4c14-ba8e-ba868d46443a.jpg)

choose to report the highest likelihood parameter set with the monotonic constraint y loss 2 , ≤ y loss 4 , ≤ y loss 6 , in Table V. We compute the Bayesian evidence without and with the monotonic constraint and obtain a Bayes factor

of ln ( B w o constr . . ./ w constr . . ) = 0 85 . ± 0 02, which falls in the . 'inconclusive' category based on the Jeffreys' Scale [78] to distinguish the two setups.

Figure 5 compares the posterior distribution for the averaged rapidity loss ⟨ y loss ⟩ as a function of y init with and without the monotonic constraints, y loss 2 , ≤ y loss 4 , ≤ y loss 6 , . The 90% posterior band shrinks almost by half for y init ∈ [3 , 6]. Because the values for y loss 2 , and y loss 4 , are also anti-correlated in the posterior distribution, the ⟨ y loss ⟩ at small y init also receive some minor effects from the monotonic constraint.

Figure 6 shows the marginalized posterior distributions for individual rapidity loss parameters on the diagonal and the correlations between the different parameters in the off-diagonal corners.

The diagonal panels show that imposing the monotonic constraint on the rapidity loss parameters leads to more peaked marginalized posterior distributions because the prior space is smaller with the constraint. While the peak position of y loss 2 , is barely affected by the constraint, the peak position of y loss 4 , is shifted to a slightly smaller value and the position of y loss 6 , to a slightly larger value. The marginal distribution for the rapidity loss fluctuations σ y loss is almost unaffected. We have further checked that the marginalized posterior distributions of the other 16 parameters are not affected by imposing the monotonic constraint.

The upper right off-diagonal panels show the correlation between the rapidity loss parameters without the monotonic constraint. We find strong anti-correlations between y loss 2 , and y loss 4 , and between y loss 4 , and y loss 6 , . Such anti-correlation is a consequence of the piece-wise parameterization in Eq. (1) when the experimental measurements at 200 and 19.6 GeV set constraints on the average rapidity loss near their beam rapidity 5.36 and 3.04, respectively. Looking at the constrained case in the lower left part of the off-diagonal panels, we also observe a sharp cut in the y loss 6 , vs. y loss 4 , distribution, which is a direct consequence of the y loss 4 , ≤ y loss 6 , constraint in the prior.

## VII. SENSITIVITY ANALYSIS

The fast model emulator lets us quantitatively study how individual model parameters affect different experimental observables.

Starting with a global sensitivity analysis, we compute the first-order Sobol' indices [85] for individual experimental observables over the entire prior region. Similar analyses were performed for heavy-ion collisions with (2+1)D simulations assuming longitudinal boost invariance [43, 86]. The first-order Sobol' indices quantify the importance of each parameter for a given observable O ( θ ) by performing a decomposition of the variance of O ( ) over the parameter space · θ ≡ ( θ , . . . θ 1 m ). The firstorder Sobol' index for the model parameter θ i is defined

$$S _ { i } ^ { \mathbb { O } } \equiv \frac { \text{Var} _ { \theta _ { i } } \left ( \mathbb { E } _ { \theta _ { - i } } ( \mathbb { O } ( \theta ) | \theta _ { i } ) \right ) } { \text{Var} _ { \theta } ( \mathbb { O } ( \theta ) ) }, \ \ i = 1, \dots, m. \quad ( 2 0 )$$

Here, the index i runs over the dimension of the model parameter space ( m = 20). We sample each model parameter θ i independently from a uniform distribution over its corresponding prior range in Table I. The expression E θ -i ( O ( θ ) | θ i ) performs an average of the observable O over a parameter subspace θ -i ≡ θ \ θ i by holding the i -th parameter θ i at a fixed value. Then Sobol' index S O i measures the fraction of the variance of the observable O from the parameter θ i relative to its total variance.

Figure 7 shows the first-order Sobol' indices for identified particle yields d N/ y d and ⟨ p T ⟩ in the most central events and charged hadron v ch 2 { 2 } at 0-5% and 20-30% centrality bins. The particle yields of charged pions have strong sensitivities to the rapidity loss parameters y loss ,n . The measurements at 200 GeV are sensitive to the rapidity loss parameters at incoming rapidity y init = 4 and y init = 6, while the strongest sensitivity shifts to y loss 2 , for the lower collision energy measurements. Comparing the Sobol' indices between pion and proton yields, we find that the proton yields show additional sensitivity to the baryon stopping parameter λ B and the string's longitudinal profile parameter σ string η .

The averaged transverse momenta of pions and protons are mostly sensitive to the initial-state hotspot size in the transverse plane σ string x and the magnitude of preequilibrium flow α preFlow . Both model parameters have strong effects on the development of hydrodynamic radial flow. The hotspot's transverse size controls the magnitude of the early-stage pressure gradients. A stronger pre-equilibrium flow leads to faster fireball expansion at the early time. We also find the pion's mean p T at 200 GeV shows some sensitivity to the specific bulk viscosity of the QGP. However, the magnitudes of their Sobol' indices are dwarfed compared to those model parameters in the initial-state and pre-equilibrium stages. A similar situation is observed for mid-rapidity anisotropic flow. The charged hadron elliptic flow coefficients at central and semi-peripheral collisions are mostly sensitive to the hotspot's transverse size and pre-equilibrium flow over the entire prior range. The sensitivity to the specific shear viscosity of the QGP is very small. We checked that the charged hadron triangular flow shows very similar global sensitivity as the elliptic flow observables.

Figure 8 further represents the first-order Sobol' indices for rapidity-dependent observables. We find that the charged hadron yields at forward rapidity show a stronger sensitivity to the specific bulk viscosity of the QGP and also the amount of rapidity loss in the beam remnant α rem . We find a similar global sensitivity for charged hadron elliptic flow at central and forward rapidity regions.

Although it is informative to quantify the global sensitivity of experimental observables to the model parameter over the entire prior parameter space shown in Figs. 7

FIG. 7. The first-order Sobol' indices for selected midrapidity observables at different collision energies √ s NN and centrality bins. The black bars indicate the 95% confidence interval.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-0dfff1a8-39f0-4565-b60a-3073c2211879.jpg)

and 8, the responses could be different in different local regions of the parameter space. Therefore, we further compute the response coefficients between experimental observables and model parameters using the posterior

FIG. 8. The first-order Sobol' indices for selected observables at central and forward rapidity regions. The black bars indicate the 95% confidence interval.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-aeebe1d3-25e7-4b4c-b579-7895e01bffb6.jpg)

distribution. We follow the approach in Ref. [87] and define the ensemble-averaged response coefficient matrix,

$$\mathcal { R } _ { a i } \equiv \left \langle \frac { \partial y _ { a } } { \partial \theta _ { i } } \right \rangle _ { \text{post} } = \sum _ { j } \langle \delta y _ { a } \delta \theta _ { j } \rangle _ { \text{post} } \langle \delta \theta _ { j } \delta \theta _ { i } \rangle _ { \text{post} } ^ { - 1 }. \ \ ( 2 1 ) \quad \text{cosity} \\ \quad \text{sensit} \\ \quad \text{is rela}$$

Here, the index a runs over individual experimental observables in Table II, and the index i runs over the model parameters in Table I. The ensemble average ⟨·⟩ post is performed on the posterior distribution to extract the model's response to the experimental observables near the real measurements. To make the response coefficients unitless, we define the normalized experimental observables δy a = ( y a ( θ j ) -¯ ) y a /σ y a and model parameters δθ i = ( θ i -¯ ) θ i /σ θ i , where ¯ ( ¯ ) and y a θ i σ y a ( σ θ i ) are the mean and standard deviation of the observable y a (parameter θ i ) over the posterior distribution, respectively. The right-hand side of Eq. (21) should be interpreted as a matrix multiplication.

Figure 9 shows the response matrix between the experimental observables and model parameters over the posterior distribution. We observe that identified particle yields show strong sensitivity to the rapidity loss parameters, consistent with the results from the global analysis above. The signs of the response coefficients provide additional information compared to Sobol' indices. Interestingly, the mid-rapidity charged hadron yields at 19.6 GeV show positive correlations with the rapidity loss parameters while the forward d N ch / d η negatively correlates with these parameters. This pattern results from the energy-momentum conservation imposed in the 3D-Glauber model. More particle production near mid-rapidity would require more energy loss during the initial state, which leaves less available energy at forward rapidity regions. Particle yields show significant anti-correlation with the parameter α shadowing because a smaller value of α shadowing produces more string sources from the binary collisions and increases mid-rapidity particle production.

Different from the global sensitivity analysis above, we find the observables' response coefficients to the hotspot transverse size σ string x and pre-equilibrium flow α preFlow are small in the posterior distribution. This difference is because the marginalized posterior distributions for these two parameters are sharply peaked around the edge of their prior range, as shown in Fig. 1. The correlations in the local region are small and different from the global averages.

We find that the observables' responses to the QGP transport properties in Fig. 9 are more in line with our physics intuition.

The identified particles' mean p T shows sizable negative correlations with the parameters related to the specific bulk viscosity of the QGP. This result aligns with our intuition that the bulk viscosity acts as a resistance to reduce the development of hydrodynamic radial flow.

The parameters related to the specific shear viscosity of the QGP have a significant negative correlation with the charged hadrons' anisotropic flow coefficients at 200 GeV. The pseudo-rapidity dependent v 2 ( η ) measurements also show strong sensitivity to the QGP shear viscosity, especially the v 2 at forward rapidity shows some sensitivity to the η 2 parameter. Although the correlation is relatively weak, we find a negative correlation between charged hadron v n at 19.6 GeV and the shear viscosity at µ B = 200 MeV. And the η 4 parameter, which controls the QGP shear viscosity at µ B = 400 MeV, is negatively correlated with the charged hadron v n at 7.7 GeV.

## VIII. CONCLUSION

This work presents a systematic Bayesian analysis of particle production and flow measurements at the RHIC BES program using an event-by-event (3+1)D hydrodynamics + hadronic transport framework. We obtain robust statistical constraints on the QGP transport properties and various aspects of relativistic nuclear dynamics. This work expands upon Ref. [47] and provides the complete 20-dimensional model posterior distribution using an improved Gaussian Process model emulator. Utilizing the statistical tools developed in Ref. [48], we systemati-

FIG. 9. Response matrix R ai averaged over the posterior distribution. The observables in the left column correspond to observables at √ s NN = 200 GeV collisions. The right column shows the response to observables at √ s NN = 19 6 GeV and . 7 7 GeV. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00537v1-figure-59fe293f-4a9d-4585-946f-80182ec3b8bb.jpg)

cally analyze the model's posterior distributions obtained from three model emulators with different accuracy. Using the statistical measures, namely the Kullback-Leibler divergence and Bayes factor, we demonstrate that the accuracy of the model emulator plays a crucial role in getting robust posterior distribution from the experimental measurements. With an independent Bayesian analysis framework [88], we verify that the posterior distributions of the initial-state rapidity loss and the specific shear and bulk viscosity of QGP obtained in this work are consistent with the results in Ref. [47], demonstrating the reproducibility of our results.

We sample 100 parameter sets from the posterior distributions and perform large-scale event-by-event simu- lations at each parameter set using our (3+1)D hybrid model. These results verify that the posterior distribution obtained with the model emulator is reliable for the real model. We further make predictions for a set of p T -differential observables and estimate their systematical errors based on the variations among the 100 parameter sets. Good agreements with the experimental data are found. Based on this approach, we will perform more systematic model predictions in future work. For community users who do not have sufficient computing resources to perform simulations at multiple posterior parameter sets, we provide a high-likelihood parameter set in Table V.

Lastly, we report a sensitivity analysis of how individ-

ual experimental observables respond to different model parameters, enabled by the fast model emulator. We find commons and differences between the global sensitivity analysis over the entire prior region and the response matrix from the posterior distribution. Both sensitivity analyses provide useful physics insights into the phenomenological modeling of relativistic heavy-ion collisions. We will extend the sensitivity analysis to a more systematic experimental design study in future work.

We provide the following results from this work for community use:

- · Training datasets in Table II [89]
- · Trained model emulator PCSK (LHD + HPP) [89]
- · Posterior parameter samples from the MCMC [89]
- · An updated Bayesian analysis package with the pocoMC integration [88]
- [1] P. Achenbach et al. , The present and future of QCD, Nucl. Phys. A 1047 , 122874 (2024), arXiv:2303.02579 [hep-ph].
- [2] M. Arslandok et al. , Hot QCD White Paper, (2023), arXiv:2303.17254 [nucl-ex].
- [3] H. Caines (STAR), The RHIC Beam Energy Scan: STAR'S Perspective, in 44th Rencontres de Moriond on QCD and High Energy Interactions (2009) pp. 375-378, arXiv:0906.0305 [nucl-ex].
- [4] B. Mohanty (STAR), STAR experiment results from the beam energy scan program at RHIC, J. Phys. G 38 , 124023 (2011), arXiv:1106.5902 [nucl-ex].
- [5] J. T. Mitchell (PHENIX), The RHIC Beam Energy Scan Program: Results from the PHENIX Experiment, Nucl. Phys. A 904-905 , 903c (2013), arXiv:1211.6139 [nucl-ex].
- [6] G. Odyniec, Future of the beam energy scan program at RHIC, EPJ Web Conf. 95 , 03027 (2015).
- [7] M. M. Aggarwal et al. (STAR), An Experimental Exploration of the QCD Phase Diagram: The Search for the Critical Point and the Onset of De-confinement, (2010), arXiv:1007.2613 [nucl-ex].
- [8] X. Luo and N. Xu, Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview, Nucl. Sci. Tech. 28 , 112 (2017), arXiv:1701.02105 [nucl-ex].
- [9] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan, Phys. Rept. 853 , 1 (2020), arXiv:1906.00936 [nucl-th].
- [10] C. Shen and L. Yan, Recent development of hydrodynamic modeling in heavy-ion collisions, Nucl. Sci. Tech. 31 , 122 (2020), arXiv:2010.12377 [nucl-th].
- [11] X. An et al. , The BEST framework for the search for the QCD critical point and the chiral magnetic effect, Nucl. Phys. A 1017 , 122343 (2022), arXiv:2108.13867 [nucl-th].
- [12] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, The iEBE-VISHNU code package for relativis-

## ACKNOWLEDGMENTS

We thank Bj¨ orn Schenke and Wenbin Zhao for their involvement during the early stage of the project and useful discussions. This work is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under DOE Award No. DE-SC0021969 and DE-SC0024232. H. R. and C. S. were supported in part by the National Science Foundation (NSF) within the framework of the JETSCAPE collaboration (OAC2004571). C. S. acknowledges a DOE Office of Science Early Career Award. S.-A. J. acknowledges support from the BAND Fellowship funded by the National Science Foundation CSSI program under award number OAC2004601. Numerical simulations presented in this work were partly performed at the Wayne State Grid, and we gratefully acknowledge their support. This research was done using resources provided by the Open Science Grid (OSG) [90-93], which is supported by the National Science Foundation award #2030508 and #1836650.

- tic heavy-ion collisions, Comput. Phys. Commun. 199 , 61 (2016), arXiv:1409.8164 [nucl-th].
- [13] J. H. Putschke et al. , The JETSCAPE framework, (2019), arXiv:1903.07706 [nucl-th].
- [14] B. Schenke, C. Shen, and P. Tribedy, Running the gamut of high energy nuclear collisions, Phys. Rev. C 102 , 044905 (2020), arXiv:2005.14682 [nucl-th].
- [15] G. Nijs, W. van der Schee, U. G¨rsoy, and R. Snellings, u Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum, Phys. Rev. C 103 , 054909 (2021), arXiv:2010.15134 [nucl-th].
- [16] L.-G. Pang, H. Petersen, and X.-N. Wang, Pseudorapidity distribution and decorrelation of anisotropic flow within the open-computing-language implementation CLVisc hydrodynamics, Phys. Rev. C 97 , 064918 (2018), arXiv:1802.04449 [nucl-th].
- [17] C. Shen and B. Schenke, Longitudinal dynamics and particle production in relativistic nuclear collisions, Phys. Rev. C 105 , 064905 (2022), arXiv:2203.04685 [nucl-th].
- [18] G. D. Moore, Shear viscosity in QCD and why it's hard to calculate, in Criticality in QCD and the Hadron Resonance Gas (2020) arXiv:2010.15704 [hep-ph].
- [19] A. Nakamura and S. Sakai, Transport coefficients of gluon plasma, Phys. Rev. Lett. 94 , 072305 (2005), arXiv:heplat/0406009.
- [20] H. B. Meyer, A Calculation of the shear viscosity in SU(3) gluodynamics, Phys. Rev. D 76 , 101701 (2007), arXiv:0704.1801 [hep-lat].
- [21] N. Astrakhantsev, V. Braguta, and A. Kotov, Temperature dependence of shear viscosity of SU (3)gluodynamics within lattice simulation, JHEP 04 , 101, arXiv:1701.02266 [hep-lat].
- [22] L. Altenkort, A. M. Eller, A. Francis, O. Kaczmarek, L. Mazur, G. D. Moore, and H.-T. Shu, Viscosity of pureglue QCD from the lattice, Phys. Rev. D 108 , 014503 (2023), arXiv:2211.08230 [hep-lat].
- [23] H. Song, S. A. Bass, U. Heinz, T. Hirano, and C. Shen,

- 200 A GeV Au+Au collisions serve a nearly perfect quark-gluon liquid, Phys. Rev. Lett. 106 , 192301 (2011), [Erratum: Phys.Rev.Lett. 109, 139904 (2012)], arXiv:1011.2783 [nucl-th].
- [24] I. A. Karpenko, P. Huovinen, H. Petersen, and M. Bleicher, Estimation of the shear viscosity at finite netbaryon density from A + A collision data at √ s NN = 7 7 . -200 GeV, Phys. Rev. C 91 , 064901 (2015), arXiv:1502.01978 [nucl-th].
- [25] C. Shen and U. Heinz, The road to precision: Extraction of the specific shear viscosity of the quark-gluon plasma, Nucl. Phys. News 25 , 6 (2015), arXiv:1507.01558 [nuclth].
- [26] S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions, Phys. Rev. Lett. 115 , 132301 (2015), arXiv:1502.01675 [nuclth].
- [27] B. Schenke, C. Shen, and P. Tribedy, Multi-particle and charge-dependent azimuthal correlations in heavy-ion collisions at the Relativistic Heavy-Ion Collider, Phys. Rev. C 99 , 044908 (2019), arXiv:1901.04378 [nucl-th].
- [28] S. Ryu, J.-F. Paquet, C. Shen, G. Denicol, B. Schenke, S. Jeon, and C. Gale, Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider, Phys. Rev. C 97 , 034910 (2018), arXiv:1704.04216 [nucl-th].
- [29] H. Niemi, K. J. Eskola, and R. Paatelainen, Event-byevent fluctuations in a perturbative QCD + saturation + hydrodynamics model: Determining QCD matter shear viscosity in ultrarelativistic heavy-ion collisions, Phys. Rev. C 93 , 024907 (2016), arXiv:1505.02677 [hep-ph].
- [30] C. Shen and S. Alzhrani, Collision-geometry-based 3D initial condition for relativistic heavy-ion collisions, Phys. Rev. C 102 , 014909 (2020), arXiv:2003.05852 [nucl-th].
- [31] C. Shen, Studying QGP with flow: A theory overview, Nucl. Phys. A 1005 , 121788 (2021), arXiv:2001.11858 [nucl-th].
- [32] S. Pratt, E. Sangaline, P. Sorensen, and H. Wang, Constraining the Eq. of State of Super-Hadronic Matter from Heavy-Ion Collisions, Phys. Rev. Lett. 114 , 202301 (2015), arXiv:1501.04042 [nucl-th].
- [33] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Applying Bayesian parameter estimation to relativistic heavy-ion collisions: simultaneous characterization of the initial state and quark-gluon plasma medium, Phys. Rev. C 94 , 024907 (2016), arXiv:1605.03954 [nuclth].
- [34] J. Auvinen, J. E. Bernhard, S. A. Bass, and I. Karpenko, Investigating the collision energy dependence of η /s in the beam energy scan at the BNL Relativistic Heavy Ion Collider using Bayesian statistics, Phys. Rev. C 97 , 044905 (2018), arXiv:1706.03666 [hep-ph].
- [35] J. E. Bernhard, J. S. Moreland, and S. A. Bass, Bayesian estimation of the specific shear and bulk viscosity of quark-gluon plasma, Nature Phys. 15 , 1113 (2019).
- [36] G. Nijs, W. van der Schee, U. G¨rsoy, and R. Snellings, u Transverse Momentum Differential Global Analysis of Heavy-Ion Collisions, Phys. Rev. Lett. 126 , 202301 (2021), arXiv:2010.15130 [nucl-th].
- [37] D. Everett et al. (JETSCAPE), Multisystem Bayesian constraints on the transport coefficients of QCD matter, Phys. Rev. C 103 , 054904 (2021), arXiv:2011.01430 [hep-

ph].

- [38] D. Everett et al. (JETSCAPE), Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging, Phys. Rev. Lett. 126 , 242301 (2021), arXiv:2010.03928 [hep-ph].
- [39] J. E. Parkkila, A. Onnerstad, and D. J. Kim, Bayesian estimation of the specific shear and bulk viscosity of the quark-gluon plasma with additional flow harmonic observables, Phys. Rev. C 104 , 054904 (2021), arXiv:2106.05019 [hep-ph].
- [40] J. E. Parkkila, A. Onnerstad, S. F. Taghavi, C. Mordasini, A. Bilandzic, M. Virta, and D. J. Kim, New constraints for QCD matter from improved Bayesian parameter estimation in heavy-ion collisions at LHC, Phys. Lett. B 835 , 137485 (2022), arXiv:2111.08145 [hep-ph].
- [41] D. R. Phillips et al. , Get on the BAND Wagon: A Bayesian Framework for Quantifying Model Uncertainties in Nuclear Dynamics, J. Phys. G 48 , 072001 (2021), arXiv:2012.07704 [nucl-th].
- [42] M. R. Heffernan, C. Gale, S. Jeon, and J.-F. Paquet, Early-Times Yang-Mills Dynamics and the Characterization of Strongly Interacting Matter with Statistical Learning, Phys. Rev. Lett. 132 , 252301 (2024), arXiv:2306.09619 [nucl-th].
- [43] M. R. Heffernan, C. Gale, S. Jeon, and J.-F. Paquet, Bayesian quantification of strongly interacting matter with color glass condensate initial conditions, Phys. Rev. C 109 , 065207 (2024), arXiv:2302.09478 [nucl-th].
- [44] A. Mankolli et al. (JETSCAPE), 3D Multi-system Bayesian Calibration with Energy Conservation to Study Rapidity-dependent Dynamics of Nuclear Collisions, (2023), arXiv:2401.00402 [nucl-th].
- [45] C. Shen and B. Schenke, Initial state and hydrodynamic modeling of heavy-ion collisions at RHIC BES energies, PoS CPOD2017 , 006 (2018), arXiv:1711.10544 [nuclth].
- [46] C. Shen, Dynamic modeling for heavy-ion collisions, EPJ Web Conf. 259 , 02001 (2022), arXiv:2108.04987 [nuclth].
- [47] C. Shen, B. Schenke, and W. Zhao, Viscosities of the Baryon-Rich Quark-Gluon Plasma from Beam Energy Scan Data, Phys. Rev. Lett. 132 , 072301 (2024), arXiv:2310.10787 [nucl-th].
- [48] H. Roch, S. A. Jahan, and C. Shen, On model emulation and closure tests for 3+1D relativistic heavy-ion collisions, (2024), arXiv:2405.12019 [nucl-th].
- [49] B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, Phys. Rev. C 82 , 014903 (2010), arXiv:1004.1408 [hep-ph].
- [50] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Production of photons in relativistic heavy-ion collisions, Phys. Rev. C 93 , 044906 (2016), arXiv:1509.06738 [hep-ph].
- [51] G. S. Denicol, C. Gale, S. Jeon, A. Monnai, B. Schenke, and C. Shen, Net baryon diffusion in fluid dynamic simulations of relativistic heavy-ion collisions, Phys. Rev. C 98 , 034916 (2018), arXiv:1804.10557 [nucl-th].
- [52] S. A. Bass et al. , Microscopic models for ultrarelativistic heavy ion collisions, Prog. Part. Nucl. Phys. 41 , 255 (1998), arXiv:nucl-th/9803035.
- [53] M. Bleicher et al. , Relativistic hadron hadron collisions in the ultrarelativistic quantum molecular dynamics model, J. Phys. G 25 , 1859 (1999), arXiv:hep-ph/9909407.
- [54] C. Shen, G. Denicol, C. Gale, S. Jeon, A. Monnai, and

- B. Schenke, A hybrid approach to relativistic heavy-ion collisions at the RHIC BES energies, Nucl. Phys. A 967 , 796 (2017), arXiv:1704.04109 [nucl-th].
- [55] C. Shen and B. Schenke, Dynamical initial state model for relativistic heavy-ion collisions, Phys. Rev. C 97 , 024907 (2018), arXiv:1710.00881 [nucl-th].
- [56] C. Shen, A. Noble, J.-F. Paquet, B. Schenke, and C. Gale, Illuminating early-stage dynamics of heavy-ion collisions through photons at RHIC BES energies, PoS HardProbes2023 , 042 (2024), arXiv:2307.08498 [nucl-th].
- [57] D. Kharzeev, Can gluons trace baryon number?, Phys. Lett. B 378 , 238 (1996), arXiv:nucl-th/9602027.
- [58] W. Zhao, S. Ryu, C. Shen, and B. Schenke, 3D structure of anisotropic flow in small collision systems at energies available at the BNL Relativistic Heavy Ion Collider, Phys. Rev. C 107 , 014904 (2023), arXiv:2211.16376 [nucl-th].
- [59] A. Monnai, B. Schenke, and C. Shen, Equation of state at finite densities for QCD matter in nuclear collisions, Phys. Rev. C 100 , 024907 (2019), arXiv:1902.05095 [nucl-th].
- [60] G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, Derivation of transient relativistic fluid dynamics from the Boltzmann equation, Phys. Rev. D 85 , 114047 (2012), [Erratum: Phys.Rev.D 91, 039902 (2015)], arXiv:1202.4551 [nucl-th].
- [61] A. Bazavov et al. (HotQCD), Chiral crossover in QCD at zero and non-zero chemical potentials, Phys. Lett. B 795 , 15 (2019), arXiv:1812.08235 [hep-lat].
- [62] S. Borsanyi, Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz, P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo, QCD Crossover at Finite Chemical Potential from Lattice Simulations, Phys. Rev. Lett. 125 , 052001 (2020), arXiv:2002.02821 [hep-lat].
- [63] W. A. Hiscock and L. Lindblom, Stability and causality in dissipative relativistic fluids, Annals Phys. 151 , 466 (1983).
- [64] C. Chiu and C. Shen, Exploring theoretical uncertainties in the hydrodynamic description of relativistic heavy-ion collisions, Phys. Rev. C 103 , 064901 (2021), arXiv:2103.09848 [nucl-th].
- [65] P. Huovinen and H. Petersen, Particlization in hybrid models, Eur. Phys. J. A 48 , 171 (2012), arXiv:1206.3371 [nucl-th].
- [66] M. Plumlee, O. S¨rer, S. M. Wild, and M. Y.-H. Chan, u surmise 0.2.1 Users Manual , Tech. Rep. Version 0.2.1 (NAISE, 2023).
- [67] S. Brooks, A. Gelman, G. Jones, and X. Meng, Handbook of Markov Chain Monte Carlo , Chapman &amp; Hall/CRC Handbooks of Modern Statistical Methods (CRC Press, 2011).
- [68] R. Trotta, Bayes in the sky: Bayesian inference and model selection in cosmology, Contemporary Physics 49 , 71-104 (2008).
- [69] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The mcmc hammer, Publications of the Astronomical Society of the Pacific 125 , 306-312 (2013).
- [70] M. Karamanis, F. Beutler, J. A. Peacock, D. Nabergoj, and U. Seljak, Accelerating astronomical and cosmological inference with preconditioned monte carlo, Monthly Notices of the Royal Astronomical Society 516 , 1644 (2022).
- [71] M. Karamanis, D. Nabergoj, F. Beutler, J. A. Peacock, and U. Seljak, pocomc: A python package for accelerated
- bayesian inference in astronomy and cosmology, arXiv preprint arXiv:2207.05660 (2022).
- [72] B. I. Abelev et al. (STAR), Systematic Measurements of Identified Particle Spectra in pp, d + Au and Au+Au Collisions from STAR, Phys. Rev. C 79 , 034909 (2009), arXiv:0808.2041 [nucl-ex].
- [73] B. B. Back et al. (PHOBOS), Charged-particle pseudorapidity distributions in Au+Au collisions at s NN ( ) 1 / 2 = 62.4-GeV, Phys. Rev. C 74 , 021901 (2006), arXiv:nuclex/0509034.
- [74] B. Alver et al. (PHOBOS), System size, energy, pseudorapidity, and centrality dependence of elliptic flow, Phys. Rev. Lett. 98 , 242302 (2007), arXiv:nucl-ex/0610037.
- [75] L. Adamczyk et al. (STAR), Harmonic decomposition of three-particle azimuthal correlations at energies available at the BNL Relativistic Heavy Ion Collider, Phys. Rev. C 98 , 034918 (2018), arXiv:1701.06496 [nucl-ex].
- [76] L. Adamczyk et al. (STAR), Beam Energy Dependence of the Third Harmonic of Azimuthal Correlations in Au+AuCollisions at RHIC, Phys. Rev. Lett. 116 , 112302 (2016), arXiv:1601.01999 [nucl-ex].
- [77] L. Adamczyk et al. (STAR), Bulk Properties of the Medium Produced in Relativistic Heavy-Ion Collisions from the Beam Energy Scan Program, Phys. Rev. C 96 , 044904 (2017), arXiv:1701.07065 [nucl-ex].
- [78] R. Trotta, Bayes in the sky: Bayesian inference and model selection in cosmology, Contemporary Physics 49 , 71 (2008), https://doi.org/10.1080/00107510802066753.
- [79] G. Pihan, A. Monnai, B. Schenke, and C. Shen, Unveiling baryon charge carriers through charge stopping in isobar collisions, (2024), arXiv:2405.19439 [nucl-th].
- [80] C. Shen, B. Schenke, and W. Zhao, The effects of pseudorapidity-dependent observables on (3+1)D Bayesian Inference of relativistic heavy-ion collisions, EPJ Web Conf. 296 , 14001 (2024), arXiv:2312.09325 [nucl-th].
- [81] L. Adamczyk et al. (STAR), Inclusive charged hadron elliptic flow in Au + Au collisions at √ s NN = 7.7 - 39 GeV, Phys. Rev. C 86 , 054908 (2012), arXiv:1206.5528 [nucl-ex].
- [82] B. I. Abelev et al. (STAR), Centrality dependence of charged hadron and strange hadron elliptic flow from s(NN)**(1/2) = 200-GeV Au + Au collisions, Phys. Rev. C 77 , 054901 (2008), arXiv:0801.3466 [nucl-ex].
- [83] A. Sch¨fer, I. Karpenko, X.-Y. Wu, J. Hammelmann, and a H. Elfner (SMASH), Particle production in a hybrid approach for a beam energy scan of Au+Au/Pb+Pb collisions between √ s NN = 4.3 GeV and √ s NN = 200.0 GeV, Eur. Phys. J. A 58 , 230 (2022), arXiv:2112.08724 [hepph].
- [84] G. Denicol, A. Monnai, and B. Schenke, Moving forward to constrain the shear viscosity of QCD matter, Phys. Rev. Lett. 116 , 212301 (2016), arXiv:1512.01538 [nuclth].
- [85] I. M. Sobol', On sensitivity estimation for nonlinear mathematical models, Matematicheskoe modelirovanie 2 , 112 (1990).
- [86] D. Liyanage, Y. Ji, D. Everett, M. Heffernan, U. Heinz, S. Mak, and J.-F. Paquet, Efficient emulation of relativistic heavy ion collisions with transfer learning, Phys. Rev. C 105 , 034910 (2022), arXiv:2201.07302 [nucl-th].
- [87] E. Sangaline and S. Pratt, Toward a deeper understanding of how experiments constrain the underlying physics

- of heavy-ion collisions, Phys. Rev. C 93 , 024908 (2016), arXiv:1508.07017 [nucl-th].
- [88] H. Roch and S. A. Jahan, Hendrik1704/gpbayestools-hic: v1.1.0 (2024).
- [89] S. A. Jahan, H. Roch, and C. Shen, Bayesian analysis of (3+1)D relativistic nuclear dynamics with the RHIC beam energy scan data, 10.5281/zenodo.12807556 (2024).
- [90] R. Pordes et al. , The Open Science Grid, J. Phys. Conf.
- Ser. 78 , 012057 (2007).
- [91] I. Sfiligoi, D. C. Bradley, B. Holzman, P. Mhashilkar, S. Padhi, and F. Wurthwrin, The pilot way to Grid resources using glideinWMS, WRI World Congress 2 , 428 (2009).
- [92] OSG, Ospool (2006).
- [93] OSG, Open science data federation (2015).
Syed Afrid Jahan Hendrik Roch Chun Shen
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.