A transverse momentum differential global analysis of Heavy Ion Collisions

Govert Nijs, 2, ∗ Wilke van der Schee, † Umut Gürsoy, ‡ and Raimond Snellings 5, § Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Institute for Theoretical Physics and Center for Extreme Matter and Emergent Phenomena, Utrecht University, 3584 CC Utrecht, The Netherlands Theoretical Physics Department, CERN, CH-1211 Genève 23, Switzerland Institute for Gravitational and Subatomic Physics (GRASP), Utrecht University, 3584 CC Utrecht, The Netherlands Nikhef, 1098 XG Amsterdam, The Netherlands

Introduction -Heavy ion collisions (HIC) performed at RHIC and LHC have led to a wealth of data from which the formation of a quark-gluon plasma (QGP) can be inferred. This existence can be deduced by having a model of initial conditions directly after the collision, a hydrodynamic phase with certain transport properties and lastly a hadronic phase of which the results can be compared with experimental results. Even though robust conclusions on the existence of a QGP can be reached from qualitative features of the data, such as the anisotropy of the soft particles or the quenching of high momentum partons, for a quantitative understanding it is paramount to have a careful understanding of all parameters involved in all stages.
Early studies performing such a comprehensive analysis where all parameters in all stages can be varied simultaneously include [1][2][3][4][5][6][7]. This is done in similarity to modelling in cosmology [8], where cosmological parameters have to be inferred from the Cosmic Microwave Background as well as Large Scale Structure analysis. Full simulations themselves are computationally expensive, and hence an emulator trained on few carefully selected design points is used to evaluate the likelihood of parameters using a Markov Chain Monte Carlo (mcmc). This, together with typically flat prior probability distributions, leads to final (Bayesian) posterior distributions for the chosen parameters.
For a precision study of QGP properties the scope of the full model is important, as artificially restricting e.g. the range of initial conditions may pose unphysical restrictions on hydrodynamic transport. In this Letter we will present the widest set of initial conditions studied to-date combined with hydrodynamics with varying second order transport coefficients, containing a total of 21 varying parameters (boldface in this Letter). Most importantly, we perform a global analysis including experimental data with transverse momentum dependence as well as particle identification for spectra and anisotropic flow coefficients for PbPb collisions at 2.76 and 5.02 TeV, together with identified spectra for pPb collisions at 5.02 TeV.
Model -For the initial conditions we use the T R ENTo model parametrization [9,10]. In this model nucleons of Gaussian width w are positioned by a fluctuating Glauber model separated by a distance of at least d min , and all nucleons consist of n c randomly placed constituents having a Gaussian width of v min + χ struct (w − v min ), with v min = 0.2 fm. Nucleons are wounded depending on their overlap such that the cross section matches the proton-proton result. Constituents of wounded nucleons contribute to the left and right thickness functions T A and T B with norm Nγ/n c , where γ fluctuates according to a gamma distribution of width σ fluct √ n c . These function are finally combined to a final parton density by T = ( 1 2 T p A + 1 2 T p B ) 1/p with p a free parameter.
The prehydrodynamic evolution consists of a freestreaming phase lasting for a time τ fs , with the new feature of introducing an effective velocity v fs (see also the discussion).
For the hydrodynamic evolution, we solve the conservation equations for the stress-energy tensor, with the stress-energy tensor given by the hydrodynamic constitutive relation: T µν = ρu µ u ν − (P + Π)∆ µν + π µν , where ∆ µν = g µν − u µ u ν , and we use the mostly minus convention for the metric. Here ρ, P , Π and π µν are the energy density, pressure, bulk viscous pressure and the traceless transverse shear stress respectively. The equations of motion for the bulk pressure Π and the shear stress π µν are given by the 14-moment approximation [11]   keep only the transport coefficients used in [5]: The pressure is given in terms of the energy density by the hybrid HotQCD/HRG equation of state [12,13]. We parameterize the first order transport coefficients η and ζ in terms of the dimensionless ratios η/s and ζ/s. In particular, η/s = a+b(T −T c )(T /T c ) c , with a minimal value a = (η/s) min at T c = 154 MeV, a slope b = (η/s) slope and a curvature c = (η/s) crv . The bulk viscosity ζ/s is described by an unnormalized Cauchy distribution with height (ζ/s) max , width (ζ/s) max and peak temperature (ζ/s) T0 . The second order transport coefficients τ Π , δ ΠΠ , λ Ππ , τ π , δ ππ , φ 7 , τ ππ and λ πΠ are also given in terms of dimensionless ratios. Of these, we fix to the values from kinetic theory [11], with δ = 1/3 − c 2 s , while we vary the shear and bulk relaxation times τ π and τ Π as well as one other second order coefficient τ ππ . We vary these according to the ratios Finally, the hydrodynamic fluid undergoes particlization at a temperature T switch , whereby viscous contributions as well as resonances are included according to the algorithms presented in [13,14]. These hadrons are then evolved using the SMASH hadronic cascade code [15][16][17]. Experimental data -To compare our model to experiment we start with the dataset used in [4]: PbPb charged particle multiplicity dN ch /dη at 2.76 [18] and 5.02 TeV [19], transverse energy dE T /dη at 2.76 TeV [20], identified yields dN/dy and mean p T for pions, kaons and protons at 2.76 TeV [21], integrated anisotropic flow v n {k} for both 2.76 and 5.02 TeV [22] and p T fluctuations [23] at 2.76 TeV. On top of this we added identified transverse momentum spectra using six coarse grained p T -bins separated at (0.5, 0.75, 1.0, 1.4, 1.8, 2.2, 3.0) GeV both for PbPb at 2.76 [21] and pPb at 5.02 TeV [24], anisotropic identified flow coefficients using the same p T bins (statistics allowing) at 2.76 [25] and 5.02 TeV [26]. As in [27] we useṽ n {k} anisotropic flow coefficients for pPb at 5.02 TeV [28] [29] as well as mean p T for pions, kaons and protons at 5.02 TeV [30]. All of these use representative centrality classes, whereby we also specifically included high multiplicity pPb classes for its anisotropic flow coefficients, giving a total of 418 and 96 datapoints for PbPb and pPb collisions respectively. Posterior distribution -In order to estimate the likelihood of all 21 parameters (bold in the model) we ran our full PbPb (pPb) model at 1000 (2000) design points located on a Latin Hypercube in the parameter space using 6k (40k) events per design point (the parameter ranges can be found in the posterior distributions later) [31]. For each system we apply a transformation to 25 principal components (PCs), for which we train Gaussian emulators [13,27,32]. Crucially, the emulator also estimates its own uncertainty (which we validated) and through the PCA this includes correlations among the datapoints. Full details as well as emulator results can be found in our companion paper [33].
Using either PbPb only or both PbPb and pPb emulators we ran a Markov Chain Monte Carlo (mcmc) employing the EMCEE2.2 code [13,27,34], using 600 walkers for approximately 15k steps. This led to the converged posterior distributions in Fig. 1, shown with (solid) and without (dashed) the pPb data. Fig. 2 shows results from 100 random samples of the posterior distribution for a representative selection of our datapoints. In general these compare well, even for p T -differentiated identified v n {2} distributions for both central and peripheral collisions.  For pPb the posterior distributions are significantly wider than the experimental uncertainties, since even for 2000 design points the model is sufficiently complicated that a significant emulator uncertainty remains. It is for this reason that including pPb for the posterior (blue solid versus green dashed in Fig. 1) does not change the probabilities as much as perhaps expected, though for parameters especially sensitive to small and short-lived systems better constraints are obtained (n c , τ fs , w, d min and σ fluct ).
Perhaps the most striking feature in Fig. 1 is that the posterior for the maximum of ζ/s peaks at zero. This is in contrast to previous work [4,35] that prefers a positive bulk viscosity in order to reduce the mean p T . This, however, makes it hard to describe the p T identified spectra ( Fig. 2 (middle,top)).
Given the scope of our 21 parameter model it is perhaps not surprising that constraints on the parameters and in particular the second order transport coefficients are not that strong. We do however see interesting correlations, as shown in Fig. 3. As expected (η/s) min and (η/s) slope are negatively correlated. Perhaps the most interesting correlation is between τ fs and τ π sT/η: it is possible to have a rather long free streaming time, but only if τ π is relatively small. This correlation indeed guarantees FIG. 5. We show an interesting physical correlation between τππ/τπ and (η/s)min with posterior distributions fitted to experimental data (left) or six sets of generated data from random parameter settings (shown in red). All of these show a negative Pearson's correlation averaging -0. 24 6. We show posteriors for (ζ/s)max in variations using the model presented, or limiting to a subset with fewer observables (as in [4], labelled Duke), fewer parameters (as in [4], labelled Duke), or without the pPb system. The new pTdifferential observables are the most significant addition that led to the small bulk viscosity in Fig. 4. the quick applicability of (viscous) hydrodynamics. In the pre-equilibrium stage, the larger τ fs the larger the deviation from hydrodynamics, see also [33]. On the other hand, a shorter shear relaxation time dampens large deviations from viscous hydrodynamics more quickly. In this way, we can interpret the joint constraint on τ π sT/η and τ fs from the posterior distribution as a preference for the fluid to quickly hydrodynamize [36]. Another strong negative correlation is between v fs and χ struct , and indeed for v fs = 1 our χ struct distribution is in agreement with [10]. This highlights the importance of gaining a better understanding of the initial stages of the collision. We also find that n c and χ struct are negatively correlated.
We obtain tight constraints on our prehydrodynamic flow parameter v fs . Indeed the preferred value is close to unity, with perhaps an unlikely option of a small velocity combined with a tiny χ struct . From Fig. 3 we see that either there is a short τ fs with v fs ≈ 0.85, or a longer τ fs with v fs = 1. The first option is consistent with an equation of state that is slightly below the conformal limit; indeed at T = 0.4 GeV the pressure equals 85% of its conformal value [12,13]. In this scenario the fluid is initialized with a bulk pressure close to its ideal hydrodynamic value at an early time (see also [33]). In the second scenario the fluid will start further from equilibrium, but here the influence of v fs on the initial geometry is dominant, as is also clear from the correlation with χ struct .
A closure test is crucial for any model that attempts parameter extractions as ours. For this we extracted posterior distributions for generated 'data' sets at six random points in our parameter space. Apart from confirming the probability distributions this can also lead to physical insights of wider applicability than by just using the true experimental data. An example is shown in Fig 5 showing that often (η/s) min correlates negatively with τ ππ /τ π . The interpretation is that both have similar effects on the elliptic flow. Both the shear viscosity and τ ππ are dissipative, and hence tend to isotropize the plasma, which suggests an explanation for this correlation (see also the correlation between (η/s) min and τ π sT/η in Fig. 3). Discussion -From the posterior distributions it is possible to obtain the 90% confidence limit of the viscosities, as shown in Fig. 4. For lower temperatures the viscosity is consistent with the canonical string theory value of 1/4π [37] (note also that stringy models exist with a lower viscosity [38]). There is a clear tendency for η/s to increase, as expected from the running of the coupling constant. As already shown the bulk viscosity is found to be small, which is consistent with an approximately conformal theory. We also obtain mild constraints on τ π sT/η 7 and τ ππ /τ π 1.5. The value for τ π sT/η compares well with the holographic (4 − log(4) ≈ 2.61, [39]) and weak coupling (5, [11]) results. The τ ππ /τ π value is consistent with the weak coupling result (10/7 ≈ 1.43, [40]) and agrees well with the holographic result (88/35(2 − log(2)) ≈ 1.92, [41]).
We performed several extra mcmc analyses in order to better understand our small bulk viscosity, shown in Fig. 6. There we varied the observables used (all versus those in [4] labelled Duke) and similarly our varying parameters and lastly including pPb or not. Setting both observables, parameters and systems to those in [4] reproduces their bulk viscosity. Including pPb, more parameters and, most importantly, including our p T -differential observables all reduce the bulk viscosity, explaining the result in Fig. 4.
It is still debated whether matter formed in pPb collisions can be described by hydrodynamics [42][43][44]. Indeed, one of our main motivations of this study was to shed light on this question, by computing posterior probabilities with and without pPb collisions. If our framework manages to fit pPb well this gives further evidence for a hydrodynamic picture. In general we agree well with pPb observables, but the mean kaon p T seems to deviate significantly (see also [33]). This can either imply a deviation from the hydrodynamic picture, but given that our pPb observables are much harder to emulate it could also indicate a more advanced analysis within hydrodynamics or a more advanced initial stage is needed.
Our model can be improved in two directions. Firstly, our initial state, prehydrodynamic phase and particlization are not based on a microscopic theory and in particular the transition to hydrodynamics is not smooth [45]. It would be interesting to investigate this point by including a more physically motivated initial stage. Secondly, we were only able to use data that can be reliably estimated using about 6k events, whereas much more sophisticated data is available. Our data includes the widest set available for a global analysis so far, but nevertheless only roughly 20 principal components are non-trivial. This is much more than in e.g. [7,10,13] where up to 8 PCs contain over 99% of the non-trivial information. Nevertheless the question remains whether it is possible to estimate so many parameters with relatively limited experimental data (see also [46] where it is found that using only the limited dataset it is difficult to obtain much stronger constraints than the given prior range). In the future it will hence be important to incorporate more non-trivial data, perhaps using some approximations to reduce computation time (see also [33,47]).