Determination of Quark-Gluon-Plasma Parameters from a Global Bayesian Analysis

TL;DR

This paper presents a Bayesian analysis method to quantitatively extract quark-gluon plasma properties from heavy-ion collision data, improving parameter constraints and model calibration.

Contribution

It introduces a novel Bayesian methodology for simultaneous calibration and parameter extraction of quark-gluon plasma properties from collision data.

Findings

Quantitative constraints on the temperature dependence of shear viscosity.

Successful calibration of the collision model to LHC data.

Demonstration of a universal, extensible analysis framework.

Abstract

The quality of data taken at RHIC and LHC as well as the success and sophistication of computational models for the description of ultra-relativistic heavy-ion collisions have advanced to a level that allows for the quantitative extraction of the transport properties of the Quark-Gluon-Plasma. However, the complexity of this task as well as the computational effort associated with it can only be overcome by developing novel methodologies: in this paper we outline such an analysis based on Bayesian Statistics and systematically compare an event-by-event heavy-ion collision model to data from the Large Hadron Collider. We simultaneously probe multiple model parameters including fundamental quark-gluon plasma properties such as the temperature-dependence of the specific shear viscosity $\eta/s$, calibrate the model to optimally reproduce experimental data, and extract quantitative constraints for all parameters simultaneously. The method is universal and easily extensible to other data and collision models.

Equations

$$P({\mathbf{x}}_{\star }|X,Y,{\mathbf{y}}_{\text{exp}})\propto P(X,Y,{\mathbf{y}}_{\text{exp}}|{\mathbf{x}}_{\star })P({\mathbf{x}}_{\star }).$$

$$P(X,Y,{\mathbf{y}}_{\text{exp}}|{\mathbf{x}}_{\star })\propto \exp \left[-\frac{1}{2}(\mathbf{y}({\mathbf{x}}_{\star })-{\mathbf{y}}_{\text{exp}}{)}^{⊺}{\Sigma }_y^{-1}(\mathbf{y}({\mathbf{x}}_{\star })-{\mathbf{y}}_{\text{exp}})\right],$$

$$(\eta /s)(T)=(\eta /s{)}_{\text{min}}+(\eta /s{)}_{\text{slope}}(T-T_c)(T-T_c{)}^{(\eta /s{)}_{\text{curvature}}}$$

Full text

1 Introduction

Relativistic heavy-ion collisions produce a hot, dense phase of strongly-interacting matter commonly known as the quark-gluon plasma (QGP), which rapidly expands and freezes into hadrons [1, 2, 3, 4, 5, 6, 7]. Since the QGP is not directly observable – only final-state hadrons are detected – present research seeks to quantify the fundamental properties of the QGP, such as its transport coefficients and the nature of the initial state, through comparisons of experimental measurements to computational model calculations.

Computational models must take a set of input parameters including the physical properties of interest, simulate the full time-evolution of heavy-ion collisions, and produce outputs analogous to experimental measurements. The true values of the physical properties are then extracted by calibrating the input parameters so that the model output optimally reproduces the experimental data. This generic recipe is called “model-to-data comparison”. Challenges faced by this type of analysis are the amount of computational effort required to scan the parameter space of the model and correlations among the input parameters which may affect multiple observables, so that they cannot be constrained independently.

This article introduces the use of a Bayesian model-to-data analysis for the extraction of QGP properties such as the temperature-dependence of its specific shear viscosity. It is designed to serve as an example of many different possible applications. Due to space constraints, only a broad sketch of the analysis can be given – for a rigorous description we refer the reader to [8, 9].

2 Bayesian Model to Data Comparison

Figure 1 provides a schematic overview of the of the different components of a Baysian model-to-data comparison. Starting point is a computational physics model with a set of model parameters encoding the physics that one wishes to extract from the data. Here, we use a full event-by-event heavy-ion collision model, VISHNU, based on relativistic viscous fluid dynamics with a microscopic hadronic afterburner [10, 11, 12]. The model has 12 input parameters that encode the initial condition, temperature dependent shear- and bulk viscosities and a couple of additional quantities, such as the thermalization time of the QGP and the transition temperature from the hydrodynamic evolution to the microscopic evolution. We calibrate to multiplicity and flow data from the Large Hadron Collider (LHC) – for the sake of simplicity we focus here on data taken by ALICE at 2.76 and 5.02 TeV beam energy respectively [13, 14].

In principle the model can be evaluated on a fine grid of points in its 12-dimensional parameter space. One can then utilize algorithms such as Markov Chain Monte Carlo (MCMC) to rigorously explore the complex high-dimensional parameter space. However, performing the MCMC analysis requires a very large number of model evaluations in parameter space – often thousands or millions, depending on the problem at hand. Heavy-ion collision models may run for several hours, so a direct MCMC approach may require in excess of $10^{14}$ CPU hours and is thus intractable. The situation is exacerbated when studying event-by-event fluctuations as opposed to average quantities: while event-averaged models save computation time by using a smooth initial condition and single hydrodynamic calculation, event-by-event models have realistic, fluctuated initial conditions, each of which requires its own hydrodynamic treatment. Many thousands of complete events are necessary at each point in parameter space to capture event-by-event fluctuations.

These limitations may be overcome through a modern Bayesian method for analyzing computationally expensive models [15, 16, 17]. A set of salient model parameters is chosen for calibration – the set should include any fundamental physical properties of interest – and the model is evaluated at a relatively small number of points that cover the full parameter ranges of interest. The selection of the points at which to evaluate the model is done via a Latin hypercube algorithm that ensures a full and smooth coverage of the entire parameter space. Once these design points are chosen, the full physics model is run at these points and Gaussian process emulators are trained [18] on the observables calculated from the model output. Once trained, the Gaussian process emulators provide a continuous picture of the parameter space. and act as a fast surrogate to the full model: they predicts model output at arbitrary points in parameter space with negligible computational cost. This effectively removes most practical barriers and enables parameter calibration through standard techniques such as MCMC. Emulators have been successfully used to study a wide range of physical systems, including galaxy formation [19] and heavy-ion collisions [20, 21, 22, 8, 9].

The final step in the parameter estimation method is to use MCMC to calibrate the model parameters to optimally reproduce experimental observables, thereby extracting probability distributions for the true values of the parameters. According to Bayes’ theorem, the probability for the true parameters $\mathbf{x}_{\star}$ is

[TABLE]

The left-hand side is the posterior: the probability of $\mathbf{x}_{\star}$ given the design $X$, computed observables $Y$, and experimental data $\mathbf{y}_{\text{exp}}$. On the right-hand side, $P(\mathbf{x}_{\star})$ is the prior probability – encapsulating initial knowledge of $\mathbf{x}_{\star}$ – and $P(X,Y,\mathbf{y}_{\text{exp}}|\mathbf{x}_{\star})$ is the likelihood: the probability of observing $(X,Y,\mathbf{y}_{\text{exp}})$ given a proposal $\mathbf{x}_{\star}$:

[TABLE]

where $\mathbf{y}_{\text{exp}}$ is the experimental data and $\Sigma$ is the covariance matrix, which is the total of experimental statistical and systematic uncertainty, model statistical uncertainty and GP predictive uncertainty. For simplicity, we place a uniform prior on the model parameters, i.e. the prior is constant within the design range and zero outside.

Figure 2 compares simulated observables to experimental data and demonstrates the MCMC calibration. The top two rows have explicit model calculations at each of the 300 design points (chosen by a Latin hypercube algorithm) at which the actual physics model was evaluated. All model parameters vary across their full ranges, leading to the large spread in computed observables. The results of the MCMC calibration can be seen in the bottom two rows: these show emulator predictions of 100 random samples from the posterior distribution. Here, the model has been calibrated to the experiment, so its calculations are clustered tightly around the data – although some uncertainty remains since the samples are drawn from a posterior distribution of finite width. Overall, the calibrated model provides an excellent simultaneous fit to all observables except for the mean $p_{T}$ of anti-protons (and to a lesser extent kaons), which is systematically overpredicted by 10%.

3 Determination of Initial Conditions and the Temperature-Dependence of the Shear Viscosity

The primary result of the analysis is the posterior distribution for the model parameters, Fig. 3. Here, the diagonal elements are marginal distributions for each model parameter (all other parameters integrated out), while the off-diagonals are joint distributions showing correlations among pairs of parameters. Operationally, these are all histograms of MCMC samples and the diagonal represents the individual probability distributions for all parameters of the model. One should note that the posterior distribution provides far more information than a simple least square fit for the extraction of optimum parameter values. Full probability distributions for all parameters as well as their pairwise correlations are given, enabling a rigorous assessment on how meaningful a particular parameter is for the physics model and how well the data actually constrains the values of the parameters.

The first result we wish to highlight is the TRENTo entropy deposition parameter $p$, which has a remarkably narrow distribution peaked at essentially zero with approximate 90% uncertainty of $\pm 0.1$. This implies that initial state entropy deposition is roughly proportional to the geometric mean of participant nuclear thickness functions, $s\sim\sqrt{\tilde{T}_{A}\tilde{T}_{B}}$. This confirms our previous analysis of the TRENTo model which demonstrated that $p\approx 0$ simultaneously produces the correct ratio between initial state ellipticity and triangularity and fits multiplicity distributions for a variety of collision systems [11]. We observe little correlation between $p$ and any other parameters, suggesting that its optimal value is mostly factorized from the rest of the model. Further, recall that the $p$ parameter smoothly interpolates among different classes of initial condition models; Fig. 4 shows an expanded view of the posterior distribution along with the approximate $p$-values for the other models: The EKRT and IP-Glasma models [23, 24] lie squarely in the peak – this helps explain their success – while the KLN and wounded nucleon models are considerably outside.

The shear viscosity parameters $(\eta/s)_{\text{min,slope,curvature}}$ set the temperature dependence of $\eta/s$ according to the ansatz

[TABLE]

for $T>T_{c}$. The full parametrization also includes a constant $(\eta/s)_{\text{hrg}}$ for $T<T_{c}$; this parameter was included in the calibration but yielded an essentially flat posterior distribution, implying that it has little to no effect. This is not surprising, since hadronic viscosity is largely handled by UrQMD, not the hydrodynamic model. Therefore, we omit $(\eta/s)_{\text{hrg}}$ from the posterior distribution visualizations.

Examining the marginal distributions for $\eta/s$ min, slope and curvature, we observe the distribution for $(\eta/s)_{\text{min}}$ (the value of $\eta/s$ at the transition temperature $T_{c}$ = 0.154 GeV) has a narrow peak at 0.06, below the KSS bound $1/4\pi\approx 0.08$ but consistent with it within 90% uncertainty. On the other hand, zero $\eta/s$ is excluded at the 90% level. The slope parameter has a broad peak, although zero slope is excluded, thus confirming a temperature-dependent (non-constant) $\eta/s$ at the 90% level. The curvature parameter is not constrained, but exhibits a strong correlation with the slope.

Figure 5 visualizes the estimated temperature dependence of $\eta/s$, including uncertainty, by inserting the posterior samples for the $\eta/s$ parameters back into Eq. 3. The figure reveals that the uncertainty is smallest at temperatures less than $T\approx{}$225 MeV. We hypothesize that this is the most important temperature range for the present observables at LHC energies – perhaps it is where the system spends most of its time and hence where most anisotropic flow develops, for instance – and thus the data provide a “handle” for $\eta/s$ around 200 MeV. Data at other beam energies and other, more sensitive observables could provide additional handles at different temperatures, enabling a more precise estimate of the temperature dependence of $\eta/s$ and possibly constrain the curvature in addition to the minimum and slope.

4 Other Applications and Outlook

Bayesian model-to-data comparisons have been well established and powerful tools in many different areas of science, e.g. in cosmology, particle physics and climate sciences. In relativistic heavy-ion physics they are fairly novel, having been pioneered by the MADAI collaboration [28, 20, 22, 8]. One application of particular note is the extraction of the QCD equation of state probed in relativistic heavy-ion collisions using these techniques [21]. While the outcome of the analysis – being compatible with Lattice field theory predictions – may not come as a big surprise, it is nevertheless non-trivial, since it confirms the applicability of a Lattice equation of state calculated in the full equilibrium and infinite time limit for the short-lived mesoscopic system that constitutes a relativistic heavy-ion collision. At Quark Matter 2017, roughly half a dozen different presentations utilized Bayesian model-to-data comparisons for the extraction of a variety of different QGP properties – ranging from initial conditions in pA systems [29], the sub-nucleonic structure of the initial state in AA collisions [30], the temperature-dependence of the QGP shear- and bulk-viscosity [31, 32] to the extraction of the heavy-quark transport coefficient [33].

Looking into the future, an analysis framework that can handle arbitrary numbers of inputs and outputs, systematically calculates quantitative constraints on all inputs simultaneously, and quickly evaluates the efficacy of physical models, will be of significant use for the field and enable multiple different quantitative studies and extractions of QGP properties.

The Duke QCD group acknowledges support by grants no. NSF-ACI-1550225 (NSF) and DE-FG02-05ER41367 (DOE). CPU time was provided by the Open Science Grid, supported by DOE and NSF, as well as the DOE funded National Energy Research Scientific Computing Center (NERSC).