Global fluid fits to identified particle transverse momentum spectra from heavy-ion collisions at the Large Hadron Collider

Transverse momentum spectra of identified particles produced in heavy-ion collisions at the Large Hadron Collider are described with relativistic fluid dynamics. We perform a systematic comparison of experimental data for pions, kaons and protons up to a transverse momentum of 3 GeV$/c$ with calculations using the FluiduM code package to solve the evolution equations of fluid dynamics, the TrENTo model to describe the initial state and the FastReso code to take resonance decays into account. Using data in five centrality classes at the center-of-mass collision energy per nucleon pair $\sqrt{s_\text{NN}}=2.76\,\text{TeV}$, we determine systematically the most likely parameters of our theoretical model including the shear and bulk viscosity to entropy ratios, the initialization time, initial density and freeze-out temperature through a global search and quantify their posterior probability. This is facilitated by the very efficient numerical implementation of FluiduM and FastReso. Based on the most likely model parameters we present predictions for the transverse momentum spectra of multi-strange hadrons as well as identified particle spectra from Pb-Pb collisions at $\sqrt{s_\text{NN}}=5.02\,\text{TeV}$.


I. INTRODUCTION
High-energy heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) produce a fluid consisting of quarks and gluons, the fundamental constituents of Quantum Chromodynamics (QCD) [1][2][3][4]. The produced fluid is particularly interesting because it is described on microscopic level by a renormalizable and fundamental quantum field theory. While first principle calculations of the macroscopic fluid properties are challenging, phenomenological and theoretical studies are motivated by an increasing amount of experimental results. Remarkably, data and models suggest that a fluid dynamic expansion might be behind some of the recent results of collective behaviour in yet smaller proton-nucleus and proton-proton systems [5][6][7]. Alternative descriptions in terms of initialstate physics and medium-less hadron production are also being developed [5,[8][9][10][11], all of which questions the uniqueness of a fluid-like response of the quark-gluon plasma (QGP). Therefore the resolution to the origins of collective behaviour will likely rely on quantitative rather than just qualitative agreement between data and model. To this end, we present a new framework for systematic studies of soft hadronic observables based on up-to-date and efficient modelling of heavy-ion collisions.
We combine the successful initial condition model TrENTo [12], with the recent viscous relativistic fluid dynamics implementation FluiduM [13] and the novel resonance decay procedure FastReso [14]. The mode splitting implemented in FluiduM allows for a very fast evolution with a single event taking mere seconds to compute. In our work we use an equation of state p(T ) based on recent Lattice-QCD calculations [15,16], see [13], and include both shear and bulk viscous corrections in the evolution and particle freeze-out. In addition we use an enlarged set of resonance decays [17][18][19] based on the 2016 edition of the Particle Data Group book [20].
In the absence of precise first principle calculations, the phenomenological description of heavy-ion collisions has a number of open parameters at different stages of the evolution. They can be estimated indirectly from the comparison of simulations to experimental results. Obviously, a too large number of such parameters can limit the precision of this estimate. Moreover, covering a multidimensional parameter space is computationally expensive and requires efficient implementation of the model. In the present work, we determine the specific shear and bulk viscosities of the QGP, as well as the freeze-out temperature T fo , the starting time of a fluid description τ 0 and the initial entropy profile normalization.
Previous multi-observable model-to-data fits focused on the centrality dependence of momentum integrated quantities, like particle multiplicity, mean transverse momentum or flow harmonics [21]. In this work we perform a systematic study of more differential data, namely, transverse momentum spectra with p T < 3 GeV/c of π, K, and p in five centrality classes of Pb-Pb collisions at the center-of-mass energy per nucleon pair √ s NN = 2.76 TeV at the LHC, and compare them with fluid dynamic simulations.
We summarize the details of initial condition, evolution and hadronization procedures in Sec. II. We discuss the fit procedure and determination of its uncertainties in arXiv:1909.10485v1 [hep-ph] 23 Sep 2019 Sec. III. We then provide the best fit results and predictions for additional observables in Sec. IV and Sec. IV B. Finally, we discuss the analysis and future directions in Sec. V.

II. SETUP
In this section we briefly describe the different components of our theoretical model. We start with the time evolution as implemented in FluiduM [13] solving the equations of relativistic fluid dynamics with shear and bulk viscosity and corresponding relaxation times. Subsequently we turn to the initial conditions, specifically the shape of the energy density in the transverse plane for which we use the TrENTo model [12]. Finally, kinetic freeze-out and the implementation of strong resonance decays is done using FastReso [22].

A. Hydrodynamic evolution: FluiduM
To solve the relativistic fluid equations of motion, we use the code package FluiduM [13]. It is based on the theoretical framework of relativistic fluid dynamics with mode expansion [23][24][25], where the fluid dynamic fields are decomposed in terms of a backgroundfluctuation splitting, similar to what is done for example in cosmology. Schematically, we write the fluid fields Φ(τ, r, φ, η) = Φ 0 (τ, r) + Φ 1 (τ, r, φ, η). The nonlinear evolution equations for an azimuthally and Bjorken boost symmetric background Φ 0 (τ, r) are solved first, while azimuthally and rapidity dependent perturbations Φ 1 (τ, r, φ, η) around this are then studied separately. The evolution equations for both the background and the perturbations around them can be implemented with very accurate and highly efficient numerical algorithms [13].
For the present paper we are interested in azimuthally averaged transverse momentum spectra of identified particles in the mid-rapidity region and do not consider azimuthally and rapidity-dependent perturbations. Neglecting terms that are of quadratic or higher order in perturbation amplitudes, we need only the background solution to the fluid evolution equations as calculated from FluiduM. The corresponding equations of motion have been analyzed from a mathematical perspective, and with an emphasis on their causality structure in ref. [26].
We note here that the current implementation of FluiduM features a flow at vanishing net baryon number chemical potential, based on a state-of-the-art thermodynamic equation of state [15,16], as well as shear and bulk viscous dissipation. For the present paper we assume the shear viscosity to entropy ratio η/s to be independent of temperature. The bulk viscosity to entropy ratio ζ/s is taken to be temperature dependent, however. Specifi-cally, we assume it to be of the Lorentzian form with the peak temperature T peak = 175 MeV and ∆T = 24 MeV [27]. The maximum value (ζ/s) max is taken as a fit parameter. Shear and bulk relaxation times are assumed to be determined by the relations [28] τ shear η/( + p) where is the energy density, p is the pressure, c s is the (temperature dependent) velocity of sound, and a = 0.1 fm/c is a small offset such that a causal evolution of the radial expansion is indeed ensured [26]. For more details on the implementation we refer to [13].

B. Initial conditions: TrENTo
In general terms, a characterization of the initial conditions for Israel-Stewart type fluid dynamics with azimuthal rotation and longitudinal boost symmetry as used for the background in FluiduM consists of the temperature T , radial fluid velocity u r , two independent components of shear stress π φ φ and π η η as well as bulk viscous pressure π bulk on some initial Cauchy surface, such as τ = τ 0 . In the present work we neglect initial radial flow and assume initially π φ φ = π η η = π bulk = 0. This choice is respecting relativistic causality [26].
The shape of the initial entropy density distribution in the transverse plane (which determines the temperature through the thermodynamic equation of state) is taken from the initial state model TrENTo [12], with an overall normalization factor that we take as a fit parameter. The parameters of TrENTo have been taken as in ref. [12], in particular we selected the reduced thickness parameter p = 0, the fluctuation parameter k = 1.4, the nucleon width σ = 0.6 fm and the inelastic nucleon-nucleon cross section σ NN inel = 6.4 fm 2 . Using this set of parameters we have generated the transverse density T R (x, y) for 10 7 events with impact parameter sampled from the range b ∈ [0 fm, 20 fm] and randomized event plane angle. It is convenient to shift the events in the transverse plane such that d 2 x T R ( x) = 0.
The integrated transverse density d 2 xT R ( x) is expected to be monotonously related to the total final charged particle multiplicity, therefore we used this quantity to divide the generated events into narrow multiplicity classes of one percent. Each of these centrality classes can be seen as an ensemble of events with random orientation in the transverse plane.
For each centrality class we calculate the averaged or expected entropy density profile as We introduce here a normalization constant Norm i for each centrality class i. Ideally, the initial state model should take care of centrality dependence and all centrality classes would have one identical normalization, however we observe some tension with the data for the choice p = 0 which we lift by allowing the normalization to be centrality class dependent. We have taken out the initialization time τ 0 to already take into account the main effect of the longitudinal expansion (Bjorken flow) at early times. The initial temperature as a function of radius is then obtained using the equation of state.
While it is convenient for the theoretical description to work with rather narrow centrality classes, they are typically somewhat larger in the experimental results. There are now two possible strategies to deal with this. The first would be to calculate particle spectra for each of the narrow classes and to combine (average) them in a convenient way afterwards. The second strategy is to produce averaged entropy densities for the larger centrality classes by averaging the corresponding distributions from the more narrow classes and to propagate those. The difference in experimental observables between both procedures can be taken as an estimate for the importance of fluctuations. We have compared both strategies and found the difference for transverse momentum spectra to be rather small, of the order of 1% for central collisions. Because of the advantage with respect to computational costs, we follow therefore mainly the second strategy.
C. Freeze-out and resonance decays: FastReso As the system cools down and dilutes, it crosses from a quark-gluon plasma to a fluid dominated by hadronic degrees of freedom. The fluid dynamic description of the latter breaks down eventually, because particle scatterings are no longer efficient in maintaining (first chemical and then kinetic) equilibrium. This necessitates the conversion of fluid fields, such as temperature and flow velocity, to the distribution of hadronic degrees of freedom.
The dynamics of hadronization is not completely understood, but lattice QCD calculations show that below the QCD pseudo-critical temperature T pc = 156 ± 1.5 MeV [29,30], color neutral hadrons become the dominant degrees of freedom of the plasma. In particular the equation of state approaches that of a hadron resonance gas (HRG) [18].
Around or somewhat below 155 MeV in temperature, fluid fields are customary converted to particle distributions using Cooper-Frye procedure [31]. The spectrum of hadron species a on the freeze-out hypersurface Σ is given by the following integral where ν a is the degeneracy factor of spin/polarization states and f a is a particle distribution function, which, in addition to particle energy in fluid rest-frameĒ p ≡ −u ν p ν , may also depend on the local temperature T (x), fluid velocity u µ (x), chemical potential µ(x), viscous shear-stress π µν (x) and bulk viscous pressure π bulk (x). Chemical freeze-out takes place when particle species changing interactions are no longer able to keep up with the expansion rate. However, in practice, a simpler criterion based solely on the freeze-out temperature is used and the freeze-out surface Σ is assumed to be a surface of constant temperature T = T fo . One sometimes includes after chemical freeze-out and before kinetic freeze-out a phase described by fluid dynamics but for a fluid in partial chemical equilibrium. Such a fluid is governed by a number of conservation laws in addition to the ones for energy and momentum. We have implemented this in our theoretical model but found eventually that the improvement of transverse momentum spectra of the studied particles species is not significant. For this reason we use in the present work a simpler prescription with only a single, chemical and kinetic freeze-out.
On the freeze-out surface we take the particle distribution function to be given by the equilibrium Bose-Einstein or Fermi-Dirac distribution (depending on the species), modified by additional corrections due to bulk and shear viscous dissipation, For the viscous corrections we use the commonly employed parametrizations [32,33] δf δf shear = f eq (1 ± f eq ) π ρν p ρ p ν 2( + p)T 2 .
Here m is the mass of the primary resonance.
After chemical freeze-out the populations of unstable resonances decrease as a consequence of their decays and feed the spectra of long lived particles. This large modification of the pion, kaon and proton spectra can be calculated by decaying all (sufficiently unstable) resonances. An efficient procedure to calculate these direct decays was recently introduced by some of us in ref. [14]. The main idea is to apply the decay maps to the primary distributions in Eq. (4) before doing the surface integral. The resulting distribution function of decay products can be decomposed into irreducible components (with respect to rotations in the fluid rest frame) that are pre-computed and stored [22]. Furthermore, for the case of azimuthally symmetric and boost-invariant surface, the freeze-out integrals over space-time rapidity and azimuthal angle can also be pre-computed. Parametrizing the remaining 1+1 dimensional freeze-out surface in radial coordinates by (τ (α), r(α)) where α ∈ (0, 1) is some parameter, the Cooper-Frye freeze-out integral simplifies to one-dimensional integral over α, Here K eq i (p T , u r ), K shear i (p T , u r ) and K bulk i (p T , u r ) are rapidity and azimuthal angle integrated decay kernels [14]. The kernels have implicit dependence on scalars like freeze-out temperature or decay constants which do not vary on the freeze-out surface. The spectra of pions, kaons and protons as calculated with Eq. (8) can then be compared to the experimentally measured p T differential spectra of identified hadrons.
For the calculation of freeze-out kernels in Eq. (8), we use the publicly available code FastReso to perform strong and electroweak decays of unstable hadrons up to mass m ≈ 3 GeV 1 . We use the list of ∼ 700 resonances from refs. [17][18][19], which is based on all listed states (also less well established * states) in the Particle Data Group 2016 publication [20]. This is approximately twice the number of resonances used previously [14]. To perform a scan in chemical freeze-out temperature, we varied it in the range T chem ∈ [130, 0.180] MeV with 0.5 MeV increments and zero baryon chemical potential. The transverse momentum p T (in GeV) and the radial fluid velocity u r have been discretized each on a 81 point non-linear grid in the range of [0, 3.5].

A. Global fit procedure
To summarize, our theoretical description has currently the free parameters η/s, (ζ/s) max , the initialization time τ 0 , the freeze-out temperature T fo and they are assumed to be independent of the centrality class. In addition, we have the normalization constants Norm i for the initial entropy profile which depend on the centrality classes. In order to find the most likely model parameters, in this work we aim at fitting the p T -differential spectra of pions, kaons, and protons in five centrality intervals: 0-5%, 5-10%, 10-20%, 20-30% and 30-40% for Pb-Pb collisions at √ s NN = 2.76 TeV measured by the ALICE Collaboration [34]. We choose to restrict to the soft particle momentum range p T < 3 GeV/c, a region which is believed to be described by a fluid dynamic approximation to QCD dynamics. It is sensitive to radial flow, the viscous transport coefficients and the initial conditions of the plasma [35][36][37][38]. Nine model parameters are left free to vary simultaneously in specific intervals (see Table I), in which the physical values are expected to be located based on physical considerations and previous work [12,38,39]. Of course it is important to check a posteriori that the best fit values are indeed inside these intervals.  In order to determine which combination of the parameters provides the best description of the experimental data we search for the global minimum of where x i is the experimental value of the transverse momentum spectrum at some p T bin for a particular hadron species and centrality class, y i is the corresponding model prediction (for a given set of model parameters) and σ i = σ 2 i,sys + σ 2 i,stat is the sum (in quadrature) of the systematic and statistical uncertainties of the corresponding experimental data point. The sum goes over the five centrality classes (0-5%, 5-10%, 10-20%, 20-30% and 30-40%), three particle species (π, K, p) and the number of p T intervals in the fit range up to 3 GeV/c (N π pT = 41, N K pT = 36 and N p pT = 34). The total number of degrees of freedom is accordingly N dof = 555 − 9.
Note that the degree of correlation in the systematic uncertainties as a function of p T in the experimental measurements is not reported and we do not consider such correlations in the fit.
Furthermore, the simulations themselves might have considerable systematic uncertainties. For example, our model assumes a rather simple freeze-out picture without a detailed modeling of hadronic scatterings and the dissipative corrections to the single-particle distribution functions on the freeze-out surface are arguably somewhat uncertain. Also, our model neglects currently a possible temperature dependence of the shear viscosity to entropy ratio. Independently from this, also completely new physics might affect the experimental data in the low transverse momentum regime, for example pion condensation, see [40,41] and references therein. It is hard to predict a priori the change to model results due to such effects but for the interpretation of results it is important to keep in mind that theoretical uncertainties exist.
It remains to find the values of the model parameters which correspond to the fit of the experimental measurements with minimum χ 2 . We have explored here different strategies. What works best eventually is to discretize the model parameters on a hypercubic lattice and to use numerical interpolation between the lattice points. This allows to determine the χ 2 landscape systematically.
We discretized the parameter ranges by 10 equidistant values for each parameter, which correspond to 10 5 different model calculations for each centrality class. Let us note here that, thanks to streamlined fluid dynamic evolution and resonance decay procedures in our framework, one model simulation for a particular set of parameters takes only a few tens of seconds per centrality on a single core and even in the exhaustive search with 10 5 simulations, the entire fit can be performed with a rate of 1 day/centrality class using a ∼ 100 core machine.
Once all 10 5 simulations have been computed we use an order-5 spline interpolation between them and apply a numerical minimization technique to find the lowest value of χ 2 and its position. For this minimization we used a gradient descent algorithm, to find a local minimum, given a starting point. The global minimum is found by repeating the local minimization for a set of random starting points -in practise we used 100 randomly generated starting points. The best fit found in this way gives χ 2 /N dof = 1.40. As a check of the numerical interpolation, we have also calculated χ 2 directly for this configuration with a compatible result.
The best fit parameters obtained in this way are reported in Table III.

B. Uncertainties and correlations of model parameters
In order to study correlations between pairs of model parameters according to their posterior probability distribution we use two methods. Firstly, two-dimensional slices of the nine dimensional χ 2 landscape are computed with the remaining parameters kept at their global best fit (minimum χ 2 ) value. The results are shown in Fig. 1.
The first four panels on the top show the correlations of the initial entropy profile normalization Norm i with respect to the other four parameters. As an example we report the normalization for the centrality interval 0-5%, Norm 1 . Thanks to the factored out scaling with initialization time τ 0 in Eq. (3), the different Norm i are observed to be almost independent from the other parameters. Only a rather weak correlation is observed between the initial normalization and η/s as well as τ 0 . This could be due to the combined effects of viscous entropy production at early times and the delayed generation of radial flow for larger τ 0 values. In the other six panels of Fig. 1 the correlations of the remaining parameter pairs are shown. We see positive correlation between (ζ/s) max and η/s, between τ 0 and η/s as well as between τ 0 and T fo . On the other hand, negative correlations are instead observed between T fo and η/s as well as between τ 0 and (ζ/s) max . Finally no strong and clear correlation is observed between T fo and (ζ/s) max .
In order to quantify the information that is visually available in Fig. 1, we also determine numerically the form of the χ 2 as a function of the nine model parameters in the vicinity of the minimum. In terms of deviations from the best fit value ∆r = (∆Norm i , ∆τ 0 /(fm/c), ∆(η/s), ∆(ζ/s), ∆T fo /MeV) we find The diagonal values of the inverse matrix A −1 can formally be understood as variances in the fit parameters in a Gaussian approximation to the posterior probability ∼ e −χ 2 /2 , Moreover, the matrix quantifies correlations between the fitted parameters, again in a Gaussian approximation to the posterior distribution. Note that these considerations assume that the uncertainties that enter Eq. (9) are independent and normally distributed. We show the resulting matrix ρ ij in Table II. The uncertainties on the model parameters according to eq. (11) are shown in Table III as uncertainties from the χ 2 landscape. One remarks here that the latter are actually rather small. On the one side, this illustrates nicely that the experimental data are of high quality and have a high power to constrain theoretical models. On the other side, some remarks of caution about a too straight-forward interpretation are in order: It is known that estimating model parameters including their uncertainty is difficult for situations with large N dof and when the minimum χ 2 deviates substantially from its statistical expectation value (for a complete theoretical model) χ 2 = N dof . This problem arises indeed for us when we attempt a global fit for the full range of transverse momenta and all five centrality classes with a single set of parameters. As a characteristic one may calculate the "goodness of fit" Q = 1 − F χ 2 (χ 2 , N dof ) where F χ 2 (x, ν) gives the cumulative χ 2 distribution function with ν degrees of freedom. For the full global fit we find Q = 1.7 × 10 −9 , which is indeed very small. This can be understood as the probability for the observed minimal χ 2 given the data are correctly described by the model and all deviations from it arise indeed due to independent Gaussian fluctuations of the experimental data points. The fact that the goodness of fit is so small means that the theoretical model as it is currently implemented is in fact incomplete. As we will see in the next section, the situation is not as bad, and our fluid model is certainly competitive with other attempts for theoretical descriptions, at least by visual inspection. Certain physics features might be missing, specifically in the low transverse momentum region for pions. Nevertheless, we should take the experimental data and the goodness of fit seriously. This leaves us with the problem to estimate the uncertainty of model parameters for an incomplete theoretical model.

C. Estimation of systematic uncertainties
In order to quantify how well the model parameters of the fluid description can actually be constrained from transverse momentum spectra, we cannot rely purely on the fit uncertainties, which are unrealistically small. As discussed above, the underlying reason is that the theoretical model is not complete. This can be seen directly from the goodness of fit estimate, but also indirectly from the fact that the outcome for the most likely model parameters depends on how the fits are being done in detail. In this subsection we will discuss this latter point, and estimate systematic uncertainties of the model parameters through variations of the fitting scheme.
The first check consists in fitting the five centrality classes separately and estimating the model parameters as a function of centrality. In addition to quantifying uncertainties, this test might also reflect possible temperature dependence of transport coefficients (specifically η/s). On the left hand side of Fig. 2 we show the result for the most likely model parameters when they are determined separately for the different centrality classes (full circles). The error bars illustrate the corresponding uncertainties according to eq. (11), determined from the χ 2 landscape. One finds that the variation arising from the centrality dependence is somewhat larger than the Best fit parameters and their uncertainties determined from the χ 2 landscape through eq. (11), and from the variation of the fitting procedure as reported in Fig. 2. For the global fit we find χ 2 /N dof = 1.40.
calculated fit uncertainties.
In a similar way, we also perform the fit restricting to kaons (without pions and protons) or to kaons and pions (without protons). This is done globally with respect to centrality. The results for the most likely model parameters obtained in this way are shown on the right hand side of Fig. 2 (open stars). In Fig. 2 the red lines represent the values obtained from the global fit reported in Table  III. From the variations shown in Fig. 2 we determine systematic uncertainties of the fitted model parameters and report them in Table III in the right-most column. Specifically, we take this uncertainty to be the maximal deviation seen in Fig. 2 from the best fit parameter.

IV. RESULTS AND DISCUSSION
The final step in the modeling workflow is to compute observables with best fit parameters, Table III, and to make predictions for observables not used in the fit. In this work, the simulations are performed for Pb-Pb collisions at √ s NN = 2.76 TeV for the centrality intervals 0-5%, 5-10%, 10-20%, 20-30% and 30-40%. First, we compare the fitted p T -differential spectra of identified hadrons to experimental measurements. Then we study the derived quantities, like the total multiplicities and mean-p T for different hadron species. Finally, we make model calculations for observables not used in the fit. Namely, we compute the p T spectra for strange and multi-strange baryons at the same collision energy and centrality classes at √ s NN = 2.76 TeV and we make predictions for the pion, kaons, and protons spectra in Pb-Pb collisions at √ s NN = 5.02 TeV.  Table III. A. Fitted particle spectra of π, K, p In Fig. 3 (top panels) we show the transverse momentum differential spectra of identified light hadrons π, K and p using our best fit parameters (χ 2 /N dof = 1.40) listed in Table III (lines) and we compare our results with the ALICE measurement (symbols). The bottom panels show the data to model ratio with shaded areas representing the combined experimental uncertainties.
Simulations are in overall good quantitative agreement with the experimental measurements. Kaon and proton spectra are reproduced within 10%-20% accuracy and within 3σ of experimental errors from the data or the entirep T < 3 (GeV/c) momentum range in all the centrality classes. The pion spectra is reproduced well in a narrower 0.5 < p T < 2.5 (GeV/c) momentum range, while low-p T pions are systematically underpredicted and make major contributions to the relatively large χ 2 /N dof = 1.40 in the fit. We checked that exclud- ing soft pions from the fit results in significantly smaller χ 2 /N dof 0.6 and the minimum moves out from the parameter ranges given in Table I. Such discrepancies in the pion spectra are well known and have been observed both in hydrodynamic simulations [38,39,42,43] and blast-wave fits with resonance decays [44,45].
The enhancement of low-p T pion spectra is typically attributed to the feed-down of resonance decays [46]. However, even after we included a considerably larger set of primary resonances [17][18][19], the agreement of soft pion spectra improved only marginally. Additional physics effects like finite width of resonance decays [47], the presence of pion condensation in heavy-ion collisions [40,41] or going beyond linearised viscous corrections to the freeze-out spectra [43] are being studied.
We would like to note here that our simulations show flat data to model ratio for protons for the considered momentum range and in all centrality classes. In studies simulating hadronic phase after the chemical freezeout [38,39], which is not included in our framework, protons receive additional boost, resulting in a harder spectrum. At early stages of this work, when using a single normalization for initial entropy profile, we indeed observed different proton slope with momenta compared to the experimental measurements. However by allowing centrality dependent normalizations, which exhibit a linear centrality dependence (see Fig. 2), the agreement for proton spectra improved. This indicates that some of the effects attributed to hadronic rescatering might be also reproduced by the change in initial conditions.
In addition to the data to model comparison of partice spectra, we can compute other derived observables: particle multiplicity and mean p T . In the top panel of Fig. 4 we compare our results of total charged and identified particle multiplicities at mid-rapidity as a function of collision centrality for pions, kaons, and protons with the ALICE measurements [34]. Our simulations give a reasonably good description of the centrality dependence of the charged hadron multiplicity. However, also in this case we see a tension with the pion and total charged hadron yields, especially in most central collisions, which is a clear reflection of the underestimation of the low p T pion spectra observed in Fig. 3. In the bottom panel of Fig. 4 we compare the mean transverse momentum p T , for pions, kaons, and protons as a function of centrality between our simulations and the experiment [34]. While p T of kaons agrees very well with the experimental measurements, the p T of pions and protons is systematically different. We note that similar discrepancies are observed in other hydrodynamic simulations [21,39,43] and none appears able to reproduce data within the very small experimental uncertainties.
To our best knowledge no recent heavy-ion simulations are able to produce a uniformly good description of identified particle spectra in central to mid-central nucleusnucleus collisions if experimental uncertainties are taken seriously. The pioneering studies of [49] showed excellent agreement of identified particle spectra for ideal hydrodynamic simulations, but the agreement worsened when effects of viscosity were included. In the EKRT model [50], pion spectra are described well at the expense of overpredicted kaon and proton yields. In Ref. [39] where the effect of both bulk viscosity and hadronic rescattering were studied, the data to model agreement is arguably on the same level as in our work, although we employ a single freeze-out approximation. We note here that the extensive Bayesian analyses of refs. [21,51] have concentrated on momentum integrated observables. In summary, the excellent quality of experimental data of identified particle spectra indicates the need of including additional physics in hydrodynamic simulations of heavy-ion collisions.
B. Strange, multi-strange and energy dependence of particle spectra The benefit of finding the optimal parameters of our hybrid model is that many other observables, not used in the fit, can be directly predicted. This is an important step in validating the physics picture behind the model. Therefore we use the fluid dynamic evolution with the best fit parameters to compute the p T spectra of strange and multi-strange hadrons (Λ, Ξ, Ω) and compare the results with the ALICE measurements [52,53].
The comparison can be seen in Fig. 5 for the 10-20% (left panel) and 20-40% (right panel) centrality intervals. From the comparison one can see that if the value of T fo ≈ 136.9 MeV is kept the same as obtained from the best fit (solid lines), the experimental p T -differential spectra of strange and multi-strange baryons are underestimated by the simulation. This effect is more pronounced for the Λ baryons, which shows a ∼ 40% − 50% discrepancy, while for the Ξ and Ω the simulation and data tend to agree for p T > 2 GeV/c. In previous work [39] it was observed that strange and multi-strange baryons are much more sensitive to the freeze-out temperature than pions, kaons and protons. Indeed, if we increase the value of T fo to 145 MeV, while keeping all other parameters fixed, the simulation shows better agreement with data at low momentum, see in Fig. 5 (dashed lines), but Ξ is then over-predicted for p T > 2 GeV/c. The tendency of strange and multi-strange hadrons preferring higher freeze-out temperatures [54], is sometimes used as an evidence for the scenario of sequential hadronization where the switching from quark to hadron degrees of freedom occurs at different temperatures for different particle flavours [55,56]. However, one should not discount the possibility that additional resonance feed-down might improve the agreement with data. Indeed, by approximately doubling the list of primary hadrons [17][18][19], we observed a significant feed-down for Λ baryons. Further extensions of decay channels and global fits including the strange particles would certainly reduce the apparent discrepancy.
Finally, we can use our model to make predictions for the p T -differential spectra of pions, kaons and protons in Pb-Pb collisions at √ s NN = 5.02 TeV. At higher collision energies, nuclei have more energy to deposit in the collision area, which ultimately results in an increased final particle multiplicity and higher initial QGP energy density. However as the increase of multiplicity is fractional, the fundamental properties of the QGP are not expected to change substantially and we can use the same best fit model parameters to predict particle spectra at higher energies. The only change made is the overall normalization Norm i of the initial entropy density profile. The normalization at √ s NN = 5.02 TeV Pb-Pb collisions is fixed by doing a fit to the published unidentified charged hadron multiplicity as a function of the collision central-  ity [57] 2 . We report the result in Fig. 6, together with the model calculations for integrated yields of pions, kaons and protons as a function of centrality. The corresponding plots for the p T -differential spectra of pions, kaons and protons in Pb-Pb collisions at √ s NN = 5.02 TeV are reported in Fig. 7 for the centrality intervals 0-5%, 5-10%, 10-20%, 20-30% and 30-40%. The p T -spectra at √ s NN = 5.02 TeV are higher and flatter than the ones at √ s NN = 2.76 TeV, which illustrates that stronger radial flow has been developed in the systems with larger final multiplicities at the higher collision energy.

V. SUMMARY AND CONCLUSION
In summary, we have performed a global fit of transverse momentum particle spectra for identified pions, kaons and protons in five centrality classes based on a relativistic fluid approximation to QCD dynamics including a realistic thermodynamic equation of state as well as shear and bulk viscous dissipation, see ref. [ [57]. Prediction for the mid-rapidity densities dN/dy (|y| <0.5) of pions, kaon and protons are also reported.
One immediate result is the outcome for the most likely model parameters. They are summarized in Table III. In particular, the initialization time of the fluid description comes out relatively low, τ 0 = 0.27 fm/c. For the shear viscosity to entropy ratio we find η/s = 0.22, and for the peak value of the bulk viscosty to entropy ratio (ζ/s) max = 0.05. The combined chemical and kinetic freeze-out temperature is determined to be 136.9 MeV.
Moreover, from a quadratic expansion of χ 2 , corre- sponding to a Gaussian approximation to the posterior probability of the model parameters, we determine also their uncertainties as well as their correlation matrix, see Tables II and III. It becomes apparent that these uncertainties extracted from the χ 2 landscape are rather small. This underlines the quality of the experimental data and their high power to constrain theoretical models. However, one must also say that the best fitting model parameters lead to χ 2 /N dof = 1.40 with N dof = 546. The deviation from the expectation value χ 2 = N dof is actually relatively large, which strictly speaking, implies that it is very unlikely that the current theoretical model correctly describes all of the observed physics. In other words, the residual deviations in Fig. 3 are statistically significant. The tension concerns in particular pions in the region of low transverse momenta. We may speculate which physics effect our model is missing.
One possibility that comes to mind is that contributions from the feed down of decaying resonances have for some reason been underestimated. We have checked this possibility by doing our calculation with two different sets of hadronic resonances. While the current implementation uses the rather large set of ∼ 700 resonances of ref. [58], we have also tried a smaller set based on an earlier listing [14] and found the difference for the low-p T pions to be rather small. Of course, it cannot be fully excluded that an even larger set, or a more detailed resolution of the resonant widths, could remedy the problem.
Another interesting possibility is a non-thermal production mechanism for low-momentum pions such as from evolving coherent fields or condensates. An idea how this can happen in an out-of-equilibrium scenario is the one of a disoriented chiral condensate, see [59] for a review. Further work is needed to see whether such contributions from coherent fields and fluid dynamics can be reconciled.
Given that the current theoretical model is incomplete, it is rather difficult to determine its model parameters and the corresponding uncertainty. In particular, although straight-forward to calculate, the uncertainty from the χ 2 landscape as quoted in Table III can not be taken as a complete estimate of uncertainty in a situation where the theoretical description is itself not yet complete. For this reason we have also studied how our best fit parameters change when the procedure for their determination is varied. Specifically, in Fig. 2 we show how the best fit parameters change if the fit is not done globally, i. e. for all centralities and all three particle species, but rather separately for individual centrality classes (and all three species), separately for kaons or for kaons and pions together (but including all centralities). One observes that this leads indeed to a sizeable variation of the model parameters and we estimate on this basis the uncertainties from fit variations in Table III.
In conclusion we find that a fluid dynamic description of transverse momentum spectra for identified pions, kaons and protons works reasonably but with statistically significant residuals. The experimental data are now of a rather high quality and we expect that they will indeed allow to find a more complete theoretical description in the future.