Bayesian estimation of the specific shear and bulk viscosity of the quark-gluon plasma with additional flow harmonic observables

The transport properties of the strongly-coupled quark-gluon plasma created in ultra-relativistic heavy-ion collisions are extracted by Bayesian parameter estimate methods with the latest collision beam energy data from LHC. This Bayesian analysis includes sophisticated flow harmonic observables for the first time. We found that the temperature dependence of specific shear viscosity shows weaker than the previous studies. The results prefer a lower value of specific bulk viscosity and a higher switching temperature to reproduce additional observables. However, the improved statistical uncertainties both on the experimental data and hydrodynamic calculations with additional observables do not help to reduce the final credibility ranges much, indicating a need for improving the dynamical collision model before the hydrodynamic takes place. In addition, the sensitivities of experimental observables to the parameters in hydrodynamic model calculations are quantified. It is found that the analysis benefits most from the symmetric cumulants and non-linear flow modes, which mostly reflect non-linear hydrodynamic responses, in constraining the temperature dependence of the specific shear and bulk viscosity in addition to the previously used flow coefficients.

The transport properties of the strongly-coupled quark-gluon plasma created in ultra-relativistic heavy-ion collisions are extracted by Bayesian parameter estimate methods with the latest collision beam energy data from the CERN Large Hadron Collider. This Bayesian analysis includes sophisticated flow harmonic observables for the first time. We found that the temperature dependence of specific shear viscosity appears weaker than in the previous studies. The results prefer a lower value of specific bulk viscosity and a higher switching temperature to reproduce additional observables. However, the improved statistical uncertainties both on the experimental data and hydrodynamic calculations with additional observables do not help to reduce the final credibility ranges much, indicating a need for improving the dynamical collision model before the hydrodynamic takes place. In addition, the sensitivities of experimental observables to the parameters in hydrodynamic model calculations are quantified. It is found that the analysis benefits most from the symmetric cumulants and non-linear flow modes, which mostly reflect non-linear hydrodynamic responses, in constraining the temperature dependence of the specific shear and bulk viscosity in addition to the previously used flow coefficients.

I. INTRODUCTION
The primary goal of heavy-ion physics is to investigate and understand the strongly coupled color-deconfined matter, quark-gluon plasma (QGP), which is produced in ultrarelativistic collisions between heavy ions. The QGP is believed to be the predominant form of matter during the first phases of the early universe. This matter behaves like a near-perfect fluid with the smallest specific shear viscosity, the ratio of the shear viscosity to the entropy density (η/s), of any known substance in nature [1].
The most important remaining open questions in the field are the location of the critical point (T c ) in the QCD phase diagram and temperature dependence of specific shear (η) and bulk (ζ) viscosities of the QGP. The flow analysis in the Large Hadron Collider (LHC) has been very successful and provides valuable information to the field [2][3][4][5][6]. For example, the main constraints for the QGP properties using the Bayesian analysis [7] in the theory came from the ALICE measurements [2,8,9] with both low and high beam energy data. Even though the Bayesian analyses [7,[10][11][12][13] were successful, the current uncertainties from these works are large because of statistical limitations of the data, limited observables used for the analysis, and computational constraints. In addition to the aforementioned limitations, pinning down the absolute value of η/s at T c has a few challenges. First, a principle calculation to describe the initial conditions (IC) is still under development. Second, extracting the temperature dependence of η/s(T ) has been complicated with the existence of the bulk viscosity [14,15]. However, * jasper.parkkila@cern.ch large flow found in small systems like proton-proton (pp) collisions was striking and opened up the importance of gluon fluctuations within protons and certainly, the experimental data would help to improve the understanding of IC both for small and large systems [16]. There are newer observables that give much stronger constraints to the theory [3,4], showing good sensitivities especially to η/s(T ) and ζ/s(T ). The correlation strength measured in [3,4] was experimentally decomposed into two components of linear and non-linear flow modes in [5,17] for the first time in the field, which gives a better understanding of our harmonic analysis and its origin with both LHC Run 1 (2009Run 1 ( -2013 and Run 2 (2015Run 2 ( -2018 data.
In this work, we extend the Bayesian parameter estimation methods employed in [7] with larger statistic LHC Run 2 results [18,19] as well as a few additional observables [5,6] for the first time which require substantial computational power. This work also allows us to quantify the sensitivity of each observable to the hydrodynamic model parameters in a controlled way. In Sec. II, we present a brief overview of Bayesian analysis methods and model setups. The experimental observables are described in Sec. III. Model parameters and calibrations are explained in Sec. IV. The results are presented in Sec. V, after which Sec. VI summarizes our results and findings.

II. BAYESIAN ANALYSIS
There have been a number of studies that utilized Bayesian methods for heavy ion collisions [20][21][22][23][24]. We employ the recent state-of-the-art development in [7] for our present study. We define a vector of model parameters x, and a set of experimental data y that will be compared with model calculations. Bayes's theorem gives the posterior distributions for the model parameters as P (x|y) ∝ P (y|x)P (x). (1) Here P (y|x) is the likelihood, which quantifies the model agreement with the data. The prior P (x) encapsulates initial knowledge on the parameters. The model parameters are then extracted from the posterior distributions. We follow the same procedures as [7], where the model is first evaluated at a small O(10 2 ) number of 'design' parameter points. The resulting discrete set of model predictions is then made continuous by the use of a Gaussian process (GP) emulator, which thereby can be used to systematically probe the parameter space with Markov chain Monte-Carlo (MCMC) methods.

A. Hydrodynamic model and parameters
The model used in this analysis consists of the T R ENTo model [25] for the initial condition, which is connected with free streaming to a 2+1 dimensional causal hydrodynamic model VISH(2+1) [26,27]. The evolution continues after particlization via the UrQMD model [28,29].
A hydrodynamic modeling relies on the energy and momentum conservation laws of the fluid dynamics. The conservation is expressed in terms of where T µν (x) is the energy-momentum tensor. In the case of viscous hydrodynamics, the energy-momentum tensor becomes where is the energy density, P is the local pressure given by the equation of state, and ∆ µν = g µν − u µ u ν is a projector onto the transverse four-velocity. The shear-and bulk viscosities are encoded into π µν and Π, respectively. Free parameters of this model include the initial conditions, η/s(T ) and ζ/s(T ), characterized by a total of 14 model parameters, which together control the prominent features of the model. The parameter set, described in detail in later sections, will enable simultaneous characterization of the initial state and medium response, including any correlations.
Each event consists of a single initial condition given energy density profile and a hydrodynamic simulation followed by multiple samples of the freeze-out hypersurface. The parameter estimation is conducted using 500 parameter design points, sampled evenly from the parameter space using the Latin hypercube scheme [30,31]. At each design point, the model is used to generate around 3×10 5 events with the corresponding parametrization, with each event surface sampled ten times to produce a total of 3 × 10 6 events for 0-60% centrality ranges. A large number of events is generated to ensure a better accuracy for high harmonic observables. A GP emulator is then trained to produce predictions for the observables in between the design points, after which the predictions are validated against a validation set. See [7] for details of the emulator. Using the emulator to produce predictions in continuous parameter space, the final posterior distribution is created using MCMC sampling.

B. Calibrating the model parameters
The parameter estimation attempts to calibrate the model parameters for the model to optimally reproduce experimental observables. With Bayesian methods, the optimal parameters are characterized by probability distributions for their true values. As given by Bayes's theorem, the probability for the true parameters x is The left-hand side is the posterior : the probability of x given the design X, computed observables Y , and the experimental data y exp . On the right-hand side, P (x ) is the prior probability, encapsulating the initial knowledge of x , and P (X, Y, y exp |x ) is the likelihood: the probability of observing (X, Y, y exp ) given a proposal x . The likelihood may be computed using the principal component GP emulators as where z = z (x ) are the principal components predicted by the emulators, z exp is the principal component transform of the experimental data y exp , and Σ z is the covariance (uncertainty) matrix. The covariance matrix encodes all the experimental and model uncertainties [32].
In the principal component space, the covariance matrix can be expressed as where Σ exp z is the matrix for experimental errors and Σ GP z = diag(σ z,1 (z ) 2 , σ z,2 (z ) 2 , . . . , σ z,k (z ) 2 ) is the diagonal GP emulator covariance matrix, representing the model statistical and GP predictive uncertainty. Additionally, σ sys m is a free parameter ranging from zero to one, with the purpose of including all remaining uncertainties arising from the model imperfections. All model parameters are given a uniform prior. Together with the likelihood (5) and Bayes' theorem (4), the posterior probability can evaluated at an arbitrary point in the parameter space. To construct the posterior distribution, an MCMC method can be used, which generates random walks through parameter space by accepting or rejecting proposal points based on the posterior probability.

III. EXPERIMENTAL OBSERVABLES
In the previous studies, the centrality dependence of identified particle yields dN/dy and mean transverse momenta p T for charged pions, kaons, and protons as well as two-particle anisotropic flow coefficients v n for n = 2, 3, 4 were used. The observables are measured by the ALICE Collaboration in Pb-Pb collisions at √ s NN = 2.76 TeV [2,8,9].   In this work, we mainly focus on the larger statistic higher beam energy collisions at √ s NN = 5.02 TeV, which give better precision. In addition to the above mentioned observables, we include higher harmonic flow coefficients v n [5] (up to n = 9), the normalized symmetric cumulants, NSC(m, n) [6], and the non linear flow mode coefficients χ n,mk [5]. The anisotropic pressure-driven expansion of the QGP, commonly referred to as anisotropic flow, can be characterized by a Fourier decomposition of the azimuthal particle distributions as where v n quantifies the magnitude of the n th harmonic flow, and ψ n its direction. NSC(m, n) quantifies the correlations between event-by-event fluctuations of flow harmonics of different orders, NSC(m, [3,6] and χ n,mk measures the contribution of lower order harmonic flows to higher order harmonics (i.e., χ 4,22 is the non linear contribution of v 4 originating from v 2 ; see the details in [5]). These additional observables give better sensitivity to the medium properties and initial conditions as demonstrated in Refs. [3][4][5][6]. As shown in Fig. 1 and 2, a model calculation with the best-fit parametrization given by maximum a posteriori (MAP) from the previous Bayesian analysis [7] shows deviations of the measurements for the flow coefficients from n = 5 and the non-linear flow mode coefficients from χ n,mk (n = 4). The black filled and open circles represent the higher ( √ s NN = 5.02 TeV) and lower ( √ s NN = 2.76 TeV) energy data points, respectively, whereas the red and orange bands represent the higher and lower energy model calculations. The v 2 − v 4 values calculated from data were used in model calibration, and as seen in Fig. 1, the calculations agrees well for v 4 . However, a discrepancy is seen for v 3 with an underestimation of the calculations for the centrality up to ≈ 45% for both energies, and an even larger discrepancy in v 2 for the higher energy calculation while the lower energy calculation agrees well except for the low centrality of 0 − 10%. For higher harmonics (≥ v 5 ) the deviation is still visible. The model calculations for the non linear flow mode coefficients in Fig. 2 agree within ±15% for χ 4,22 and +15% for the higher energy model calculation of χ 5,23 , while the lower energy model calculation goes to -30% in central collisions and even larger than 50% at high centrality ranges. The discrepancies between data and model calculations are significantly larger from χ 6,222 ; however, for χ 6,24 , most of the higher energy data points agree with the calculations within systematic uncertainties.
Model calculations reproduce the value for NSC(3, 2) up to the 40 − 50% centrality class, which is shown in Fig. 3. Both model predictions underestimate the values of NSC(4, 2) for all centrality classes presented. The model calculations overestimate NSC(4, 3) for the lower energy data and give similar results for the higher energy. However, the results show clear differences between the two beam energies. The differences get larger toward the central collisions. While it is negative for the lower energy data and the model calculations, the measurement at 5.02 TeV shows the change of the signs in central collisions. Also, the magnitudes are smaller in lower energy collisions, which is attributed to the increasing contribution from the nonlinear hydrodynamic response in v 4 [6].
In order to quantify the agreement of the models with the data, the χ 2 test was performed in the same way as in Eq. (5) in Ref [4] for the centrality range 5-50%. The results are shown in Fig. 4 for the flow coefficients, nonlinear flow mode coefficients, and the normalized symmetric cumulants. A significant difference is observed between the χ 2 values for v n of higher and lower energies at n ≤ 4. The χ 2 values for v 5 are larger for both beam energies with similar magnitudes. The higher energy χ 2 value for χ 4,22 is significantly larger than the one from the lower energy as shown in Fig. 4. The disagreement is still significant for χ 5,23 and χ 6,222 . For NSC(m, n), the χ 2 values are larger for higher harmonics at both beam energies. The χ 2 is especially large for the higher beam energy NSC(4, 3). In our calculations of the observables, we used the same methods also used in experimental analysis in Refs. [3][4][5][6]. Our centrality classes in this study are chosen to match those used for the experimental data. We define the multiplicity range for each centrality class by simulating events using the MAP-parametrization from [7], and sorting the resulting minimum-bias events by chargedparticle multiplicity dN/dη at midrapidity (|η| < 0.5). The identified dN/dη and p T were evaluated by counting and averaging the particle species at midrapidity (|y| < 0.5). The experimental data are readily corrected and extrapolated to zero p T [9], and therefore no additional processing is required while preparing the comparison. For the identified dN/dη, only protons were used in model calibration, as the model did not reproduce the spectra of the other species with any of the parametrizations. Finally, we calculated flow coefficients and other observables for charged particles within the kinematic range of the ALICE detector using the same methods as in the data analyses [5,6]. A summary of all the observables that are included in the Bayesian analysis is given in Table I. The table presents the particle species, kinematic cuts, and centrality classes for each observable.

IV. PARAMETER ESTIMATION USING NEW LHC MEASUREMENTS
The model to be evaluated in this analysis consists of multiple stages, of which a brief overview will be given next. Altogether, the model setup includes the parametric T R ENTo initial conditions, free-streaming Minimum η/s(T ) 0 -0.2 (η/s) slope Slope of η/s(T ) above Tc 0 -4 (η/s)crv Curvature of η/s(T ) above Tc pre-equilibrium dynamics, and the VISH(2+1) hydrodynamic model for medium evolution. Furthermore, the model performs the hadronization and includes UrQMD hadronic cascade. The model setup used is identical to the one developed and used in [7], except for the number of hypersurface samples taken after evolution. In this work, exactly ten events are sampled from the hypersurface regardless of the cumulative number of particles. The centrality definition is shared for all parametrizations. With close to fixed initial stage parameters, the possible effects of a shared centrality definition should be negligible.
Our main focus will be to investigate the effects of the higher harmonic observables on the temperature dependence of the transport coefficients. The parametrizations of the transport coefficients are [7] (η/s)(T ) = (η/s)(T c ) + (η/s) slope (T − T c ) T T c (η/s)crv for the ratio of shear viscosity and bulk viscosity over entropy, respectively. Based on previous work, it is known that the lowest value of η/s(T ) is around the critical temperature T c , close to the universal minimum 1/(4π). The temperature dependence of η/s(T ) is moderate, and increasing with higher values of temperature. Within close proximity of 150 to 500 MeV, the slope of η/s(T ) is approximately linear. The bulk viscosity over entropy ratio ζ/s(T ) is expected to peak around T c , and to decrease at higher values of temperature.
With this knowledge, we may construct our priori, and assume the initial parameter ranges. The chosen parameter ranges are loosely based on the optimal parameters found in [7]. It was found that in most cases, by taking the optimal parameters in [7] as the center points of the prior range and expanding the range slightly based on a reasonable σ value, those parameters could be further optimized with the additional observables. In this study, we have kept the initial stage parameter ranges narrow around the MAP values found in [7] with the assumption that the additional observables affect mostly the transport coefficients. Very small variation was allowed to give the parameters space to adjust for minor differences.
The included and varied parameters, of which there are 14 in total, are summarized in Table II. The parametric T R ENTo initial conditions comprise an ansatz in terms of five parameters: a normalization factor Norm, entropy deposition parameter p, standard deviation of the nuclear multiplicity fluctuations σ fluct , Gaussian-shaped nucleon width w, and minimum allowed distance between nucleons d min . The initial conditions are assumed to be already well constrained and presumably not affected by the addition of medium effect sensitive observables. The range for free-streaming time τ fs characterizing the allotted time for pre-equilibrium dynamics was kept relatively large.
The rest of the parameters are the components of the transport coefficient parametrizations, and the switching temperature T switch describing the temperature at which the hadronization begins to take place. The initial ranges given for these parameters are more generous, although large deviations in the final parameters compared to the previous study are not expected. The prior range for the transport coefficients is plotted in Fig. 5 among some parametrizations from other related studies [7,13,33,34]. The parametrizations are valid only up to the corresponding limits of the model; 100 MeV in the case of EKRT and 150 MeV for JETSCAPE. We note that the parametrizations EKRT+param0 and EKRT+param1 were not obtained through Bayesian analysis and we do not consider the slightly higher η/s at around T = 100 MeV in our prior. Furthermore, we do not consider the large ζ/s reported with the PTB particlization by the JETSCAPE Collaboration [13]. Nevertheless, the ζ/s obtained using the Grad or CE particlization distributions are within our prior, considering the temperature limit T switch > 150 MeV.
The model is calibrated to the latest Pb-Pb collision data at √ s NN = 5.02 TeV from the ALICE Collabora-

T(GeV)
FIG. 5. The specific shear (η/s) and bulk (ζ/s) viscosity ratios as a function of temperature. The region plotted in red visualizes the prior range used in this study. Other curves represent some of the parametrizations found in previous studies: the best-fit η/s(T ) with EKRT initial conditions [33,34], parametrizations from the JETSCAPE Collaboration with three different particlization distributions [13], and a recent parametrization found in [7].
tion [5,6,18,19]. Figures 6-9 show the calculations of each observable using the design parametrizations obtained from the prior distribution. The yellow curves represent the calculations corresponding to each design point parametrization, which are used in training the GP emulator, whereas the red curves represent emulator predictions corresponding to random points sampled from the posterior distribution.   V. RESULTS Figure 10 highlights the posterior and marginal distributions for select components of the transport parameters. The primary components, η/s slope, η/s(T c ), (ζ/s) max in the transport parametrizations are well constrained. The initial condition parameters are well constrained within the narrow prior range. Figure 11 presents the estimated temperature dependence of η/s(T ) and ζ/s(T ) according to the parametrizations from Eqs. (8) and (9), respectively. The shaded region around the curves represents the 90%credibility region. This region reflects all uncertainties coming from the finite width of the posterior distribution, experimental statistical and systematic uncertainties, statistical uncertainties in model calculations, predictive uncertainty from the GP emulator, and systematic model bias. With high probability, the true curve is located within this region. Table III presents the best-fit MAP parameters from our analysis. We list here the important findings: (i) While the temperature dependence of η/s(T ) is similar to what was obtained in [7], the curvature of η/s(T ) is slightly stronger, resulting in lower values at higher temperatures above T c .
(ii) A notable change is the lower (ζ/s(T )) max in order to reproduce the additional observables. The obtained ζ/s(T ) is smaller than those found in the previous Bayesian analyses [7,13] where the additional observables were not included. A similar value was reported in Ref. [11]. On average this represents a value an order of magnitude lower compared to the lattice QCD calculation [35] and the parametrizations used in [36,37], where the parametrizations were tuned to simultaneously reproduce lower harmonic v n as well as the charged particle multiplicity and the low-p T region of the charged hadron spectra.
(iii) The switching temperature on the other hand is higher than the one found in the aforementioned studies, where on average T switch is located around ≈ 0.150 GeV. As discussed in [4,5,17], the additional observables, the non linear response modes and the correlations between flow harmonics are sensitive to viscous corrections to the equilibrium distribution at hadronic freeze-out [38][39][40][41] and seem to prefer the higher switching temperature.

T(GeV)
Calibrated to PbPb s NN = 5.02 TeV   FIG. 11. The 90%-credibility region for the shear (top) and bulk (bottom) viscosity to entropy ratio is given as a blue band. The blue line represents the median of the credible range. The MAP parametrization from [7] as well as the corresponding 90%-credibility range are plotted as green dashed curves.
We performed high-statistics hydrodynamic calculations with the new parametrization. Figures 12 and 13 present the calculations for the flow coefficients v n and non linear flow mode coefficients, respectively. The v n is reproduced within 10% agreement for n = 2 up to n = 4. For the fifth harmonic, the calculations underestimate the data. The new parametrization estimates the data better in central and peripheral collisions but deviates significantly in the peripheral region. The magnitude of the successive harmonics from v 6 is not quite captured by the calculations within the statistical uncertainty. Furthermore, with our new parametrization, the   FIG. 14. Normalized symmetric cumulants [NSC(m,n)) from two hydrodynamical calculations are compared to the experimental data [6] at center of mass energy of 5.02 TeV. The blue band is calculated with the MAP-parametrization from this work, whereas the red band uses the parametrization from [7].
also improved compared to the parametrization from [7], as indicated by the ratio plots. In this case, the lower harmonic non linear flow mode coefficients are no longer overestimated, and the magnitude and centrality dependence are correctly captured. We note that the nonlinear flow mode coefficients have not been included in the model calibration in [7], whereas, coefficients up to χ 6,33 and χ 6,222 were used in this analysis. Figure 14 presents the calculation of the normalized SC using our obtained parametrization. The performance of the new parametrization and the one from [7] are comparable for NSC (3,2). For NSC (4,2) and NSC(4,3), the centrality dependence is better described by the new parametrization. However, both parametrizations are unable to reproduce the strong centrality dependence of NSC(4,2), underestimating the most data points in the most peripheral collisions. The multiplicity and the mean p T calculations are compared to the results from [7] in Fig. 15. Our parametriza- tion improves the estimate of the proton multiplicity and gives the same charged particle multiplicity for 5.02 TeV collisions, while the pion and kaon multiplicities are not in good agreement with the experimental data, as similarly found in [7] for 2.76 TeV calculations. Interestingly, the parametrizations from [7] mainly utilizing 2.76 TeV data give better agreement with pions and kaons in 5.02 TeV collisions than our results while overestimating the proton yields approximately 10%. Agreement of the calculated mean p T with the experimental data is good for all particle species as well as with the results from [7] for both beam energies. Refining this analysis by including low beam energy data in the future will help us to understand the beam energy dependence on various observables.
Finally, Fig. 16 shows the χ 2 values with the bestfit MAP parameters extracted from this study for each observable. They are compared to the ones from [7]. The χ 2 values for our new calculation only seem to improve v 3 and v 5 for the v n observable. For the non linear flow mode coefficients the χ 2 values are improved up to χ 6,222 with our new parametrization, while the χ 2 values for the higher harmonics are worse than in [7]. For NSC, the χ 2 from the new calculation is worse for NSC (3,2) and NSC(4,3), but is improved for NSC (4,2). We note that the larger statistical error in the calculations using the parametrization of Ref. [7] lowers the corresponding χ 2 values, slightly affecting the direct comparison between the two parametrizations. The sign change of NSC (4,3) in most central collisions is not reproduced by the models while the beam energy-dependent magnitudes are better described with new parametrizations. We leave those differences for future research work where the present results should be refined by including experimental data from the lower energy beam data.
As a final study in this analysis, we conduct a simple sensitivity analysis of the included observables to the model transport parameters. The sensitivity of each observable is evaluated using the GP emulator by observing the relative difference in the magnitude of the observable between two parameter points in the parameter space. The difference can be formulated as whereÔ(x) andÔ(x ) represent the values of an observable at parameter points x and x , respectively [13]. In this study, we choose a reference parameter point x to be the one representing the MAP values obtained in this analysis (see Table III). To probe the sensitivity of a parameter j, another point is defined as x = (x 1 , x 2 , . . . , (1 + δ)x j , . . . , x p ), where δ is a small value representing a percentile change in the parameter space. We have used a value δ = 0.1, although larger values were observed to yield similar results.
We then calculate a final sensitivity index for each observable and parameter pair in various centrality classes as Figure 17 presents the evaluated sensitivity for each observable against the transport parameters. The sensitivity was evaluated over four centrality classes from 5% to 40% and averaged for the final plot. We did not observe large differences in the sensitivity between the individual centrality classes.
For v n , we can verify a known fact that the sensitivity of the flow coefficients is generally very limited to the temperature dependence of η/s(T ) [33], although as expected, the sensitivity to the average η/s , in this case, represented by η/s(T c ), is very strong, and increasing at higher harmonics. The sensitivity of the v n to the (ζ/s) peak is visible, and also in this case the higher harmonics provide stronger constraints. Based on previous studies, the non-linear flow mode coefficients χ n,mk are known to be sensitive to η/s(T ) at the freeze-out temperature. This is reflected by the observed sensitivity to η/s(T c ) as well as T c . By far, the normalized symmetric cumulants provide the strongest constraints to the temperature dependence of η/s(T ). This is confirmed by higher sensitivity for the other components of η/s(T ), and not only η/s(T c ), which is also higher.
Two other parameters have also been included in this study: the free-streaming time scale τ fs and the switching temperature T switch . On average, the observables are reported to be generally weakly sensitive to τ fs , apart from the symmetric cumulants and χ 6,33 . Furthermore, most of the observables, such as v n , χ 6,mk and the NSC(m,n), are seen to be highly sensitive to the switching temperature T switch . In both cases, the results reported here regarding τ fs and T switch are not compatible with what has been observed in [13]. 18. Sensitivity of the mean multiplicity yields and mean transverse momenta pT to the model parameters. Figure 18 presents the sensitivity of the multiplicity and the p T to the model parametrizations. Most prominently, the switching temperature affects the proton multiplicity. Furthermore, we observe a comparatively large sensitivity of p T to the free-streaming time scale τ fs . In the case of the transport parameters, the effect on the observables is relatively small. It is observed that p T acts as a subtle constraint to the parameters describing the specific bulk viscosity. The posterior distribution for all parameters is shown in Fig. 19.

VI. DISCUSSIONS
In summary, we performed a Bayesian analysis with the recently available data from ALICE Collaboration [5,6,18,19] as an extension of the work [7]. We found that the temperature dependence of η/s(T ) is similar to what was obtained in [7] and that the curvature of η/s(T ) above T c is slightly lower at higher temperatures, showing weak temperature dependence of η/s. Notable changes include the lower (ζ/s(T )) max and the higher switching temperature T switch to reproduce additional observables such as symmetric cumulants and non-linear flow coefficients. However, the improved statistical uncertainties on both the experimental data and hydrodynamic calculations do not help to reduce the final credibility ranges. It is also noticeable that v 5 is still underestimated as observed in [7]. It is worthwhile to mention that the differences for v 2 , v 3 and NSC(4,2) still remain about 5-10% for 5.02 TeV. The sign change of NSC (4,3) in most central collisions is not reproduced by the models while the beam energy-dependent magnitudes are better described with new parametrizations. We leave those differences for future research work in which the present results should be refined by including the lower energy beam data. The parameter sensitivity analysis for the observables conducted in this study indicates that observables such as the symmetric cumulants and non linear flow modes provide a strong constraining power which, however, is still underutilized in [7] as well as the other Baysian analyses [10,11,13]. In our study, we confirm that the flow coefficients alongside the symmetric cumulants and non linear flow mode can provide some of the strongest constraints for the temperature dependence of η/s(T ) and T switch . Improving aspects of the collision model, for example by replacing the initial state model with others like EKRT [33,34], IP-Glasma [42], and AMPT [43,44] with incorporation of nucleon substructure [45] in the initial conditions through an improved dynamical collision model before the hydrodynamic takes place [46,47], might help to improve the understanding of the uncertainties of the extracted QGP properties and/or the model building blocks.