we can therefore select only the first $m^{'}$ elements in the principal component space, and use the new training data set $\mathcal{D^{'}} = \{(\mathbf{x_i}, \mathbf{z_i})\}$ to train $m^{'}$ Gaussian Process emulators then make predictions of $\mathbf{z^{\*}}$ at each arbitrary parameter $\mathbf{x}$. The dimension of the output space is reduced by $m-m^{'}$.
Once the output $\mathbf{z^{\*}}$ of the principal components are predicted by the GP emulators, one would perform a inverse transformation to predict the original output $\mathbf{y^{\*}}$.

\begin{equation}
\mathbf{y^{\*}} = \frac{1}{\sqrt{n}} \mathbf{z^{\*}} V.
\end{equation}

§ LATIN HYPERCUBE SAMPLING

As discussed previously, the Gaussian process emulator serves as a fast surrogate model of the actual model. The idea is to use a few training data points while still being able to cover the parameter space with controlled uncertainty. The training data points play a vital role here. However, if a fixed number of values are sampled in each dimension, the number of total sampled points grows exponentially with the dimension of the input parameter space.

[A 5-point latin hypercube sampling] A 5-point latin hypercube sampling. Each parameter is sampled only once in every interval.

On the other hand, Latin hypercube sampling (LHS) doesn't depend on the dimension of the input. Instead, one specifies the number of samples $N$ one would like to draw in advance. The range of each variable is partitioned into $N$ non-overlapping intervals based on equal probability size. One value from each interval is selected at random with respect to the probability density in the interval. The $N$ values obtained for $x_1$ are paired in a random manner with the $N$ values for $x_2$. These $N$ pairs of $(x_1, x_2)$ are combined in a random manner with the $N$ values of $x_3$, and so on, until a set of $N$ $n$-tuples is formed.

§ BAYESIAN PARAMETER INFERENCE IN ACTION

§.§ A simple model for bulk physics

In this section I would present the Bayesian analysis results for a simple bulk model [151]. In this model, the final-state charged particle multiplicity and elliptic flow are related to an “evolution history averaged” viscosity:

\begin{equation}
N_{ch}=N_0 \cdot \Tilde{E}_T \left[1+(\eta/s)_{eff}\right], \ v_2=\epsilon_2 \cdot \exp\left[-50 \frac{(\eta/s)_{eff}}{N^{1/3}_{ch}}\right],
\end{equation}

where $(\eta/s)_{eff}$ is the temperature average of $\eta/s$:

\begin{equation}\label{eqn:eta_s_def}
(\eta/s)_{eff}=\frac{\int_{T_{min}}^{T_{max}}T^p(\eta/s) dT}{\int_{T_{min}}^{T_{max}}T^p dT}.
\end{equation}

The parameters $N_0=9, \ T_{min}=0.13, \ p=-1$ are assumed to be fixed. $T_{max}, \ \Tilde{E}_T, \ \epsilon_2$ are inputs depending on the centrality.

We propose a simple parameterization of $\eta/s$ with respect to temperature $T$:

\begin{equation}
\frac{\eta}{s}(T) = a +
\begin{cases}
b (T - T_s), & T>T_s,\\
c (T_s - T), & T<T_s,\\
\end{cases}
\end{equation}

$T_s=0.18$ is also fixed in this study but can be varied. And the parameters of our model are $a,b,c$.

One of the goals of this study is to go over the process of Bayesian analysis. Another goal is to test the validity of Bayesian analysis. The “experimental data” will be generated by the model, meaning we know the true value of the parameters and thus $\eta/s(T)$. Therefore we can also test whether we can constrain the model parameters given data with different levels of uncertainty. To achieve this goal, I need to minimize the uncertainties introduced by the model and the Gaussian process emulator.

§.§ Train surrogate model with Gaussian process

The first step in the Bayesian analysis is to generate the Latin hypercube sampling. In order to reduce the emulator uncertainty, 1000 sampling points are generated:

Latin hypercube sampling of 1000 design points.

As seen from Fig. <ref>, LHS yields almost uniform distribution on each of the parameters, as expected.

[Variance of the principal components] Left: Variance of the principal components by their index. Right: The ratio of the total variance of the first $n$ PC to the total variance of all the PC.

[Comparison between model calculation and emulator prediction at 5 randomly sampled parameter points] Comparison between model calculation and emulator prediction at 5 randomly sampled parameter points. The band of the GP represents $95\%$ CL.

The next step is to perform principal component analysis (PCA). From the model definition, one can see that $N_{ch}$ and $v_2$ are correlated. The total dimension of the “experiment” data $\mathbf{y}$ is $16$. But training 16 Gaussian process emulators would be redundant. In fact, after performing the PCA, one can see the variance is dominated by the first few components, and the first five will be used in this study.

Now the Gaussian process emulators on the selected principal components can be trained. Before proceeding further, it is important to check the performance of these emulators (emulator validation).

From Fig. <ref>, one can see the emulators can match with model calculations at those randomly generated new parameter points.

[Pseudo data generated with $15\%$, $5\%$, and $1\%$ statistical uncertainty] Pseudo data generated with $15\%$, $5\%$, and $1\%$ statistical uncertainty.

Now I would like to revisit the definition of the covariance matrix in Eq. <ref>. $\Sigma_{exp}$ can be further decomposed as $\Sigma_{exp}=\Sigma_{sys}+\Sigma_{stat}$. One often doesn't know the off-diagonal values in $\Sigma_{sys}$ from experiment publications but the diagonal values are generally available. As for the model covariance matrix $\Sigma_{M}$, the use of Gaussian process and principal component analysis both introduce uncertainties ($\Sigma_{emulator},\Sigma_{truncation}$) in addition to the statistical and systematic uncertainties from model calculations: $\Sigma_{M}=\Sigma_{emulator}+\Sigma_{truncation}+\Sigma_{model}$. In this example, the model calculation is exact. So only the experimental data uncertainty and emulator uncertainty are varied.

§.§ Inference with different levels of experimental data

The “experimental data” are prepared with three different levels of statistical uncertainties. If the data are sensitive to the parameters, we should be able to constrain the parameters even when sizable uncertainties are introduced. And the constraining power should increase as the uncertainties decrease.

Posterior distributions of the parameters $(a,b,c)$ inferred from data with $15\%$ of statistical uncertainties.

Posterior distributions of the parameters $(a,b,c)$ inferred from data with $5\%$ of statistical uncertainties.

Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties.

First, the posterior distribution of the parameters with $15\%$ of experimental data uncertainty is shown. In Fig. <ref>, the diagonal plots are the posterior distribution of each of the parameters. The off-diagonal plots display the joint distributions between pairs of the parameters. The truth values are indicated by black lines, the center of the posterior distribution are indicated by red lines and the values that maximize the posterior distribution (also referred as the maximum a posteriori value) are labeled by the blue lines. In this case, only $b$ is constrained while the MAP estimation is still off. Next, the result with $5\%$ of data uncertainty is shown in Fig. <ref>. $b$ is well constrained by the data and the prediction of $a$ and $c$ are much closer to their truth values. Lastly, the posterior distribution inferred from data with only $1\%$ of uncertainty is shown in Fig. <ref>. In this case, both $a$ and $b$ are constrained very well. However, parameter $c$ shows some discrepancy between the truth value and the two estimations.

Even though we failed to constrain the parameters for the case with $15\%$ of data uncertainty, the emulator prediction still matches the data pretty well (see Fig. <ref>).

[Comparison between data with $15\%$ uncertainty and emulator prediction with $95\%$ CL] Top row: Comparison between data with $15\%$ uncertainty and emulator prediction with $95\%$ CL.  Bottom row: Ratio between emulator prediction and data.

The posterior distribution of $\eta/s(T)$ is what people are really interested in. One can see that as the experimental uncertainty reduces, the posterior estimation gradually recovers the truth value (see Fig. <ref>).

[Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$]Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$. Top left: Inferred from data with $15\%$ uncertainty. Top right: Inferred from data with $5\%$ uncertainty. Bottom: Inferred from data with $1\%$ uncertainty.

As seen from Fig. <ref>, <ref>, <ref> and <ref>, parameter $b$ is always the most constrained among the three parameters as varying this parameter will affect the slope of $\eta/s(T)$ curve above $T_s$. Parameter $a$ is also constrained for the last two cases because varying $a$ will shift the entire $\eta/s(T)$ curve up and down. $c$ is least constrained since the temperature range it affects ($T_s-T_{min}=0.05$ GeV) is much smaller than the temperature range $a$ affects ($\approx 0.5$ GeV).

§.§ Inference with different levels of emulator uncertainty

[Comparison between model calculation and emulator prediction at 5 randomly sampled parameter points] Comparison between model calculation and emulator prediction at 5 randomly sampled parameter points. The band of the GP represents $95\%$ CL. The emulator only uses 20 design points and 3 principal components.

[Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties] Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties. The Gaussian process emulator is only using 50 design points and 5 principal components in this case.

[Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties] Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties. The Gaussian process emulator is only using 20 design points and 3 principal components in this case.

Next, I want to check the performance of Bayesian analysis when the uncertainties of the Gaussian process emulator are much larger than the experimental data uncertainties. The data uncertainty is fixed to be $1\%$ from now on. First, I only use $20$ design points to train the emulator and pick just the first $3$ of the principal components (as opposed to $1000$ and $5$ in the previous setup). While it is quite remarkable that Gaussian process emulator can still mimic the model (see Fig. <ref>), we are no longer able to constrain the parameters as shown in Fig. <ref>. However, the emulator is still able to predict the data (see Fig. <ref>). If we increase the number of design points to $50$ and pick the first $5$ principal components, the performance is surprisingly almost the same as the setup in the previous section (see Fig. <ref> and Fig. <ref>).

[Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$]Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$. Top Left: Inferred from data with $1\%$ uncertainty and Gaussian process emulators with minimal uncertainty (1000 design points and 5 principal components). Top Right: Inferred from data with $1\%$ uncertainty and Gaussian process emulators with only 50 design points and 5 principal components. Bottom: Inferred from data with $1\%$ uncertainty and Gaussian process emulators with only 20 design points and 3 principal components.

[Comparison between data with $1\%$ uncertainty and emulator prediction with $95\%$ CL] Top row: Comparison between data with $1\%$ uncertainty and emulator prediction with $95\%$ CL. The emulator is just using 20 design points and 3 principal components. Bottom row: Ratio between emulator prediction and data.

§.§ Inference with systematic model uncertainty

[Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties] Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties. The data are generated with the linear parameterization of $\eta/s(T)$ while the model is assuming a quadratic form.

[Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$]Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$. Left: Inferred from data with $1\%$ uncertainty and Gaussian process emulators with minimal uncertainty. Right: Inferred from data generated with a linear form and Gaussian process emulators trained with a quadratic form.

[Comparison between data with $1\%$ uncertainty and emulator prediction with $95\%$ CL] Top row: Comparison between data with $1\%$ uncertainty and emulator prediction with $95\%$ CL. The emulator is trained with a different model. Bottom row: Ratio between emulator prediction and data.

[Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties] Posterior distributions of the parameters $(a,b,c)$ inferred from data with $1\%$ of statistical uncertainties. The data are generated with the quadratic parameterization of $\eta/s(T)$ while the model is assuming a linear form.

[Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$]Comparison between the $95\%$ CL posterior distribution, the central estimation, the MAP estimation and the truth value of $\eta/s(T)$. Left: Inferred from data with $1\%$ uncertainty and Gaussian process emulators with minimal uncertainty. Right: Inferred from data generated with a quadratic form and Gaussian process emulators trained with a linear form.

[Comparison between data with $1\%$ uncertainty and emulator prediction with $95\%$ CL] Top row: Comparison between data with $1\%$ uncertainty and emulator prediction with $95\%$ CL. The emulator is trained with a different model. Bottom row: Ratio between emulator prediction and data.

In the previous two sections, I explored the scenarios with large data uncertainty or emulator uncertainty. The model, on the other hand, is used to generate the mock experimental data. In reality, our current modeling of nature is most likely incomplete or not precise, which translates to systematic model uncertainty. In this section, I introduce a different parameterization of $\eta/s(T)$ by changing its linear form into a quadratic form:

\begin{equation}
\frac{\eta}{s}(T) = a +
\begin{cases}
b (T - T_s)^2+c (T-T_s), & T>T_s\\
d (T_s - T)^2+e(T_s - T), & T<T_s\\
\end{cases}
\end{equation}

The mock data are still generated using the linear parameterization of $\eta/s(T)$ and with $1\%$ uncertainty. There are no truth values in this case, as we are using different models for training the emulator and generating the data. The quadratic parameterization will fall back to the linear one if we set the two second-order coefficients to zero. And this is what we see in the posterior distribution (see Fig. <ref>). The central values for the second-order coefficients are close to $0$, which is what the data suggests. The reason that we are not able to constrain the two second-order coefficients is probably due to their small effect on the final observables. Since $|T-T_s|<0.5$, the second-order term has a smaller impact on $\eta/s_{eff}$ compared to the linear term.

It is, however, still meaningful to compare the posterior distribution of $\eta/s(T)$ (see Fig. <ref>). We can see that by introducing two extra parameters, the posterior distribution with $95\%$ CL is much wider at low temperature and at higher temperature since we are constraining the two second order coefficients. There is no systematic model uncertainty in this case as the parameterization in the model can reduce to the one used for generating the data. But still much constraining power is lost.

The other direction can be explored as well, meaning the linear parameterization of $\eta/s$ is used in the model, but the data are generated with a quadratic parameterization. The posterior distribution is shown in Fig. <ref> and Fig. <ref>. The parameters seem to be well constrained although the parameterization used in the model is inherently wrong. The model with the linear parameterization is still able to describe the experiment data with very little statistical uncertainty (see Fig. <ref>). Because it is really difficult to constrain $\eta/s(T)$ if only integrated quantities over $\eta/s(T)$ are observed. If one probes even higher temperatures or more differential observables, maybe the systematic bias can show up in Fig. <ref>.

§ SUMMARY

This Chapter laid out the fundamentals for performing Bayesian analysis. We utilize techniques like Latin hypercube sampling, Gaussian process emulator, and principal component analysis to be able to quickly predict model calculations at arbitrary parameter points. Markov Chain Monte Carlo is then used to draw the posterior distributions of the parameters.

I then performed Bayesian analysis using a simple model for the bulk medium evolution with different scenarios for the covariance matrix and the $\eta/s(T)$ parameterization. The lessons learned from this exploration are the following:

* Reducing experimental data and emulator uncertainty can improve the constraining power. Even if the model has no systematic uncertainty, one may be unable to constrain the parameters if too much experimental or emulator uncertainty is introduced. The improvement of experimental uncertainty may require years in the field of heavy ion collision. The emulator uncertainty can be reduced by using more design points which also requires more computation time (by the scale of millions of CPU hours for our current project). Nevertheless, even with significant data uncertainty or emulator uncertainty, we are still able to recover the truth values of the parameters within our posterior distribution, meaning we are gaining knowledge about the truth values of the parameters by doing Bayesian analysis.

* Some parameters may be difficult to constrain if their impact on the final observables is small. In the model we have used, the slope for the lower part of the $\eta/s(T)$ curve (parameter $c$) is very challenging to constrain because it only governs a small range of temperature, and the observables we looked at are both integrated over $T$. If we want to constrain it, we should look at more sensitive observables (i.e., differential observables measured at low temperatures). In practice, we should include as many as experimental observables in our calibration to help constrain the parameters.

* Using a more complicated parameterization will likely reduce the constraining power while using a parameterization with insufficient flexibility will introduce systematic model uncertainty. It is essential to assume a reasonable parameterization with our prior knowledge about the underlying physics. A very good example is for the model where we know that the density deposition function $f(T_A, T_B)$ is scale invariant ($f(cT_A, cT_B)=cf(T_A, T_B)$). It is therefore natural to assume that the parameterization to be a generalized mean $f(T_A, T_B)=(\frac{T_A^p+T_B^p}{2})^{1/p}$[60]. If the form of the function is completely unknown, we can assume it to be in the form of a Taylor expansion or Gaussian random field [152].

CHAPTER: BAYESIAN MODEL-TO-DATA COMPARISON RESULTS

In this chapter, I will apply Bayesian analysis along with other techniques discussed in Chapter <ref> to estimate the various parameters explored in Chapter <ref>.

§ PREVIOUS ENERGY LOSS STUDIES USING BAYESIAN ANALYSIS

[Calculation of the D-meson $R_{AA}$ and $v_2$ in PbPb collisions at 5.02 TeV] Calculation of the D-meson $R_{AA}$ and $v_2$ in PbPb collisions at 5.02 TeV, taking the parameters randomly drawn from the posterior distributions. The improved Langevin model is used for parton energy loss. [153]

[Posterior predictive distributions of $R_{AA}$ using MATTER+LBT] Posterior predictive distributions of $R_{AA}$ using MATTER+LBT compared to data. Dashed lines show model calculation using median values of parameters. A virtuality dependent parameterization of $\hat{q}$ is used. [126]

[Emulated calculations of inclusive jet and single hadron suppression] Emulated calculations of inclusive jet and single hadron suppression after model parameters are systematically calibrated to these observables. In each set of plots, the $R_{AA}$ is shown in the upper panel, and the ratio to experiment is shown in the lower panel. The blue bands correspond to 95% credible limit given by the model emulator. The gray bands in the ratio plot denote $\pm10\%$ and $\pm20\%$ level of discrepancy respectively. [111]

The application of Bayesian analysis to jet energy loss in heavy ion collisions is inspired by the application of the same technique to soft observables in heavy ion collisions [43, 139]. Started in 2018, early attempts focused on open heavy flavor observables (see Fig. <ref>)[110, 78]. As only the leading heavy parton was tracked, oversampling (sample an ensemble of heavy flavor partons in one event) was possible. Therefore enough statistics could be achieved for calculating even the flow coefficients with high precision. Next, people started to study charged hadron and jet observables with Bayesian analysis [154, 126]. Specifically, Ref. [126] uses the JETSCAPE framework with the multi-stage MATTER+LBT approach for parton energy loss (see Fig. <ref>). More recently, Ref. [111] performed the calibration on charged hadron, D meson and inclusive jet $R_{AA}$ (see Fig. <ref>).

§ CALIBRATION SETUP

§.§ Parameter design

As we have seen in Chapter <ref>, parameters like $\tau_0, T_c$ have minor influence on $R_{AA}$ so they will be fixed in this analysis. What will be explored are the coupling constant $\alpha_s$, the switching scale $Q_s$ and $c_1,c_2$ in the parameterization of $\hat{q}$. The prior range for these parameters are considered to be:

Prior ranges for the parameters in our Bayesian calibration.

Parameter     Range

$\alpha_s$      0.1 - 0.5

$Q_s$      1.5 - 4

$c_1$       1 - 10

$c_2$       50 - 300

Distribution of input parameter $\alpha_s$ and $Q_s$ from 50 sampled design points within the prior range.

Those ranges are selected based on our simple exploration in Chapter <ref> and also in Ref. [137]. Due to constrained computation budget, I use 50 design points drawn from Latin hypercube sampling (see Fig. <ref>).

§.§ Experimental observables

There exist a wide range of experimental observables from different collision systems and at different collision energies and centralities. Due to a constrained computation budget, I will focus only on PbPb collision at 5.02 TeV and 0-10% centrality. A previous study [110] shows that calibrating to different collision energies independently or at the same time may yield slightly different posteriors or one may need to use different values for the same parameter in other collision systems [111]. I will leave such an exploration to a future study.

As for the experimental observables that will be in this calibration, I choose $R_{AA}$ for charged hadrons, D mesons, and inclusive jets. Calibrating to these observables together should give a better constraint on the transport properties of the QGP. One thing to notice is that data points for the charged hadron and D meson $R_{AA}$ below $7$ GeV are not considered in the calibration as they are affected by non-perturbative effects, medium response, recombination contributions to hadronization, etc. We could also consider flow observables like $v_2,v_3$ for the D mesons, but that would require around $10$ times more statistics to get meaningful results in central collisions. Another class of observables is the jet substructure observable which will be ignored as well since contributions from medium response and background subtraction are still being investigated.

§.§ Model calculation at the design points

Now I want to show the model calculation from all the design points as an indication of the prior range of our prediction in the observable space. If the model calculations from all design points are all above or below the data, it is unlikely we will be able to describe that observable. Fortunately, the model calculations span a wide range in $R_{AA}$ and can cover the data for all five observables (see Fig. <ref>).

For each design point, roughly 400k events are generated and are distributed among 400 fluid simulations which give rise to about 1000 events per fluid simulation. This number is only around $4\%$ of what was used in Chapter <ref>, so we observe pretty big fluctuations in the calculations. Those fluctuations will impact our calibration. I will verify the validity of our Bayesian analysis against model calculation fluctuations in Chapter <ref> and Chapter <ref>.

[Model calculation using parameters from all the design points] Model calculation using parameters from all the design points. From top to bottom, left to right, the measurements are charged hadron $R_{AA}$ from CMS [142], D meson $R_{AA}$ from ALICE [143], inclusive jet $R_{AA}$ from ATLAS [155] and two inclusive jet $R_{AA}$ with different jet radius from ALICE [155]. To better see how different calculations are distributed among these observables, each design point is assigned with a different color.

Comparison between emulator predictions and model calculations at 5 random design points.

[Comparison between emulator predictions and model calculations at all design points] Comparison between emulator predictions and model calculations at all design points. The first five plots are scatter plots that plots the model calculation and emulator prediction for each observable in pairs. If the model calculation and emulator prediction are the same, they should be positioned on the $y=x$ line (the black solid line). The last plot shows the histograms of the relative difference between model calculation and emulator prediction for different observables.

[Relation between variances and principal components] Left: The variance explained by each principal component. Right: The cumulative variance explained by the first $n$ principal components.

§ EMULATOR TRAINING AND VALIDATION

The emulator is the essential part that makes the Bayesian analysis possible. One can view it as a black box that accepts arbitrary input parameters and outputs predictions of the observables. One should first ensure that the emulator is working correctly before running the MCMC to get the posterior distribution of the parameters. Otherwise, one might be interpreting a biased posterior distribution. The key ingredients when training the emulator are the kernel to use, the number of principal components, and the treatment of uncertainties (covariance matrix). Below are the reasonings for choosing the specific settings for the emulator in this work:

* The kernel I use is a linear combination of the radial based function (RBF) kernel (Eq. <ref>) and the white noise kernel. The rationale for this particular choice is that the various $R_{AA}$'s are generally smooth and monotonic increasing functions (at least in the $p_T$ range we are interested in). So the smooth RBF kernel should be able to fit them well. The fluctuations in our simulation can be absorbed by the white noise kernel. I have also tried the Matérn kernel which is a generalization of the RBF kernel [150]. It has an additional variable $\nu$ that controls its smoothness (Matérn kernel with $\nu\rightarrow \infty$ is equivalent to the RBF kernel). A more detailed comparison between different kernels will be carried out in Section <ref> and Appendix. <ref>. Unless specifically stated, the kernel used in the emulator will be the RBF+white noise kernel.

* The number of principal components can be determined by examining the hyper parameters of the emulator for each PC (like the length scale $l$ in the RBF kernel or the noise level in the white noise kernel). One can expect that using too few PC will under-fit the training data. Using more PC can capture finer details of the model calculations but also risk of treating statistical fluctuations as real input. In this study, 5 PCs that explain around $95\%$ of the total variance show the best performance for closure tests and will be used hereafter (see Fig. <ref>). The effects of the number of PCs on closure tests will be studied in more detail in Section <ref>.

* As discussed in Chapter <ref>, the covariance matrix takes the form:

\begin{equation}
\Sigma=\Sigma_{sys}+\Sigma_{stat}+\Sigma_{emulator}+\Sigma_{truncation}+\Sigma_{model}.
\end{equation}

The $\Sigma_{emulator}$ contains the total covariance of all the Gaussian processes. The $\Sigma_{truncation}$ contains the total covariance of all the remaining principal components not considered by the emulator. $\Sigma_{model}$ is not considered in this analysis as it is challenging to quantify. As for the experimental covariance $\Sigma_{sys}+\Sigma_{stat}$, it is assumed to be diagonal in this analysis. The effect of off-diagonal terms in $\Sigma_{sys}$ will be discussed in Appendix. <ref>.

Some of the arguments listed above are more or less qualitative. In the next section, I will introduce a quantitative way to determine the optimal settings for the Gaussian process emulator.

§.§ Emulator validation

The mean $\mu$ standard deviation $\sigma$ when fitting Gaussian to the distribution of the relative difference between model calculation and emulator prediction.

Gaussian fit     Mean $\mu$     Standard deviation $\sigma$

$R_{AA}^{h^\pm}$, CMS      -0.0126     0.124

$R_{AA}^{D}$, ALICE      0.00360     0.297

$R_{AA}^{jet}$, ATLAS      0.000912     0.0677

$R_{AA}^{jet}$, ALICE      -0.00433     0.0892

$R_{AA}^{jet}$, ALICE      0.0000105     0.0835

Before running MCMC with the emulator, the emulator's performance on the training data should be validated. If the emulator can not reproduce the training data, it can not serve as a surrogate for the model. A direct comparison between emulator predictions and model calculations at $5$ random design points can be seen in Fig. <ref>. The emulator predictions fit the model calculations pretty well and cut off some statistical fluctuations in our model calculations. One can also notice that the emulator prediction gives larger uncertainty bands for model calculations with larger fluctuations.

To see the performance of the emulator on all design points, the model calculation and emulator prediction for each observable from all design points are plotted in pairs in the first five plots in Fig. <ref>. The emulator seems to perform the best predicting inclusive jet $R_{AA}$, followed by charged hadron $R_{AA}$ and finally the $D$ meson $R_{AA}$. In the last plot in Fig. <ref>, the histograms of the relative difference between model calculations and emulator predictions for different observables are shown. The distributions are all be fitted by a Gaussian distribution centered near the origin, meaning there is little systematic bias introduced into the emulator. The mean and standard deviation are listed in Table. <ref>. The uncertainty is around $6-9\%$ when predicting inclusive jet $R_{AA}$, $12\%$ when predicting charged hadron $R_{AA}$, and $30\%$ when predicting D meson $R_{AA}$. When calculating the emulator prediction at one design point, the training data will exclude data from that specific design point. This means we are training different emulators for each design point and the emulator does not know the truth values when making predictions. If all the design points are used to train a single emulator, the uncertainty of the relative difference will be slightly reduced (around $5-7\%$ when predicting inclusive jet $R_{AA}$, $10\%$ when predicting charged hadron $R_{AA}$, and $25\%$ when predicting D meson $R_{AA}$).

§.§ Closure test

[Mock data using simulation statistical uncertainty and using experimental uncertainty] Left Column: mock data using statistical fluctuations from simulation and no systematic fluctuations. Right Column: mock data using statistical fluctuations and systematic fluctuations from experiments.

Closure tests of $\alpha_s$ and $Q_s$ using 8 sets of mock data with statistical fluctuations from model calculations.

Closure tests of $c_1$ and $c_2$ plotted as the ratio of $H(Q^2)$ using 8 sets of mock data with statistical fluctuations from model calculations.

Closure tests of $\alpha_s$ and $Q_s$ using 8 sets of mock data with statistical and systematic fluctuations from experiments.

Closure tests of $c_1$ and $c_2$ plotted as the ratio of $H(Q^2)$ using 8 sets of mock data with statistical and systematic fluctuations from experiments.

Being able to predict the training data does not guarantee that our emulator can constrain the model parameters well. If the data are not sensitive to some parameters, those parameters may end up with a wide posterior. If the data are degenerate (can be described with multiple sets of parameters), the posterior distributions may end up with multiple peaks. If the emulator captures too much fluctuation in the model calculation, the distance in the parameter space may be distorted in the observable space. These scenarios can be checked by performing a cross-validation closure test, where one design point is taken out from the training process and treated as the truth. The emulator then calibrates to model calculations (mock data) at that design point. Since we know the truth values for the parameters in this case, we can make a comparison between the posterior distribution of the parameters and the truth.

There are two cases to consider for the fluctuations of the mock data. First, the statistical fluctuations from the model simulation can be used. In our study, this means small fluctuations for inclusive jet $R_{AA}$ and charged hadron $R_{AA}$ and large fluctuations for D meson $R_{AA}$. Second, we can use the statistical and systematic fluctuations from experiments. This amounts to large fluctuations for the inclusive jet $R_{AA}$ data at ALICE, and large fluctuations at high $p_T$ for the rest of the three $R_{AA}$ results. These two cases are shown in Fig. <ref> for a single design point. For each case, I will show the results for the closure test at $8$ random design points. Specifically, when showing the closure test results for $c_1$ and $c_2$, their posterior distributions are not directly shown. Instead, I introduce the modification to $\hat{q}^{HTL}$ in Eq. <ref> as:

\begin{equation}
H(Q^2, c_1, c_2) = H(t, c_1, c_2) =\frac{c_0}{1+c_1\ln^2(t)+c_2\ln^4(t)}=\frac{1+c_1\ln^2(t_0)+c_2\ln^4(t_0)}{1+c_1\ln^2(t)+c_2\ln^4(t)}.
\label{eq:Hq}
\end{equation}

The ratio between the posterior distributions and the truth value of $H(Q^2)$ is plotted. This is because $c_1$ and $c_2$ does not affect the $\hat{q}$ independently. The posterior of $H(Q^2)$ is what we are actually looking for.

In Fig. <ref>, the posterior distribution of $\alpha_s$ and $Q_s$ using $8$ random sets of mock data are shown. The truth values of $\alpha_s$ and $Q_s$ are shown by black lines. In Fig. <ref>, the ratio of $H(Q^2)$ between the posterior distribution and the truth are plotted. The results demonstrate that the $H(Q^2)$ are well enclosed by the inferred 60% and
95% confidence regions.

Next, the closure tests are done with experimental uncertainties as shown in Fig. <ref> and Fig. <ref>. A similar agreement with the truth values is observed. The two closure tests indicate that the emulator is stable with both simulation and experimental level of uncertainty and can generally infer the parameters from mock data.

The closure tests, however, rely on the fact that the mock data and the training data use the same underlying physics model. We still don't know if our model can describe the experiment observables given the appropriate choice of parameters.

§ A QUANTITATIVE METHOD TO DETERMINE THE OPTIMAL SETTINGS FOR THE EMULATOR

In the previous section, closure tests at $8$ random design points were performed with model statistical uncertainty or real experimental uncertainty for the covariance matrix $\Sigma_{exp}$. However, it is difficult to tell how well the emulator is doing by just looking at those posteriors, let alone compare the performance of one type of kernel with another type. I would like to find a quantitative measure of the emulator's performance. By performance, I mean how well the emulator can recover the truth from the mock data. Naturally, one can define this quantity $\Delta_d$, which measures the deviation from the truth value for a specific parameter while using the $d$-th design point as the truth:

\begin{equation}
\Delta_d= \int \left(\frac{|p-p_{truth}|}{|p_{max}-p_{min}|}\right)^l P(p)dp = \frac{1}{N}\sum_{j=1}^N \left(\frac{|p_j-p_{truth}|}{|p_{max}-p_{min}|}\right)^l,
\end{equation}

where the summation denotes a Monte Carlo sampling of the integration, and $l$ can be $1,2,3$ or just any positive real number. The larger $l$ is, the more penalty is put for values further away from the truth.

$\Delta_d$ is for a cross-validation test at one design point, we can then average over it over all design points and get the overall performance of the emulator:

\begin{equation}
\langle\Delta\rangle=\frac{1}{N_{design}}\sum_d \Delta_d.
\end{equation}

Suppose there are infinitely many design points uniformly distributed among the prior range, one can calculate the values of $\langle\Delta\rangle$ analytically with a uniform posterior or a Gaussian posterior distribution centered at the truth with different variances (see Table. <ref>). As one would expect, $\langle\Delta\rangle$ gets smaller when $l$ increases or the posterior centers more around the truth. Fig. <ref> shows the change in $\langle\Delta\rangle$ versus $\sigma$ in a Gaussian if a Gaussian posterior centered at the truth is used. $\langle\Delta\rangle$ approaches the value calculated with a uniform posterior when $\sigma\rightarrow\infty$, as expected.

Values of $\langle\Delta\rangle$ assuming a uniform posterior or a Gaussian posterior and infinite number of uniformly distributed design points.

$\langle\Delta\rangle$     uniform posterior     Gaussian ($\sigma=0.5$)     Gaussian ($\sigma=0.2$)     Gaussian ($\sigma=0.1$)

$l=1$            0.333                    0.267                          0.141                          0.0753

$l=2$            0.167                    0.112                         0.0323                         0.00904

$l=3$            0.100                    0.0588                        0.00954                       0.00139

[$\langle\Delta\rangle$ calculated with Gaussian posterior with different $l$ and $\sigma$] $\langle\Delta\rangle$ calculated with Gaussian posterior centered at the truth versus the variance $\sigma$ of the Gaussian posterior. The analytical result with a uniform posterior are shown with dashed lines.

When calculated with the $50$ design points used in this analysis, the $\langle\Delta\rangle$ values are shown in Table. <ref>. The results are very close to those in Table. <ref>, which confirms that Latin hypercube sampling (LHS) indeed samples the design points uniformly from the assumed uniform prior.

Values of $\langle\Delta\rangle$ for different parameters using a uniform posterior and the sampled design points.

$\langle\Delta\rangle(\alpha_s)$     $\langle\Delta\rangle(c_1)$     $\langle\Delta\rangle(c_2)$     $\langle\Delta\rangle(Q_s)$

$l=1$          0.333          0.333     0.333     0.333

$l=2$          0.167          0.166     0.166     0.167

$l=3$          0.100          0.0996     0.0994     0.100

§.§ Comparison with other performance metrics

[Comparison between different closure test metrics with a Gaussian posterior] Change in $P(p_{truth})$ (top left), $AIC$ (top right), $\langle\Delta\rangle \ (l=1)$ (middle left),  $\langle\Delta\rangle \ (l=2)$ (middle right),  $\langle\Delta\rangle \ (l=10)$ (bottom left) and  $\langle\Delta\rangle \ (l=1/5)$ (bottom right) with different bias and $\sigma$ of the Gaussian posterior.

There are other metrics for measuring the performance of the closure test, namely the posterior at the truth $P(p_{truth})$ and the Akaike Information Criterion (AIC) [156] which is defined as:

\begin{equation}
\end{equation}

where $\mathcal{L}_{max}$ is the maximum likelihood in the parameter space and $k$ is the number of parameters.

We can compare the variation of these metrics with a simple $1$-dimensional Gaussian posterior distribution centered away from the truth value $p_{truth}$ and some finite variance as shown in Fig. <ref>. The posterior at the truth and the $AIC$ flatten in some regions of the plane, which does not seem ideal. $AIC$ behaves particularly bad in this example as it is entirely indifferent to the change in the bias. $\langle\Delta\rangle$ looks more reasonable, yielding a larger penalty to both larger bias and larger $\sigma$. When the posterior becomes more complex, it becomes more complicated to compare one metric with another. I believe the story here is similar to the choice of the kernel in the Gaussian process emulator: there is no single kernel that can handle every possible situation. In this study, I will use $\langle\Delta\rangle$ with $l=2$ as the metric for calculating the deviation of the posterior from the truth in closure tests.

§.§ Connection to the Kullback-Leibler Divergence

The Kullback-Leibler Divergence is defined as:

\begin{equation}
D_{KL}(P|Q)=\int p(x)\ln\left(\frac{p(x)}{q(x)}\right)dx
\end{equation}

between two probability distribution $p(x)$ and $q(x)$. It is a measure of how one probability distribution $P$ is different from a second, reference probability distribution, $Q$. In Bayesian statistics, it can be seen as a measure of the information gain in moving from a prior distribution to a posterior distribution. We will calculate the $D_{KL}$ with the extracted posterior distribution in Section <ref>.

However, $D_{KL}$ is ill-defined (diverges) when the prior distribution $q(x)$ is a Dirac delta function centered at $x_0$ (here we use the limit representation $\delta(x-x_0)=\lim_{\sigma\rightarrow 0} \frac{1}{\sqrt{\pi}\sigma}e^{-(x-x_0)^2/\sigma^2}$):

\begin{equation}\label{eqn:DKL_equivalence}
\begin{split}
D_{KL}(P|Q)&=\int \lim_{\sigma\rightarrow 0} p(x)\ln\left(\frac{p(x)}{\frac{1}{\sqrt{\pi}\sigma}e^{-(x-x_0)^2/\sigma^2}}\right)dx \\
&=\int p(x)\ln p(x) dx + \lim_{\sigma\rightarrow 0} \ln{\sqrt{\pi}\sigma} + \lim_{\sigma\rightarrow 0}\frac{1}{\sigma^2} \int  p(x)(x-x_0)^2dx \\
&=\lim_{\sigma\rightarrow 0}\frac{1}{\sigma^2} \int  p(x)(x-x_0)^2dx + \mathcal{O}(\frac{1}{\sigma}).
\end{split}
\end{equation}

However, when we take the ratio between two $D_{KL}$, or argue that the divergence caused by $\frac{1}{\sigma^2}$ is independent from the posterior $p(x)$, the relevant part is:

\begin{equation}
\int  p(x)(x-x_0)^2dx,
\end{equation}

which is equivalent to the definition of $\Delta$ when $l=2$. If we use other limit representations of the Dirac delta function, for example:

\begin{equation}
\delta(x-x_0)=\lim_{\sigma\rightarrow 0} \frac{1}{2\sigma\Gamma(1+\frac{1}{l})}e^{-|x-x_0|^l/\sigma^l},
\end{equation}

where $l\ge 1$, we should get to the same equivalence relation with the specific $l$. When $l<1$, the derivation would be the same except that the relevant part is no longer the dominant contribution in Eq. <ref>.

The above derivations tell us that $\Delta$ is really measuring the finite relevant part in $D_{KL}$ when the prior is a Dirac delta function. $D_{KL}$ was not used before as a metric in closure test because of the divergence, now we have “regulated” it.

One should notice that we are actually measuring information loss in closure tests, since the prior is a Dirac delta function which is the most certain probability distribution. In this case, we want the $D_{KL}$ (information loss) to be as small as possible. When we employ a uniform prior distribution in Bayesian analysis, we are calculating information gain because a uniform prior means as least information as  possible. This time larger $D_{KL}$ means we are constraining the parameter better.

§.§ Comparison between different kernels

Six types of kernel are compared in this section:

* the RBF kernel (equivalent to the Matérn kernel with $\nu \rightarrow \infty$):

\begin{equation}
\end{equation}

* the Matérn ($\nu=3/2$) kernel:

\begin{equation}
\end{equation}

* the Matérn ($\nu=5/2$) kernel:

\begin{equation}
\end{equation}

* the RBF + white noise kernel.

* the Matérn ($\nu=3/2$) + white noise kernel.

* the Matérn ($\nu=5/2$) + white noise kernel.

A lesser common choice of $\nu$ in the Matérn kernel is $\nu=1/2$ which gives the following kernel function:

\begin{equation}
\end{equation}

Comparison between Matérn kernels with different $\nu$.

[The posterior function fitted by  Matérn kernels with different $\nu$.] The posterior function fitted by  Matérn kernels with different $\nu$. Red points represent the training data, the slid black line is the mean of the posterior, and the grey band represents the $1\sigma$ credible region.

The reason for choosing these half-integer values $p+1/2$ for $\nu$ is that the Matérn kernel can be simplified into a product of an exponential and a polynomial of order $p$. The Matérn kernels with different $\nu$ are plotted in Fig. <ref>. As we can see, a larger $\nu$ means a stronger correlation at small distances and a weaker correlation at large distances. This is also observed in Fig. <ref> where the posterior is trained with only one or two data point(s). The solid line is the mean and the grey band represents the $1\sigma$ credible region. The Matérn ($\nu=1/2$) kernel is used only in Fig. <ref> and Fig. <ref> for comparison purposes, since it predicts non differentiable functions which do not fit in our use case. In our use case where we use the Gaussian process to train a model that maps input parameters to observables, it is not clear how this mapping would look like for an arbitrary parameterization. In our brief parameter exploration in Section <ref>, it seems like the observables have monotonic relations with the model parameters. However, no conclusive observation can be drawn. This is another motivation for why we need to compare the performance of different kernels.

Experimental uncertainty is used in the covariance matrix $\Sigma_{exp}$. $\langle\Delta\rangle$ are calculated with different kernels and varied number of principal components. The results are shown in Fig. <ref>. Reading from the numbers in the plots, $\langle\Delta\rangle(\alpha_s)$ are generally the smallest, followed by $\langle\Delta\rangle(Q_s)$. $\langle\Delta\rangle(c_1)$ and $\langle\Delta\rangle(c_2)$ fluctuates around $0.167$ (indicated by the black dashed line) which is the analytical result assuming a uniform posterior. This means the emulators are having trouble recovering these two parameters given the current level of uncertainties.

Comparison of $\langle\Delta\rangle$ for different parameters with different kernels and number of principal components.

In order to pick the optimal settings for the emulator, the product of all the $\langle\Delta\rangle$ will be compared (see Fig. <ref>). In this test, the Matérn ($\nu=5/2$) + white noise kernel with $6$ principal components performs the best. However, the RBF + white noise kernel with $5$ principal components performs at a similar level, and we will stick to it for consistency with previous studies. Closure test result using the Matérn ($\nu=5/2$) + white noise kernel with $6$ principal components is shown in Appendix. <ref>. Another observation is that the white noise kernel improves the overall performance, especially for the Matérn kernel.

Another quantity we can look at is the variance of $\Delta$ ($\sigma^\Delta$). $\sigma^\Delta$ is not an indication of the emulator's performance, as it is possible to have large mean and small variance (or small mean and large variance) at the same time. Nevertheless, the variance limits the distribution of $\Delta$ to the right side of the mean $\langle\Delta\rangle$. In Fig. <ref>, the product of all the variance of $\Delta$ are shown. Interestingly, the smallest values are achieved at around $5$ or $6$ of PC and with the kernels that includes the white noise kernel.

The method discussed in this section is not limited to the particular Bayesian analysis in this work. It can be easily applied to other Bayesian analysis projects utilizing the Gaussian process emulator as a surrogate model. Key ingredients for training the emulator, including choice of the kernel and the number of principal components can all be determined via this method. The level of constraint on each parameter, can also be reflected by the magnitude of $\langle\Delta\rangle$. The caveat, however, is that these observations assume that we are using the correct model. No systematic uncertainty of the model is taken into account here.

Comparison of the product of all $\langle\Delta\rangle$ for different parameters with different kernels, training data selection, and number of principal components.

Comparison of the product of all $\langle\Delta\rangle$ for different parameters with different kernels, training data selection, and number of principal components.

§ BAYESIAN ANALYSIS RESULTS

The posterior distribution of the model parameters. The emulator is using 5 PC and the RBF + white noise kernel and trained from 50 design points.

[The prior and posterior distribution of $H(Q^2)$] Top: The prior (grey), 95% credible region of the posterior (blue), and 60% credible region of the posterior (deep blue) of $H(Q^2)$ define in Eq. <ref>. Bottom: The corresponding information gain (Kullback-Leibler divergence $D_{KL}$).

Comparison between the posterior distribution of the observables and experiment data.

In this section, I will first show the Bayesian analysis result calibrating to all five experiment data sets. This is the best estimation of the model parameters we can get at the current stage (see Fig. <ref>, Fig. <ref>). The posterior distribution of the observables compared to data are shown in Fig. <ref>. Compared to the parameters we used in Chapter <ref> (where $\alpha_s=0.3, \ c_1=10, \ c_2=100, \  Q_s=2$), the posterior distribution in Fig. <ref> suggests similar value for $\alpha_s$ but a slightly smaller $Q_s$. From previous variation of $Q_s$ we know a smaller $Q_s$ would shift the charged hadron and D meson $R_{AA}$ up, possibly making up for the disagreement in charged hadron $R_{AA}$ at high $p_T$ seen in Fig. <ref>. Meanwhile, the constraint on $c_1$ is much wider and a peak over large $c_2$ is observed. However, from the tests done in Sec. <ref>, the constraints on $c_1$ and $c_2$ should be not be taken too seriously as on average the emulator's performance is just slightly better than using a uniform posterior. One can also see why $c_1$ and $c_2$ are hard to constrain by looking at Fig. <ref>. Even with much higher precision, the D meson and charged hadron $R_{AA}$ results are difficult to distinguish between different $\hat{q}$ parameterization except at high $p_T$ where the experimental data uncertainty are still large.

In Fig. <ref>, the top plot shows the prior, 95% credible region, and 60% credible region of $H(Q^2)$ as a function of $\ln(Q^2)$. The Kullback-Leibler divergence (defined in Section <ref>) peaked near the lower limit of the prior range of $Q_s$ which is $1.5$. This is likely due to the fact that the second-order term and the fourth-order term in $H(Q^2)$ are comparable in magnitude in this region. This shows that the constraint is more on the joint distribution of $c_1$ and $c_2$ at low $Q$ than their individual values. At large $Q$, the fourth-order term dominates $H(Q^2)$ and the small $D_{KL}$ values here suggest that not much information is gained for the posterior of $c_2$.

§.§ Sensitivity to different observables

How different observables help with constraining the parameters is another interesting topic. There are essentially three types of observables we are calibrating on (charged hadron, D meson, and inclusive jet $R_{AA}$). In this section, I start with only calibrating to the three inclusive jet $R_{AA}$ results (see Fig. <ref>). Compared to Fig. <ref>, the major difference is a strong preference of large $c_2$ and weaker constraint on $Q_s$. Next, the calibration is done with both the inclusive jet $R_{AA}$ results and the $D$ meson $R_{AA}$ result as shown in Fig. <ref>. This time a combination of small $c_1$ and large $c_2$ is slightly preferred. The constraint on $Q_s$ is still weaker compared to Fig. <ref>. We can further divide the data into different $p_T$ ranges and test the sensitivity of the posterior to different combinations of the segmented data. I'll leave that to a future study.

The posterior distribution of the model parameters. The emulator is only calibrating to the three inclusive jet $R_{AA}$ results.

The posterior distribution of the model parameters. The emulator is calibrating to the three inclusive jet $R_{AA}$ results and the $D$ meson $R_{AA}$ result.

§ SUMMARY

In this chapter, I have applied Bayesian analysis to the multi-stage energy loss model discussed in Chapter <ref>, and Chapter <ref>. I first performed closure tests on the emulator to validate its performance. With a quantitative measure $\langle\Delta\rangle$, the optimal settings for the Gaussian process emulator are determined, and the emulator's constraining power on each parameter is studied. Then I showed the sensitivity of the posterior distribution to different observables. In the end, I have successfully constrained the $\alpha_s$ and $Q_s$ and explained why $c_1$ and $c_2$ are difficult to constrain within the current setup. Continued discussion of the stability of the posterior distribution to additional fluctuations and the effect of correlated experimental systematic uncertainty or using the Matérn kernel in the emulator can be found in Appendix. <ref>. This analysis is the first result of a simultaneous description of $R_{AA}$ data for multiple hadron species and inclusive jets with a multi-stage energy loss description.

However, there are still many things to improve in our analysis:

* Include more realistic modeling into the JETSCAPE framework. For example, the MATTER model could use different $Q_s$ and $\hat{q}$ for different quark species. The $\hat{q}$ parameterization could be more flexible and dependent on quark mass. We should also apply a hybrid hadronization model for the heavy quarks and hadronic rescattering which are absent in current simulations. A more complex pre-equilibrium modeling and a (3+1)D hydrodynamic background are also things to improve upon.

* Include more observables into our calibration. Due to the current limitations, the flow coefficients are hard to simulate to a reasonable precision, but are helpful for studying the geometry and flow of the medium. Considering data from more centrality and collision systems is also important.

* Increase the statistics for each design point and the number of design points. The current level of model statistical uncertainty is not enough to constrain $c_1$ and $c_2$ (see Fig. <ref>).

CHAPTER: CONCLUSION

Ultra-relativistic heavy ion collisions provide the opportunity for studying the extremely hot and dense form of matter, namely the quark gluon plasma (QGP). The degrees of freedom in the QGP, as its name suggests, are quarks and gluons. The dynamics of the QGP is believed to be dominated by quantum chromodynamics (QCD).

In this thesis, I have adopted a multi-stage approach for simulating the evolution of the QGP during the collision. Specifically, the description is divided into two parts depending on the energy scale we are looking at: the bulk medium and the hard probes. The evolution of the bulk medium which produces large number of soft hadrons in the final state has been successfully described by relativistic viscous hydrodynamics with phenomenology modeling of the initial condition and stopping condition.
As for the hard probes (leading partons with $p_T>10$ GeV, jets), they are created at the early stage of the collision and acting as probes to the evolution of the medium. In this study, their creation is simulated by PYTHIA and their interaction with the QGP medium described by transport models like MATTER and LBT. The MATTER model takes care of the in-medium DGLAP evolution of the highly virtual partons. LBT describes the evolution of on-shell partons inside the medium and contains both the elastic and inelastic scattering kernel. A virtuality dependent parameterization of the transport coefficient $\hat{q}$ is proposed in order to explain the smaller $\hat{q}$ extracted from experiments with higher collision energy.

All these evolution models have been integrated into the JETSCAPE framework. This allows a systematic comparison between different energy loss mechanisms which is carried out in Section <ref>. The effects of employing MATTER, LBT, or the multi-stage MATTER+LBT approach on the charged hadron and D meson $R_{AA}$ have been studied. We found that the MATTER model with the virtuality dependent $\hat{q}$, is approaching the vacuum DGLAP at high $p_T$. However, it effectively reduces the path length the hard parton spent in the LBT phase, thus enhances the $R_{AA}$. Using the current “optimal” parameters and the MATTER+LBT approach, we have predicted a slightly lower $R_{AA}$ for D meson compared to charged hadron at high $p_T$ across multiple centralities. The experimental data is yet to confirm this observation due to the lack of D meson $R_{AA}$ beyond $100$ GeV.

Despite a good description of the experimental data with the MATTER+LBT approach, it is difficult to determine the optimal values for the model parameters and how far we can get in terms of describing the experimental data. The uncertainty in experimental data and simulation should also be taken into account. For that, I have utilized a state-of-the-art Bayesian model-to-data comparison framework which has been proven successful when applied to the soft sector of heavy ion collisions. In fact, the bulk medium evolution is simulated with the optimal extracted parameters from a previous Bayesian model-to-data comparison.
The effects of emulator uncertainty, experiment uncertainty and model uncertainty, have been investigated in Section <ref> and Section <ref>.
I have proposed a new metric for measuring the performance of the Gaussian process emulator against the choice of kernel and number of principal components in closure tests. The sensitivity of the posterior distribution of the parameters to different observables and uncertainty levels are investigated before we draw a conclusion from the final posterior distribution calibrated to charged hadron, D meson and inclusive jet $R_{AA}$ simultaneously. We found a similar peak of $\alpha_s$ and a slightly lower peak of $Q_s$ compared to our previous manual parameter fitting. And the current level of uncertainty from both experiment and simulation failed to constrain the $c_1$ and $c_2$ parameter in the parametrization of $\hat{q}$.

Our study is a step forward towards a simultaneous description of energy loss for multiple hadron species and jet types in heavy ion collisions, with proper treatment of the uncertainties from different sources and a qualitative method for choosing the optimal Gaussian process emulator. However, there are still many places that can be improved. On the theoretical side, the recombination mechanism and hadronic rescattering should be included in the energy loss calculation. A more flexible parameterization of the transport coefficient could also be considered. When calibrating to data, the Bayesian analysis should include more observables, and the simulation statistics for a single design point as well as the total number of design points could be increased.

With better and expanded experimental measurements from LHC (Run 3) and upcoming electron ion collider (EIC), and advancement in both theoretical modeling and computing power, we are looking forward to a better understanding of the properties of QGP in the coming years.

CHAPTER: DERIVING THE LANGEVIN EQUATION FROM THE FOKKER-PLANCK EQUATION

One way to show the equivalence between the Fokker-Planck and the Langevin equation is to covert the Langevin equation to a path integral expression and recognize the Fokker-Planck equation as an Euclidean Schrodinger equation, which also has a path integral representation [157]. Another way is first to derive the equation of motion for the probability distribution $\rho(\vec{x}, \vec{p}, t)$ of finding a particle in the interval $(x,x+dx), (p,p+dp)$ at time $t$ for one realization of the random kick $\xi(t)$, then average $\rho(\vec{x}, \vec{p}, t)$ over a ensemble of $\xi(t)$.

The normalization constraint on $\rho(\vec{x}, \vec{p}, t)$ reads:

\begin{equation}
\int_{-\infty}^{\infty} d^3x_i \int_{-\infty}^{\infty} d^3p_i \rho(\vec{x}, \vec{p}, t)  = 1,
\end{equation}

which gives the continuity equation in the phase space:

\begin{equation}
\begin{split}
\frac{\partial}{\partial t} \rho(\vec{x}, \vec{p}, t) & = - \frac{\partial }{\partial x_i} (\frac{\partial x_i}{\partial t} \cdot \rho(\vec{x}, \vec{p}, t)) - \frac{\partial }{\partial p_i} (\frac{\partial p_i}{\partial t} \rho(\vec{x}, \vec{p}, t))\\
& = - \frac{\partial }{\partial x_i} \left(\frac{p_i}{E} \cdot \rho(\vec{x}, \vec{p}, t)\right) - \frac{\partial}{\partial p_i} \left((-\eta_D p_i + \xi_i))\cdot \rho(\vec{x}, \vec{p}, t)\right) \\
& = \left(-\frac{p_i}{E} \frac{\partial}{\partial x_i} + p_i \frac{\partial}{\partial p_i} \eta_D p_i\right) \rho(\vec{x}, \vec{p}, t) - \xi_i(t) \frac{\partial}{\partial p_i} \rho(\vec{x}, \vec{p}, t) \\
& = - L_0 \rho(\vec{x}, \vec{p}, t) - L_1(t) \rho(\vec{x}, \vec{p}, t),
\label{eqn:probability_flow_equation}
\end{split}
\end{equation}

\begin{align}
& L_0 = \frac{p_i}{E} \frac{\partial}{\partial x_i} - \frac{\partial}{\partial p_i} \eta_D p_i\\
& L_1 = \xi_i(t) \frac{\partial}{\partial p_i}.
\end{align}

Performing a change of variables $\rho(\vec{x}, \vec{p}, t) = e^{-L_0 t} \sigma(\vec{x}, \vec{p}, t)$ will transform Eq. <ref> into:

\begin{equation}
\frac{\partial}{\partial t} \sigma(\vec{x}, \vec{p}, t) = - e^{L_0 t} L_1(t) e^{-L_0 t} \sigma(\vec{x}, \vec{p}, t) = - V(t) \sigma(\vec{x}, \vec{p}, t),
\end{equation}

which is now a first order linear partial differential equation (PDE) and has the following solution:

\begin{equation}
\sigma(t) = \exp \left[- \int_0^t dt' V(t')\right]\sigma(0).
\end{equation}

Performing an ensemble average over the Gaussian random kick $\xi(t)$ leads to:

\begin{equation}
\left<\sigma(t)\right>_{\xi} = \exp\left[\frac{1}{2} \int_0^t dt_1 \int_0^t dt_2 \left<V(t_1) V(t_2)\right>_{\xi}\right] \sigma(0).
\label{eqn:characteristic function}
\end{equation}

And the integral in Eqn. <ref> can be calculated as:

\begin{align*}
\frac{1}{2} \int_0^t dt_1 \int_0^t dt_2 \left<V(t_1) V(t_2)\right>_{\xi} &= \frac{1}{2} \int_0^t dt_1 \int_0^t dt_2 \left<e^{L_0 t_1} L_1(t_1) e^{-L_0 t_1} e^{L_0 t_2} L_1(t_1) e^{-L_0 t_1}\right> \\
& = \frac{1}{2} \int_0^t dt_1 \int_0^t d t_2 \left<e^{L_0 t_1} \xi_i(t_1) \frac{\partial}{\partial p_i} e^{-L_0 t_1} e^{L_0 t_2} \xi_j(t_2) \frac{\partial}{\partial p_j} e^{-L_0 t_1}\right> \\
& = \frac{1}{2} \int_0^t dt_1 \int_0^t d t_2 \left<e^{L_0 t_1}  \frac{\partial}{\partial p_i} e^{-L_0(t_1-t_2)}  \xi_i(t_1) \xi_j(t_2) \frac{\partial}{\partial p_j} e^{-L_0 t_1}\right> \\
& = \frac{1 }{2} \int_0^t dt_1 e^{L_0t_1}\frac{\partial}{\partial p_i} \kappa^{ij} \frac{\partial}{\partial p_j} e^{-L_0 t_1},
\end{align*}

which means:

\begin{equation}
\frac{\partial}{\partial t} \left<\sigma(t)\right>_{\xi} = \frac{1}{2} e^{L_0 t} \frac{\partial}{\partial p_i} \kappa^{ij} \frac{\partial}{\partial p_j} \left<\sigma(t)\right>_{\xi},
\end{equation}

If one defines $p(\vec{x}, \vec{p}, t)=\left<\rho(\vec{x}, \vec{p}, t)\right>_{\xi}=e^{-L_0 t}\left<\sigma(\vec{x}, \vec{p}, t)\right>_{\xi} $, Eq. <ref> becomes:

\begin{equation}\label{eqn:fokker_plank_equation_full}
\frac{\partial}{\partial t} p(\vec{x}, \vec{p}, t) =  - \frac{p_i}{E} \frac{\partial}{\partial x_i} p(\vec{x}, \vec{p}, t) - \frac{\partial}{\partial p_i}(\eta_D p_i p(\vec{x}, \vec{p}, t)) + \frac{1}{2} \frac{\partial}{\partial p_i} \kappa^{ij} \frac{\partial}{\partial p_j} p(\vec{x}, \vec{p}, t).
\end{equation}

This is the Fokker-Planck equation for the full phase distribution $p(\vec{x}, \vec{p}, t)$. The equivalence between the Fokker-Plank and Langevin equation under a Gaussian random force is now proved.

CHAPTER: MORE ON BAYESIAN ANALYSIS

§ CLOSURE TEST RESULTS WITH A SIMPLE ANALYTICAL BULK MODEL

[Comparison of $\langle\Delta\rangle$ for different parameters with different kernels and number of principal components] Comparison of $\langle\Delta\rangle$ for different parameters with different kernels and number of principal components. The last plot shows the product of all $\langle\Delta\rangle$. $1\%$ of statistical model uncertainty is introduced. The Matérn kernel uses $\nu=3/2$.

[Comparison of $\langle\Delta\rangle$ for different parameters with different kernels and number of principal components] Comparison of $\langle\Delta\rangle$ for different parameters with different kernels and number of principal components. The last plot shows the product of all $\langle\Delta\rangle$. $5\%$ of statistical model uncertainty is introduced. The Matérn kernel uses $\nu=3/2$.

[Comparison of $\langle\Delta\rangle$ for different parameters with different kernels, number of principal components and model statistical uncertainties] Comparison of $\langle\Delta\rangle$ for different parameters with different kernels, number of principal components and model statistical uncertainties. The last plot shows the product of all $\langle\Delta\rangle$. The Matérn kernel uses $\nu=3/2$.

We can also calculate $\langle\Delta\rangle$ using the analytical bulk model discussed in Section <ref> which has three parameters $a$, $b$, and $c$. The benefit of using this simple model is that the level of statistical model uncertainty can be tuned very easily. $50$ design points are sampled and four levels of statistical model uncertainty $1\%,5\%,10\%,15\%$ are investigated.
In Fig. <ref> and Fig. <ref>, the effect of additional statistical model uncertainty on the choice of kernel can be seen. The inclusion of the white noise kernel helps reducing $\langle\Delta\rangle$.

Fig. <ref> compares the $\langle\Delta\rangle$ results between different levels of statistical model uncertainty. As expected, the smaller the statistical model uncertainty, the smaller $\langle\Delta\rangle$ for all parameters. However, one should notice that even with just $1\%$ of statistical model uncertainty, $\langle\Delta\rangle(c)$ is around $0.13$ which translates to a Gaussian posterior with $\sigma\approx 0.2$ on average. This is consistent with the findings from Section <ref> that the parameter $c$ is the most difficult to constrain as it only affects a small temperature range in the $\eta/s$ parameterization.

§ STABILITY OF THE POSTERIOR TO FLUCTUATIONS

The closure test offers a excellent check that the emulator can reproduce mock data (model calculation) at many design points. But the relation between the posterior distribution and the uncertainty level of the training data is still not explored. It is difficult to reduce the fluctuations of the training data as it requires running more events. The other direction is easier to explore and can be investigated in two ways.

[Relation between variances and principal components] Left: The variance explained by each principal component. Right: The cumulative variance explained by the first $n$ principal components. Only $1/3$

First, we can reduce the statistics when generating the training data. Only $1/3$ of events are now used for each design point. This way the statistical fluctuations for each observable are scaled up proportionally. In Fig. <ref>, we can see that now $5$ principal components only explains around $86\%$ of the total variance compared to $95\%$ before.  Fig. <ref> shows the posterior distribution of the parameters. Compared to whats shown in Fig. <ref>, the posterior distribution is not altered much in this case, except for less constraint on $c_1$ and $c_2$.

The posterior distribution of the model parameters. Only $1/3$ of events are used for each design point.

One question can be answered is how much does going from using $1/3$ of the statistics to using all the statistics improve the constraint on the parameters. In Fig. <ref>, we can see that $\langle\Delta\rangle(\alpha_s)$ and $\langle\Delta\rangle(Q_s)$ does improve. However, $\langle\Delta\rangle$ for $c_1$ and $c_2$ still float around the dashed line representing result calculated with a uniform posterior, showing that $c_1$ and $c_2$ are indeed the most difficult to constrain. It is still unclear how much using more design points and more statistics for each design point can help constrain these two parameters at this moment.

Comparison of $\langle\Delta\rangle$ for different parameters with the same settings for the emulator but different statistics for the training data.

Another way is to add Gaussian noise to all the model calculations for each design point. This time every observable gets an equal amount of additional statistical fluctuation. In Fig. <ref>, Fig. <ref>, Fig. <ref>, and Fig. <ref>, we can see the constraints on each observable gradually reduce as more fluctuations are introduced. One should notice that additional noise with $0.02$ standard deviation is already larger than the model statistical fluctuation for low $p_T$ $R_{AA}$ data.

The posterior distribution of the model parameters. A Gaussian noise with zero $\mu=0$ and $\sigma=0.02$ is added to all the training data.

The posterior distribution of the model parameters. A Gaussian noise with zero $\mu=0$ and $\sigma=0.05$ is added to all the training data.

The posterior distribution of the model parameters. A Gaussian noise with $\mu=0$ and $\sigma=0.1$ is added to all the training data.

The posterior distribution of the model parameters. A Gaussian noise with $\mu=0$ and $\sigma=0.2$ is added to all the training data.

§ EXPERIMENTAL COVARIANCE

At the moment the experimental covariance matrix is assumed to be diagonal since the off diagonal correlations are not reported by experiments. However there are possible systematic correlations among different observables. A simple way to account for such correlation is to assume a Pearson correlation only among observables in the same class:

\begin{equation}
\Sigma_{exp}^{sys}=\frac{1}{\sigma_i \sigma_j} exp[-\frac{1}{2}(\frac{x_{i}-x_j}{l})^2]
\end{equation}

where $\sigma_i, \sigma_j$ are the standard deviation of the two observable, $x_i-x_j$ measures the distance between the two observable in the space they are defined (for example, the distance in $p_T$ when measuring $R_{AA}$). The hyper-parameter $l$ controls the correlation length and can be varied. In principle, a larger $l$ means more correlation across the measurements and less constraint on the parameters. The posterior distribution using this experimental covariance matrix is shown in Fig. <ref> and similar qualitative behavior is observed as in Fig. <ref>.

[The posterior distribution of the model parameters] The posterior distribution of the model parameters. A off diagonal systematic correlation with $l=0.1(p_{T,max}-p_{T,min})$ is assumed to be present within each observed $R_{AA}$.

§ POSTERIOR DISTRIBUTION USING THE MATÉRN KERNEL

The choice of kernel in the emulator also plays a important role in Bayesian analysis. The Matérn kernel [150] encodes prior assumption of the smoothness of the underlying function. The posterior distribution using the Matérn($\nu=5/2$)+white noise kernel with $6$ principal components is shown in Fig. <ref>. This is the settings found to give the smallest product of all $\langle\Delta\rangle$ in Section <ref>. Similar posterior distribution is seen for $\alpha_s$ and $Q_s$ compared to the RBF + white noise kernel in Fig. <ref>. The main difference is that now there is a also strong peak for $c_1$. As mentioned in Section <ref>, we shouldn't expect tight constraint on $c_1$ and $c_2$ as the average $\langle\Delta\rangle$ for these two parameters is close to the one calculated with a uniform posterior.

The result of closure test using 8 random design points is shown in Fig. <ref> and Fig. <ref>. No qualitative difference is found compared to Fig. <ref> and Fig. <ref> which uses the RBF+white noise kernel. This is why $\langle\Delta\rangle$ or other metrics is needed for a quantitative comparison between different kernels.

[The posterior distribution of the model parameters] The posterior distribution of the model parameters. The emulator is using the Matérn kernel with $\nu=5/2$ and $N_{PC}=6$.

[Closure tests of $\alpha_s$ and $Q_s$ using 8 sets of mock data] Closure tests of $\alpha_s$ and $Q_s$ using 8 sets of mock data. The emulator is using the Matérn kernel with $\nu=5/2$ and $N_{PC}=6$.

[Closure tests of $c_1$ and $c_2$ plotted as the ratio of $H(Q^2)$ using 8 sets of mock data] Closure tests of $c_1$ and $c_2$ plotted as the ratio of $H(Q^2)$ using 8 sets of mock data. The emulator is using the Matérn kernel with $\nu=5/2$ and $N_{PC}=6$.

[4]
P. D. Group et al., “Review of particle physics,” The European
Physical Journal C-Particles and Fields, vol. 3, no. 1-4, pp. 1–783, 1998.

[5]
A. Bazavov et al., “Equation of state in ( 2+1 )-flavor QCD,” Phys. Rev., vol. D90, p. 094503, 2014.

[6]
. N. S. A. Committee et al., “Reaching for the horizon: The 2015 long
range plan for nuclear science,” 2017.

[7]
A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann,
S. Mukherjee, P. Petreczky, C. Schmidt, D. Smith, et al., “Freeze-out
conditions in heavy ion collisions from qcd thermodynamics,” Physical
review letters, vol. 109, no. 19, p. 192302, 2012.

[8]
C. Gattringer and K. Langfeld, “Approaches to the sign problem in lattice
field theory,” International Journal of Modern Physics A, vol. 31,
no. 22, p. 1643007, 2016.

[9]
M. Buballa, “Njl-model analysis of dense quark matter,” Physics
Reports, vol. 407, no. 4, pp. 205–376, 2005.

[10]
A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, “Mapping the
phases of quantum chromodynamics with beam energy scan,” Physics
Reports, vol. 853, pp. 1–87, 2020.

[11]
M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, “Glauber modeling
in high energy nuclear collisions,” Ann. Rev. Nucl. Part. Sci.,
vol. 57, pp. 205–243, 2007.

[12]
Z. Qiu, Event-by-event hydrodynamic simulations for relativistic heavy-ion
The Ohio State University, 2013.

[13]
J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, “Hydrodynamic predictions
for 5.02 TeV Pb-Pb collisions,” Phys. Rev., vol. C93, no. 3,
p. 034912, 2016.

[14]
H. Niemi, K. J. Eskola, R. Paatelainen, and K. Tuominen, “Predictions for
5.023 TeV Pb + Pb collisions at the CERN Large Hadron Collider,” Phys.
Rev., vol. C93, no. 1, p. 014912, 2016.

[15]
K. Aamodt et al., “Higher harmonic anisotropic flow measurements of
charged particles in Pb-Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV,” Phys. Rev. Lett., vol. 107, p. 032301, 2011.

[16]
J. E. Bernhard, J. S. Moreland, and S. A. Bass, “Bayesian estimation of the
specific shear and bulk viscosity of quark–gluon plasma,” Nature
Physics, vol. 15, no. 11, pp. 1113–1117, 2019.

[17]
R. D. Weller and P. Romatschke, “One fluid to rule them all: viscous
hydrodynamic description of event-by-event central p+ p, p+ pb and pb+ pb
collisions at s= 5.02 tev,” Physics Letters B, vol. 774, pp. 351–356,

[18]
S. Chatrchyan et al., “Study of high-pT charged particle suppression in
PbPb compared to $pp$ collisions at $\sqrt{s_{NN}}=2.76$ TeV,” Eur.
Phys. J., vol. C72, p. 1945, 2012.

[19]
M. Aaboud, G. Aad, B. Abbott, B. Abeloos, D. Abhayasinghe, S. Abidi,
O. AbouZeid, N. Abraham, H. Abramowicz, H. Abreu, et al., “Measurement
of the nuclear modification factor for inclusive jets in pb+ pb collisions at
$\sqrt{s_{NN}}$= 5.02 tev with the atlas detector,” Physics Letters B,
vol. 790, pp. 108–128, 2019.

[20]
V. Khachatryan et al., “Charged-particle nuclear modification factors
in PbPb and pPb collisions at $ \sqrt{s_{NN}}=5.02 $ TeV,” JHEP,
vol. 04, p. 039, 2017.

[21]
A. M. Sirunyan et al., “Measurement of the ${B}^{\pm}$ Meson Nuclear
Modification Factor in Pb-Pb Collisions at $\sqrt{{s}_{NN}}=5.02\text{
}\text{ }\mathrm{TeV}$,” Phys. Rev. Lett., vol. 119, no. 15,
p. 152301, 2017.

[22]
A. M. Sirunyan et al., “Nuclear modification factor of D$^0$ mesons in
PbPb collisions at $\sqrt{s_\mathrm{NN}} = 5.02$ TeV,” Phys. Lett.,
vol. B782, pp. 474–496, 2018.

[23]
X. Yao, W. Ke, Y. Xu, S. A. Bass, and B. Müller, “Coupled boltzmann
transport equations of heavy quarks and quarkonia in quark-gluon plasma,”
Journal of High Energy Physics, vol. 2021, no. 1, pp. 1–38, 2021.

[24]
B. Müller, The physics of the quark-gluon plasma, vol. 225.

[25]
S. Chatrchyan, V. Khachatryan, A. Sirunyan, A. Tumasyan, W. Adam, T. Bergauer,
M. Dragicevic, J. Erö, C. Fabjan, M. Friedl, et al., “Multiplicity
and transverse momentum dependence of two-and four-particle correlations in
ppb and pbpb collisions,” Physics Letters B, vol. 724, no. 4-5,
pp. 213–240, 2013.

[26]
S. Zhang, Y. Ma, G. Ma, J. Chen, Q. Shou, W. He, and C. Zhong, “Collision
system size scan of collective flows in relativistic heavy-ion collisions,”
Physics Letters B, vol. 804, p. 135366, 2020.

[27]
M. Gyulassy, “The qgp discovered at rhic,” in Structure and dynamics of
elementary matter, pp. 159–182, Springer, 2004.

[28]
B. Back, M. Baker, M. Ballintijn, D. Barton, B. Becker, R. Betts, A. Bickley,
R. Bindel, A. Budzanowski, W. Busza, et al., “Centrality dependence of
charged-hadron transverse-momentum spectra in d+au collisions at
$\sqrt{s_{NN}}$= 200 gev,” Physical review letters, vol. 91, no. 7,
p. 072302, 2003.

[29]
S. S. Adler, S. Afanasiev, C. Aidala, N. Ajitanand, Y. Akiba, A. Al-Jamel,
J. Alexander, K. Aoki, L. Aphecetche, R. Armendariz, et al., “Absence
of suppression in particle production at large transverse momentum in
$\sqrt{s_{NN}}$= 200 gev d+au collisions,” Physical review letters,
vol. 91, no. 7, p. 072303, 2003.

[30]
B. Abelev, J. Adam, D. Adamova, A. Adare, M. Aggarwal, G. A. Rinella,
M. Agnello, A. Agocs, A. Agostinelli, Z. Ahammed, et al., “Transverse
momentum distribution and nuclear modification factor of charged particles in
p+ pb collisions at $\sqrt{s_{NN}}$= 5.02 tev,” Physical review
letters, vol. 110, no. 8, p. 082302, 2013.

[31]
V. Khachatryan, A. M. Sirunyan, A. Tumasyan, W. Adam, T. Bergauer,
M. Dragicevic, J. Erö, M. Friedl, R. Fruehwirth, V. M. Ghete, et al., “Nuclear effects on the transverse momentum spectra of charged
particles in ppb collisions at $\sqrt{s_{NN}}=5.02$ tev,” The European
Physical Journal C, vol. 75, no. 5, pp. 1–25, 2015.

[32]
C. collaboration et al., “Measurement of inclusive jet production and
nuclear modifications in ppb collisions at $\sqrt{s_{NN}}$= 5.02 tev,” arXiv preprint arXiv:1601.02001, 2016.

[33]
B. Abelev, J. Adam, D. Adamová, A. Adare, M. Aggarwal, G. Aglieri Rinella,
M. Agnello, A. Agocs, A. Agostinelli, Z. Ahammed, et al., “J/$\psi$
production and nuclear effects in p-pb collisions at $\sqrt{s_{NN}} $= 5.02
tev,” Journal of High Energy Physics, vol. 2014, no. 2, pp. 1–26,

[34]
B. Abelev, J. Adam, D. Adamová, M. M. Aggarwal, G. A. Rinella, M. Agnello,
A. Agostinelli, N. Agrawal, Z. Ahammed, N. Ahmad, et al., “Measurement
of prompt d-meson production in p- pb collisions at $\sqrt{s_{NN}}$= 5.02
tev,” Physical review letters, vol. 113, no. 23, p. 232301, 2014.

[35]
C. collaboration et al., “Study of b meson production in ppb collisions
at $\sqrt{s_{NN}}$= 5.02 tev using exclusive hadronic decays,” arXiv
preprint arXiv:1508.06678, 2015.

[36]
S. Chatrchyan, V. Khachatryan, A. M. Sirunyan, A. Tumasyan, W. Adam, E. Aguilo,
T. Bergauer, M. Dragicevic, J. Erö, C. Fabjan, et al.,
“Observation of long-range, near-side angular correlations in ppb collisions
at the lhc,” Physics Letters B, vol. 718, no. 3, pp. 795–814, 2013.

[37]
C. F. Perdrisat, V. Punjabi, and M. Vanderhaeghen, “Nucleon Electromagnetic
Form Factors,” Prog. Part. Nucl. Phys., vol. 59, pp. 694–764, 2007.

[38]
D. Kharzeev, E. Levin, and M. Nardi, “The Onset of classical QCD dynamics in
relativistic heavy ion collisions,” Phys. Rev., vol. C71, p. 054903,

[39]
H.-J. Drescher, A. Dumitru, A. Hayashigaki, and Y. Nara, “The Eccentricity in
heavy-ion collisions from color glass condensate initial conditions,” Phys. Rev., vol. C74, p. 044905, 2006.

[40]
B. Schenke, P. Tribedy, and R. Venugopalan, “Fluctuating Glasma initial
conditions and flow in heavy ion collisions,” Phys. Rev. Lett.,
vol. 108, p. 252301, 2012.

[41]
B. Schenke, P. Tribedy, and R. Venugopalan, “Event-by-event gluon
multiplicity, energy density, and eccentricities in ultrarelativistic
heavy-ion collisions,” Phys. Rev., vol. C86, p. 034908, 2012.

[42]
H. Song and U. W. Heinz, “Causal viscous hydrodynamics in 2+1 dimensions for
relativistic heavy-ion collisions,” Phys. Rev., vol. C77, p. 064901,

[43]
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., vol. C94, no. 2, p. 024907, 2016.

[44]
S. Pratt and G. Torrieri, “Coupling relativistic viscous hydrodynamics to
boltzmann descriptions,” Physical Review C, vol. 82, no. 4, p. 044901,

[45]
M. McNelis, D. Everett, and U. Heinz, “Particlization in fluid dynamical
simulations of heavy-ion collisions: The is3d module,” Computer Physics
Communications, vol. 258, p. 107604, 2021.

[46]
S. A. Bass et al., “Microscopic models for ultrarelativistic heavy ion
collisions,” Prog. Part. Nucl. Phys., vol. 41, pp. 255–369, 1998.
[Prog. Part. Nucl. Phys.41,225(1998)].

[47]
M. Cacciari, M. Greco, and P. Nason, “The P(T) spectrum in heavy flavor
hadroproduction,” JHEP, vol. 05, p. 007, 1998.

[48]
M. Cacciari, S. Frixione, and P. Nason, “The p(T) spectrum in heavy flavor
photoproduction,” JHEP, vol. 03, p. 006, 2001.

[49]
A. Majumder, “Incorporating space-time within medium-modified jet-event
generators,” Physical Review C, vol. 88, no. 1, p. 014909, 2013.

[50]
R. J. Fries, B. Müller, C. Nonaka, and S. Bass, “Hadron production in
heavy ion collisions: Fragmentation and recombination from a dense parton
phase,” Physical Review C, vol. 68, no. 4, p. 044902, 2003.

[51]
J. Putschke, K. Kauder, E. Khalaj, A. Angerami, S. Bass, S. Cao, J. Coleman,
L. Cunqueiro, T. Dai, L. Du, et al., “The jetscape framework,” arXiv preprint arXiv:1903.07706, 2019.

[52]
J. S. Moreland, J. E. Bernhard, and S. A. Bass, “Estimating initial state and
quark-gluon plasma medium properties using a hybrid model with nucleon
substructure calibrated to $p$-Pb and Pb-Pb collisions at
$\sqrt{s_\mathrm{NN}}=5.02$ TeV,” 2018.

[53]
“An introduction to pythia 8.2,” Computer physics communications,
vol. 191, pp. 159–177, 2015.

[54]
B. L. Combridge, “Associated Production of Heavy Flavor States in p p and
anti-p p Interactions: Some QCD Estimates,” Nucl. Phys., vol. B151,
pp. 429–456, 1979.

[55]
D. Everett, W. Ke, J.-F. Paquet, G. Vujanovic, S. Bass, L. Du, C. Gale,
M. Heffernan, U. Heinz, D. Liyanage, et al., “Phenomenological
constraints on the transport properties of qcd matter with data-driven model
averaging,” Physical review letters, vol. 126, no. 24, p. 242301,

[56]
J. S. Moreland and R. A. Soltz, “Hydrodynamic simulations of relativistic
heavy-ion collisions with different lattice quantum chromodynamics
calculations of the equation of state,” Phys. Rev., vol. C93, no. 4,
p. 044913, 2016.

[57]
C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, “The iebe-vishnu
code package for relativistic heavy-ion collisions,” Computer Physics
Communications, vol. 199, pp. 61–85, 2016.

[58]
W. Israel and J. M. Stewart, “Transient relativistic thermodynamics and
kinetic theory,” Annals Phys., vol. 118, pp. 341–372, 1979.

[59]
J. D. Bjorken, “Highly relativistic nucleus-nucleus collisions: The central
rapidity region,” Physical review D, vol. 27, no. 1, p. 140, 1983.

[60]
J. S. Moreland, J. E. Bernhard, and S. A. Bass, “Alternative ansatz to
wounded nucleon and binary collision scaling in high-energy nuclear
collisions,” Phys. Rev., vol. C92, no. 1, p. 011901, 2015.

[61]
S. S. Adler et al., “Centrality dependence of pi0 and eta production at
large transverse momentum in $\sqrt{s_{NN}}$ = 200-GeV d+Au collisions,”
Phys. Rev. Lett., vol. 98, p. 172302, 2007.

[62]
B. Abelev, J. Adam, D. Adamová, M. M. Aggarwal, G. A. Rinella, M. Agnello,
A. Agostinelli, N. Agrawal, Z. Ahammed, N. Ahmad, et al.,
“Multiparticle azimuthal correlations in p-pb and pb-pb collisions at the
cern large hadron collider,” Physical Review C, vol. 90, no. 5,
p. 054901, 2014.

[63]
G. Aad, T. Abajyan, B. Abbott, J. Abdallah, S. A. Khalek, A. Abdelalim,
R. Aben, B. Abi, M. Abolins, O. AbouZeid, et al., “Measurement with
the atlas detector of multi-particle azimuthal correlations in p+ pb
collisions at $\sqrt{s_{NN}}$= 5.02 tev,” Physics Letters B, vol. 725,
no. 1-3, pp. 60–78, 2013.

[64]
G. Aad, B. Abbott, J. Abdallah, S. A. Khalek, R. Aben, B. Abi, M. Abolins,
O. AbouZeid, H. Abramowicz, H. Abreu, et al., “Measurement of
long-range pseudorapidity correlations and azimuthal harmonics in
$\sqrt{s_{NN}}$= 5.02 tev proton-lead collisions with the atlas detector,”
Physical Review C, vol. 90, no. 4, p. 044906, 2014.

[65]
A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, and D. Teaney,
“Matching the nonequilibrium initial stage of heavy ion collisions to
hydrodynamics with qcd kinetic theory,” Physical Review Letters,
vol. 122, no. 12, p. 122302, 2019.

[66]
O. Linnyk, E. Bratkovskaya, V. Ozvenchuk, W. Cassing, and C. Ko, “Dilepton
production in nucleus-nucleus collisions at top sps energy within the
parton-hadron-string dynamics (phsd) transport approach,” arXiv
preprint arXiv:1107.3402, 2011.

[67]
D. Everett et al., “Multi-system Bayesian constraints on the transport
coefficients of QCD matter,” Phys. Rev. C, vol. 103, no. 5,
p. 054904, 2021.

[68]
C. Peterson, D. Schlatter, I. Schmitt, and P. M. Zerwas, “Scaling violations
in inclusive e+ e- annihilation spectra,” Physical Review D, vol. 27,
no. 1, p. 105, 1983.

[69]
S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, “Linearized Boltzmann transport
model for jet propagation in the quark-gluon plasma: Heavy quark
evolution,” Phys. Rev., vol. C94, no. 1, p. 014909, 2016.

[70]
Y. Oh, C. M. Ko, S. H. Lee, and S. Yasui, “Heavy baryon/meson ratios in
relativistic heavy ion collisions,” arXiv preprint arXiv:0901.1382,

[71]
C.-W. Hwang, “Charge radii of light and heavy mesons,” The European
Physical Journal C-Particles and Fields, vol. 23, no. 3, pp. 585–590, 2002.

[72]
Z.-w. Lin, T. G. Di, and C. M. Ko, “Charm meson scattering cross-sections by
pion and rho meson,” Nucl. Phys., vol. A689, pp. 965–979, 2001.

[73]
S. Cao, G.-Y. Qin, and S. A. Bass, “Heavy-quark dynamics and hadronization in
ultrarelativistic heavy-ion collisions: Collisional versus radiative energy
loss,” Phys. Rev., vol. C88, p. 044907, 2013.

[74]
Y. Xu, M. Nahrgang, J. E. Bernhard, S. Cao, and S. A. Bass, “A data-driven
analysis of the heavy quark transport coefficient,” Nucl. Phys.,
vol. A967, pp. 668–671, 2017.

[75]
M. He, R. J. Fries, and R. Rapp, “$\mathbf{D_s}$-Meson as Quantitative Probe
of Diffusion and Hadronization in Nuclear Collisions,” Phys. Rev.
Lett., vol. 110, no. 11, p. 112301, 2013.

[76]
M. He, R. J. Fries, and R. Rapp, “Thermal Relaxation of Charm in Hadronic
Matter,” Phys. Lett., vol. B701, pp. 445–450, 2011.

[77]
M. He, R. J. Fries, and R. Rapp, “Heavy Flavor at the Large Hadron Collider
in a Strong Coupling Approach,” Phys. Lett., vol. B735, pp. 445–450,

[78]
W. Ke, Y. Xu, and S. A. Bass, “Linearized Boltzmann-Langevin model for heavy
quark transport in hot and dense QCD matter,” Phys. Rev., vol. C98,
no. 6, p. 064901, 2018.

[79]
S. Plumari, W. M. Alberico, V. Greco, and C. Ratti, “Recent thermodynamic
results from lattice QCD analyzed within a quasi-particle model,” Phys. Rev., vol. D84, p. 094004, 2011.

[80]
F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco, “Estimating
the charm quark diffusion coefficient and thermalization time from D meson
spectra at energies available at the BNL Relativistic Heavy Ion Collider and
the CERN Large Hadron Collider,” Phys. Rev., vol. C96, no. 4,
p. 044905, 2017.

[81]
J. Uphoff, O. Fochler, Z. Xu, and C. Greiner, “Heavy quark production at RHIC
and LHC within a partonic transport model,” Phys. Rev., vol. C82,
p. 044906, 2010.

[82]
J. Uphoff, O. Fochler, Z. Xu, and C. Greiner, “Elliptic Flow and Energy Loss
of Heavy Quarks in Ultra-Relativistic heavy Ion Collisions,” Phys.
Rev., vol. C84, p. 024908, 2011.

[83]
J. Uphoff, F. Senzel, Z. Xu, and C. Greiner, “Momentum imbalance of D mesons
in ultra-relativistic heavy-ion collisions at LHC,” Phys. Rev.,
vol. C89, no. 6, p. 064906, 2014.

[84]
J. Uphoff, O. Fochler, Z. Xu, and C. Greiner, “Elastic and radiative heavy
quark interactions in ultra-relativistic heavy-ion collisions,” J.
Phys., vol. G42, no. 11, p. 115106, 2015.

[85]
H. Song and U. W. Heinz, “Suppression of elliptic flow in a minimally viscous
quark-gluon plasma,” Phys. Lett., vol. B658, pp. 279–283, 2008.

[86]
T. Song, H. Berrehrah, D. Cabrera, J. M. Torres-Rincon, L. Tolos, W. Cassing,
and E. Bratkovskaya, “Tomography of the Quark-Gluon-Plasma by Charm
Quarks,” Phys. Rev., vol. C92, no. 1, p. 014910, 2015.

[87]
W. A. Horowitz and M. Gyulassy, “Heavy quark jet tomography of Pb + Pb at
LHC: AdS/CFT drag or pQCD energy loss?,” Phys. Lett., vol. B666,
pp. 320–323, 2008.

[88]
W. A. Horowitz, “Testing pQCD and AdS/CFT Energy Loss at RHIC and LHC,”
AIP Conf. Proc., vol. 1441, no. 1, pp. 889–891, 2012.

[89]
Z. Xu and C. Greiner, “Thermalization of gluons in ultrarelativistic heavy
ion collisions by including three-body interactions in a parton cascade,”
Phys. Rev., vol. C71, p. 064901, 2005.

[90]
A. Majumder, “The In-medium scale evolution in jet modification,” 2009.

[91]
A. Majumder, R. J. Fries, and B. Muller, “Photon bremsstrahlung and diffusive
broadening of a hard jet,” Phys. Rev., vol. C77, p. 065209, 2008.

[92]
A. Majumder and B. Muller, “Higher twist jet broadening and classical
propagation,” Phys. Rev., vol. C77, p. 054903, 2008.

[93]
X.-f. Guo and X.-N. Wang, “Multiple scattering, parton energy loss and
modified fragmentation functions in deeply inelastic e A scattering,” Phys. Rev. Lett., vol. 85, pp. 3591–3594, 2000.

[94]
A. Majumder, “Hard collinear gluon radiation and multiple scattering in a
medium,” Phys. Rev., vol. D85, p. 014023, 2012.

[95]
R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff, “The
Landau-Pomeranchuk-Migdal effect in QED,” Nucl. Phys., vol. B478,
pp. 577–597, 1996.

[96]
R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff,
“Radiative energy loss of high-energy quarks and gluons in a finite volume
quark - gluon plasma,” Nucl. Phys., vol. B483, pp. 291–320, 1997.

[97]
R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff,
“Radiative energy loss and p(T) broadening of high-energy partons in
nuclei,” Nucl. Phys., vol. B484, pp. 265–282, 1997.

[98]
M. Gyulassy, P. Levai, and I. Vitev, “Jet quenching in thin quark gluon
plasmas. 1. Formalism,” Nucl. Phys., vol. B571, pp. 197–233, 2000.

[99]
M. Gyulassy, P. Levai, and I. Vitev, “Jet quenching in thin plasmas,” Nucl. Phys., vol. A661, pp. 637–640, 1999.

[100]
M. Gyulassy, P. Levai, and I. Vitev, “Reaction operator approach to
nonAbelian energy loss,” Nucl. Phys., vol. B594, pp. 371–419, 2001.

[101]
M. Gyulassy, P. Levai, and I. Vitev, “NonAbelian energy loss at finite
opacity,” Phys. Rev. Lett., vol. 85, pp. 5535–5538, 2000.

[102]
M. Gyulassy, P. Levai, and I. Vitev, “Jet tomography of Au+Au reactions
including multigluon fluctuations,” Phys. Lett., vol. B538,
pp. 282–288, 2002.

[103]
P. B. Arnold, G. D. Moore, and L. G. Yaffe, “Photon emission from
ultrarelativistic plasmas,” JHEP, vol. 11, p. 057, 2001.

[104]
P. B. Arnold, G. D. Moore, and L. G. Yaffe, “Photon emission from quark gluon
plasma: Complete leading order results,” JHEP, vol. 12, p. 009, 2001.

[105]
P. B. Arnold, G. D. Moore, and L. G. Yaffe, “Photon and gluon emission in
relativistic plasmas,” JHEP, vol. 06, p. 030, 2002.

[106]
P. B. Arnold, G. D. Moore, and L. G. Yaffe, “Effective kinetic theory for
high temperature gauge theories,” JHEP, vol. 01, p. 030, 2003.

[107]
N. Armesto et al., “Comparison of Jet Quenching Formalisms for a
Quark-Gluon Plasma 'Brick',” Phys. Rev., vol. C86, p. 064904, 2012.

[108]
S. A. Bass, C. Gale, A. Majumder, C. Nonaka, G.-Y. Qin, T. Renk, and
J. Ruppert, “Systematic Comparison of Jet Energy-Loss Schemes in a
realistic hydrodynamic medium,” Phys. Rev., vol. C79, p. 024901,

[109]
B.-W. Zhang, E. Wang, and X.-N. Wang, “Heavy quark energy loss in a nuclear
medium,” Physical review letters, vol. 93, no. 7, p. 072301, 2004.

[110]
Y. Xu, J. E. Bernhard, S. A. Bass, M. Nahrgang, and S. Cao, “Data-driven
analysis for the temperature and momentum dependence of the heavy-quark
diffusion coefficient in relativistic heavy-ion collisions,” Physical
Review C, vol. 97, no. 1, p. 014907, 2018.

[111]
W. Ke and X.-N. Wang, “Qgp modification to single inclusive jets in a
calibrated transport model,” Journal of High Energy Physics,
vol. 2021, no. 5, pp. 1–51, 2021.

[112]
J. C. Collins, D. E. Soper, and G. Sterman, “Factorization of hard processes
in qcd,” in Perturbative QCD, pp. 1–91, World Scientific, 1989.

[113]
R. S. Thorne, S. Bailey, T. Cridge, L. A. Harland-Lang, A. Martin, and
R. Nathvani, “Updates of pdfs in the mmht framework,” arXiv preprint
arXiv:1907.08147, 2019.

[114]
S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin,
C. Schmidt, D. Stump, and C.-P. Yuan, “New parton distribution functions
from a global analysis of quantum chromodynamics,” Physical Review D,
vol. 93, no. 3, p. 033006, 2016.

[115]
S. Alekhin, J. Blümlein, S. Moch, and R. Plačakytė, “Parton
distribution functions, $\alpha$ s, and heavy-quark masses for lhc run ii,”
Physical Review D, vol. 96, no. 1, p. 014011, 2017.

[116]
R. D. Field, Applications of perturbative QCD.
Basic Books, 1989.

[117]
S. Cao, C. Sirimanna, and A. Majumder, “The medium modification of
high-virtuality partons,” arXiv preprint arXiv:2101.03681, 2021.

[118]
N. Armesto, L. Cunqueiro, and C. A. Salgado, “Q-pythia: a medium-modified
implementation of final state radiation,” The European Physical Journal
C, vol. 63, no. 4, pp. 679–690, 2009.

[119]
S. Caron-Huot, “Transverse diffusion (“$\hat{q}$”) at weak coupling,”
Nuclear Physics A, vol. 820, no. 1-4, pp. 115c–118c, 2009.

[120]
K. M. Burke, A. Buzzatti, N. Chang, C. Gale, M. Gyulassy, U. Heinz, S. Jeon,
A. Majumder, B. Müller, G.-Y. Qin, et al., “Extracting the jet
transport coefficient from jet quenching in high-energy heavy-ion
collisions,” Physical Review C, vol. 90, no. 1, p. 014909, 2014.

[121]
M. Gyulassy, P. Lévai, and I. Vitev, “Reaction operator approach to
non-abelian energy loss,” Nuclear Physics B, vol. 594, no. 1-2,
pp. 371–419, 2001.

[122]
A. Majumder and C. Shen, “Suppression of the high-p t charged-hadron r a a at
the lhc,” Physical Review Letters, vol. 109, no. 20, p. 202301, 2012.

[123]
X.-F. Chen, T. Hirano, E. Wang, X.-N. Wang, and H. Zhang, “Suppression of
high-p t hadrons in pb+ pb collisions at energies available at the cern large
hadron collider,” Physical Review C, vol. 84, no. 3, p. 034902, 2011.

[124]
B. Schenke, C. Gale, and S. Jeon, “Martini: An event generator for
relativistic heavy-ion collisions,” Physical Review C, vol. 80, no. 5,
p. 054913, 2009.

[125]
G.-Y. Qin, J. Ruppert, C. Gale, S. Jeon, G. D. Moore, and M. G. Mustafa,
“Radiative and collisional jet energy loss in the quark-gluon plasma at the
bnl relativistic heavy ion collider,” Physical review letters,
vol. 100, no. 7, p. 072301, 2008.

[126]
S. Cao, Y. Chen, J. Coleman, J. Mulligan, P. Jacobs, R. Soltz, A. Angerami,
R. Arora, S. Bass, L. Cunqueiro, et al., “Determining the jet
transport coefficient $\hat{q}$ from inclusive hadron suppression
measurements using bayesian parameter estimation,” Physical Review C,
vol. 104, no. 2, p. 024905, 2021.

[127]
Y. Mehtar-Tani, C. A. Salgado, and K. Tywoniuk, “Jets in qcd media: from color
coherence to decoherence,” Physics Letters B, vol. 707, no. 1,
pp. 156–159, 2012.

[128]
J. Casalderrey-Solana and E. Iancu, “Interference effects in medium-induced
gluon radiation,” Journal of High Energy Physics, vol. 2011, no. 8,
pp. 1–49, 2011.

[129]
Y. Mehtar-Tani, C. A. Salgado, and K. Tywoniuk, “Antiangular ordering of gluon
radiation in qcd media,” Physical Review Letters, vol. 106, no. 12,
p. 122002, 2011.

[130]
A. Kumar, A. Majumder, and C. Shen, “Energy and scale dependence of $\hat{q}$
and the “jet puzzle”,” Physical Review C, vol. 101, no. 3,
p. 034908, 2020.

[131]
Y. Tachibana, “Medium response to jet-induced excitation: theory overview,”
Nuclear Physics A, vol. 982, pp. 156–162, 2019.

[132]
T. Luo, “Medium response in jet quenching,” Nuclear Physics A,
vol. 1005, p. 121992, 2021.

[133]
Y. Tachibana, N.-B. Chang, and G.-Y. Qin, “Full jet in quark-gluon plasma with
hydrodynamic medium response,” Physical Review C, vol. 95, no. 4,
p. 044909, 2017.

[134]
W. Chen, S. Cao, T. Luo, L.-G. Pang, and X.-N. Wang, “Effects of jet-induced
medium excitation in $\gamma$-hadron correlation in a+ a collisions,” Physics Letters B, vol. 777, pp. 86–90, 2018.

[135]
Y. Tachibana, J. Collaboration, et al., “Hydrodynamic response to jets
with a source based on causal diffusion,” Nuclear Physics A,
vol. 1005, p. 121920, 2021.

[136]
J. Casalderrey-Solana, J. G. Milhano, D. Pablos, K. Rajagopal, and X. Yao,
“Jet wake from linearized hydrodynamics,” Journal of High Energy
Physics, vol. 2021, no. 5, pp. 1–50, 2021.

[137]
A. Kumar, Y. Tachibana, C. Sirimanna, G. Vujanovic, S. Cao, A. Majumder,
Y. Chen, L. Du, R. Ehlers, D. Everett, et al., “Inclusive jet and
hadron suppression in a multi-stage approach,” arXiv preprint
arXiv:2204.01163, 2022.

[138]
A. Kumar, Y. Tachibana, D. Pablos, C. Sirimanna, R. Fries, A. Majumder,
A. Angerami, S. Bass, S. Cao, Y. Chen, et al., “Jetscape framework: p+
p results,” Physical Review C, vol. 102, no. 5, p. 054906, 2020.

[139]
J. E. Bernhard, Bayesian parameter estimation for relativistic heavy-ion
PhD thesis, Duke U., 2018-04-19.

[140]
S. A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst,
L. Gerland, M. Hofmann, S. Hofmann, J. Konopka, et al., “Microscopic
models for ultrarelativistic heavy ion collisions,” Progress in
Particle and Nuclear Physics, vol. 41, pp. 255–369, 1998.

[141]
A. M. Sirunyan, A. Tumasyan, W. Adam, F. Ambrogi, E. Asilar, T. Bergauer,
J. Brandstetter, E. Brondolin, M. Dragicevic, J. Erö, et al.,
“Nuclear modification factor of d0 mesons in pbpb collisions at snn= 5.02
tev,” Physics Letters B, vol. 782, pp. 474–496, 2018.

[142]
V. Khachatryan, A. M. Sirunyan, A. Tumasyan, W. Adam, E. Asilar, T. Bergauer,
J. Brandstetter, E. Brondolin, M. Dragicevic, J. Erö, et al.,
“Charged-particle nuclear modification factors in pbpb and ppb collisions at
$\sqrt{s_{NN}} = 5.02$ tev,” Journal of High Energy Physics,
vol. 2017, no. 4, pp. 1–45, 2017.

[143]
A. collaboration et al., “Measurement of d0, d+, d*+ and ds+ production
in pb-pb collisions at $\sqrt{s_{NN}}$= 5.02 tev,” Journal of High
Energy Physics, vol. 2018, no. 10, 2018.

[144]
A. Kumar et al., “JETSCAPE framework: $p+p$ results,” Phys. Rev.
C, vol. 102, no. 5, p. 054906, 2020.

[145]
E. Norrbin and T. Sjostrand, “Production and hadronization of heavy quarks,”
The European Physical Journal C-Particles and Fields, vol. 17, no. 1,
pp. 137–161, 2000.

[146]
G. Vujanovic et al., “Multi-stage evolution of heavy quarks in the
quark-gluon plasma,” Nucl. Phys. A, vol. 1005, p. 121965, 2021.

[147]
C. Park, A. Angerami, S. Bass, S. Cao, J. Coleman, L. Cunqueiro, T. Dai, L. Du,
H. Elfner, D. Everett, et al., “Multi-stage jet evolution through qgp
using the jetscape framework: inclusive jets, correlations and leading
hadrons,” arXiv preprint arXiv:1902.05934, 2019.

[148]
K. M. Burke et al., “Extracting the jet transport coefficient from jet
quenching in high-energy heavy-ion collisions,” Phys. Rev. C,
vol. 90, no. 1, p. 014909, 2014.

[149]
D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: The
mcmc hammer,” PASP, vol. 125, pp. 306–312, 2013.

[150]
E. Schulz, M. Speekenbrink, and A. Krause, “A tutorial on gaussian process
regression: Modelling, exploring, and exploiting functions,” Journal of
Mathematical Psychology, vol. 85, pp. 1–16, 2018.

[151]
J.-F. Paquet and S. A. Bass, “Effective viscosities in a hydrodynamically
expanding boost-invariant qcd plasma,” Physical Review C, vol. 102,
no. 1, p. 014903, 2020.

[152]
M. Xie, W. Ke, H. Zhang, and X.-N. Wang, “Information field based global
bayesian inference of the jet transport coefficient,” arXiv preprint
arXiv:2206.01340, 2022.

[153]
Y. Xu, Data-driven Analysis of Heavy Quark Transport in Ultra-relativistic
Heavy-ion Collisions.
PhD thesis, Duke University, 2019.

[154]
F.-L. Liu, W.-J. Xing, X.-Y. Wu, G.-Y. Qin, S. Cao, and X.-N. Wang, “Qlbt: A
linear boltzmann transport model for heavy quarks in a quark-gluon plasma of
quasi-particles,” arXiv preprint arXiv:2107.11713, 2021.

[155]
S. Acharya, D. Adamova, A. Adler, J. Adolfsson, M. M. Aggarwal, G. A. Rinella,
M. Agnello, N. Agrawal, Z. Ahammed, S. Ahmad, et al., “Measurements of
inclusive jet spectra in p p and central pb-pb collisions at $\sqrt{s_{NN}}$=
5.02 tev,” Physical Review C, vol. 101, no. 3, p. 034911, 2020.

[156]
H. Akaike, “A new look at the statistical model identification,” IEEE
transactions on automatic control, vol. 19, no. 6, pp. 716–723, 1974.

[157]
J. Zinn-Justin, Quantum field theory and critical phenomena, vol. 171.
Oxford university press, 2021.

[158]
S. Odaka et al., “Measurement of the running of effective QED coupling
at large momentum transfer in the space - like region,” Phys. Rev.
Lett., vol. 81, p. 2428, 1998.

[159]
P. Braun-Munzinger and J. Wambach, “The Phase Diagram of Strongly-Interacting
Matter,” Rev. Mod. Phys., vol. 81, pp. 1031–1050, 2009.

[160]
K. Rajagopal and F. Wilczek, “The Condensed matter physics of QCD,” in At the frontier of particle physics. Handbook of QCD. Vol. 1-3 (M. Shifman
and B. Ioffe, eds.), pp. 2061–2151, 2000.

[161]
M. A. Stephanov, “QCD phase diagram: An Overview,” PoS, vol. LAT2006,
p. 024, 2006.

[162]
D. Higdon, J. D. McDonnell, N. Schunck, J. Sarich, and S. M. Wild, “A
Bayesian Approach for Parameter Estimation and Prediction using a
Computationally Intensive Model,” J.Phys., vol. G42, no. 3,
p. 034009, 2015.

[163]
D. J. Gross and F. Wilczek, “Ultraviolet Behavior of Nonabelian Gauge
Theories,” Phys. Rev. Lett., vol. 30, pp. 1343–1346, 1973.

[164]
H. D. Politzer, “Reliable Perturbative Results for Strong Interactions?,”
Phys. Rev. Lett., vol. 30, pp. 1346–1349, 1973.

[165]
Y. Akiba et al., “The Hot QCD White Paper: Exploring the Phases of QCD
at RHIC and the LHC,” 2015.

[166]
A. O’Hagan, “Bayesian analysis of computer code outputs: a tutorial,”
Rel.Engin.Sys.Safety, vol. 91, no. 10–11, pp. 1290–1300, 2006.

[167]
C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine
Cambridge, MA: MIT Press, 2006.

[168]
G. L. J. Steve Brooks, Andrew Gelman and X.-L. Meng, Handbook of Markov
Chain Monte Carlo.
Cambridge, MA: Chapman Hall/CRC., 2006.

[169]
O. Philipsen, “Lattice QCD at non-zero temperature and baryon density,” in
Modern perspectives in lattice QCD: Quantum field theory and high
performance computing. Proceedings, International School, 93rd Session, Les
Houches, France, August 3-28, 2009, pp. 273–330, 2010.

[170]
R. Pasechnik and M. Šumbera, “Phenomenological Review on Quark–Gluon
Plasma: Concepts vs. Observations,” Universe, vol. 3, no. 1, p. 7,

[171]
R. Rapp and E. V. Shuryak, “Resolving the anti-baryon production puzzle in
high-energy heavy ion collisions,” Phys. Rev. Lett., vol. 86,
pp. 2980–2983, 2001.

[172]
B. I. Abelev et al., “Energy dependence of pi+-, p and anti-p
transverse momentum spectra for Au+Au collisions at $\sqrt{s_{NN}}$ = 62.4
and 200-GeV,” Phys. Lett., vol. B655, pp. 104–113, 2007.

[173]
F. Riek and R. Rapp, “Quarkonia and Heavy-Quark Relaxation Times in the
Quark-Gluon Plasma,” Phys. Rev., vol. C82, p. 035201, 2010.

[174]
P. Braun-Munzinger, J. Stachel, and C. Wetterich, “Chemical freezeout and the
QCD phase transition temperature,” Phys. Lett., vol. B596,
pp. 61–69, 2004.

[175]
S. S. Adler et al., “Common suppression pattern of eta and pi0 mesons
at high transverse momentum in Au+Au collisions at $\sqrt{s_{NN}}$ =
200-GeV,” Phys. Rev. Lett., vol. 96, p. 202301, 2006.

[176]
S. S. Adler et al., “Single electrons from heavy flavor decays in p+p
collisions at $\sqrt{s_{NN}}$ = 200-GeV,” Phys. Rev. Lett., vol. 96,
p. 032001, 2006.

[177]
S. S. Adler et al., “Nuclear modification of electron spectra and
implications for heavy quark energy loss in Au+Au collisions at
$\sqrt{s_{NN}}$ - 200-GeV,” Phys. Rev. Lett., vol. 96, p. 032301,

[178]
R. Nouicer, “Formation of Dense Partonic Matter in High Energy Heavy-Ion
Collisions: Highlights of RHIC Results,” in Advanced Studies
Institute on Symmetries and Spin (SPIN-Praha-2008) Prague, Czech Republic,
July 20-26, 2008, 2009.

[179]
N. Armesto, C. A. Salgado, and U. A. Wiedemann, “Medium induced gluon
radiation off massive quarks fills the dead cone,” Phys. Rev.,
vol. D69, p. 114003, 2004.

[180]
M. G. Mustafa, “Energy loss of charm quarks in the quark-gluon plasma:
Collisional versus radiative,” Phys. Rev., vol. C72, p. 014905, 2005.

[181]
R. Stock, “The Parton to hadron phase transition observed in Pb + Pb
collisions at 158-GeV per nucleon,” Phys. Lett., vol. B456,
pp. 277–282, 1999.

[182]
K. Adcox et al., “Formation of dense partonic matter in relativistic
nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX
collaboration,” Nucl. Phys., vol. A757, pp. 184–283, 2005.

[183]
Z. Qiu and U. W. Heinz, “Event-by-event shape and flow fluctuations of
relativistic heavy-ion collision fireballs,” Phys. Rev., vol. C84,
p. 024911, 2011.

[184]
Wikipedia contributors, “Jet (particle physics),” 2004.
[Online; accessed 22-July-2004].

[185]
L. L. Frankfurt and M. I. Strikman, “Hard Nuclear Processes and Microscopic
Nuclear Structure,” Phys. Rept., vol. 160, pp. 235–427, 1988.

[186]
K. J. Eskola, H. Paukkunen, and C. A. Salgado, “EPS09: A New Generation of
NLO and LO Nuclear Parton Distribution Functions,” JHEP, vol. 04,
p. 065, 2009.

[187]
P. Romatschke and U. Romatschke, “Viscosity Information from Relativistic
Nuclear Collisions: How Perfect is the Fluid Observed at RHIC?,” Phys.
Rev. Lett., vol. 99, p. 172301, 2007.

[188]
M. He, R. J. Fries, and R. Rapp, “Heavy-Quark Diffusion and Hadronization in
Quark-Gluon Plasma,” Phys. Rev., vol. C86, p. 014903, 2012.

[189]
H. van Hees, M. Mannarelli, V. Greco, and R. Rapp, “Nonperturbative
heavy-quark diffusion in the quark-gluon plasma,” Phys. Rev. Lett.,
vol. 100, p. 192301, 2008.

[190]
A. Beraudo, A. De Pace, M. Monteno, M. Nardi, and F. Prino, “Heavy flavors in
heavy-ion collisions: quenching, flow and correlations,” Eur. Phys.
J., vol. C75, no. 3, p. 121, 2015.

[191]
A. Beraudo, A. De Pace, M. Monteno, M. Nardi, and F. Prino, “Heavy-flavour
production in high-energy d-Au and p-Pb collisions,” JHEP, vol. 03,
p. 123, 2016.

[192]
T. Lang, H. van Hees, J. Steinheimer, and M. Bleicher, “Elliptic flow and
nuclear modification factors of $D$-mesons at FAIR in a Hybrid-Langevin
approach,” 2013.

[193]
E. V. Shuryak, “Quantum Chromodynamics and the Theory of Superdense
Matter,” Phys. Rept., vol. 61, pp. 71–158, 1980.

[194]
D. H. Rischke, “The Quark gluon plasma in equilibrium,” Prog. Part.
Nucl. Phys., vol. 52, pp. 197–296, 2004.

[195]
J.-Y. Ollitrault, “Anisotropy as a signature of transverse collective
flow,” Phys. Rev., vol. D46, pp. 229–245, 1992.

[196]
K. H. Ackermann et al., “Elliptic flow in Au + Au collisions at
$\sqrt{s_{NN}}$ = 130 GeV,” Phys. Rev. Lett., vol. 86, pp. 402–407,

[197]
T. Epelbaum and F. Gelis, “Pressure isotropization in high energy heavy ion
collisions,” Phys. Rev. Lett., vol. 111, p. 232301, 2013.

[198]
F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, “The Color Glass
Condensate,” Ann. Rev. Nucl. Part. Sci., vol. 60, pp. 463–489, 2010.

[199]
Z. Xu and C. Greiner, “Transport rates and momentum isotropization of gluon
matter in ultrarelativistic heavy-ion collisions,” Phys. Rev.,
vol. C76, p. 024911, 2007.

[200]
W. Cassing and E. L. Bratkovskaya, “Parton transport and hadronization from
the dynamical quasiparticle point of view,” Phys. Rev., vol. C78,
p. 034919, 2008.

[201]
W. Cassing and E. L. Bratkovskaya, “Parton-Hadron-String Dynamics: an
off-shell transport approach for relativistic energies,” Nucl. Phys.,
vol. A831, pp. 215–242, 2009.

[202]
H. Niemi, K. J. Eskola, and R. Paatelainen, “Event-by-event fluctuations in a
perturbative QCD + saturation + hydrodynamics model: Determining QCD matter
shear viscosity in ultrarelativistic heavy-ion collisions,” Phys.
Rev., vol. C93, no. 2, p. 024907, 2016.

[203]
M. Alvioli, H. J. Drescher, and M. Strikman, “A Monte Carlo generator of
nucleon configurations in complex nuclei including Nucleon-Nucleon
correlations,” Phys. Lett., vol. B680, pp. 225–230, 2009.

[204]
P. Huovinen and H. Petersen, “Particlization in hybrid models,” Eur.
Phys. J., vol. A48, p. 171, 2012.

[205]
R. J. Fries, B. Muller, C. Nonaka, and S. A. Bass, “Hadronization in heavy
ion collisions: Recombination and fragmentation of partons,” Phys.
Rev. Lett., vol. 90, p. 202303, 2003.

[206]
X.-N. Wang, “Role of multiple mini - jets in high-energy hadronic
reactions,” Phys. Rev., vol. D43, pp. 104–112, 1991.

[207]
M. Bleicher et al., “Relativistic hadron hadron collisions in the
ultrarelativistic quantum molecular dynamics model,” J. Phys.,
vol. G25, pp. 1859–1896, 1999.

[208]
H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz, and W. Soeldner,
“Charmonium properties in hot quenched lattice QCD,” Phys. Rev.,
vol. D86, p. 014509, 2012.

[209]
J. Novak, K. Novak, S. Pratt, J. Vredevoogd, C. Coleman-Smith, and R. Wolpert,
“Determining Fundamental Properties of Matter Created in Ultrarelativistic
Heavy-Ion Collisions,” Phys. Rev., vol. C89, no. 3, p. 034917, 2014.

[210]
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., vol. 114, p. 202301, 2015.

[211]
J. Auvinen, J. E. Bernhard, S. A. Bass, and I. Karpenko, “Investigating the
collision energy dependence of $\eta$/s in the beam energy scan at the BNL
Relativistic Heavy Ion Collider using Bayesian statistics,” Phys.
Rev., vol. C97, no. 4, p. 044905, 2018.

[212]
N. Benitez, “Bayesian photometric redshift estimation,” Astrophys.
J., vol. 536, pp. 571–583, 2000.

[213]
A. Hocker, H. Lacker, S. Laplace, and F. Le Diberder, “A New approach to a
global fit of the CKM matrix,” Eur. Phys. J., vol. C21, pp. 225–259,

[214]
B. P. Abbott et al., “Properties of the Binary Black Hole Merger
GW150914,” Phys. Rev. Lett., vol. 116, no. 24, p. 241102, 2016.

[215]
J. Goodman and J. Weare, “Ensemble samplers with affine invariance,” Comm.App.Math.Comp.Sc., vol. 5, no. 1, pp. 65–80, 2010.

[216]
C. Peterson, D. Schlatter, I. Schmitt, and P. M. Zerwas, “Scaling Violations
in Inclusive $e^+$ $e^-$ Annihilation Spectra,” Phys. Rev., vol. D27,
p. 105, 1983.

[217]
S. K. Das, F. Scardina, S. Plumari, and V. Greco, “Heavy-flavor in-medium
momentum evolution: Langevin versus Boltzmann approach,” Phys. Rev.,
vol. C90, p. 044901, 2014.

[218]
J. Ghiglieri, G. D. Moore, and D. Teaney, “Jet-Medium Interactions at NLO in
a Weakly-Coupled Quark-Gluon Plasma,” JHEP, vol. 03, p. 095, 2016.

[219]
E. Braaten and R. D. Pisarski, “Simple effective Lagrangian for hard thermal
loops,” Phys. Rev., vol. D45, no. 6, p. R1827, 1992.

[220]
G. D. Moore and D. Teaney, “How much do heavy quarks thermalize in a heavy
ion collision?,” Phys. Rev., vol. C71, p. 064904, 2005.

[221]
P. B. Arnold, “Symmetric path integrals for stochastic equations with
multiplicative noise,” Phys. Rev., vol. E61, pp. 6099–6102, 2000.

[222]
J. F. Gunion and G. Bertsch, “HADRONIZATION BY COLOR BREMSSTRAHLUNG,” Phys. Rev., vol. D25, p. 746, 1982.

[223]
M. G. Mustafa, D. Pal, D. K. Srivastava, and M. Thoma, “Radiative energy loss
of heavy quarks in a quark gluon plasma,” Phys. Lett., vol. B428,
pp. 234–240, 1998.

[224]
Y. L. Dokshitzer and D. E. Kharzeev, “Heavy quark colorimetry of QCD
matter,” Phys. Lett., vol. B519, pp. 199–206, 2001.

[225]
R. Abir, C. Greiner, M. Martinez, M. G. Mustafa, and J. Uphoff, “Soft gluon
emission off a heavy quark revisited,” Phys. Rev., vol. D85,
p. 054012, 2012.

[226]
Wikipedia contributors, “Plagiarism — Wikipedia, the free
encyclopedia,” 2004.
[Online; accessed 22-July-2004].

[227]
S. Caron-Huot and G. D. Moore, “Heavy quark diffusion in QCD and N=4 SYM at
next-to-leading order,” JHEP, vol. 02, p. 081, 2008.

[228]
S. Caron-Huot and G. D. Moore, “Heavy quark diffusion in perturbative QCD at
next-to-leading order,” Phys. Rev. Lett., vol. 100, p. 052301, 2008.

[229]
R. Rapp and H. van Hees, “Heavy Quark Diffusion as a Probe of the Quark-Gluon
Plasma,” 2008.

[230]
R. Baier, Y. L. Dokshitzer, A. H. Mueller, and D. Schiff, “Radiative energy
loss of high-energy partons traversing an expanding QCD plasma,” Phys.
Rev., vol. C58, pp. 1706–1713, 1998.

[231]
D. K. Srivastava, R. Chatterjee, and S. A. Bass, “Transport dynamics of
parton interactions in pp collisions at energies available at the CERN Large
Hadron Collider,” Phys. Rev., vol. C97, no. 6, p. 064910, 2018.

[232]
D. K. Srivastava, R. Chatterjee, and S. A. Bass, “Landau-Pomeranchuk-Midgal
effect and charm production in $pp$ collisions at energies available at the
CERN Large Hadron Collider using the parton cascade model,” Phys.
Rev., vol. C98, no. 5, p. 054910, 2018.

[233]
B. B. Abelev et al., “Azimuthal anisotropy of D meson production in
Pb-Pb collisions at $\sqrt{s_{\rm NN}} = 2.76$ TeV,” Phys. Rev.,
vol. C90, no. 3, p. 034904, 2014.

[234]
J. Liu, C. Shen, and U. Heinz, “Pre-equilibrium evolution effects on
heavy-ion collision observables,” Phys. Rev., vol. C91, no. 6,
p. 064906, 2015.
[Erratum: Phys. Rev.C92,no.4,049904(2015)].

[235]
S. K. Das, M. Ruggieri, F. Scardina, S. Plumari, and V. Greco, “Effect of
pre-equilibrium phase on $R_{AA}$ and $v_2$ of heavy quarks in heavy ion
collisions,” J. Phys., vol. G44, no. 9, p. 095102, 2017.

[236]
L. Keegan, A. Kurkela, A. Mazeliauskas, and D. Teaney, “Initial conditions
for hydrodynamics from weakly coupled pre-equilibrium evolution,” JHEP, vol. 08, p. 171, 2016.

[237]
W. van der Schee, P. Romatschke, and S. Pratt, “Fully Dynamical Simulation of
Central Nuclear Collisions,” Phys. Rev. Lett., vol. 111, no. 22,
p. 222302, 2013.

[238]
S. Acharya et al., “Event-shape engineering for the D-meson elliptic
flow in mid-central Pb-Pb collisions at $\sqrt{s_{\rm NN}} =5.02$ TeV,”
JHEP, vol. 02, p. 150, 2019.

[239]
C. A. G. Prado, J. Noronha-Hostler, R. Katz, A. A. P. Suaide, J. Noronha, M. G.
Munhoz, and M. R. Cosentino, “Event-by-event correlations between soft
hadrons and $D^0$ mesons in 5.02 TeV PbPb collisions at the CERN Large Hadron
Collider,” Phys. Rev., vol. C96, no. 6, p. 064903, 2017.

[240]
R. Derradi de Souza, T. Koide, and T. Kodama, “Hydrodynamic Approaches in
Relativistic Heavy Ion Reactions,” Prog. Part. Nucl. Phys., vol. 86,
pp. 35–85, 2016.

[241]
R. Derradi de Souza, J. Takahashi, T. Kodama, and P. Sorensen, “Effects of
initial state fluctuations in the final state elliptic flow measurements
using the NeXSPheRIO model,” Phys. Rev., vol. C85, p. 054909, 2012.

[242]
H. Niemi, G. S. Denicol, H. Holopainen, and P. Huovinen, “Event-by-event
distributions of azimuthal asymmetries in ultrarelativistic heavy-ion
collisions,” Phys. Rev., vol. C87, no. 5, p. 054901, 2013.

[243]
S. Juchem, W. Cassing, and C. Greiner, “Quantum dynamics and thermalization
for out-of-equilibrium $\phi^4$ theory,” Phys. Rev., vol. D69,
p. 025006, 2004.

[244]
S. Juchem, W. Cassing, and C. Greiner, “Nonequilibrium quantum field dynamics
and off-shell transport for phi**4 theory in (2+1)-dimensions,” Nucl.
Phys., vol. A743, pp. 92–126, 2004.

[245]
H. Berrehrah, E. Bratkovskaya, T. Steinert, and W. Cassing, “A dynamical
quasiparticle approach for the QGP bulk and transport properties,” Int. J. Mod. Phys., vol. E25, no. 07, p. 1642003, 2016.

[246]
C. Sasaki and K. Redlich, “Bulk viscosity in quasi particle models,” Phys. Rev., vol. C79, p. 055207, 2009.

[247]
M. Bluhm, B. Kampfer, and K. Redlich, “Bulk and shear viscosities of the
gluon plasma in a quasiparticle description,” Phys. Rev., vol. C84,
p. 025201, 2011.

[248]
S. Cao et al., “Towards the extraction of heavy-quark transport
coefficients in quark-gluon plasma,” 2018.

[249]
A. Beraudo et al., “Extraction of Heavy-Flavor Transport Coefficients
in QCD Matter,” Nucl. Phys., vol. A979, pp. 21–86, 2018.

[250]
V. Ozvenchuk, O. Linnyk, M. I. Gorenstein, E. L. Bratkovskaya, and W. Cassing,
“Shear and bulk viscosities of strongly interacting “infinite”
parton-hadron matter within the parton-hadron-string dynamics transport
approach,” Phys. Rev., vol. C87, no. 6, p. 064903, 2013.

[251]
G. Policastro, D. T. Son, and A. O. Starinets, “The Shear viscosity of
strongly coupled N=4 supersymmetric Yang-Mills plasma,” Phys. Rev.
Lett., vol. 87, p. 081601, 2001.

[252]
P. Kovtun, D. T. Son, and A. O. Starinets, “Viscosity in strongly interacting
quantum field theories from black hole physics,” Phys. Rev. Lett.,
vol. 94, p. 111601, 2005.

[253]
S. Mattiello and W. Cassing, “Shear viscosity of the Quark-Gluon Plasma from
a virial expansion,” Eur. Phys. J., vol. C70, pp. 243–249, 2010.

[254]
S. S. Adler et al., “Centrality dependence of direct photon production
in $\sqrt{s_{NN}}$ = 200-GeV Au + Au collisions,” Phys. Rev. Lett.,
vol. 94, p. 232301, 2005.

[255]
B. de la Cruz, “W, Z and photon production in CMS,” Nucl. Phys.,
vol. A910-911, pp. 139–146, 2013.

[256]
J. Adam et al., “Anisotropic flow of charged particles in Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV,” Phys. Rev. Lett.,
vol. 116, no. 13, p. 132302, 2016.

[257]
K. J. Eskola, H. Paukkunen, and C. A. Salgado, “Nuclear PDFs at NLO - status
report and review of the EPS09 results,” Nucl. Phys., vol. A855,
pp. 150–157, 2011.

[258]
T. Sjostrand, S. Mrenna, and P. Z. Skands, “A Brief Introduction to PYTHIA
8.1,” Comput. Phys. Commun., vol. 178, pp. 852–867, 2008.

[259]
P. Nason, S. Dawson, and R. K. Ellis, “The Total Cross-Section for the
Production of Heavy Quarks in Hadronic Collisions,” Nucl. Phys.,
vol. B303, pp. 607–633, 1988.

[260]
C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass, and U. Heinz, “The
iEBE-VISHNU code package for relativistic heavy-ion collisions,” Comput. Phys. Commun., vol. 199, pp. 61–85, 2016.

[261]
G. S. Denicol, S. Jeon, and C. Gale, “Transport Coefficients of Bulk Viscous
Pressure in the 14-moment approximation,” Phys. Rev., vol. C90,
no. 2, p. 024912, 2014.

[262]
J. Adam et al., “Centrality dependence of the charged-particle
multiplicity density at midrapidity in Pb-Pb collisions at $\sqrt{s_{\rm
NN}}$ = 5.02 TeV,” Phys. Rev. Lett., vol. 116, no. 22, p. 222302,

[263]
J. Stachel, A. Andronic, P. Braun-Munzinger, and K. Redlich, “Confronting LHC
data with the statistical hadronization model,” J. Phys. Conf. Ser.,
vol. 509, p. 012019, 2014.

[264]
N. Jacazio, “Production of identified charged hadrons in Pb–Pb collisions
at $\sqrt{{s}_{NN}}=$ 5.02 TeV,” Nucl. Phys., vol. A967,
pp. 421–424, 2017.

[265]
N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, “A New method for measuring
azimuthal distributions in nucleus-nucleus collisions,” Phys. Rev.,
vol. C63, p. 054906, 2001.

[266]
N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, “Flow analysis from
multiparticle azimuthal correlations,” Phys. Rev., vol. C64,
p. 054901, 2001.

[267]
A. Bilandzic, R. Snellings, and S. Voloshin, “Flow analysis with cumulants:
Direct calculations,” Phys. Rev., vol. C83, p. 044913, 2011.

[268]
S. Chatrchyan et al., “Multiplicity and Transverse Momentum Dependence
of Two- and Four-Particle Correlations in pPb and PbPb Collisions,” Phys. Lett., vol. B724, pp. 213–240, 2013.

[269]
L. Adamczyk et al., “Observation of $D^0$ Meson Nuclear Modifications
in Au+Au Collisions at $\sqrt{s_{NN}}=200$ GeV,” Phys. Rev. Lett.,
vol. 113, no. 14, p. 142301, 2014.
[Erratum: Phys. Rev. Lett.121,no.22,229901(2018)].

[270]
D. Caffarri, “Open heavy-flavour and quarkonium production in Pb-Pb and p-Pb
collisions measured by the ALICE detector at the LHC,” in Proceedings, 51st Rencontres de Moriond on QCD and High Energy Interactions:
La Thuile, Italy, March 19-26, 2016, pp. 273–276, 2016.

[271]
A. Andronic et al., “Heavy-flavour and quarkonium production in the LHC
era: from proton–proton to heavy-ion collisions,” Eur. Phys. J.,
vol. C76, no. 3, p. 107, 2016.

[272]
S. Wicks, W. Horowitz, M. Djordjevic, and M. Gyulassy, “Elastic, inelastic,
and path length fluctuations in jet tomography,” Nucl. Phys.,
vol. A784, pp. 426–442, 2007.

[273]
S. Wicks, W. Horowitz, M. Djordjevic, and M. Gyulassy, “Heavy quark jet
quenching with collisional plus radiative energy loss and path length
fluctuations,” Nucl. Phys., vol. A783, pp. 493–496, 2007.

[274]
P. B. Gossiaux and J. Aichelin, “Towards an understanding of the RHIC single
electron data,” Phys. Rev., vol. C78, p. 014904, 2008.

[275]
P. B. Gossiaux, R. Bierkandt, and J. Aichelin, “Tomography of a quark gluon
plasma at RHIC and LHC energies,” Phys. Rev., vol. C79, p. 044906,

[276]
M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, “Azimuthal
correlations of heavy quarks in Pb + Pb collisions at $\sqrt{s}=2.76$ TeV at
the CERN Large Hadron Collider,” Phys. Rev., vol. C90, no. 2,
p. 024907, 2014.

[277]
Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, “A Multi-phase transport
model for relativistic heavy ion collisions,” Phys. Rev., vol. C72,
p. 064901, 2005.

[278]
S. et al., “Hot and Dense QCD matter, Unravealing the Mysteries of the
Strongly Interacting Quark-Gluon-Plasma,” White Paper (Brookhaven
National lab), 2012.

[279]
C. Collaboration, “$D^{0}$ meson $v_{n}$ harmonics in PbPb collisions at
$5.02~\mathrm{TeV}$,” 2016.

[280]
S. Jaelani, “Measurement of the D-meson Nuclear Modification Factor and
Elliptic Flow in Pb–Pb Collisions at $\sqrt{s_{NN}}=$ 5.02 TeV with ALICE
at the LHC,” Int. J. Mod. Phys. Conf. Ser., vol. 46, p. 1860018,

[281]
T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya, “Charm
production in Pb + Pb collisions at energies available at the CERN Large
Hadron Collider,” Phys. Rev., vol. C93, no. 3, p. 034906, 2016.

[282]
R. Sharma, I. Vitev, and B.-W. Zhang, “Light-cone wave function approach to
open heavy flavor dynamics in QCD matter,” Phys. Rev., vol. C80,
p. 054902, 2009.

[283]
A. Majumder, “Calculating the jet quenching parameter $\hat{q}$ in lattice
gauge theory,” Physical Review C, vol. 87, no. 3, p. 034905, 2013.

[284]
M. Laine and A. Rothkopf, “Towards understanding thermal jet quenching via
lattice simulations,” arXiv preprint arXiv:1310.2413, 2013.

[285]
H. Liu, K. Rajagopal, and U. Wiedemann, “Calculating the jet quenching
parameter from ads,” CFT [hep-ph/0605178] Preprint.

[286]
F.-L. Lin and T. Matsuo, “Jet quenching parameter in medium with chemical
potential from ads/cft,” Physics Letters B, vol. 641, no. 1,
pp. 45–49, 2006.

[287]
R. Abir and A. Majumder, “Drag-induced radiative energy loss from semihard
heavy quarks,” Phys. Rev., vol. C94, no. 5, p. 054902, 2016.

[288]
A. Majumder, “Incorporating space-time within medium-modified jet event
generators,” Phys. Rev., vol. C88, p. 014909, 2013.

[289]
X.-N. Wang and X.-f. Guo, “Multiple parton scattering in nuclei: Parton
energy loss,” Nucl. Phys., vol. A696, pp. 788–832, 2001.

[290]
T. Sjöstrand, “The PYTHIA Event Generator: Past, Present and Future,”
Comput. Phys. Commun., vol. 246, p. 106910, 2020.

[291]
A. Kumar, A. Majumder, and C. Shen, “Energy and scale dependence of $\hat{q}$
and the “JET puzzle”,” Phys. Rev.
C, vol. 101, no. 3, p. 034908, 2020.

[292]
S. Cao et al., “Determining the jet transport coefficient $\hat{q}$
from inclusive hadron suppression measurements using Bayesian parameter
estimation,” Phys. Rev. C, vol. 104, no. 2, p. 024905, 2021.

[293]
S. Cao et al., “Multistage Monte-Carlo simulation of jet modification
in a static medium,” Phys. Rev. C, vol. 96, no. 2, p. 024909, 2017.

[294]
D. Everett et al., “Phenomenological constraints on the transport
properties of QCD matter with data-driven model averaging,” Phys. Rev.
Lett., vol. 126, no. 24, p. 242301, 2021.

[295]
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., vol. 15, no. 11, pp. 1113–1117, 2019.

[296]
G. S. Denicol, H. Niemi, E. Molnar, and D. H. Rischke, “Derivation of
transient relativistic fluid dynamics from the Boltzmann equation,” Phys. Rev., vol. D85, p. 114047, 2012.
[Erratum: Phys. Rev.D91,no.3,039902(2015)].

[297]
W. Israel, “Nonstationary irreversible thermodynamics: A causal relativistic
theory,” Annals of Physics, vol. 100, no. 1–2, pp. 310 – 331, 1976.

[298]
G.-Y. Qin and A. Majumder, “A pQCD-based description of heavy and light
flavor jet quenching,” Phys. Rev. Lett., vol. 105, p. 262301, 2010.

[299]
T. Luo, S. Cao, Y. He, and X.-N. Wang, “Multiple jets and $\gamma$-jet
correlation in high-energy heavy-ion collisions,” Phys. Lett.,
vol. B782, pp. 707–716, 2018.

[300]
Y. He, T. Luo, X.-N. Wang, and Y. Zhu, “Linear Boltzmann Transport for Jet
Propagation in the Quark-Gluon Plasma: Elastic Processes and Medium
Recoil,” Phys. Rev. C, vol. 91, p. 054908, 2015.
[Erratum: Phys.Rev.C 97, 019902 (2018)].

[301]
J. Noronha-Hostler, B. Betz, J. Noronha, and M. Gyulassy, “Event-by-event
hydrodynamics $+$ jet energy loss: A solution to the $R_{AA} \otimes v_2$
puzzle,” Phys. Rev. Lett., vol. 116, no. 25, p. 252301, 2016.

[302]
W.-t. Deng and X.-N. Wang, “Multiple Parton Scattering in Nuclei: Modified
DGLAP Evolution for Fragmentation Functions,” Phys. Rev., vol. C81,
p. 024902, 2010.

[303]
X.-N. Wang, ed., Quark-Gluon Plasma 5.
New Jersey: World Scientific, 2016.

[304]
G. Vujanovic, J.-F. Paquet, C. Shen, G. S. Denicol, S. Jeon, C. Gale, and
U. Heinz, “Exploring the influence of bulk viscosity of QCD on dilepton