Temperature dependence of $\eta/s$ of strongly interacting matter: effects of the equation of state and the parametric form of $(\eta/s)(T)$

We investigate the temperature dependence of the shear viscosity to entropy density ratio $\eta/s$ using a piecewise linear parametrization. To determine the optimal values of the parameters and the associated uncertainties, we perform a global Bayesian model-to-data comparison on Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV and Pb+Pb collisions at $2.76$ TeV and $5.02$ TeV, using a 2+1D hydrodynamical model with the EKRT initial state. We provide three new parametrizations of the equation of state (EoS) based on contemporary lattice results and hadron resonance gas, and use them and the widely used $s95p$ parametrization to explore the uncertainty in the analysis due to the choice of the equation of state. We found that $\eta/s$ is most constrained in the temperature range $T\approx 150$--$220$ MeV, where, for all EoSs, $0.08<\eta/s<0.23$ when taking into account the 90\% credible intervals. In this temperature range the EoS parametrization has only a small $\sim 10\%$ effect on the favored $\eta/s$ value, which is less than the $\sim 30\%$ uncertainty of the analysis using a single EoS parametrization. Our parametrization of $(\eta/s)(T)$ leads to a slightly larger minimum value of $\eta/s$ than the previously used parametrizations. When we constrain our parametrization to mimic the previously used parametrizations, our favored value is reduced, and the difference becomes statistically insignificant.


I. INTRODUCTION
The main goal of the ultrarelativistic heavy-ion collisions at the Large Hadron Collider (LHC) and the Relativistic Heavy-Ion Collider (RHIC) is to understand the properties of the strongly interacting matter produced in these collisions. In recent years the main interest has been in extracting the dissipative properties of this QCD matter from the experimental data (e.g. [1][2][3][4][5]), in particular its specific shear viscosity: the ratio of shear viscosity to entropy density η/s (for a review, see Refs. [6][7][8][9]). The field has matured to a level where a global Bayesian analysis of the parameters can provide statistically meaningful credibility ranges to the temperature dependence of η/s [10][11][12]. These credibility ranges agree with earlier results like those obtained using the EKRT model [13].
However, with the exception of papers like Refs. [14][15][16], the equation of state (EoS) is taken as given in the models used to extract the η/s ratio from the data. In particular, the EoS parametrization s95p [17] was used in many studies in the literature. This parametrization is based on by now outdated lattice data [18], and recent studies have reported an approximate 60% [16] or 30% increases [19] in the extracted value of η/s when switching from s95p to a contemporary lattice-based EoS. Further- * auvinen@ipb.ac.rs more, even if the errors of the contemporary lattice QCD calculations overlap, there is still a small tension between the trace anomalies obtained using the HISQ [20,21] and stout [22,23] discretization schemes. Consequently the EoSs differ, and if the procedure to extract η/s from the data is as sensitive to the details of the EoS as Refs. [16,19] claim, this tension may lead to additional uncertainties in the η/s values extracted from the heavyion collision data.
In the Bayesian analysis mentioned previously [10,11], the temperature dependence of η/s was assumed to be monotonously increasing above the QCD transition temperature T c . In a Bayesian analysis the slope parameter of such parametrization is always constrained to be non-negative, and limiting the final slope parameter to zero would require extremely strong constraints from the experimental data. Therefore, by construction, the analysis leads to an η/s increasing with temperature above T c , even if there is no physical reason to exclude a scenario where η/s is constant in a broad temperature range above T c . A more flexible parametrization, which does not impose such constraints, is thus needed to determine the temperature dependence of η/s.
In this work we address both the sensitivity of the extracted η/s to the EoS used in the model calculation, and the temperature dependence of η/s in the vicinity of the QCD transition temperature. We perform a Bayesian analysis of the results of EKRT + hydrodynamics calculations [13,24], and the data obtained in √ s NN = 200 GeV Au+Au collisions [25][26][27], and Pb+Pb collisions at 2.76 TeV [28][29][30] and 5.02 TeV [30,31]. To study the temperature dependence of η/s we use a piecewise linear parametrization in three parts: linearly decreasing and increasing regions at low and high temperatures are connected by a constant-value plateau of variable range. With this parametrization, data favoring a strong temperature dependence will lead to large slopes and a narrow plateau; conversely, an approximately constant η/s can be obtained with small slope parameter values and a wide plateau. To explore the sensitivity to the EoS, we use four different parametrizations: the well-known s95p parametrization, and three new parametrizations based on contemporary lattice QCD results. A comparison of the final probability distributions of the parameters will tell whether the most probable parameter values depend on the EoS used, and whether that difference is significant when the overall uncertainty in the fitting procedure is taken into account.

II. EQUATION OF STATE
In lattice QCD the calculation of the equation of state (EoS) usually proceeds through the calculation of the trace anomaly, Θ(T ) = (T ) − 3p(T ), where and p are energy density and pressure, respectively. Thermodynamical variables are subsequently derived from it using so-called integral method [32]. Therefore we base our EoS parametrizations on the trace anomaly and obtain pressure from the integral Once the pressure is known, the energy and entropy densities can be calculated, (T ) = Θ(T ) + 3p(T ), and s(T ) = [ (T ) + p(T )]/T , respectively. To make a construction of chemical freeze-out at T ≈ 150 MeV temperature possible, we use the hadron resonance gas (HRG) trace anomaly at low temperatures instead of the lattice QCD result. Equally important is that this choice allows for energy and momentum conserving switch from fluid degrees of freedom to particle degrees of freedom without any non-physical discontinuities in temperature and/or flow velocity 1 . Furthermore, it gives us a consistent value for the pressure at T low required for the evaluation of pressure (see Eq. (1)). 1 Energy and momentum conservation require that the fluid EoS is that of free particles, and that the degrees of freedom are the same in the fluid and particles [33]. If the dissipative corrections are small, switch from fluid consistent with the contemporary lattice QCD results [34] to particles in the UrQMD [35] or SMASH [36] hadron cascades at T = 150 MeV temperature leads to roughly 9-10% or 6-7% loss in both total energy and entropy, respectively.
As a baseline, we use the s95p parametrization [17], where HRG containing hadrons and resonances below M < 2 GeV mass from the 2004 PDG summary tables [37] is connected to the parametrized hotQCD data from Ref. [18]. To explore the effects of various developments during the last decade, we first connect the HRG based on the PDG 2004 particle list [37] to parametrized contemporary lattice data obtained using the HISQ discretization scheme [20,21]. The lattice spacing, a, is related to the temperature and temporal lattice extent, N t , as a = 1/(N t T ). Since the lattice spacing (N t ) dependence is small for this action, we use these results at fixed lattice spacing N t = 8, 10 and 12. We name our parametrizations according to the convention used to name s95p, and label this parametrization s87h 04 . 's87' signifies entropy density reaching 87% of its ideal gas value at T = 800 MeV, the letter 'h' refers to the HISQ action, and the subscript '04' to the vintage of the PDG particle list (2004). Note that even if our parametrization differs from the lattice trace anomaly in the hadronic phase, it agrees with the contemporary lattice calculations which show that at T = 800 MeV the entropy density reaches 87-88% of the ideal gas value (c.f. Fig. 8 of Ref. [21]).
The number of well-established resonances has increased since 2004, so we base our parametrization s88h 18 on HRG containing all 2 strange and non-strange hadrons and resonances in the PDG 2018 summary tables 3 [40], and on the same HISQ lattice data [20,21] we used for s87h 04 . Furthermore, there is a slight tension in the trace anomaly between the HISQ and stout discretization schemes. To explore whether this difference has any effect on hydrodynamical modeling, we construct the parametrization s83s 18 using PDG 2018 resonances, and the continuum extrapolated lattice data obtained using the stout discretization [22,23]. The second letter 's' in the label refers now to the stout action, and the subscript '18' to the vintage of the particle list. The details of these parametrizations are shown in Appendix A.
In the top and middle panels of Fig. 1, we show the parametrized trace anomalies, and the lattice data as used to make them: continuum extrapolated for the stout action, and at fixed lattice spacing for the HISQ action, since its lattice spacing (N t ) dependence is small. As seen in the topmost panel, the most noticable change in the lattice results during the last decade is the reduction of the peak of the trace anomaly (cf. s95p to others). Also, as mentioned, the lattice results obtained using the 2 With the exception of f 0 (500). See Refs. [38,39]. 3 Note that PDG Meson Summary Table and Baryon Summary  Table contain (almost) all states listed by the PDG, and are different from the PDG Meson Summary Tables and Baryon Summary Tables we use [ FIG. 1. The trace anomaly (top and middle) and the speed of sound squared (bottom) as functions of temperature in the four parametrizations of the EoS compared to the lattice data obtained using the HISQ [20,21] and stout [22,23] discretization schemes.
HISQ and stout actions slightly differ around the peak, and consequently s83s 18 differs from s87h 04 and s88h 18 . The higher peak does not, however, mean a lower speed of sound. As shown in the lowest panel of Fig. 1, the speed of sound in the s95p parametrization is not significantly lower than in the other parametrizations, but the temperature region where it is low is broader than in the other parametrizations. Thus we expect s95p to be effectively softer than the other EoSs. On the other hand, the speed of sound in s88h 18 depicts a characteristic dip below the speed of sound in the other parametrizations. This is a consequence of the parametrization of the trace anomaly in that temperature region.
As known, the HRG trace anomaly is below the lattice results [20,22,23] at low temperatures. This difference has been interpreted to indicate the existence of yet unobserved resonance states [34,42]. The need for further states has also been seen in the study of the strangeness baryon correlations on the lattice [43], and confirmed by the S-matrix based virial expansion [44]. However, we do not include predicted states from any model 4 in this work, since we do not know how they would decay, but use the states from the PDG summary tables only. Consequently the parametrized trace anomaly is slightly below even the most generous error bars of the lattice results around T ≈ 150-160 MeV temperature, as shown in the middle panel of Fig. 1.
On the other hand, whether we use the PDG 2004 or 2018 particle list causes only a tiny difference in the trace anomaly. The main difference between the s87h 04 and s88h 18 parametrizations arises from the connection of the HRG to the lattice parametrization. When parametrizing s88h 18 we wanted the trace anomaly to reach its lattice values soon above T c = 155 MeV, whereas we allowed s87h 04 to agree with lattice at larger temperature where the lattice trace anomaly drops below the HRG trace anomaly-for details, see Appendix A. Consequently the s88h 18 parametrization rises above the HRG values leading to the characteristic dip in the speed of sound (lowest panel in Fig. 1). Note that the s83s 18 parametrization does not depict a similar dip in the speed of sound, since the lower peak and larger errors of the continuum extrapolated stout action result allow the parametrized trace anomaly to drop below the HRG values immediately.

III. HYDRODYNAMICAL MODEL
We employ a fluid dynamical model used previously in Refs. [13,24,[46][47][48]. The spacetime evolution is computed numerically in (2+1) dimensions [49], and the longitudinal expansion is accounted for by assuming longitudinal boost invariance. We also neglect here the bulk viscosity and the small net-baryon number. The evolution of the shear-stress tensor π µν is described by the second-order Israel-Stewart formalism [50], with the coefficients of the non-linear second-order terms obtained by using the 14-moment approximation in the ultrarelativistic limit [51,52]. The shear relaxation time is related to the shear viscosity by τ π = 5η/( + p), where is energy density in the local rest frame, and p is thermodynamic pressure.
Transverse momentum spectra of hadrons are computed by using the Cooper-Frye freeze-out formalism at a constant-temperature surface, followed by all 2-and 3-body decays of unstable hadrons. The chemical freezeout is encoded into the EoS as described in Ref. [53], and the fluid evolves from chemical to kinetic freeze-out in partial chemical equilibrium (PCE) [54]. The kinetic and chemical freeze-out temperatures T dec and T chem are left as free parameters to be determined from the experimental data through the Bayesian analysis. The dissipative corrections δf to the momentum distribution at the freeze-out are computed according to the usual 14moment approximation δf k ∝ f 0k k µ k ν π µν , where f 0k is the equilibrium distribution function, and k µ is the fourmomentum of the hadron.
The remaining input to fluid dynamics are the EoS, initial conditions, and the shear viscosity. The different options for EoS were discussed in the previous section, and the initial conditions will be detailed in the next section. The temperature dependence of the shear viscosity η/s is parametrized in three parts, controlled by T H , the lower bound of the temperature range where η/s has its minimum value, (η/s) min , and the width of this temperature range, W min : where the additional parameters are the linear slopes below T H and above T Q = T H + W min , denoted by S HG and S QGP , respectively.
We note that bulk viscosity and chemical nonequilibrium are related [55,56]. Even if we ignore the bulk viscosity, some of its effects are accounted for by the fugacities in a chemically frozen fluid: At temperatures below T chem the isotropic pressure is reduced compared to the equilibrium pressure due to the different chemical composition. Thus introducing the chemical freeze-out changes not only the particle yields w.r.t. evolution in equilibrium, but similarly to the bulk viscosity, reduces the average transverse momentum of hadrons too. However, this affects the evolution only when temperature is below T chem , and in contrast to the bulk viscosity, there is e.g. no entropy production associated with the chemical freeze-out and subsequent chemical non-equilibrium [58].
Finally, we emphasize that we solve the spacetime evolution from the hot QGP all the way to the kinetic freezeout as a single continuous fluid dynamical evolution. This is different from the hybrid models used e.g. in Refs. [10][11][12] where the evolution below some switching temperature is solved with a microscopic hadron cascade. The advantage of the fluid dynamical evolution without a cascade stage is that the transport properties are continuous in the whole temperature range. Note that in the hybrid models the switching from fluid dynamics to hadron cascade introduces an unphysical discontinuity in e.g. η/s that is O(1) in the cascade [57], but O(0.1) in fluid dynamical simulations at switching. Another advantage of our approach is that we can freely parametrize the viscosity in the hadronic matter too, and constrain it using the experimental data.

IV. INITIAL CONDITIONS
The initial energy density profiles are determined using the EKRT model [59][60][61] based on the NLO perturbative QCD computation of the transverse energy, and a gluon saturation conjecture. The latter controls the transverse energy production through a local semi-hard clear thickness function at transverse location (x, y). The essential free parameters in the EKRT model are the proportionality constant K sat in the saturation condition, and the constant β controlling the exact definition of the minijet transverse energy at NLO [60]. The setup used here is identical to the one used in Refs. [13,24,48], where β = 0.8, and K sat is left as a free parameter to be determined from the data. We note that K sat is independent of the collision energy √ s NN and nuclear mass number A, so that once K sat is fixed the √ s NN and A dependence of the initial conditions is entirely determined from the QCD dynamics of the EKRT model. With a given p sat the local energy density at the formation time τ p = 1/p sat can be written as This we further evolve to the same proper time τ 0 = 1/p min , where p min = 1 GeV, at every point in the transverse plane where p sat > p min by using 0+1 dimensional Bjorken hydrodynamics with the assumption = 3p. In the EKRT model, fluctuations in the product of the nuclear thickness functions, T A T A , give rise to the eventby-event fluctuations in the energy density through p sat in Eq. (3). Moreover, the centrality dependence of the initial conditions arises from the centrality dependence of T A T A . A full treatment of the dynamics in heavy-ion collisions would take the event-by-event fluctuations into account by evolving each event separately. However, to make the present study computationally feasible, we omit the evolution of such fluctuations here; instead, for each centrality class, we average a large number of these fluctuating initial states, and compute the fluid dynamical evolution only for the averaged initial distributions.
The computed energy densities are not linear in K sat nor in T A T A , and different averaging procedures can lead to significantly different event-averaged initial conditions. In the previous event-by-event EKRT studies [13,24,48] a fair agreement was obtained between the data and the computed √ s NN , A, and centrality dependencies of the charged hadron multiplicity. To preserve as much as possible of this agreement, we average the initial conditions by averaging the initial entropy distributions: We compute first a large set of initial energy density profiles using the procedure detailed in Ref. [13]. Each of the generated energy density profiles is converted to an entropy density profile by using the EoS which will be used later during the evolution. The entropy density profiles are then averaged, and the average entropy density profile is converted to an average energy density profile using the same EoS. In the event-by-event framework the centrality classes were determined from the final multiplicity distribution. However, this way of classifying events is not available here, as it would require fluid dynamical evolution of each of the fluctuating initial conditions. Instead, we predetermine the centrality classes according to the number of wounded nucleons in the sampled Monte-Carlo nuclear configurations, which were used to construct the eventby-event initial conditions. The number of wounded nucleons are computed using the nucleon-nucleon cross section σ NN = 42, 64, and 70 mb for 200 GeV, 2.76 TeV, and 5.02 TeV collisions, respectively. We note that the nucleon-nucleon cross section does not enter in the computation of the initial conditions, but they are used here only in the centrality classification. In the context of the full event-by-event modeling we have tested that the final results are only weakly sensitive to the precise way of the centrality classification.
Let us consider each combination of the free parameters as a point x in the 8-dimensional input (parameter) space, the model output y( x) as a corresponding point in the 90-dimensional output space (space of observables), and the experimental data y exp as the target point in the 5 Charged particle multiplicities at RHIC are averages over two adjacent PHENIX centrality classes; for example, at (10-20)% centrality N ch is an average over (10-15)% and (15)(16)(17)(18)(19)(20)% classes, (20-30)% is an average over (20-25)% and (25-30)% classes, and so on. This applies also for RHIC identified particle data at (10-20)% centrality. 6 We consider an average of measured protons and antiprotons as the target value for the proton multiplicity at RHIC. space of observables. With these definitions we can formulate the posterior probability distribution P ( x| y exp ) of the best-fit parameter values by utilizing Bayes' theorem: where P ( x) is the prior probability distribution of input parameters and P ( y exp | x) is the likelihood function Here Σ is the covariance matrix representing the uncertainties related to the model-to-data comparison. As a function with an eight-dimensional domain, the posterior probability distribution P ( x| y exp ) is too complicated to evaluate and analyze fully. Instead, we produce samples of it with a parallel tempered Markov chain Monte Carlo [62] based on the emcee sampler [63]. An ensemble of random walkers is initialized in the input parameter space based on the prior probability 7 and each proposed step in parameter space is accepted or rejected based on the change in the value of the likelihood function. At a large number of steps, the distribution of the taken steps is expected to converge to the posterior distribution.
Evaluating the output y( x) of the fluid dynamical model at every point x where the random walker might enter is a computationally impossible task. Therefore we approximate the output using Gaussian process (GP) emulators [64] (see Appendix B). Each GP is able to provide estimates for only one observable, so to keep the number of required emulators manageable, we perform a principal component analysis (PCA) to reduce the dimension of the output space from 90 observables into k = 6 most important principal components. Further details about the PCA are described in Appendix C. We utilize the scikit-learn Python module [65] and in particular the submodules sklearn.gaussian process and sklearn.decomposition.PCA in the model emulation.
Thus, in our likelihood function (5), we replace y( x) with the GP estimate in the principal component space z GP ( x) (likewise y exp is transformed to z exp ), and include the emulator estimation error into the covariance matrix: where Σ exp z is the (originally diagonal) experimental error matrix transformed to principal component space, and is the GP emulator covariance matrix obtained from the emulator (see Appendix B).
To work, the GP emulators must be conditioned with a set of training points, { z( x i )}, created by running the fluid dynamical model with several different parameter combinations { x i }. For the present investigation, we have produced 170 training points for each EoS, distributed evenly in the input parameter space 8 using minmax Latin hypercube sampling [66]. The emulation quality was then checked by using the trained emulator to predict the results at 30 additional test points, which were not part of the training data. An example of the results of this confirmation process is shown in Fig. 2

for 2.76
TeV Pb+Pb collisions using the s95p parametrization.

VI. RESULTS
The marginal posterior probability distributions for each parameter are obtained from the full 8-dimensional probability distribution (see Section V) by integrating over the other seven parameters. The resulting distributions when using the four investigated EoSs are shown in Figs. 3 and 4. In these figures the range of the x-axis illustrates the prior range of values, with the exception of T chem which range depends on the EoS 9 . The median values of these distributions provide a good approximation for the most probable values, and are listed both in the legends of the figures, and in Table I. The 90% credible intervals-i.e. the range which covers 90% of the distribution around the median-are shown as errors in Table I. Two dimensional projections of the probability distribution depicting correlations between parameter pairs are shown in Appendix D. 8 The restriction T dec < T chem does not apply to the training points. 9 For s83s 18 , s87h 04 , and s88h 18 , the prior range is 120 MeV < T chem < T 0 , where T 0 is the temperature where the parametrization deviates from the HRG (see Appendix A). For s95p the range is 120 < T chem /MeV < 180.

A. Nuisance parameters
The analysis involves three parameters which are not directly related to the transport properties of produced matter: K sat , T dec , and T chem . The probability distributions for these three "nuisance" parameters, shown in Fig. 3, are nicely peaked, and the parameters have well defined constraints. For the chemical freeze-out temperature, the median is T chem = 153-155 MeV, which is compatible with the values obtained using the statistical hadronization model [67]. Note that the difference in the median is not due to the resonance content of the EoS, but due to a complicated interplay of the softness of the EoS, shear, and build-up of the flow. Nevertheless, the particle ratios are the dominant factor in constraining T chem .
For K sat and T dec , we see a common trend where s95p gives a distribution which peaks at the lowest value of the four EoSs, followed by s87h 04 , and the highest peak values are shared by s83s 18 and s88h 18 with almost identical distributions. The obtained values for the EKRT normalization parameter, K sat ≈ 0.5, are compatible with the values found previously [13], and the small differences between different EoS parametrizations result from slightly different entropy production during the evolution. Differences seen in the kinetic freeze-out temperature T dec = 126-132 MeV are also small, and seem to follow the conventional rule of thumb: a softer EoS requires a lower freeze-out temperature to create hard enough proton p T distributions. On the other hand, differences in the median values of all these three parameters are smaller than the credibility limits, and thus not statistically meaningful.  At first sight (η/s) min depicts the behavior described in Refs. [16,19]: the favored value is lower for s95p than for the newer parametrizations (see Fig. 3 and Table I). However, the effect is noticeably smaller than seen in those studies-only ∼ 10-20%-and well within the 90% credible intervals (± ∼ 30%) of the analysis. The comparison of η/s for different EoSs is further complicated by the large number of parameters controlling the temperature dependence of η/s. The probability distributions of parameters T H , W min , S HG , and S QGP , shown in Fig. 4, are very broad extending to the whole prior range in most cases, and thus do not possess any clearly favored values. However, the wide posterior distributions of the (η/s)(T ) parameters are partly caused by the inherent ambiguity in the chosen parametrization: for a given temperature T , multiple parameter combinations can generate similar values of η/s. For example, at low temperatures (η/s)(T ) is mostly determined by S HG and T H , but it is better constrained than either of these parameters. The rea-son is that S HG and T H are not independent, but slightly anti-correlated-the correlations between the pairs of parameters are shown in Appendix D. Thus it is more illustrative to construct the probability distribution for η/s values w.r.t. temperature, and plot the median and credibility intervals of this distribution as shown in Figs. 5 and 6. In the upper panel of Fig. 5 we show the median of (η/s)(T ) for each EoS parametrization, and the union and intersection of the 90% credible intervals of all four distributions. The union of the credibility intervals provides insight on the total uncertainty in the analysis including the uncertainty from the EoS parametrization, whereas the difference between the union and intersection illustrates how much of the uncertainty comes from the EoS parametrizations. To emphasize the result using state-of-the-art EoSs, the lower panel of Fig. 5 depicts the median and credibility intervals for the parametriza- tions s83s 18 and s88h 18 only. In the same panel two older results from Ref. [13], and a recent theoretical prediction from Ref. [68] are shown as well. To make it possible to distinguish the credibility intervals for each EoS separately, η/s for each EoS at various temperatures with associated uncertainties is shown in Fig. 6.
We obtain well constrained η/s in a temperature range 150 < ∼ T /MeV < ∼ 220, where the median values of η/s are practically constant for all contemporary EoSs, and s95p leads to modest temperature dependence well within the credibility intervals. Within this temperature range η/s is constrained between 0.08 and 0.23 by the 90% credible intervals. In particular, for the state-of-the-art EoSs (s83s 18 and s88h 18 ), we obtain even tighter limits 0.12 < η/s < 0.23 within this range, and the well constrained region extends to slightly higher temperature. For further details see Fig. 5 and Table II. Interestingly η/s at 130 MeV (or at 150 MeV in case of s95p) temperature differs from the favored value (median) of the (η/s) min parameter (compare Tables I and II)  MeV) (see Fig. 4 and Table I). This seemingly counterintuitive behavior is due to the fat tails of T H distributions extending to larger temperatures, and thus broadening the region where S HG affects the η/s values. Consequently we see the lowest η/s values at T ≈ 200 MeV temperature ( Fig. 5 and Table II), where the effect of the lower (η/s) min value of the s95p parametrization is also visible. It is not surprising that we get the best constraints on η/s in the temperature range 150 < ∼ T /MeV < ∼ 220. As was shown in Ref. [46], the temperature range where v 2 is most sensitive to the shear viscosity is only slightly broader than this, and higher order anisotropies are sensitive to shear at even narrower temperature ranges 10 .
Even if the uncertainties remain large, we can see qualitative differences in the high temperature behavior of η/s, where s95p seems to favor earlier and more rapid rise of η/s with increasing temperature (Figs. 5 and 6), a difference which is visible in the S QGP parameter as well (Fig. 4).
Considering earlier results in the literature this is intriguing. Alba et al. [16] used an EoS based on contemporary stout action data called PDG16+/WB2+1, and observed that the reproduction of the LHC data ( √ s NN = 5.02 TeV) required larger constant η/s for this EoS than for s95p. On the other hand, they were able to use the same value of constant η/s for both EoSs to reproduce the RHIC data. They interpreted this to mean that at large temperatures s95p would necessitate lower values of η/s, but we see an opposite behavior. In a similar fashion we see a difference between the high temperature behavior obtained using the HISQ (s88h 18 and s87h 04 ) and stout action based EoSs (s83s 18 ), but the differences are way smaller than the credibility intervals, and thus cannot be considered meaningful.
At temperatures below 150 MeV we again see expanding credibility intervals, and a tendency of η/s to increase with decreasing temperature, but hardly any sensitivity to the EoS. Anisotropies measured at RHIC en-ergy are sensitive to the shear viscosity in the hadronic phase [46,47], and since Schenke et al. in Ref. [19] saw sensitivity to the EoS using RHIC data only, we would have expected some sensitivity to the EoS at low temperatures. The difference may arise from the bulk viscosity which depended on the EoS as well in Ref. [19], or from a different EoS in the hadronic phase. As mentioned, our EoSs are based on known resonance states, whereas the EoSs in Refs. [16,19] follow the lattice results closely. Better fit to lattice results can be obtained by including predicted but unobserved resonance states in the HRG. We plan to study how the inclusion of these states might affect the results, once we have concocted a plausible scheme for their decays, so that we can evaluate their contribution to the EoS after chemical freeze-out in a consistent manner.
Furthermore, unlike in Ref. [19] where a hadron cascade was used to describe the evolution in the hadronic phase, in our approach the change in the EoS can also be partly compensated by a change in the freeze-out temperature instead of shear viscosity. As shown in Appendix D, there is indeed an anti-correlation between T dec and (η/s) min . Therefore forcing the system to freeze out at the same temperature, independent of the EoS, would increase the difference in (η/s) min . However, the anti-correlation is rather weak ∼ −0.4(−0.2) for s88h 18 (s95p), and thus requiring EoS independent T dec would not change (η/s) min a lot.
Our result of a very slowly rising η/s with decreasing temperature in the hadronic phase (i.e., below T ≈ 150 MeV) may look inconsistent with microscopic calculations predicting relatively large η/s ∼ 1 in the hadronic phase [57,69,70]. However, our result is for a chemically frozen HRG, while the microscopic calculations usually give η/s in chemical equilibrium. At a given temperature the entropy density s PCE in a chemically frozen HRG can be significantly larger than the entropy density in chemical equilibrium s CE , and as a consequence η/s PCE can be way smaller than η/s CE . We may obtain an approximation for the η/s in a chemically equilibrated system as η/s CE = (η/s PCE )(s PCE /s CE ) [13], since, as a first order approximation, η depends only weakly on the chemical non-equilibrium [71]. In our case, where T chem = 154 MeV, the ratio of entropies in a chemically frozen to a chemically equilibrated system is ∼ 3.5 at T = 100 MeV (∼ 1.8 at T = 130 MeV) which is sufficient to bring our results to the level described in Ref. [70].
In Fig. 5 we also made a comparison to the earlier results of Ref. [13] and the recent quasiparticle model prediction from Ref. [68]. As expected, the earlier results from Ref. [13,24] are not far from the present analysis, and param1 is practically within the 90% credible interval in the whole temperature range. On the other hand, constant η/s = 0.2 is below the s95p limits at high temperatures, but as discussed, the overall sensitivity to η/s at high temperatures is low. Interestingly the prediction of the quasiparticle model of Ref. [68] comes very close to our values for η/s around T c , although the region where η/s is low is narrower than what we found here. This is intriguing, since the quasiparticle model was tuned to reproduce the stout action EoS, i.e., our EoS s83s 18 , which in our analysis leads to the broadest region where η/s is almost constant.
The small value of η/s and its weak temperature dependence in the temperature range 150 < ∼ T /MeV < ∼ 220 may indicate that the QGP is strongly coupled not only in the immediate vicinity of T c , but in a broader temperature region. This was first proposed in Ref. [72], and agrees with the lattice QCD calculations that indicate the presence of hadronlike resonances in QGP in a similar or slightly broader temperature interval [73][74][75][76]. The strongly coupled nature of QGP can also be seen in the large value of the coupling constant defined in terms of the free energy of static quark anti-quark pairs [77]. In any case, our result for (η/s)(T ) is compatible with the lattice QCD calculations, which indicate that weakly coupled QGP picture may be applicable only for T > 350 MeV [77][78][79][80].

C. The effect of the parametric form
When we use the state-of-the-art EoSs (s88h 18 and s83s 18 ), our result for the minimum value of η/s is higher than the result obtained in an earlier Bayesian analysis of Ref. [10]: 0.12 < η/s < 0.23 vs. η/s = 0.07 +0.05 −0.04 . An important difference in these analyses is that Ref. [10] assumed the minimum of η/s to occur at fixed T = 154 MeV temperature, and η/s to rise linearly above that temperature. Moreover, below T = 154 MeV they used a hadron cascade to model the evolution, and the transport properties of the hadronic phase were thus fixed.
To explore how much the results depend on the form of the (η/s)(T ) parametrization, we mimic the parametrization used in Ref. [10] by constraining the plateau in our parametrization to be very small (0 < W min /MeV < 2), and the minimum to appear close to T c (150 < T H /MeV < 160). The resulting temperature dependence of η/s for the s88h 18 and s95p parametrizations is shown in Fig. 7, and compared to our full result (the behavior of the s87h 04 and s83s 18 parametrizations is similar to s88h 18 ).
The change in parametrization reduces the minimum value of η/s to 0.12 +0.03 −0.03 for s88h 18 and 0.06 +0.04 −0.04 for s95p, which are consistent with the value obtained in Ref. [10]. Another interesting change is seen in the high-temperature behavior. In the full analysis the s95p parametrization leads to the largest η/s at large temperatures, but the restricted parametrization causes s95p to favor the lowest η/s at large temperatures. As seen previously, s95p favors the lowest η/s at 200 < T /MeV < 250 temperatures (see Fig. 6 and Table II), which in the restricted parametrization dictates the behavior at much higher temperatures as well.
Nevertheless, even if the results depend on the form of the parametrization, the credibility intervals overlap and the results are consistent. The only deviation from this rule is for the s95p parametrization around T = 160 MeV temperature where the difference is statistically significant (see Fig. 7). We have also checked that when we use the favored parameter values, the typical differences in the fit to the data due to different parametric forms are only ∼ (1 − 3)%.
Similarly, we can mimic temperature independent η/s by constraining the priors of the S HG and S QGP parameters close to zero. We have checked that such a choice does not increase the sensitivity of η/s to the EoS parametrization, and that the median values of the constant η/s = (η/s) min were only ≈ 10% larger than the median values for η/s at T = 200 MeV for the full parametrization. Again, a sign of v 2 being most sensitive to shear viscosity in the 150 < ∼ T /MeV < ∼ 220 temperature range [46].
Thus, in the Bayesian analysis the parametric form of η/s does affect the results, and is therefore a kind of prior whose effects are difficult to quantify. On the other hand, the credibility intervals overlap in all the cases, which emphasizes their importance: The "true" value could be anywhere within the credibility interval, and there is still a 10% chance it is outside of it. GeV compared to PHENIX data [25]. Middle panel: Pb+Pb at √ sNN = 2.76 TeV compared to ALICE data [28]. Bottom panel: Pb+Pb at √ sNN = 5.02 TeV compared to ALICE data [31].
Finally, as an overall quality check, we show how well the favored parameter combinations reproduce the experimental data. This is done by drawing 1000 samples from each posterior distribution and using the Gaussian process emulator to predict the simulation output for these values. The results for charged and identified particle multiplicities, identified particle p T , and the elliptic flow  [26]. Right panels: Pb+Pb at √ sNN = 2.76 TeV compared to ALICE data [29] .
tively. The overall agreement with the data is quite good for all observables, and the analysis is able to find equally good data fits for all four EoSs. As normal for thermal models, the charged particle multiplicities tend to be underestimated due to the tension between pion multiplicity on one hand, and kaon and proton multiplicities on the other hand. As the analysis makes a compromise between too few pions and too many kaons and protons, the overall charged particle multiplicity (which is dominated by pions) will remain below the data. Also the mean transverse momentum of pions is slightly too large, which may prove difficult to alleviate without the introduction of bulk viscosity [4] and/or improved treatment of resonances during the hadronic phase [81].

VII. SUMMARY
In this work, we have introduced three new parametrizations of the equation of state based on the contemporary lattice data: • s87h 04 connects the HRG based on the PDG 2004 particle list to parametrized lattice data obtained using the HISQ discretization scheme.
• s88h 18 is based on the HRG containing all strange and non-strange hadrons and resonances in the PDG 2018 summary tables, and the same HISQ lattice data as s87h 04 .
• s83s 18 is constructed using the PDG 2018 resonances, and the continuum extrapolated lattice data obtained using the stout discretization.
We used these new parametrizations and the older s95p parametrization to examine how sensitive the shear viscosity over entropy density ratio η/s is to the equation of state. We assumed a piecewise linear parametrization for (η/s)(T ), and determined the probability distributions of the best-fit parameter values within the EKRT framework using a Bayesian statistics approach. Using charged and identified particle multiplicities, identified particle mean transverse momenta, and elliptic flow at three different collision energies as calibration data, we were able to constrain the value of η/s to be between 0.08 and 0.23 with 90% credibility in the temperature range 150 < ∼ T /MeV < ∼ 220 when all EoS parametrizations are taken into account. When we constrain the EoSs to the most contemporary parametrizations s83s 18 and s88h 18 , we obtain 0.12 < η/s < 0.23 in the above mentioned temperature range. As the differences between the EoSs are well covered by the 90% credible intervals, the earlier results obtained using the s95p parametrization remain valid. The weak sensitivity to the EoS is consistent with the old ideal fluid results √ sNN = 200 GeV compared to STAR data [27]. Middle panel: Pb+Pb at √ sNN = 2.76 TeV compared to ALICE data [30]. Right panel: Pb+Pb at √ sNN = 5.02 TeV compared to ALICE data [30].
for flow and EoS: Based on flow alone, it is difficult to distinguish an EoS with a smooth crossover from an EoS without phase transition [82]. Thus when the differences between EoSs are just details in the crossover, the differences in flow, which should be compensated by different shear viscosity, are small, and consequently differences in the extracted η/s are small. The overall agreement with the data is quite good, and similar to Refs. [13,24], where event-by-event fluctuations were included to the framework of EKRT initial conditions and fluid dynamics, albeit without the Bayesian analysis. The good agreement achieved here is partly due to the EKRT initial conditions. In particular the centrality and √ s NN dependence of hadron multiplicities follow mainly from the QCD dynamics of the EKRT model. A noticeable difference to the earlier event-by-event analysis is that here we used identified hadron multiplicities as constraint, which led to the chemical freeze-out temperature T chem ∼ 154 MeV, and a slight overshoot of the pion average p T compared to the data. In the earlier analysis T chem ∼ 175 MeV was used to reproduce the average p T data, which in turn led to too large proton multiplicity. It is possible to solve this tension by introducing bulk viscosity [4], but that is left for a future work. We emphasize that compared to the (in principle) more detailed hydro + cascade models our hydro + partial chemical equilibrium approach has two major advantages: It allows us to parametrize η/s(T ) so that it is continuous in the whole temperature range, and at the same time it gives us a possibility to constrain the viscosity also in the hadronic phase. Inclusion of event-by-event fluctuations to the analysis would provide access to several new flow observables such as higher flow harmonics v n , and flow correlations, which may give tighter constraints in broader temperature interval on η/s(T ). However, within the current uncertainties of the fitting procedure, we cannot exclude the possibility that the effect of the EoS remains negligible even when η/s at T > 220 MeV becomes better under control.
Since the sensitivity of flow to shear viscosity at high temperatures is low, observables based on high p T particles may be useful to constrain, not only the preequilibrium dynamics [83][84][85], but also the properties of the fluid when it is hottest.

ACKNOWLEDGMENTS
We thank V. Mykhaylova    anomaly, and its first and second derivatives are continuous. This requirement provides constraints for three parameters, d 0 , d 1 , and d 2 , and leaves the remaining seven, d 3 , d 4 , d 5 , n 3 , n 4 , n 5 , and T 0 to be fixed by minimizing a χ 2 fit to the data. Fitting the powers n 3 -n 5 would be a highly non-linear problem, but we simplify the problem by requiring that the powers are integers, and using brute force: We make a fit with all the integer values 5 ≤ n 3 ≤ 40, n 3 < n 4 ≤ 41, and n 4 < n 5 ≤ 42, and choose the values n 3 , n 4 , and n 5 which lead to the smallest χ 2 . When the powers and T 0 are kept fixed, minimizing χ 2 requires only a simple matrix inversion. Thus to fix T 0 we are able to cast χ 2 as a function of only a single parameter, T 0 . We require that 155 ≤ T 0 /MeV ≤ 190, and search for the value of T 0 which minimizes χ 2 .
To obtain the continuum limit in the lattice calculations of the trace anomaly, one has to perform interpolation in the temperature, and then perform continuum extrapolations (see e.g. [23]). This procedure can introduce additional uncertainties when providing parametrization of the lattice results. As mentioned in the main text, the lattice spacing (N t ) dependence of the lattice results on the trace anomaly is small in the case of the HISQ discretization scheme for N t ≥ 8. In fact, for T > 230 MeV and T < 170 MeV there is no statistically significant N t dependence, so in these temperature ranges we can use the HISQ lattice results with N t = 8, 10 and 12. In the peak region, 170 < T /MeV < 230, the N t = 8 HISQ results are slightly higher than the N t = 10 and N t = 12 results, and therefore have been omitted from the fits. At temperatures above 800 MeV only lattice results with N t = 6 and 4 are available [20,21]. To take the larger discretization errors of the N t = 6 and 4 results into account, we follow Ref. [21], scale them by factors 1.4 and 1.2, and include systematic errors of 40% and 20%, respectively. Contrary to the HISQ action results, we employ the continuum extrapolated stout action results [22,23] for simplicity. The resulting parameters are shown in Table III. We find that only the parametrization s88h 18 requires the use of all six terms in Eq. A1. In the cases of s83s 18 and s87h 04 we are able to obtain equally good fits with only five terms, and thus set d 5 to zero by hand.
For the sake of completeness, we also parametrize the HRG part of the trace anomaly as − 3p T 4 = a 1 T m1 + a 2 T m2 + a 3 T m3 + a 4 T m4 . (A2) To carry out the fit we evaluate HRG trace anomaly in temperature interval 70 < T /MeV < T high , where T high depends on the parametrization, with 1 MeV steps assuming that each point has equal "error". The limits have entirely utilitarian origin: in hydrodynamical applications the system decouples well above 70 MeV temperature and only a rough approximation of the EoS, p = p( ), is needed at lower temperatures. On the other hand we expect to switch to the lattice parametrization below T high , and the HRG EoS above that temperature is not needed either. We fix the powers in Eq. (A2) again using brute force. We require them to be integers, go through all the combinations 0 ≤ l 1 < l 2 < l 3 < l 4 ≤ 10, fit the parameters a 1 , a 2 , a 3 , a 4 to the HRG trace anomaly evaluated with 1 MeV intervals, and choose the values l 1 , l 2 , l 3 and l 4 which minimize the χ 2 . We end up with parameters shown in Table IV. To obtain the EoS, one also needs the pressure at the lower limit of the integration (see Eq.(1)) T low = 0.07 GeV: p(T low )/T 4 low = 0.1661. Our EoSs are available in a tabulated form at arXiv as ancillary files for this paper, and at Ref. [86]. These tables also include the option of a chemically frozen hadronic stage, and a list of resonances included in the hadronic stage with their properties and decay channels. Let us assume that we do not know exactly what the model's output y for a particular input parameter x is, but we know its most probable value µ( x). We postulate that the probability distribution for the output value P (y) is a normal distribution with mean µ( x) and so far unknown width σ. Thus the probability distribution for a set Y a of N model outputs for observable a, corresponding to a set X of N points in the parameter space, is a multivariate normal distribution: where µ = µ(X) = {µ( x 1 ), . . . , µ( x N )} is the mean of the distribution, and C is the covariance matrix defined by the covariance function c( x, x ): As we are only interested in interpolating within the training data, we may set µ(X) ≡ 0, and construct the covariance function c( x, x ) in such a way that the probability distribution is narrow at the training points nevertheless. This way we minimize our a priori assumptions about the model behavior in regions of parameter space not covered by the training data 11 . Our chosen covariance function is a radial-basis function (RBF) with a noise term The hyperparameters θ = (θ 0 , θ 1 , . . . , θ n , θ noise ), where n is the dimension of the input parameter space, are not known a priori and must be estimated from training data, consisting of simulation output U computed at training points T , by maximizing the log-likelihood (see Chapter  5 Emulator prediction for the model output y 0 at a point x 0 can then be determined by writing a joint probability distribution for the output at various points in parameter space: from which we can derive the conditional predictive mean y GP ( x 0 ) and associated variance σ GP ( x 0 ) 2 as (see e.g. while for matrices (such as the covariance matrix in the likelihood function (5)) the transformation is To compare an emulator prediction z GP against physical observables, we use the inverse transformation