# On model emulation and closure tests for 3+1D relativistic heavy-ion
collisions

Hendrik Roch<EMAIL_ADDRESS>Department of Physics and Astronomy,
Wayne State University, Detroit, Michigan 48201, USA  Syed Afrid Jahan
<EMAIL_ADDRESS>Department of Physics and Astronomy, Wayne State University,
Detroit, Michigan 48201, USA  Chun Shen<EMAIL_ADDRESS>Department of
Physics and Astronomy, Wayne State University, Detroit, Michigan 48201, USA

###### Abstract

In nuclear and particle physics, reconciling sophisticated simulations with
experimental data is vital for understanding complex systems like the Quark
Gluon Plasma (QGP) generated in heavy-ion collisions. However, computational
demands pose challenges, motivating using Gaussian Process emulators for
efficient parameter extraction via Bayesian calibration. We conduct a
comparative analysis of Gaussian Process emulators in heavy-ion physics to
identify the most adept emulator for parameter extraction with minimal
uncertainty. Our study contributes to advancing computational techniques in
heavy-ion physics, enhancing our ability to interpret experimental data and
understand QGP properties.

##  I Introduction

In nuclear and particle physics, comparing complex simulations with multiple
experimental observables is important for elucidating the properties of
systems under study. A quintessential example lies in nuclear physics, which
investigates the Quark Gluon Plasma (QGP) via relativistic heavy-ion
collisions. Experiments at the Relativistic Heavy Ion Collider (RHIC) and the
Large Hadron Collider (LHC) aim to recreate the extreme conditions akin to
those prevailing in the early Universe shortly after the Big Bang Busza  _et
al._ (2018).

However, direct measurement of the resultant hot and dense QGP state remains
elusive due to its transient existence and confinement to small spatial
dimensions within collider experiments. Understanding the QGP’s properties
necessitates analyzing experimental measurements, which imprint the entire
collision history as the QGP is created, expands, and fragments into a
collection of hadrons Busza  _et al._ (2018). Theoretical frameworks for
modeling such collisions have undergone substantial advancement over recent
decades, encompassing various stages from the initial state Miller  _et al._
(2007); Loizides  _et al._ (2015); Broniowski  _et al._ (2009); Drescher and
Nara (2007); Schenke  _et al._ (2012); Paatelainen  _et al._ (2013); Moreland
_et al._ (2015); Giacalone  _et al._ (2019); Shen and Schenke (2018) to
hadronic rescattering dynamics Bass  _et al._ (1998); Weil  _et al._ (2016);
Cassing and Bratkovskaya (1999); Lin  _et al._ (2005), including pre-
equilibrium Schlichting and Teaney (2019), relativistic viscous hydrodynamics
Teaney (2010); Romatschke and Romatschke (2019), and particlization Cooper and
Frye (1974).

These models, containing numerous constituent elements, demand careful
calibration against experimental data, entailing the tuning of many parameters
simultaneously to extract physics from the observations. Moreover, the
computational demands pose a substantial challenge, with long simulation
runtime required even for a single event simulation. By inferring from
experimental data, this model parameter determination embodies an inverse
problem paradigm.

Bayesian inference, leveraging Bayes’s theorem, offers a principled approach
to address inverse problems by discerning parameter sets and their associated
uncertainties systematically Sivia and Skilling (2006). However, the
exhaustive evaluation of models across diverse parameter spaces often proves
infeasible due to computational constraints. Herein lies the utility of model
emulation, which, once trained with model outputs, furnish computationally
inexpensive surrogates for the full model. These emulators, by virtue of their
efficiency, replace the full model in Bayesian analyses, thereby facilitating
parameter calibration.

This study compares various Gaussian Process (GP) emulators for (3+1)D
simulations of relativistic heavy-ion physics in the RHIC Beam Energy Scan
program. Our objective is to discern the most adept GP emulator for a given
training dataset, characterized by minimal prediction uncertainty, for
parameter extraction via Bayesian calibration. Additionally, we aim to
quantify the influence of GP emulators with different accuracies on the model
parameter extraction through closure tests.

In Sec. II, we concisely discuss the requisite components for the Bayesian
inference approach for inverse problems. Following this, Sec. III presents
fundamental theories for training the GP emulators. Subsequently, in Sec. IV,
we compare three open-source GP emulator implementations, where we look at the
uncertainty and quality of their prediction in different scenarios. Then,
Section V analyzes the posterior distributions obtained from these different
GPs in a closure test. We also perform a local sensitivity analysis for the
observables under consideration in Sec. VI. Finally, we will conclude and give
an outlook in Sec. VII.

##  II Bayesian Inference

This section provides an overview of the application of Bayesian Inference
with multi-stage Monte Carlo event simulations in relativistic heavy-ion
collisions, highlighting their utility in mitigating the computational costs
associated with such simulations.

Multi-stage event-by-event simulation frameworks, such as iEBE-VISHNU Shen
_et al._ (2016); Bernhard  _et al._ (2016), JETSCAPE Putschke  _et al._
(2019), Trajectum Nijs  _et al._ (2021a), and iEBE-MUSIC Shen and Schenke
(2022), are essential tools for studying various physics phenomena in
relativistic heavy-ion collisions. However, these simulations are often
computationally intensive. To alleviate this computational burden, emulators
are employed to predict the model output for a given parameter set Kennedy and
O’Hagan (2001); Higdon  _et al._ (2004).

Before delving into the specifics of the emulators compared in this study, we
introduce the concept of Bayesian model calibration. Let
$\bm{\theta}\equiv(\theta_{1},\dots,\theta_{m})$ denote the parameters of the
model, and let ${\bf y}_{\rm sim}(\cdot)$ represent the model function. This
function maps the parameter set $\bm{\theta}$ to simulation outputs in
$\mathbb{R}^{d}$, where $d$ is the dimension of the observable space. These
simulation outputs are then compared to experimental values, denoted by ${\bf
y}_{\rm exp}\equiv(y_{\rm exp,1},\dots,y_{{\rm exp},d})$. Thus, a statistical
model can be formulated as:

${\bf y}_{\rm exp}={\bf y}_{\rm sim}(\bm{\theta})+\bm{\epsilon},$ (1)

where $\bm{\epsilon}$ represents the residual error following a multivariate
normal (MVN) distribution with mean $\mathbf{0}$ and covariance matrix
$\bm{\Sigma}$.

The problem at hand is inherently an inverse problem, where we seek to
identify the parameters $\bm{\theta}$ that best describe the observed data. To
address this, we adopt a Bayesian approach to inference, where probability is
interpreted as the degree of belief accounting for all available information
Sivia and Skilling (2006). Bayes’ rule provides a framework for updating our
beliefs based on observed data:

$\mathcal{P}(\bm{\theta}|{\bf y}_{\rm exp})=\frac{\mathcal{P}({\bf y}_{\rm
exp}|\bm{\theta})\mathcal{P}(\bm{\theta})}{\mathcal{P}({\bf y}_{\rm exp})}.$
(2)

Here, $\mathcal{P}(\bm{\theta}|{\bf y}_{\rm exp})$ represents the posterior
probability density function, $\mathcal{P}(\bm{\theta})$ is the prior
probability density function, $\mathcal{P}({\bf y}_{\rm exp})$ stands for the
total evidence, and $\mathcal{P}({\bf y}_{\rm exp}|\bm{\theta})$ is the
likelihood function, which, given the MVN distribution of $\bm{\epsilon}$, can
be expressed as:

$\frac{1}{\sqrt{|2\pi\bm{\Sigma}|}}\exp\left[-\frac{1}{2}({\bf y}_{\rm
sim}(\bm{\theta})-{\bf y}_{\rm exp})^{\mathsf{T}}\bm{\Sigma}^{-1}({\bf y}_{\rm
sim}(\bm{\theta})-{\bf y}_{\rm exp})\right].$ (3)

To obtain the posterior distribution, which is often intractable analytically,
Markov Chain Monte Carlo (MCMC) techniques are employed Brooks  _et al._
(2011); Trotta (2008). These methods enable the estimation of posterior
distributions for $\bm{\theta}$ and quantifying their uncertainties.

In practice, evaluating the right-hand side of Eq. (2) for numerous
$\bm{\theta}$ sets in MCMC simulations can be computationally prohibitive.
Here, GP emulators offer a solution by providing computationally efficient
surrogates for the expensive simulations Gramacy (2020); Rasmussen (2004).
These emulators are trained using simulation data for various parameter sets
$\bm{\theta}$. Once trained, the emulators can predict the mean
$\bm{\mu}(\bm{\theta})$ and covariance ${\bf C}(\bm{\theta})$ of the
simulation output ${\bf y}_{\rm exp}(\bm{\theta})$. Subsequently, the
likelihood function (Eq. (3)) can be approximated as Kennedy and O’Hagan
(2001):

$\frac{1}{\sqrt{|2\pi{\bf
V}(\bm{\theta})|}}\exp\left[-\frac{1}{2}(\bm{\mu}(\bm{\theta})-{\bf y}_{\rm
exp})^{\mathsf{T}}{\bf V}(\bm{\theta})^{-1}(\bm{\mu}(\bm{\theta})-{\bf y}_{\rm
exp})\right],$ (4)

where ${\bf V}(\bm{\theta})={\bf C}(\bm{\theta})+\bm{\Sigma}$. This
approximation facilitates the MCMC sampling process, enabling efficient
exploration of the posterior probability distribution function.

##  III Gaussian Process Emulators

This section presents a comparative analysis of different implementations of
GP emulators employed in our study. Specifically, we investigate two GP
emulator methodologies, namely the PCGP and the PCSK methods, which are
integrated into the Python package surmise by the BAND collaboration Plumlee
_et al._ (2023). Additionally, we consider the standard GP emulator available
in the Scikit-learn Python package. The Scikit GP was used in a few previous
Bayesian analyses in our field Bernhard  _et al._ (2019); Everett  _et al._
(2021a); Mäntysaari  _et al._ (2022); Shen  _et al._ (2024, 2023).

All these emulators share a basis vector approach Higdon  _et al._ (2008),
tailored to address the challenges posed by the high-dimensional output space
of the model. Principal Component Analysis (PCA) Ramsay and Silverman (1997)
is utilized to project the high-dimensional output onto a lower-dimensional
space, enabling efficient emulation. Specifically, each output for a principal
component is independently trained by a GP emulator.

We begin by elucidating the fundamental concept behind GP emulation with PCA.
Let $\\{\bm{\theta}_{1}^{\text{tr}},\dots,\bm{\theta}_{n}^{\text{tr}}\\}$
represent a training parameter set comprising $n$ unique training points for
the emulators. Each training point $\bm{\theta}_{i}^{\text{tr}}$ generates a
$d$-dimensional vector
$\bar{\mathbf{y}}_{\text{sim}}(\bm{\theta}_{i}^{\text{tr}})$, where $d$
denotes the number of observables. The model simulation results
$\bar{\mathbf{y}}_{\text{sim}}(\bm{\theta}_{i}^{\text{tr}})$ are averaged over
multiple collision events. Before applying the PCA transformation, we
standardize the training data as, $\tilde{\bf y}_{\rm
sim}(\bm{\theta}_{i}^{\rm tr})\equiv\mathbf{G}^{-1}(\bar{\bf y}_{\rm
sim}(\bm{\theta}_{i}^{\rm tr})-\bf h)$, using a centering vector $\mathbf{h}$
and scaling matrix $\mathbf{G}=\mathrm{diag}(\bf c)$. The resulting
$\tilde{\bf y}_{\rm sim}(\bm{\theta}_{i}^{\rm tr})$ has a zero mean and unity
variance over the $n$ training samples. We define a $d\times n$ matrix
$\mathbf{\Xi}\equiv\\{\tilde{\mathbf{y}}_{\text{sim}}(\bm{\theta}_{1}^{\text{tr}}),\dots,\tilde{\mathbf{y}}_{\text{sim}}(\bm{\theta}_{n}^{\text{tr}})\\}$.
PCA yields a linear transformation matrix ${\bf S}=\left[{\bf
s}_{1},\dots,{\bf s}_{p}\right]$ of dimension $d\times p$, which can project
the output matrix $\bm{\Xi}$ from a $d$-dimensional space to a $p$-dimensional
principal component space via the transformation ${\bf T}=\left[{\bf
t}_{1},\dots,{\bf t}_{p}\right]=\bf S^{\mathsf{T}}\bm{\Xi}$. In most
applications, keeping only $p\ll d$ principal components to describe almost
all the variance in the original dataset is sufficient Higdon  _et al._
(2015); Sürer  _et al._ (2022); Liyanage  _et al._ (2023); Moses Y.-H. Chan
and Wild (2024).

Subsequently, an independent GP emulator is constructed for each of the
outputs $\mathbf{t}_{l}(\bm{\theta})={\bf s}_{l}^{\mathsf{T}}\tilde{\bf
y}_{\rm sim}(\bm{\theta}_{i}^{\rm tr})$. The emulator predicts the mean
$m_{l}(\bm{\theta})$ and variance $s^{2}_{l}(\bm{\theta})$ for each principal
component $t_{l}$ with $l=1,\dots,p$, given by

$t_{l}(\bm{\theta})|\bm{t}_{l}\approx\text{N}(m_{l}(\bm{\theta}),s^{2}_{l}(\bm{\theta})),$
(5)

where
$m_{l}(\bm{\theta})=\mathbf{k}_{l}^{\mathsf{T}}\mathbf{K}_{l}^{-1}\mathbf{t}_{l}$
and
$s^{2}_{l}(\bm{\theta})=k_{l}(\bm{\theta},\bm{\theta})-\mathbf{k}_{l}^{\mathsf{T}}(\bm{\theta})\mathbf{K}_{l}^{-1}\mathbf{k}_{l}(\bm{\theta})$.
Here, $\mathbf{k}_{l}(\bm{\theta})=\left[k_{l}(\bm{\theta},\bm{\theta}^{\rm
tr}_{i})\right]^{n}_{i=1}$ denotes the covariance vector between $n$ training
parameters $\\{\bm{\theta}_{1}^{\rm tr},\dots,\bm{\theta}_{n}^{\rm tr}\\}$ and
any chosen parameter $\bm{\theta}$, and
$\mathbf{K}_{l}=\left[k_{l}(\bm{\theta}_{i}^{\rm tr},\bm{\theta}_{j}^{\rm
tr})+\delta_{i,j}r_{l,i}\right]_{i,j=1}^{n}$ represents the covariance matrix
between the $n$ training parameters. The covariance function
$k_{l}(\bm{\theta},\bm{\theta}^{\prime})$ depends on the choice of the kernel
function in the GP.111The PCSK and PCGP emulators use a Matérn kernel
(smoothness parameter $3/2$), while we choose to use a Radial Basis Function
(RBF) kernel for the Scikit GP. The Scikit GP setup was used in a few of the
previous Bayesian analyses Everett  _et al._ (2021a); Shen  _et al._ (2024).
Additionally, the term $\delta_{i,j}r_{l,i}$, along with the Kronecker delta,
accounts for uncertainty arising from stochastic simulation results for the
different events. Here $r_{l,i}$ is the square of the statistical uncertainty
of the $l$-th principal component at the training point $\bm{\theta}_{i}$.

The term $\delta_{i,j}r_{l,i}$ is the main difference between the PCGP/Scikit
GP implementations of GPs used in almost all heavy-ion physics related studies
Petersen  _et al._ (2011); Novak  _et al._ (2014); Sangaline and Pratt (2016);
Bernhard  _et al._ (2015, 2016); Moreland  _et al._ (2020); Bernhard  _et al._
(2019); Everett  _et al._ (2021b); Nijs  _et al._ (2021b, a); Everett  _et
al._ (2021a); Heffernan  _et al._ (2023a, b) and the PCSK method used in
Liyanage  _et al._ (2023). Instead of only taking into account the mean values
of the training data, the PCSK method also considers the standard deviation at
the training points. This is particularly relevant as simulations of heavy-ion
collisions with different initial states for given parameters $\bm{\theta}$
yield observables with a certain variation. The PCSK emulation incorporates
this variation in the accuracy between the resulting data points, known as
heteroskedasticity, by using stochastic kriging Ankenman  _et al._ (2010). We
explore the effects of the size of statistical errors on the trained emulator
accuracy in Appendix A.

After training, the emulators predict the mean $m_{l}(\bm{\theta})$ and
variance $s^{2}_{l}(\bm{\theta})$ for each observable at a given point
$\bm{\theta}$. The outputs for each of the GP emulators $l=1,\dots,p$ are
transformed back to the original $d$-dimensional space using the inverse PCA
transformation. Thus, the $p$-dimensional emulator prediction mean
$\bm{m}(\bm{\theta})=(m_{1}(\bm{\theta}),\dots,m_{p}(\bm{\theta}))$ and a
$p\times p$ covariance matrix $\bm{\sigma}(\bm{\theta})$ with diagonal
elements $s^{2}_{l}(\bm{\theta})$ is used, such that the inverse PCA
transformation yields

$\bm{\mathcal{G}}(\bm{\theta})\approx\mathrm{MVN}(\bm{\mu}(\bm{\theta}),{\bf
C}(\bm{\theta})),$ (6)

where $\bm{\mu}(\bm{\theta})={\bf h}+{\bf G}{\bf S}\bm{m}(\bm{\theta})$ and
${\bf C}(\bm{\theta})={\bf G}{\bf S}{\bm{\sigma}(\bm{\theta})}{\bf
S}^{\mathsf{T}}{\bf G}$ are the emulator prediction mean and covariance
respectively. The likelihood at a parameter space point $\bm{\theta}$ is
obtained using Eq. (4), and Eq. (2) is subsequently utilized to derive the
posterior probability distribution
$\mathcal{P}(\bm{\theta}|\mathbf{y}_{\text{exp}})$.

##  IV GP Emulator Comparison

In Sec. IV.1 we will describe our training setups in this work, additionally
we are going to introduce a PCA transformation for functional model parameters
in Sec. IV.2 and the characteristic measures for the emulator precision in
Sec. IV.3. We will discuss the validation of the GP emulators for different
scenarios in Sec. IV.4.

###  IV.1 Training Setup

In this work, we divide the experimental observables into several groups and
apply the principal component Gaussian Process Regression for each group
Liyanage  _et al._ (2023); Shen  _et al._ (2024). In our work, we train GP
emulators for different observable groups first separated by the system’s
collision energy. We further split the particle multiplicities observables
from the anisotropic flow coefficients and average transverse momenta because
we will test the GP emulators’ accuracy for training with multiplicities or
the logarithm of the data. The logarithmic transformation was adopted by a few
previous works Everett  _et al._ (2021a); Nijs  _et al._ (2021a).

Table 1 shows the different observables and the number of data points we
grouped for the training in the GP emulators for this analysis used before in
Ref. Shen  _et al._ (2024).

Table 1: The training data from Au+Au collisions used in this study. The datasets are separated in a way, which is later on used for the training of the emulators. The last column shows the number of data points for the training set. Training Set | Observables | Points
---|---|---
1 | $\mathrm{d}N/\mathrm{d}y|_{\pi^{+},K^{+},p}^{7.7\;\rm{GeV}}$ | 21
2 |  $\langle p_{\rm T}\rangle|_{\pi^{+},K^{+},p,\bar{p}}^{7.7\;\rm{GeV}}$, $v_{2/3}^{\rm ch}\\{2\\}|^{7.7\;\rm{GeV}}$ | 40
3 | $\mathrm{d}N/\mathrm{d}y|_{\pi^{+},K^{+},p}^{19.6\;\rm{GeV}}$ | 21
4 |  $\langle p_{\rm T}\rangle|_{\pi^{+},K^{+},p,\bar{p}}^{19.6\;\rm{GeV}}$, $v_{2/3}^{\rm ch}\\{2\\}|^{19.6\;\rm{GeV}}$ | 40
5 | $\mathrm{d}N/\mathrm{d}y|_{\pi^{+},K^{+},p,\bar{p}}^{200\;\rm{GeV}}$ | 28
6 |  $\langle p_{\rm T}\rangle|_{\pi^{+},K^{+},p,\bar{p}}^{200\;\rm{GeV}}$, $v_{2/3}^{\rm ch}\\{2\\}|^{200\;\rm{GeV}}$ | 40
7 | $\mathrm{d}N_{\rm ch}/\mathrm{d}\eta|^{19.6\;\rm{GeV}}$ | 170
8 | $\mathrm{d}N_{\rm ch}/\mathrm{d}\eta|^{200\;\rm{GeV}}$ | 170
9 | $v_{2}^{\rm ch}(\eta)|^{200\;\rm{GeV}}$ | 14

The datasets contain event-by-event simulations of Au+Au collisions at three
energies: 200, 19.6, and 7.7 GeV. The training dataset contains 1000 design
points determined with a Maximum Projection Latin Hypercube Design algorithm
Roshan Joseph  _et al._ (2020, 2015) to explore the $\bm{\theta}$-space. This
algorithm orders the design points to ensure uniform coverage of model
parameter space with a subset of design points. This feature allows us to
quantify how the emulator accuracy scales with the size of training datasets.
In the following, we denote the design points from this algorithm by “LHD
points”. For each design point, 1000 minimum bias events were simulated at 200
GeV and 2000 minimum bias events at 19.6 and 7.7 GeV each.222Points with
relative statistical uncertainty larger than 10% are discarded before training
the GP emulators. Additionally, 100 parameter sets were drawn from the
resulting posterior distributions in Ref. Shen  _et al._ (2024), and
simulations were carried out for these points. This allows us to study, in an
active learning scenario, whether the uncertainty in the emulators can be
decreased when the density of training points is increased in regions where
the likelihood of the final parameters is more significant.

The simulations were performed with a 20-dimensional parameter space,
summarized in Table 2 with the ranges of the uniform priors. Their definitions
and the physics behind them can be found in Ref. Shen  _et al._ (2024). In
this work, we want to focus on the accuracy of the GP emulators and their
impact on Bayesian parameter extraction in closure tests.

Table 2: The 20 model parameters and their prior ranges. Parameter | Prior | Parameter | Prior
---|---|---|---
$B_{G}\;[\mathrm{GeV}^{-2}]$ | $[1,25]$ | $\alpha_{\text{string tilt}}$ | $[0,1]$
$\alpha_{\rm shadowing}$ | $[0,1]$ | $\alpha_{\text{preFlow}}$ | $[0,2]$
$y_{{\rm loss},2}$ | $[0,2]$ | $\eta_{0}$ | $[0.001,0.3]$
$y_{{\rm loss},4}$ | $[1,3]$ | $\eta_{2}$ | $[0.001,0.3]$
$y_{{\rm loss},6}$ | $[1,4]$ | $\eta_{4}$ | $[0.001,0.3]$
$\sigma_{y_{\rm loss}}$ | $[0.1,0.8]$ | $\zeta_{\rm max}$ | $[0,0.2]$
$\alpha_{\rm rem}$ | $[0,1]$ | $T_{\zeta,0}\;[{\rm GeV}]$ | $[0.15,0.25]$
$\lambda_{B}$ | $[0,1]$ | $\sigma_{\zeta,+}\;[{\rm GeV}]$ | $[0.01,0.15]$
$\sigma_{x}^{\rm string}\;[{\rm fm}]$ | $[0.1,0.8]$ | $\sigma_{\zeta,-}\;[{\rm GeV}]$ | $[0.005,0.1]$
$\sigma_{\eta}^{\rm string}$ | $[0.1,1]$ | $e_{\rm sw}\;[{\rm GeV}/{\rm fm}^{3}]$ | $[0.15,0.5]$

###  IV.2 PCA Transformation for Functional Model Parameters

Our model simulations take in functional inputs for some key physics
quantities of interest, namely the average initial-state rapidity loss
$\langle y_{\rm loss}\rangle(y_{\mathrm{init}})$ and QGP specific shear and
bulk viscosities $\tilde{\eta}(\mu_{B})$ and $\tilde{\zeta}(T,\mu_{B})$ Shen
_et al._ (2024). The average initial-state rapidity loss for each parton pair
is parametrized as a piece-wise function of the incoming rapidity $y_{\rm
init}$ in the collision pair rest frame,

$\langle y_{\rm loss}\rangle=\begin{cases}y_{\rm loss,2}\frac{y_{\rm
init}}{2},\;0<y_{\rm init}\leq 2\\\ y_{\rm loss,2}+(y_{\rm loss,4}-y_{\rm
loss,2})\frac{y_{\rm init}-2}{2},\;2<y_{\rm init}<4\\\ y_{\rm loss,4}+(y_{\rm
loss,6}-y_{\rm loss,4})\frac{y_{\rm init}-4}{2},\;y_{\rm init}\geq
4\end{cases}$ (7)

with the parameter $y_{{\rm loss},n}$ being the average amount of rapidity
loss for $y_{\rm init}=n$. The $\mu_{B}$-dependence of the QGP’s specific
shear viscosity is parametrized by

$\tilde{\eta}(\mu_{B})=\begin{cases}\eta_{0}+(\eta_{2}-\eta_{0})\frac{\mu_{B}}{0.2},\;0<\mu_{B}\leq
0.2\;\mathrm{GeV}\\\
\eta_{2}+(\eta_{4}-\eta_{2})\frac{(\mu_{B}-0.2)}{0.2},\;0.2<\mu_{B}<0.4\;\mathrm{GeV}\\\
\eta_{4},\;\mu_{B}\geq 0.4\;\mathrm{GeV}\end{cases}$ (8)

where $\tilde{\eta}=\eta T/(e+P)$ with the parameters $\eta_{0}$, $\eta_{2}$,
$\eta_{4}$ are the values of the specific shear viscosity at
$\mu_{B}=0,0.2,0.4\;\mathrm{GeV}$ respectively. The functional dependence for
the QGP’s specific bulk viscosity $\tilde{\zeta}$ is assumed to be an
asymmetric Gaussian of the form

$\tilde{\zeta}(T,\mu_{B})=\begin{cases}\zeta_{\rm
max}\exp\left[-\frac{(T-T_{\zeta}(\mu_{B}))^{2}}{2\sigma^{2}_{\zeta,-}}\right],\;T<T_{\zeta}(\mu_{B})\\\
\zeta_{\rm
max}\exp\left[-\frac{(T-T_{\zeta}(\mu_{B}))^{2}}{2\sigma^{2}_{\zeta,+}}\right],\;T\geq
T_{\zeta}(\mu_{B})\end{cases},$ (9)

where $\tilde{\zeta}\equiv\zeta T/(e+P)$ with the bulk pressure peak
temperature $T_{\zeta}(\mu_{B})\equiv
T_{\zeta,0}-\frac{0.15}{1\;\mathrm{GeV}}\mu_{B}^{2}$.

Although the model parameters in Table 2 provide plenty of information to the
GP emulators, the functional form of the parametrizations is not included in
the training. One way to provide the GP emulators with this information is to
enlarge the dimension of the model parameter space and feed the GP emulators
with more sampled points for $\langle y_{\rm loss}\rangle(y_{\mathrm{init}})$
and QGP specific shear and bulk viscosities $\tilde{\eta}(\mu_{B})$ and
$\tilde{\zeta}(T,\mu_{B})$. This approach was carried out in the previous
JETSCAPE analysis and effectively increased emulator accuracy Everett  _et
al._ (2021a, b). Here, we introduce a more systematic way to deal with
functional model parameters using PCA transformation.

We use the parametrizations in Eqs. (7)-(9) to generate the corresponding
functions at each training point $\bm{\theta}_{i}^{\text{tr}}$ with 100
equally spaced points in the range $T\in\left[0,0.5\right]\;\mathrm{GeV}$ for
$\tilde{\zeta}$, $\mu_{B}\in\left[0,0.6\right]\;\mathrm{GeV}$ for
$\tilde{\eta}$ and $y_{\rm init}\in\left[0,6.2\right]$ for $\langle y_{\rm
loss}\rangle$. Then, we perform the standardization and PCA similar to the one
described in Sec. III. We proceed with as many principal components as
necessary to explain 99 percent of the variance of each parametrization and
substitute the four parameters $\zeta_{\rm max}$, $T_{\zeta,0}$,
$\sigma_{\zeta,+}$ and $\sigma_{\zeta,+}$ by the needed number of principal
components. Similarly, we substitute the $\eta_{0}$, $\eta_{2}$, $\eta_{4}$
parameters for the shear viscosity and the three rapidity loss parameters
$y_{\rm loss,2}$, $y_{\rm loss,4}$ and $y_{\rm loss,6}$.333The bulk viscosity
parameterized in Eq. (9) needs six principal components, while three
components are sufficient in the QGP specific shear viscosity in Eq. (8) and
averaged initial-state rapidity loss in Eq. (7). The emulator is then trained
with the principal components instead of the parameters of the
parametrization. When the prediction with a parameter set $\bm{\theta}$ is
performed, the original parameters are then transformed into the PCA space
with the same transformation, such that the transformation is completely
encapsulated in the emulator.

###  IV.3 Test Metrics for Model Emulators

To quantify the prediction error of the GP emulators, we define

$\mathcal{E}\equiv\sqrt{\left\langle\left(\frac{\text{prediction}-\text{truth}}{\text{truth}}\right)^{2}\right\rangle}$
(10)

for each observable in the analysis. The GP emulator provides the covariance
matrix at the prediction point, which is used as an estimation for the
prediction uncertainty. To test how reliable the GP emulators’ uncertainty is,
we define

$\mathcal{H}\equiv\ln\left(\sqrt{\left\langle\left(\frac{\text{prediction}-\text{truth}}{\text{prediction
uncertainty}}\right)^{2}\right\rangle}\right).$ (11)

For an accurate prediction, we expect the values of $\mathcal{E}\rightarrow 0$
and $\mathcal{H}\rightarrow 0$. In the case where $\mathcal{H}>0$, the
emulator gives uncertainties that are too small compared to the actual error
away from the true values; when $\mathcal{H}<0$, the returned uncertainty
estimates are too conservative.

###  IV.4 Emulator Validation

We will now validate the results of the three different emulators by
investigating $\mathcal{E}$ and $\mathcal{H}$ in various setups. First, we
will look at the dependence in terms of the number of LHD training points
(Sec. IV.4.1), and then consider the values including most of the LHD points;
we compare the performance of the three types of the emulator in Sec. IV.4.2.
We further will investigate the ansatz with the logarithmic training for the
multiplicity observables (Sec. IV.4.3), followed by the inclusion of
additional PCA transformations for parametrized curves (Sec. IV.4.4). Lastly,
we will discuss the active learning approach by including 70 additional points
from the posterior distribution (Sec. IV.4.5).

####  IV.4.1 Training Size Dependence for Model Emulation Performance

Before we consider applying all design points for the training of the
emulators, we will first investigate how the emulation accuracy (Eq. (10)) and
the quality of uncertainty estimation (Eq. (11)) scale with the size of
training points. For this analysis, we will only show three examples and focus
on the general qualitative behavior of the emulators when the training set
varies in size.

Figure 1 shows for three different observables at two different energies the
emulator RMS error $\mathcal{E}$ (left) defined in Eq. (10) and the quality
measure of the uncertainty estimation $\mathcal{H}$ (right) defined in Eq.
(11).

Figure 1: (Color Online) The averaged RMS errors for emulators $\mathcal{E}$
(left) and the quality of uncertainty estimation measure $\mathcal{H}$ (right)
for three different example observables. The top row shows the first
observable from the 7.7 GeV $\langle p_{\rm T}\rangle$ dataset, and the middle
row is taken from the first observable in the 200 GeV ${\rm d}N/{\rm d}y$
dataset. The bottom row shows an observable at midrapidity from the 200 GeV
$v_{2}(\eta)$ dataset.

In the first row, we show the first observable from the 7.7 GeV $\langle
p_{\rm T}\rangle$ dataset; in the second one, the first observable from the
200 GeV ${\rm d}N/{\rm d}y$ dataset, and in the last one the $v_{2}(\eta)$ at
midrapidity from the 200 GeV dataset. We divide the 1,000 design points into
two groups: the first $n_{\rm train}$ points for training the GP emulators and
the rest for validation. Because we use the Maximum Projection Latin Hypercube
Design, the first $n_{\rm train}$ points can give a uniform coverage of the
20-dimensional model parameter space. We note that as $n_{\rm train}$
increases to close to the full dataset, the limited number of validation
points would introduce statistical fluctuations when computing the GP test
metrics $\mathcal{E}$ and $\mathcal{H}$. Therefore, we would expect good
scaling behavior for $n_{\rm train}$ up to 700 in our test.

Starting with the $\langle p_{\rm T}\rangle$ observable (top panels), we can
see that increasing the training data size decreases the emulation uncertainty
with a power-law scaling behavior, inline with previous results Nijs  _et al._
(2021b). The quantitative power law dependence varies from one observable to
the other, therefore we do not report any fit values here. The value for the
quality measure of the GP uncertainty estimation $\mathcal{H}$ decreases from
0.5 to $\approx 0$ for $n_{\rm train}\geq 300$. It shows that the GP emulators
can provide a good estimation of their prediction uncertainty. We can also see
that the PCGP and PCSK emulators give a smaller uncertainty overall and better
uncertainty estimation compared to the Scikit GP emulator.

The test results for ${\rm d}N/{\rm d}y$ (middle panels) show an even more
clear separation between the Scikit GP emulator and the other two emulators,
with almost a factor of two larger uncertainty for $n_{\rm train}\geq 700$.
The $\mathcal{H}$ metric is close to zero, but the PCGP and PCSK emulators
even give a more conservative estimate while having a smaller prediction
uncertainty.

The $v_{2}(\eta)$ observable (bottom panels) shows a similar performance of
all three emulators in the uncertainty estimate. Still, the PCGP and PCSK
emulators seem to be slightly better at estimating the prediction uncertainty
compared to the Scikit GP emulator.

####  IV.4.2 Benchmark Emulation Performance

Figure 2: (Color Online) The averaged RMS errors $\mathcal{E}$ for the three
different types of GP emulators. Different training sets are separated by
black lines. All emulators are trained with the same 970 LHD points. Figure 3:
(Color Online) The metric for emulator uncertainty estimation $\mathcal{H}$
for the three different types of GP emulators. Black lines separate different
training sets. All emulators are trained with the same 970 LHD points.

To benchmark and compare the performance of the three types of GP emulators
over all observables, we now fix the number of training points from the LHD to
970 in the following. This leaves us 30 validation points to compute the
$\mathcal{E}$ and $\mathcal{H}$ test metrics for the emulators. The absolute
values of $\mathcal{E}$ and $\mathcal{H}$ will still contain some statistical
fluctuations.

Figure 2 displays the emulators’ RMS error $\mathcal{E}$ for all 544
experimental observables – separated by black lines according to Table 1 –
trained with the three types of GP emulators. The upper panel summarizes the
mid-rapidity observables from central to 60% peripheral centrality for Au+Au
collisions at three energies. The lower panel shows the pseudo-rapidity
dependence observables in five centrality classes at 19.6 and 200 GeV. For the
${\rm d}N_{\rm ch}/{\rm d}\eta$ datasets, we can identify the different
“$w$”-shaped structures, which are the multiplicities as a function of
rapidity in different centrality classes. Each point in the plot represents a
single observable, indicated by the ticks on the $x$-axis.

We observe that the PCGP and PCSK emulators are compatible with each other for
most observables and that their overall prediction uncertainty is lower than
that of the Scikit GP. This trend was already visible in comparing the
emulators for single observables as a function of the number of training
points in Sec. IV.4.1. The particle yield $\mathrm{d}N/\mathrm{d}y$ dataset at
$7.7\;\mathrm{GeV}$ seems to be equally challenging for all three emulators,
with uncertainties up to 100% in peripheral centrality bins, while most other
points are in the region of $2-20\%$.

Next, we want to analyze how reliable the uncertainty provided by the
emulators is with the metric $\mathcal{H}$ defined in Eq. (11). The results
for the three types of GP emulators are shown in Fig. 3, using the same
conventions as in Fig. 2. When $\mathcal{H}$ is around zero, the estimated
uncertainty is close to the true uncertainty in our validation test.

For all three types of GP emulators, most observables show a small trend
toward positive values, which means that the emulator underestimates the
uncertainty compared to the true values in the validation set. Comparing the
GP emulators, we find that the PCGP and PCSK emulators perform better than the
Scikit GP emulator. Their $\mathcal{H}$ are closer to 0 than those from the
Scikit GP emulator.

Overall, we find the PCGP and PCSK emulators from the BAND surmise package
perform better than the Scikit GP emulator with standard settings. We note
that the relatively poorer performance of the Scikit GP emulator could be
related to not having optimized the GP kernel and its parameters.

####  IV.4.3 Logarithmic transformation to training datasets

Figure 4: (Color Online) The averaged RMS error $\mathcal{E}$ for the PCSK
emulator trained with (open symbols) and without (solid markers) logarithmic
transformation to particle multiplicities observables. Both emulators are
trained with the same 970 LHD points and validated with the remaining 30
points in the initial design. Black lines separate different sets of
observables. Figure 5: (Color Online) Same layout as in Fig. 4 but for the
metric for emulator uncertainty estimation $\mathcal{H}$.

Based on the results from previous sections, we now pick the PCSK emulator for
further investigation because it is one of the two more accurate ones and
incorporates the uncertainties of the training data. In Appendix A, we discuss
the emulation results of this emulator in more detail concerning the size of
the uncertainties of the training points and how it compares to the PCGP
emulator in the case no uncertainties are provided. The results discussed in
the following are general and hold for the other two emulators as well (not
shown).

In this section, we discuss the impacts of training with the logarithm of the
multiplicity observables on the accuracy of the GP emulations. Performing
logarithmic transformation on particle multiplicities was found to be an
effective trick for yielding more accurate GP emulators Nijs  _et al._
(2021b); Everett  _et al._ (2021a). It is suspected that
$\ln(\rm{d}N/\rm{d}y)$ reduces the relative variations in the training
dataset. However, the scaling transformation described in Sec. III should be a
more systematic way to adjust the relative variations in the training data.

Figure 4 shows in full squares the $\mathcal{E}$ points for PCSK trained
directly with particle multiplicities $\rm{d}N/\rm{d}y$ (from Fig. 2) and
compares them to the results from the logarithmic trained emulators shown in
open light blue points.

We find that the uncertainty is decreased for almost all five datasets or
remains comparable to the standard results. Only for very few points does the
uncertainty increase by a few percent. The largest decrease is visible for the
large rapidity data points, which are the rightmost ones in the linear rise of
the $\mathrm{d}N/\mathrm{d}y$ datasets and the outermost points in the
“w”-shaped parts of the $\mathrm{d}N_{\rm ch}/\mathrm{d}\eta$ datasets.

Concerning the metric $\mathcal{H}$ for the uncertainty estimate in Fig. 5, we
find that for most datasets the points move closer to zero or stay reasonably
close to zero.

Overall, we find that the logarithmic transformation to the training data can
be considered a valid way to decrease the uncertainty of quantities varying
over a large range.

####  IV.4.4 PCA Transformation to Functional Model Parameters

Figure 6: (Color Online) The same layout as in Fig. 2 but for comparing the
PCSK emulators with and without PCA transformation on the functional model
parameters. Figure 7: (Color Online) The same layout as in Fig. 6 but for the
metric for emulator uncertainty estimation $\mathcal{H}$.

In this section, we will assess the impacts on the emulator accuracy from
performing PCA transformations on the functional model inputs,
$\tilde{\zeta}(T,\mu_{B})$, $\tilde{\eta}(\mu_{B})$ and $\langle y_{\rm
loss}\rangle(y_{\mathrm{init}})$, and use the principal components to fit the
GP emulators instead of the curve parameters listed in Table 2.

To explain 99% of the variance in the training design, we need three principal
components for $\tilde{\eta}(\mu_{B})$ and $\langle y_{\rm
loss}\rangle(y_{\mathrm{init}})$ and six principal components for
$\tilde{\zeta}(T,\mu_{B})$. Since the functional parametrizations for
$\tilde{\eta}(\mu_{B})$ and $\langle y_{\rm loss}\rangle(y_{\mathrm{init}})$
are simple piece-wise linear functions (see Eqs. (7) and (8)), the number of
principal components needed is the same as the number of model parameters. For
the QGP’s specific bulk viscosity, the parameterization in Eq. (9) is non-
linear. The PCA transformation indicates that two more degrees of freedom are
needed in addition to the four model parameters.

Figure 6 shows the standard PCSK emulators in full points and the results (in
open markers) from the PCSK GP trained with the principal components from the
PCA transformations on the model parameters, $\tilde{\zeta}(T,\mu_{B})$,
$\tilde{\eta}(\mu_{B})$ and $\langle y_{\rm loss}\rangle(y_{\mathrm{init}})$.
We note that because of the PCA transformation on model parameters, our
training points may no longer be evenly distributed in the 22-dimensional
model parameter space.

Comparing the full and open points across the whole set of observables, we
find that the PCA transformations do not significantly alter the overall
emulation performance. There are a few observables where the uncertainty
increases by a few percent, but there are also observables with a decreased
uncertainty.

The same statement holds for the metric for emulator uncertainty estimation
shown in Fig. 7, where overall the $\mathcal{H}$ values are not significantly
affected by the PCA transformation.

Our results indicate that the relevant model parameters listed in Table 2 are
sufficient to capture all the essential information coded in Eqs (7), (8), and
(9), when training the GP emulators. The two additional degrees of freedom in
the PCA transformed model parameter space do not significantly influence the
GP performance.

However, it is noteworthy, that the PCA transformation we introduced here
could serve as a generic approach towards feeding general functional model
inputs to the GP emulators.

In the future, this approach can be used to allow for a parameter-less
emulation of the rapidity loss and the viscosities in the hydrodynamic
evolution.

####  IV.4.5 Active Learning Approach

Figure 8: (Color Online) The same layout as in Fig. 2 but for comparing the
PCSK emulators trained with and without the high probability posterior points.
Figure 9: (Color Online) The same layout as in Fig. 8 but for the metric for
emulator uncertainty estimation $\mathcal{H}$.

The dimension of our model parameters is high. It is more computationally
efficient to reduce the emulator uncertainty only in the region of model space
near where the real experimental data lie than across the entire parameter
space.

For this purpose, we explore the following two-step approach for our emulator
training. First, we train a model emulator with the given 1000 LHD points and
obtain the model’s posterior distribution by performing a Bayesian inference
analysis on the real experimental data Shen  _et al._ (2024). Then we run full
model simulations on 100 parameter sets sampled from the posterior
distribution. We denote this dataset as the High Probability Posterior (HPP)
points in the following. Finally, we train our GP emulators using the combined
LHD + HPP datasets. Because the HPP points populate a relatively dense region
near the physical posterior region, they would improve the emulator’s accuracy
within this phase space.

In the following tests, we train the GP emulators using 1000 LHD + 70 HPP
points and evaluate their performance in the model posterior region with the
rest 30 HPP points as the validation data. We compare such a setup with the
standard GP training with only the 1000 LHD points.

In this section, we only evaluate the emulator performance by computing the
$\mathcal{E}$ and $\mathcal{H}$ metrics on the leftover 30 HPP points. The
values of the $\mathcal{E}$ and $\mathcal{H}$ are different from the results
from previous sections, in which the validations were performed over the
entire model parameter space.

We denote our approach here as “active learning” although we only have one
iteration in our exercise. In principle, this technique can be applied in
multiple iterations to increase the accuracy in the interesting regions of the
high dimensional parameter space. This allows reducing the computational cost
in a Bayesian parameter extraction using a computationally expensive physics
model. A similar approach was employed in previous works Fan  _et al._ (2023);
Liyanage  _et al._ (2023).

Figure 8 shows the comparison between the fully trained LHD data set (full
symbols) and the active learning (open symbols).

First, we can say that the overall uncertainty of all quantities decreases
significantly when the full 1000 LHD points are taken into account and the
emulator is tested with the posterior points. This is probably due to a bias
in the testing in only a preferred region of the parameter space, where the
posterior is peaked. Nevertheless, we can see that the inclusion of the
posterior points improves almost all observables. The ones that do not change
have similar uncertainty compared to the standard setup, while others change
by up to factors five to ten.

Figure 9 shows that using only the posterior points to compute $\mathcal{H}$
leads to more conservative uncertainty estimates of the emulator
($\mathcal{H}<0$), which is of similar order for most of the observables and
both methods.

##  V Closure Tests

In this section, we will explore how the different levels of emulator
performance can influence the model parameter posterior distribution in the
Bayesian inference analysis.

###  V.1 MCMC Setups

Our analysis in this paper utilizes two different implementations of posterior
sampling techniques. The first employs affine invariant MCMC, leveraging the
emcee Python package Foreman-Mackey  _et al._ (2013). Secondly, we employ
Parallel Tempering Langevin Monte Carlo (PTLMC), implemented in the surmise
package Plumlee  _et al._ (2023). The parallel tempering part enhances the
exploration of the target distribution by running multiple chains at different
“temperatures” and exchanging states periodically, which scans the parameter
space more accurately and reduces the probability of getting stuck in a local
minimum. This technique is of advantage if the posterior distribution is a
multi-modal distribution. Additionally, some Langevin dynamics are
incorporated to explore complex posterior distributions more efficiently
Roberts and Stramer (2002).

To test the different scenarios, including the different emulators and methods
discussed in Sec. IV.4, we are defining some Runs in Tab. 3. The first three
setups in our analysis compare the impact of different types of GP emulators
on posterior distributions in a standard configuration. This configuration
excludes logarithmic transformation, but incorporates PCA for the functional
model inputs, and trains the emulators using both LHD and HPP points. By
examining these Runs, we aim to assess how various emulator choices influence
the accuracy and reliability of the posterior inference process.

Run 4 mirrors the Run 1’s setup but utilizing PTLMC to obtain the posterior.
This direct comparison between MCMC and the parallel tempering variant allows
us to evaluate the efficacy of temperature-enhanced exploration techniques in
improving posterior sampling efficiency and exploring high-dimensional
parameter spaces more effectively. By employing the advanced sampling method,
we seek to uncover insights into how different sampling strategies impact the
exploration of posterior distributions and the convergence properties of
Bayesian inference algorithms.

By contrasting Runs 1 with 5, we can elucidate the influence of incorporating
an additional 95 HPP points from the posterior distribution in the training
phase on posterior inference. Similarly, comparing Runs 1 with 6 and 7 enables
us to quantify the impact of PCA and logarithmic transformation in the
training design on the posterior distributions.

These comparative analyses offer valuable insights into the various factors
that contribute to the accuracy and reliability of the GP emulators that
affect Bayesian posterior sampling techniques.

###  V.2 Closure Test

In Section IV.4, we demonstrated that the emulators can predict the training
data with a quantifiable level of uncertainty, as discussed previously.
However, the ability of the given data sets to effectively constrain model
parameters may vary. If the data exhibit insensitivity to certain parameters,
their constraints may be poor, leading to broad posterior distributions.
Additionally, if the underlying physics model contains degeneracies, where
different parameter sets can equally describe the experimental data, posterior
distributions may even become multimodal. To identify such scenarios,
conducting a closure test is needed. In a closure test, one training point is
excluded from the training set and treated as the pseudo-experimental data. By
carrying out Bayesian inference analysis on the pseudo-data, we can compare
the obtained posterior distributions with the true model parameters used to
generate the pseudo-data. Ideally, the posterior distributions should exhibit
peaks close to the true value. In cases of model degeneracy, the peaks should
at least be proximate. By doing closure tests, we can quantify how well our
method works in understanding the model and the limitations of our particular
Bayesian inverse problem.

Figure 10 shows the comparison of the three different types of GP emulators
from Runs 1, 2, and 3 listed in Tab. 3.

Figure 10: (Color Online) Closure test posterior distributions for the 20
different parameters from Tab. 2 using the three different GP emulators with
runs 1, 2, and 3 in Tab. 3. The marginalized parameter distributions are
plotted in logarithmic scale in the $y$-axis. The true value is indicated by
the solid black line. The values above each plot in the respective colors of
the lines are the median values with the 90% confidence level intervals.

We find that the PCSK and PCGP emulators show clear and narrow peaks around
the true values indicated by the black vertical line and relatively a lot
broader posterior distributions for the results using the Scikit GP emulator.
This result shows that the quality of the GP emulator has a strong influence
on the posterior distribution extracted in the Bayesian inference. The PCSK
and PCGP emulators show compatible results in the closure test.

To further quantify the closure test and compare results from different
setups, introduced in Tab. 3 we will compare a metric $\Delta$ from the
information-theoretic approach similar to Ref. Fan  _et al._ (2023). The main
contribution to the GP emulator information loss, which is mainly determined
by the second moment of the posterior distribution, using

$\Delta\equiv\frac{1}{N_{\rm
param.}}\int\left|\frac{\bm{\theta}-\bm{\theta}_{\rm truth}}{\bm{\theta}_{\rm
max}-\bm{\theta}_{\rm min}}\right|^{2}p(\bm{\theta})\;\mathrm{d}\bm{\theta},$
(12)

where $\bm{\theta}_{\rm truth}$ is the closure test point, $|\bm{\theta}_{\rm
max}-\bm{\theta}_{\rm min}|$ is the allowed prior range and $p(\bm{\theta})$
is the posterior distribution, where the closure test point is not included.
Different from the definition in Ref. Fan  _et al._ (2023), we are dividing by
$N_{\rm param.}$, the number of model parameters, to normalize the values of
$\Delta$ between zero and one. The smaller the value, the more peaked is the
posterior distribution around the true value. The last column of Tab. 3
summarizes the values of $\Delta$ for the different MCMC runs.

Table 3: Closure test runs with GP emulators trained in different settings. The metric $\Delta$ in the last column quantifies the quality of the obtained posterior distributions to the true model values in the closure test. Run | Training Design | MCMC Type | Log Transformation | PCA Transformation | Emulator | $\Delta$
---|---|---|---|---|---|---
1 | LHD + HPP | MCMC | $\times$ | $\checkmark$ | PCSK | $4.5\times 10^{-4}$
2 | LHD + HPP | MCMC | $\times$ | $\checkmark$ | PCGP | $5.7\times 10^{-4}$
3 | LHD + HPP | MCMC | $\times$ | $\checkmark$ | Scikit GP | $2.5\times 10^{-2}$
4 | LHD + HPP | PTLMC | $\times$ | $\checkmark$ | PCSK | $4.4\times 10^{-4}$
5 | LHD | MCMC | $\times$ | $\checkmark$ | PCSK | $4.5\times 10^{-2}$
6 | LHD + HPP | MCMC | $\times$ | $\times$ | PCSK | $4.3\times 10^{-4}$
7 | LHD + HPP | MCMC | $\checkmark$ | $\checkmark$ | PCSK | $6.2\times 10^{-4}$
Table 4: The 90% confidence level intervals for the 20 model parameters extracted from the posterior distributions of runs one to four (see Tab. 3). The rightmost column shows the “true” value from the closure test. Model parameter | Run 1 | Run 2 | Run 3 | Run 4 | True Value
---|---|---|---|---|---
$B_{G}\;[\mathrm{GeV}^{-2}]$ | $19.956_{-0.855}^{+0.878}$ | $19.861_{-1.037}^{+0.910}$ | $20.028_{-7.182}^{+4.304}$ | $19.962_{-0.606}^{+0.794}$ | $19.928$
$\alpha_{\rm shadowing}$ | $0.027_{-0.010}^{+0.011}$ | $0.032_{-0.012}^{+0.012}$ | $0.043_{-0.039}^{+0.067}$ | $0.028_{-0.010}^{+0.011}$ | $0.021$
$y_{{\rm loss},2}$ | $1.514_{-0.019}^{+0.017}$ | $1.511_{-0.022}^{+0.021}$ | $1.464_{-0.114}^{+0.115}$ | $1.514_{-0.019}^{+0.016}$ | $1.511$
$y_{{\rm loss},4}$ | $1.702_{-0.026}^{+0.028}$ | $1.714_{-0.033}^{+0.033}$ | $1.825_{-0.168}^{+0.190}$ | $1.703_{-0.025}^{+0.028}$ | $1.690$
$y_{{\rm loss},6}$ | $2.366_{-0.037}^{+0.038}$ | $2.361_{-0.045}^{+0.042}$ | $2.216_{-0.234}^{+0.250}$ | $2.368_{-0.035}^{+0.030}$ | $2.338$
$\sigma_{y_{\rm loss}}$ | $0.298_{-0.034}^{+0.021}$ | $0.295_{-0.027}^{+0.024}$ | $0.317_{-0.168}^{+0.182}$ | $0.303_{-0.019}^{+0.018}$ | $0.317$
$\alpha_{\rm rem}$ | $0.569_{-0.013}^{+0.013}$ | $0.573_{-0.013}^{+0.013}$ | $0.627_{-0.088}^{+0.092}$ | $0.570_{-0.011}^{+0.013}$ | $0.568$
$\lambda_{B}$ | $0.285_{-0.012}^{+0.012}$ | $0.285_{-0.014}^{+0.014}$ | $0.281_{-0.054}^{+0.056}$ | $0.286_{-0.011}^{+0.011}$ | $0.284$
$\sigma_{x}^{\rm string}\;[{\rm fm}]$ | $0.181_{-0.005}^{+0.004}$ | $0.179_{-0.006}^{+0.006}$ | $0.130_{-0.025}^{+0.037}$ | $0.181_{-0.005}^{+0.004}$ | $0.182$
$\sigma_{\eta}^{\rm string}$ | $0.120_{-0.005}^{+0.005}$ | $0.120_{-0.006}^{+0.006}$ | $0.148_{-0.031}^{+0.038}$ | $0.120_{-0.005}^{+0.005}$ | $0.122$
$\alpha_{\text{string tilt}}$ | $0.977_{-0.039}^{+0.020}$ | $0.962_{-0.053}^{+0.033}$ | $0.614_{-0.535}^{+0.354}$ | $0.982_{-0.035}^{+0.016}$ | $0.952$
$\alpha_{\text{preFlow}}$ | $0.081_{-0.011}^{+0.012}$ | $0.082_{-0.013}^{+0.017}$ | $0.155_{-0.100}^{+0.139}$ | $0.081_{-0.010}^{+0.013}$ | $0.086$
$\eta_{0}$ | $0.101_{-0.007}^{+0.006}$ | $0.100_{-0.010}^{+0.010}$ | $0.086_{-0.035}^{+0.040}$ | $0.101_{-0.006}^{+0.006}$ | $0.098$
$\eta_{2}$ | $0.073_{-0.008}^{+0.009}$ | $0.077_{-0.012}^{+0.013}$ | $0.130_{-0.099}^{+0.118}$ | $0.073_{-0.007}^{+0.009}$ | $0.073$
$\eta_{4}$ | $0.209_{-0.007}^{+0.008}$ | $0.211_{-0.010}^{+0.010}$ | $0.206_{-0.098}^{+0.079}$ | $0.209_{-0.007}^{+0.007}$ | $0.204$
$\zeta_{\rm max}$ | $0.190_{-0.005}^{+0.004}$ | $0.189_{-0.006}^{+0.006}$ | $0.182_{-0.030}^{+0.016}$ | $0.189_{-0.005}^{+0.005}$ | $0.188$
$T_{\zeta,0}\;[{\rm GeV}]$ | $0.228_{-0.003}^{+0.003}$ | $0.228_{-0.003}^{+0.003}$ | $0.226_{-0.029}^{+0.019}$ | $0.228_{-0.003}^{+0.003}$ | $0.228$
$\sigma_{\zeta,+}\;[{\rm GeV}]$ | $0.077_{-0.005}^{+0.001}$ | $0.075_{-0.003}^{+0.003}$ | $0.075_{-0.024}^{+0.023}$ | $0.077_{-0.005}^{+0.001}$ | $0.082$
$\sigma_{\zeta,-}\;[{\rm GeV}]$ | $0.041_{-0.002}^{+0.003}$ | $0.041_{-0.002}^{+0.003}$ | $0.042_{-0.021}^{+0.013}$ | $0.041_{-0.002}^{+0.003}$ | $0.041$
$e_{\rm sw}\;[{\rm GeV}/{\rm fm}^{3}]$ | $0.487_{-0.006}^{+0.005}$ | $0.486_{-0.007}^{+0.006}$ | $0.462_{-0.042}^{+0.033}$ | $0.487_{-0.006}^{+0.005}$ | $0.499$

Comparing the first three runs in Tab. 3 with the posterior distributions in
Fig.10 gives us a good impression of the interpretation of the $\Delta$
quantity. The values for the PCGP and PCSK emulators are two orders of
magnitude smaller than the Scikit GP value, which means that they are more
peaked about the true value. We also verify that the two different MCMC
approaches give consistent posterior distributions and compatible values of
$\Delta$ (run 1 vs. 4).

When the HPP points are not included in the emulator training, we can see that
the value of $\Delta$ increases by two orders of magnitude between run 1 and
run 5. This result is in line with the findings we had for the emulator
uncertainties $\mathcal{E}$ discussed in Sec. IV.4.5.

Including the PCA and logarithmic transformations in the emulators’ training
does not significantly change the quality of reproducing the true parameter
set in the closure test.

Table 4 summarizes the median values of the model parameters with 90%
confidence intervals from the obtained posterior distributions for runs 1 to
4. Because the posterior distributions from runs 1, 2, and 4 are sharply
peaked around the true parameter set, the 90% confidence intervals are
significantly narrower than those from run 3.

##  VI Observable Sensitivity Analysis

Effectively constraining model parameters through Bayesian analysis requires
understanding which parameters have a strong influence on specific
observables. In our analysis, we evaluate how parameters impact observables by
calculating their first-order partial derivatives at specific points in the
parameter space. We define the response of observable $i$ w.r.t parameter $j$
by

$\mathcal{R}_{ij}\equiv\frac{\partial\ln(\bm{\mathcal{G}}_{i}(\bm{\theta}))}{\partial\ln(\bm{\theta}_{j})}.$
(13)

The elements of this response matrix are unitless and show how strongly the
observables depend on the parameters: positive values mean an increase and
negative values mean a decrease. To calculate the above derivative numerically
we use the centered, finite difference method for first-order derivatives from
the emulator predictions.

From Fig. 11, we can identify some clear trends in how some observables are
sensitive to local changes in parameter space. As an example, we use the
closure test point for the sensitivity analysis.

Figure 11: Response matrix $\mathcal{R}_{ij}$ at the point used in the closure
test. The observables in the left figure correspond to 200 GeV collisions. The
right-figure observables come from 19.6 GeV and 7.7 GeV simulations.

Let us now look at a few features in the response matrix of the observables
w.r.t. the 20 parameters used in the analysis. The rapidity distribution and
mean transverse momentum distribution of $\mathrm{d}N/\mathrm{d}y$ of $\pi$,
$K$ and $p$ at 200 GeV show a strong positive response to the rapidity loss
parameters $y_{\mathrm{loss},4}$ and $y_{\mathrm{loss},6}$, which is natural
due to the high rapidity of the incoming particles. On the other hand, the
same observables for 7.7 GeV collisions depend strongly on the parameter
$y_{\mathrm{loss},2}$, since the beams have a much smaller rapidity. We also
see a strong dependence for all observables at least on some of the bulk
viscosity parameters. When the shear viscosity parameter $\eta_{0}$ increases,
the anisotropic flow coefficients decrease. A decrease in switching energy
between hydrodynamics and the afterburner means a higher freeze-out
temperature, which produces more particles with smaller transverse momentum.
An increase in the maximum bulk viscosity parameter leads to a decrease in the
average transverse momentum and flow observables (more bulk viscosity reduces
anisotropic flow, all the bulk viscosity parameters and the corresponding
increase or decrease in the flow coefficients confirm this trend).

For two of the parameters, $\alpha_{\rm shadowing}$ and $\alpha_{\text{string
tilt}}$, there is almost no sensitivity shown for the currently considered
observables. All the other parameters have strong sensitivities to at least
one type of observable at a given energy.

In follow-up work, we will discuss the physics interpretation of the response
matrix in more detail and perform the same parameter sensitivity analysis of
the maximum a posteriori probability (MAP) point with the real experimental
measurements at the RHIC Beam Energy Scan program.

##  VII Conclusion

In this work, we have performed a systematic comparison between three
different types of open-source GP emulators – PCGP, PCSK, and the Scikit GP –
while training them to a (3+1)D heavy-ion collision model at three collision
energies. We define a 20-dimension model parameter space and solve the
Bayesian inverse problem with 544 observables related to the RHIC Beam Energy
Scan program. In addition to the difference in the emulators, we also
investigated the differences in the final posterior distributions obtained
from a Bayesian analysis using two different sampling methods – MCMC and PTLMC
– in our analysis.

In the first part of this paper, we have seen that the PCGP and PCSK emulators
generally give more accurate predictions and better uncertainty estimation
than the Scikit GP. We compared the standard training setup with the
logarithmic transformation on the multiplicity datasets, finding, that the
uncertainties are either the same or decreased in the latter case. The
difference between the approaches becomes smaller as the size of the training
dataset increases. When comparing the standard training setup with the PCA
transformation on functional model inputs, $\tilde{\zeta}(T,\mu_{B})$,
$\tilde{\eta}(\mu_{B})$ and $\langle y_{\rm loss}\rangle(y_{\mathrm{init}})$,
we found that there is overall no significant effect in the emulation
uncertainties or the metric for the uncertainty estimation of the emulators.
Even if the emulation performance does not significantly improve, this method
can incorporate generic functional model inputs. The last effect we have
investigated in the emulator validation is the inclusion of the high
probability posterior points, which is a simple version of an “active
learning” ansatz with one iteration. The inclusion of these additional points
sampled from the posterior distribution significantly decreases the emulation
uncertainty in the phase space where the true experimental measurements
remain.

We then propagate the different emulator setups in the Bayesian inference
closure test to quantify their impacts on the posterior distributions. We
observed that the PCGP and PCSK emulators provide well-defined peaked
posterior distributions, while the Scikit GP result is less constrained. The
different posterior distributions are compared using the information loss
quantity $\Delta$. We obtained compatible posterior distributions from the
standard MCMC and PTLMC algorithms. Both posterior distributions are peaked
around the true model parameter set. This result verified that the 544
observables used in the closure test can effectively constrain the
20-dimension model parameter space and result in the correct solution in the
Bayesian inverse problem. We find that including the high probability
posterior points in the training set significantly improves the closure test.
Therefore, we recommend an active learning setup for future model emulation
training. This approach will also use computational resources more effectively
and explore more details in the interesting regions in parameter space.
Lastly, we found that the logarithmic and PCA transformations do not
significantly influence the accuracy of the extracted posterior distribution.

With the trained GP emulators, we can study the response matrix between the
model parameters and experimental observables. Such sensitivity analysis can
help us identify the most important set of observables to constrain the QGP’s
physical properties.

The code developed for this work is open source and can be downloaded from
Roch and Jahan (2024).

###### Acknowledgements.

We thank Moses Chan for his helpful comments on the surmise Python package and
the manuscript of this paper. 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 (OAC-2004571). C.S. acknowledges a DOE Office of
Science Early Career Award. Numerical simulations presented in this work were
performed at the Wayne State Grid, and we gratefully acknowledge their
support.

##  Appendix A PCSK Different Training Uncertainty

The PCSK emulator is the only emulator in our analysis, which incorporates the
statistical uncertainty at the training points. In this case, the natural
question arises, what happens if the training data is, for example, more or
less precise? We study this effect by multiplying the uncertainties of the
training data by a factor of zero (label $0\sigma$) or two (label $2\sigma$)
and comparing the trained emulator to the one with the original uncertainties.

Figure 12: (Color Online) The same layout as in Fig. 2 but for comparing the
PCSK GP emulators with different sizes of uncertainties in the training data.
Figure 13: (Color Online) The same layout as in Fig. 12 but for the metric for
emulator uncertainty estimation $\mathcal{H}$.

Figures 12 and 13 show that the RMS emulator error $\mathcal{E}$ and the
metric $\mathcal{H}$ barely changes for all the observables in our analysis.
We also tried to change the uncertainty by using only half the number of
events for the computation of the observables per design point (not shown).
This effect is of similar size to the one shown in Fig. 12.

Additionally, we compared the PCSK emulator with the uncertainty set to zero
to the PCGP emulator, which gives similar results for the two emulators (not
shown). Small differences of the order of a few percent can be observed for
$\mathcal{E}$ in the $19.6\;\mathrm{GeV}$ dataset for the multiplicities in
the larger centrality classes.

## References

  * Busza  _et al._ (2018) Wit Busza, Krishna Rajagopal, and Wilke van der Schee, “Heavy Ion Collisions: The Big Picture, and the Big Questions,” Ann. Rev. Nucl. Part. Sci. 68, 339–376 (2018), arXiv:1802.04801 [hep-ph] .
  * Miller  _et al._ (2007) Michael L. Miller, Klaus Reygers, Stephen J. Sanders, and Peter Steinberg, “Glauber modeling in high energy nuclear collisions,” Ann. Rev. Nucl. Part. Sci. 57, 205–243 (2007), arXiv:nucl-ex/0701025 .
  * Loizides  _et al._ (2015) C. Loizides, J. Nagle, and P. Steinberg, “Improved version of the PHOBOS Glauber Monte Carlo,” SoftwareX 1-2, 13–18 (2015), arXiv:1408.2549 [nucl-ex] .
  * Broniowski  _et al._ (2009) Wojciech Broniowski, Maciej Rybczynski, and Piotr Bozek, “GLISSANDO: Glauber initial-state simulation and more..” Comput. Phys. Commun. 180, 69–83 (2009), arXiv:0710.5731 [nucl-th] .
  * Drescher and Nara (2007) H. J. Drescher and Y. Nara, “Effects of fluctuations on the initial eccentricity from the Color Glass Condensate in heavy ion collisions,” Phys. Rev. C 75, 034905 (2007), arXiv:nucl-th/0611017 .
  * Schenke  _et al._ (2012) Bjoern Schenke, Prithwish Tribedy, and Raju Venugopalan, “Fluctuating Glasma initial conditions and flow in heavy ion collisions,” Phys. Rev. Lett. 108, 252301 (2012), arXiv:1202.6646 [nucl-th] .
  * Paatelainen  _et al._ (2013) R. Paatelainen, K. J. Eskola, H. Holopainen, and K. Tuominen, “Multiplicities and $p_{T}$ spectra in ultrarelativistic heavy ion collisions from a next-to-leading order improved perturbative QCD + saturation + hydrodynamics model,” Phys. Rev. C 87, 044904 (2013), arXiv:1211.0461 [hep-ph] .
  * Moreland  _et al._ (2015) J. Scott Moreland, Jonah E. Bernhard, and Steffen A. Bass, “Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions,” Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Giacalone  _et al._ (2019) Giuliano Giacalone, Pablo Guerrero-Rodríguez, Matthew Luzum, Cyrille Marquet, and Jean-Yves Ollitrault, “Fluctuations in heavy-ion collisions generated by QCD interactions in the color glass condensate effective theory,” Phys. Rev. C 100, 024905 (2019), arXiv:1902.07168 [nucl-th] .
  * Shen and Schenke (2018) Chun Shen and Björn Schenke, “Dynamical initial state model for relativistic heavy-ion collisions,” Phys. Rev. C 97, 024907 (2018), arXiv:1710.00881 [nucl-th] .
  * Bass  _et al._ (1998) S. A. Bass _et al._ , “Microscopic models for ultrarelativistic heavy ion collisions,” Prog. Part. Nucl. Phys. 41, 255–369 (1998), arXiv:nucl-th/9803035 .
  * Weil  _et al._ (2016) J. Weil _et al._ (SMASH), “Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions,” Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .
  * Cassing and Bratkovskaya (1999) W. Cassing and E. L. Bratkovskaya, “Hadronic and electromagnetic probes of hot and dense nuclear matter,” Phys. Rept. 308, 65–233 (1999).
  * Lin  _et al._ (2005) Zi-Wei Lin, Che Ming Ko, Bao-An Li, Bin Zhang, and Subrata Pal, “A Multi-phase transport model for relativistic heavy ion collisions,” Phys. Rev. C 72, 064901 (2005), arXiv:nucl-th/0411110 .
  * Schlichting and Teaney (2019) Soeren Schlichting and Derek Teaney, “The First fm/c of Heavy-Ion Collisions,” Ann. Rev. Nucl. Part. Sci. 69, 447–476 (2019), arXiv:1908.02113 [nucl-th] .
  * Teaney (2010) Derek A. Teaney, “Viscous Hydrodynamics and the Quark Gluon Plasma,” in  _Quark-gluon plasma 4_, edited by Rudolph C. Hwa and Xin-Nian Wang (2010) pp. 207–266, arXiv:0905.2433 [nucl-th] .
  * Romatschke and Romatschke (2019) Paul Romatschke and Ulrike Romatschke,  _Relativistic Fluid Dynamics In and Out of Equilibrium_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2019) arXiv:1712.05815 [nucl-th] .
  * Cooper and Frye (1974) Fred Cooper and Graham Frye, “Comment on the Single Particle Distribution in the Hydrodynamic and Statistical Thermodynamic Models of Multiparticle Production,” Phys. Rev. D 10, 186 (1974).
  * Sivia and Skilling (2006) Devinderjit Sivia and John Skilling,  _Data analysis: a Bayesian tutorial_ (OUP Oxford, 2006).
  * Shen  _et al._ (2016) Chun Shen, Zhi Qiu, Huichao Song, Jonah Bernhard, Steffen Bass, and Ulrich Heinz, “The iEBE-VISHNU code package for relativistic heavy-ion collisions,” Comput. Phys. Commun. 199, 61–85 (2016), arXiv:1409.8164 [nucl-th] .
  * Bernhard  _et al._ (2016) Jonah E. Bernhard, J. Scott Moreland, Steffen A. Bass, Jia Liu, and Ulrich 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 [nucl-th] .
  * Putschke  _et al._ (2019) J. H. Putschke _et al._ , “The JETSCAPE framework,” (2019), arXiv:1903.07706 [nucl-th] .
  * Nijs  _et al._ (2021a) Govert Nijs, Wilke van der Schee, Umut Gürsoy, and Raimond Snellings, “Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum,” Phys. Rev. C 103, 054909 (2021a), arXiv:2010.15134 [nucl-th] .
  * Shen and Schenke (2022) Chun Shen and Björn Schenke, “Longitudinal dynamics and particle production in relativistic nuclear collisions,” Phys. Rev. C 105, 064905 (2022), arXiv:2203.04685 [nucl-th] .
  * Kennedy and O’Hagan (2001) M. C. Kennedy and A O’Hagan, “Bayesian calibration of computer models,” Journal of the Royal Statistical Society. Series B 63, 425–464 (2001).
  * Higdon  _et al._ (2004) Dave Higdon, Marc Kennedy, James C. Cavendish, John A. Cafeo, and Robert D. Ryne, “Combining field data and computer simulations for calibration and prediction,” SIAM Journal on Scientific Computing 26, 448–466 (2004).
  * Brooks  _et al._ (2011) S. Brooks, A. Gelman, G. Jones, and X.L. Meng,  _Handbook of Markov Chain Monte Carlo_, Chapman & Hall/CRC Handbooks of Modern Statistical Methods (CRC Press, 2011).
  * Trotta (2008) Roberto Trotta, “Bayes in the sky: Bayesian inference and model selection in cosmology,” Contemporary Physics 49, 71–104 (2008).
  * Gramacy (2020) Robert B. Gramacy, “Surrogates: Gaussian process modeling, design and optimization for the applied sciences,” (Chapman Hall/CRC, Boca Raton, Florida, 2020) http://bobby.gramacy.com/surrogates/.
  * Rasmussen (2004) Carl Edward Rasmussen, “Gaussian processes in machine learning,” in  _Advanced Lectures on Machine Learning_, edited by Olivier Bousquet, Ulrike von Luxburg, and Gunnar Rätsch (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 63–71.
  * Plumlee  _et al._ (2023) Matthew Plumlee, Özge Sürer, Stefan M. Wild, and Moses Y-H. Chan,  _surmise 0.2.1 Users Manual_, Tech. Rep. Version 0.2.1 (NAISE, 2023).
  * Bernhard  _et al._ (2019) Jonah E. Bernhard, J. Scott Moreland, and Steffen A. Bass, “Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma,” Nature Phys. 15, 1113–1117 (2019).
  * Everett  _et al._ (2021a) D. Everett _et al._ (JETSCAPE), “Multisystem Bayesian constraints on the transport coefficients of QCD matter,” Phys. Rev. C 103, 054904 (2021a), arXiv:2011.01430 [hep-ph] .
  * Mäntysaari  _et al._ (2022) Heikki Mäntysaari, Björn Schenke, Chun Shen, and Wenbin Zhao, “Bayesian inference of the fluctuating proton shape,” Phys. Lett. B 833, 137348 (2022), arXiv:2202.01998 [hep-ph] .
  * Shen  _et al._ (2024) Chun Shen, Björn Schenke, and Wenbin 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] .
  * Shen  _et al._ (2023) Chun Shen, Björn Schenke, and Wenbin Zhao, “The effects of pseudorapidity-dependent observables on (3+1)D Bayesian Inference of relativistic heavy-ion collisions,” in  _30th International Conference on Ultrarelativstic Nucleus-Nucleus Collisions_ (2023) arXiv:2312.09325 [nucl-th] .
  * Higdon  _et al._ (2008) Dave Higdon, James Gattiker, Brian Williams, and Maria Rightley, “Computer model calibration using high-dimensional output,” Journal of the American Statistical Association 103, 570–583 (2008).
  * Ramsay and Silverman (1997) J. Ramsay and B. Silverman,  _Functional Data Analysis_ (Springer, 1997).
  * Higdon  _et al._ (2015) Dave Higdon, Jordan D. McDonnell, Nicolas Schunck, Jason Sarich, and Stefan M Wild, “A Bayesian Approach for Parameter Estimation and Prediction using a Computationally Intensive Model,” J. Phys. G 42, 034009 (2015), arXiv:1407.3017 [physics.data-an] .
  * Sürer  _et al._ (2022) Özge Sürer, Filomena M. Nunes, Matthew Plumlee, and Stefan M. Wild, “Uncertainty quantification in breakup reactions,” Phys. Rev. C 106, 024607 (2022), arXiv:2205.07119 [nucl-th] .
  * Liyanage  _et al._ (2023) Dananjaya Liyanage, Özge Sürer, Matthew Plumlee, Stefan M. Wild, and Ulrich Heinz, “Bayesian calibration of viscous anisotropic hydrodynamic simulations of heavy-ion collisions,” Phys. Rev. C 108, 054905 (2023), arXiv:2302.14184 [nucl-th] .
  * Moses Y.-H. Chan and Wild (2024) Matthew Plumlee Moses Y.-H. Chan and Stefan M. Wild, “Constructing a simulation surrogate with partially observed output,” Technometrics 66, 1–13 (2024), https://doi.org/10.1080/00401706.2023.2210170 .
  * Petersen  _et al._ (2011) Hannah Petersen, Christopher Coleman-Smith, Steffen A. Bass, and Robert Wolpert, “Constraining the initial state granularity with bulk observables in Au+Au collisions at $\sqrt{s_{\rm NN}}=200$ GeV,” J. Phys. G 38, 045102 (2011), arXiv:1012.4629 [nucl-th] .
  * Novak  _et al._ (2014) John Novak, Kevin Novak, Scott Pratt, Joshua Vredevoogd, Chris Coleman-Smith, and Robert Wolpert, “Determining Fundamental Properties of Matter Created in Ultrarelativistic Heavy-Ion Collisions,” Phys. Rev. C 89, 034917 (2014), arXiv:1303.5769 [nucl-th] .
  * Sangaline and Pratt (2016) Evan Sangaline and Scott 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] .
  * Bernhard  _et al._ (2015) Jonah E. Bernhard, Peter W. Marcy, Christopher E. Coleman-Smith, Snehalata Huzurbazar, Robert L. Wolpert, and Steffen A. Bass, “Quantifying properties of hot and dense QCD matter through systematic model-to-data comparison,” Phys. Rev. C 91, 054910 (2015), arXiv:1502.00339 [nucl-th] .
  * Moreland  _et al._ (2020) J. Scott Moreland, Jonah E. Bernhard, and Steffen A. Bass, “Bayesian calibration of a hybrid nuclear collision model using p-Pb and Pb-Pb data at energies available at the CERN Large Hadron Collider,” Phys. Rev. C 101, 024911 (2020), arXiv:1808.02106 [nucl-th] .
  * Everett  _et al._ (2021b) D. Everett _et al._ (JETSCAPE), “Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging,” Phys. Rev. Lett. 126, 242301 (2021b), arXiv:2010.03928 [hep-ph] .
  * Nijs  _et al._ (2021b) Govert Nijs, Wilke van der Schee, Umut Gürsoy, and Raimond Snellings, “Transverse Momentum Differential Global Analysis of Heavy-Ion Collisions,” Phys. Rev. Lett. 126, 202301 (2021b), arXiv:2010.15130 [nucl-th] .
  * Heffernan  _et al._ (2023a) Matthew R. Heffernan, Charles Gale, Sangyong Jeon, and Jean-François Paquet, “Bayesian quantification of strongly-interacting matter with color glass condensate initial conditions,” (2023a), arXiv:2302.09478 [nucl-th] .
  * Heffernan  _et al._ (2023b) Matthew R. Heffernan, Charles Gale, Sangyong Jeon, and Jean-François Paquet, “Early-times Yang-Mills dynamics and the characterization of strongly interacting matter with statistical learning,” (2023b), arXiv:2306.09619 [nucl-th] .
  * Ankenman  _et al._ (2010) Bruce Ankenman, Barry L. Nelson, and Jeremy Staum, “Stochastic kriging for simulation metamodeling,” Operations Research 58, 371–382 (2010).
  * Roshan Joseph  _et al._ (2020) V. Roshan Joseph, Evren Gul, and Shan Ba, “Designing computer experiments with multiple types of factors: The maxpro approach,” Journal of Quality Technology 52, 343–354 (2020), https://doi.org/10.1080/00224065.2019.1611351 .
  * Roshan Joseph  _et al._ (2015) V. Roshan Joseph, Evren Gul, and Shan Ba, “Maximum projection designs for computer experiments,” Biometrika 102, 371–380 (2015), https://academic.oup.com/biomet/article-pdf/102/2/371/9642501/asv002.pdf .
  * Fan  _et al._ (2023) W. Fan _et al._ (JETSCAPE), “A new metric improving Bayesian calibration of a multistage approach studying hadron and inclusive jet suppression,” (2023), arXiv:2307.09641 [hep-ph] .
  * Foreman-Mackey  _et al._ (2013) Daniel Foreman-Mackey, David W. Hogg, Dustin Lang, and Jonathan Goodman, “emcee: The mcmc hammer,” Publications of the Astronomical Society of the Pacific 125, 306–312 (2013).
  * Roberts and Stramer (2002) G. O. Roberts and O. Stramer, “Langevin diffusions and metropolis-hastings algorithms,” Methodol. Comput. Appl. Probab. 4, 337–357 (2002).
  * Roch and Jahan (2024) Hendrik Roch and Syed Afrid Jahan, “Hendrik1704/gpbayestools-hic: v1.0.0,”  (2024).