Investigating the collision energy dependence of $\eta$/s in RHIC beam energy scan using Bayesian statistics

We determine the probability distributions of shear viscosity over the entropy density ratio $\eta/s$ in Au+Au collisions at $\sqrt{s_{NN}}=19.6, 39$, and $62.4$ GeV, using Bayesian inference and Gaussian process emulators for a model-to-data statistical analysis that probes the full input parameter space of a transport+viscous hydrodynamics hybrid model. We find the most likely value of $\eta/s$ to be larger at smaller $\sqrt{s_{NN}}$, although the uncertainties still allow for a constant value between 0.10 and 0.15 for the investigated collision energy range.


I. INTRODUCTION
After successfully verifying the existence of new state of matter, the quark-gluon plasma (QGP), experiments at RHIC and LHC are now focused on quantifying the transport properties of the QGP and mapping out the QCD phase diagram.The beam energy scan (BES) program at RHIC in particular is aimed at probing the phase diagram of QCD matter at varying values of the baryon chemical potential.It is known that at LHC and top RHIC collision energies, the phase transition between QGP and hadronic matter is a cross-over [1,2].However, it is not yet known if there exists a critical point marking the change to first-order phase transition at lower collision energies, where the QCD matter is exposed to typically lower temperatures but higher net-baryon densities.
The beam energy scan has introduced new challenges to the theoretical modeling of relativistic heavy-ion collisions, as the initial non-equilibrium evolution of the system leading up to thermalization gains in importance and has a large effect on the outcome of the calculation.In addition, the assumption of a midrapidity plateau does not hold any longer, requiring a full (3+1)-D hydrodynamical model to describe the evolution of the medium.
Finally, the equation of state needs to accommodate the non-zero baryon chemical potential µ B .To address these challenges, we utilize a hybrid model, which combines a (3+1)-D viscous hydrodynamics model for the QGP phase of the reaction with a microscopic hadronic transport model, that describes the non-equilibrium evolution of the system before QGP formation and after its subsequent decay into final state hadrons.
In recent years, a lot of progress has been made in the field regarding the quantification of the shear viscosity over entropy ratio η/s in QCD medium.In particular, the temperature dependence of η/s has been studied in detail [3][4][5][6][7].However, a recent study [8] has found that η/s might depend also on baryon chemical potential µ B , as earlier calculations had already found to be the case in the hadron resonance gas [9,10].
Inspired by the success of Bayesian analysis at LHC energies [6], we investigate in this article the possible collision energy dependence of η/s by performing a similar analysis on the input parameters of the hybrid model used in Ref. [8].We describe the main features and the input parameters of the hybrid model in Section II.A summary of the statistical analysis procedure is provided in Section III and the details of the simulation setup are found in Section IV. Results of the analysis are presented in Section V. We summarize our findings in Section VI.

II. HYBRID MODEL
The hybrid model utilized in our study is the same as in Ref. [8].In this model, the initial non-equilibrium evolution of the system is modeled with the UrQMD hadron-string cascade [11,12], which is followed by a transport-to-hydro transition ("fluidization").The earliest possible starting time for this transition is when the two colliding nuclei have passed through each other: . At the transport-to-hydro transition the microscopic particle properties (energy, baryon number) are mapped to hydrodynamic energy, momentum and charge densities using 3-D Gaussians with "smearing" parameters W trans , W long (≡ √ 2σ).
From the hydrodynamic side, cells with high enough local energy density correspond to regions of matter in the QGP phase.
The hydrodynamic evolution is modeled with (3+1)D relativistic viscous hydrodynamics [13].For simplicity we utilize a constant value of η/s throughout the QGP evolution, which is an input parameter of the model.In an event-by-event hydrodynamic calculation at the given energy range, temperature as well as baryon chemical potential vary significantly in space-time.Therefore an equation of state which covers the large-µ B sector is a necessary ingredient for the hydrodynamic stage.For the present study a chiral model equation of state [14] is used.In this equation of state the transition between hadron and quark-gluon phases is a crossover for all µ B .
The transition from hydrodynamics back to microscopic non-equilibrium transport ("particlization") is performed when the local rest frame energy density in the hydro cells falls below a specified switching value SW .The iso-energy density hypersurface is constructed using the Cornelius routine [15].At this hypersurface hadrons are sampled via Monte Carlo procedure using Cooper-Frye prescription and including viscous corrections to their distribution functions.The generated hadrons are then fed into the UrQMD cascade to simulate final state hadronic rescatterings and resonance decays.

III. STATISTICAL ANALYSIS
The basics of Bayesian analysis procedure and model emulation have been described in detail in [16,17].In this section we summarize the main points and describe the details specific to the current implementation.

A. Bayes' theorem
Given a set X = { x k } N k=1 of points in parameter space (design points), and a corresponding set Y = { y k } N k=1 of points in observable space (output points), the Bayes' theorem states that the conditional probability distribution (aka posterior) of an input parameter combination x * , given the observation (X, Y, y exp ), is where P ( x * ) is the prior probability distribution of x * , and P (X, Y, y exp | x * ) is the likelihood of the observation (X, Y, y exp ) for a given x * .
As in Ref. [17], the likelihood in our analysis will be of the form where y * is model output for the input parameter point x * , y exp is the target value and Σ is the covariance matrix.

B. Gaussian emulator
To calculate the likelihood, one must be able to determine the model output y * for an arbitrary x * .As the simulations typically take several hours to run, it is highly impractical to run the full model during statistical analysis.Instead, the model is emulated with a Gaussian process.That is, the model is taken as a mapping of a set X of points in the parameter space to a set Y a of values of the observable y a , where Y a is assumed to have a multivariate normal distribution with µ(X) = {µ( x 1 ), ..., µ( x N )} as the mean and as the covariance matrix with covariance function σ( x, x ).The covariance matrix is further reformulated in terms of a triangular matrix S using the Cholesky decomposition Σ = SS T .
A sample from a Gaussian process is then of the form where the components of u are normally distributed random variables, u i ∼ N (0, 1).
To predict the model output y 0 in an arbitrary point x 0 , one writes the joint distribution and calculates the conditional predictive mean Typically we define the mean function µ( x) ≡ 0 and the prediction is simply with the associated predictive (co)variance Clearly the crucial component of the Gaussian process is the covariance matrix Σ.The covariance function σ( x, x ) can have any form from a simple constant to an exponential or a periodic function, and the right choice depends on the features of the emulated model.For this study we have chosen the squared-exponential covariance function with a noise term The hyperparameters θ = (θ 0 , θ 1 , ...θ n , θ noise ) are estimated from the provided data (X, Y ).
This "training" of the emulator is done by maximising the multivariate Gaussian log- In principle, using maximum likelihood values for hyperparameters is only an approximative method, as in the rigorous Bayesian approach the hyperparameters are allowed to change during calibration procedure, resulting to posterior distributions also for the hyperparameters.However, in practice 100 training points creates strong enough constraints on the Gaussian processes, that the precise value of the hyperparameters has little effect on the final results.

C. Principal component analysis
Generally, the analysis requires one emulator per data point.However, the number of emulators can be reduced with a principal component analysis, which is an orthogonal linear transformation of the data onto the directions of maximal variance.
We use our N simulation points and m observables to form a N x m data matrix Y .Data is normalized with the corresponding experimental values to obtain dimensionless quantities, and centered by subtracting the mean of each observable from all points.The singular value where S is a diagonal matrix containing the singular values, and U and V T are orthogonal matrices containing the left-and right-singular vectors, respectively.Eigenvalue decomposition of the covariance matrix Y T Y becomes Singular values in S are square roots of eigenvalues λ i , and right singular vectors v i represented by columns in V T are eigenvectors of Y T Y ; these are the principal components (PCs).The eigenvalues are proportional to the total variance of the data, such that the eigenvalue of the first principal component explains the largest fraction of the variance.
We can now write our data matrix in principal component space: Since λ 1 ≥ λ 2 ≥ ... ≥ λ m , the fraction of the total variance explained by the qth principal component, λ q /( m j=1 λ j ), becomes negligible starting from some index q < m.This allows us to define a lower-rank approximation of the original Z as Z q = √ N Y V q .
In addition to dimensional reduction, we perform data whitening, i.e. divide each principal component with their singular value, such that they all have unit variance.The full prescription for transforming a vector y from observable space to vector z in whitened principal component space is thus defined as To compare an emulator prediction z * against physical observables, we use inverse transformation

D. Likelihood function
The log-likelihood function used for determining the posterior probability distribution in this analysis has the following form: where z * is the emulator prediction at the input parameter point x * and z exp represents the experimental data transformed to principal component space.
The covariance matrix includes both experimental and emulator uncertainties: The uncertainty of experimental data Σ exp z is obtained by transforming the observable space covariance into principal component space: 1 In the previous Bayesian beam energy scan analyses [18][19][20], the experimental uncertainties were treated as simple diagonal contributions in the principal component space: Σ exp z = diag((σ 0 ) 2 , .., (σ q ) 2 ), where σ = (σ 0 , ..., σ m ) Ṽq .However, this neglected the covariances between the principal components and exaggerated the variance for the first principal component.The uncertainties were overestimated as a net result.
For many of the observables, the reported values σ 0 , ..., σ m include both statistical and systematic uncertainties.In cases where the uncertainties are reported separately, they have been added quadratically: σ 2 i = σ 2 i,stat + σ 2 i,sys .In principle, systematic uncertainties introduce also off-diagonal terms into Σ exp y .However, the task of estimating these offdiagonal terms and investigating their effect on the final results requires detailed knowledge of the detector systems in the experiments and is referred to a future study.
The emulator uncertainty is included as where σ 2 zi refers to the GP predictive variance (Eq.( 10)) for ith principal component.Finally, the term containing the determinant of the total covariance matrix |2πΣ z | is included to provide more value for points which both fit the data and have small uncertainties, as opposed to points which fit the data but have large uncertainties.

E. Data calibration using random walkers
The posterior distribution is sampled using the Markov chain Monte Carlo (MCMC) method, which is a random walk in parameter space, where each step is accepted or rejected based on a relative likelihood.It converges to posterior distribution at the limit N steps → ∞.
An important measure for the quality of the resulting random walk is the acceptance fraction, which should be between 0.2 and 0.5 [21].A very low value of the acceptance fraction indicates that the walkers are stuck, while values close to 1 mean all proposals are accepted and the walk is purely random.Another measure is the autocorrelation time, which is an estimate of the number of steps required between independent samples.This can be used to evaluate if the random walk contained enough steps to produce a proper sample size.
We initialize O(1000) random walkers in random positions in the input parameter space.This is followed by the "burn-in" phase, which allows the random walkers to converge to the high-likelihood areas.At the end of the burn-in, the samples produced during this phase are removed and the proper sampling phase is initialized.
Once the MCMC is complete, we verify that we had sufficient number of steps in the burn-in phase to cover a couple autocorrelation lengths, and enough steps in the sampling phase to cover ∼ O(10) autocorrelations.Having thousands of walkers ensures a sufficiently large sample of independent observations for generating posterior distributions.

IV. SIMULATION SETUP A. Input parameters
The parameters of the computational model encode the physical properties of the system that we wish to extract from our analysis.These are the specific shear viscosity of the QGP, η/s, the thermalization time of the system τ 0 , the initial condition smearing parameters W trans and W long as well as the hydro-to-transport transition energy density that can provide information on the temperature scale at which non-equilibrium effects in the hadronic evolution become important.The following parameter ranges were used for the design points: As described in Section II, τ 0 is constrained by the condition that the two nuclei must have passed each other, and thus the lower limit depends on collision energy: The cells in the hydro grid have a transverse length ∆x = ∆y = 0.25 fm and a longitudinal length ∆η = 0.15, which sets a soft lower limit for the smearing parameters, as the limited resolution will lead to the "over-smearing" of the Gaussians much narrower than the cell size.
As in Ref. [17], the input parameter combinations were sampled using maximin Latin hypercube method, which attempts to optimize the sample by maximizing the minimum distance between the design points.

B. Normality of the data
Principal component analysis assumes that variance provides a reasonable measure of the spread of the simulation output values, i. e. the simulation data is approximately normally distributed.However, this is not generally true for the distributions of observable values calculated from the model output.We aim to normalize the distributions by scaling the output with the corresponding experimental data and applying the Box-Cox transformation [22] on these dimensionless values ỹ = y/y exp : We check the normality of the resulting distribution with a probability plot, which is a method for comparing the quantiles of the sample data against the respective quantiles of the normal distribution.We perform a least-squares fit on the points in probability plot: the goodness of the fit can be estimated by the coefficient of determination R 2 , which should be close to 1.We use R 2 c = 0.91 as a (very lenient) critical test value for too large deviations from normal distribution.Observables with R 2 < R 2 c are excluded from the analysis.

C. Event selection
Centrality bins were determined based on the number of participants N part , such that the average number of participants N part in each bin matches the respective experimental estimate [23][24][25][26][27].
To save computation time, the centrality selection was done by running the full simulation for the 50% of the initial UrQMD events with highest number of participants N part , instead of running the full simulation for all initial states and sorting the events by final charged particle multiplicity N ch .In addition, each switching energy density hypersurface has been used 50 times for particle sampling ("oversampling") to increase statistics.As the number of independent hypersurfaces in any centrality bin is at least an order of magnitude greater, this method is unlikely to produce any notable bias in the final results.
V. RESULTS

A. Emulator verification
The quality of the GP emulator was verified against 6 randomly selected test points from the simulation output data, which were not included in the analysis.As an example, Fig. 1 shows the verification results on π + yields at

B. Presentation of analysis results
The Bayesian posterior probability distributions are illustrated in Figs.The posterior distribution for √ s N N = 19.6GeV is presented in Fig. 2. While the data supports the earliest possible starting time and longitudinally narrow Gaussians as density representations of the initial particles, the transverse plane of the initial density profile is strongly smeared.An early transition from hydro to hadron transport is preferred, which is in contrast to the earlier studies where Ω baryon was given more weight in the analysis [18][19][20].
We can investigate the quality of statistical analysis in more detail by looking the results The recently published STAR data on identified particle yields and mean transverse momenta [25] finally makes it possible to rigorously test the bulk production capabilities of relativistic heavy ion collision models over the whole beam energy scan range.However, as the proton and antiproton data suffers from additional uncertainties related to feed-down corrections, we limit ourselves to pion and kaon data in this analysis.
The results for positive-charge pion and kaon multiplicities are shown in Fig. 5. Pion yields are well reproduced for the investigated centrality range.There is a tendency to overestimate kaon yields, especially at the more central collisions.The median value result is still within statistical uncertainties, however.
To convince us that the model has captured the correct dynamical evolution of the system, it is important to check not only particle yields but also the transverse momentum spectra dN/dp T .It was postulated in Ref. [16] that the combination of integrated yield and mean p T is sufficient to describe the p T spectrum.Results for mean transverse momenta of π + and K + are shown in Fig. 6.The pion transverse momenta is slightly below reported value for all centralities, but still within error bars, whereas the agreement with kaon data is very good.Overall pion and kaon spectra are well reproduced, as shown in Fig. 7, although the kaon excess in the most central collisions is evident also on this figure.
Two-pion interferometry, commonly referred as "Hanbury-Brown-Twiss analysis", is considered as a probe of the spacetime geometry of the kinetic freezeout region [28].Results for the HBT radii are summarized in Fig. 8 for the lowest k T bin, which corresponds to the observable used for the comparison over several different experiments in Ref. [29].Again following the example of Ref. [16], we also investigate mean radii averaged over several k T bins in Fig. 9.
Using the median values for input parameters, the model is able to match R out and R long quite well, but R side is systematically underestimated.The same is true for the averaged radii.Emphasizing the importance of HBT radii has thus potential to impose strong constraints on the input parameters.
The centrality dependence of v 2 {2} for charged particles is presented in Fig. 10.Agreement with the STAR data is very good for the (5-10)% and (10-20)% centralities; in the (20-30)% centrality bin, the median input parameters produce slightly less elliptic flow compared to observations, which may indicate an issue using the same value of τ 0 for all centralities.The lower limit of τ 0 , the time for the two nuclei to pass each other, becomes smaller at higher centralities, in principle enabling earlier hydro starting times.This might also effect W long which tends to have strong correlation with τ 0 .We leave the investigation of the effect of centrality-dependent τ 0 for future work.
Finally, we compare the transverse momentum spectrum of Ω against the STAR data in Fig. 11.It is evident from the figure that Ω yield is overestimated when using the median values of the posterior distribution.However, the fit here and at 39 GeV was done using the Experimental data from [30].
1-5 GeV/c, which was calculated from a Levy fit on the p T spectra [18].Based on results at 62.4 GeV (Fig. 29), using the fully p T -integrated yield and p T is expected to produce better results.Once such data becomes available, its effect on the posterior distribution needs be investigated.The experimental data available at √ s N N = 39 GeV is a subset of that at 19.6 GeV; information about the pseudorapidity dependence of charged particle yields is not available.
As the earliest possible starting time remains favored, the τ 0 distribution shifts below 1.0 fm, and the other distributions also peak at lower values, in particular the switching energy density SW .The results for the actual observables are very similar to √ s N N = 19.6GeV.Pion and kaon yields are again well reproduced (Fig. 13), as is their mean p T (Fig. 14), resulting in a good agreement with the measured transverse momentum spectra (Fig. 15).Regarding the HBT radii (Figs.16 and 17), R out is still well described with peak parameter values.However, in addition to the underestimation of R side , we find R long to be overestimated at this energy.As R long has been associated with the freezeout time and and temperature of the system [32], and should also be sensitive to shear viscosity [33], increasing GeV.Experimental data from [25].
the weight of HBT observables at this energy might have an effect on the most probable values of η/s and SW .Also the results for v 2 {2} (Fig. 18) and Ω (Fig. 19) support the conclusions drawn from √ s N N = 19.6GeV results.Here the agreement with experimental data is good also for (20-30)% centrality, as the initial allowed hydro starting time is lower and likely appropriate for wider range of centralities.
The transverse momentum spectrum of Ω is in somewhat better agreement with the STAR data at √ s N N = 39 GeV compared to √ s N N = 19.6GeV.The total yield is still overestimated, however.Experimental data from [30].
The highest collision energy investigated in this analysis has also the most comprehensive set of observables, including information about the p T spectra of identified particles for several centralities.Instead of calibrating the model on the yield of pions and kaons separately, the K + /π + ratio was used instead.As the majority of charged particles consists of pions, calibrating the model on the centrality dependence of N ch was deemed sufficient to constrain the total number of pions when combined with the kaon/pion ratio.
The posterior probability distribution for 62.4 GeV, shown in Fig. 20, shows similar trends as the distributions at the lower energies.Early starting time for the hydrodynamic evolution is preferred, together with a longitudinally narrow Gaussian smoothing for the initial density.
The transverse width of the smearing functions has narrowed further, leading to W trans and W long being roughly equivalent for this energy.The shear viscosity over entropy density ratio has a clear peak at 0.1; overall the posterior probability distribution of η/s does not change much over the investigated beam energies.Like at 19.6 GeV, and unlike at 39 GeV, a high value of the switching energy density is preferred.As the pseudorapidity distribution is the one observable which is included in the data calibration at 19.6 GeV and 62.4 GeV, but not included in the analysis at 39 GeV, it seems that the early transition from hydrodynamics to hadron transport is mainly driven by the longitudinal expansion.
As illustrated in Fig. 21, the simulation using the peak parameter values gives a good description of the centrality dependence of charged particle multiplicity at midrapidity, especially considering the wide range of N ch in the prior range.
The full charged particle pseudorapidity distributions for (0-3)% and (20-25)% centrality at √ s N N = 62.4 GeV are shown in Fig. 22.The systematic difference between emulator prediction and actual simulation, observed at √ s N N = 19.6GeV, is present also at this higher collision energy.However, here the emulator prediction is systematically overestimating the midrapidity yield, in contrast to the systematic underestimation seen at 19.6 GeV.
The simulations run with median values agree with PHOBOS data overall quite well at (0-3)% centrality.At (20-25)% centrality, the charged particle yields are somewhat below the data both at η ≈ 0 and again at |η| > 2.0.The mean pseudorapidities remain below the measurement, as shown in Fig. 23.At this energy, it appears possible to match PHOBOS result within the parameter prior.Based on the shape of dN ch /dη from Fig. 22, matching η is likely to lead to a notable depletion of charged particles at η ≈ 0, however.
Figure 24 shows that there is some tendency for the model to overestimate the abundance of kaons over pions.Overall there is very little room variation even in the prior range of the particle ratio, so the shape of posterior will have little effect on this observable in any case.Despite the observed slight excess of kaons compared to pions, the mean transverse momenta of π − and K + is reproduced very well with the median parameter values and the p T spectra for π − , K + agrees well with PHOBOS results, as shown in Figs. 25 and 26.
The HBT analysis results presented in Figs. of HBT radii at 39 GeV and 62.4 GeV is very similar, despite the large difference in the applied value of the switching energy density, suggests that R long is more sensitive to η/s than SW .
Finally, we show the transverse momentum spectra of Ω and proton in Fig.  STAR data from [25] and PHOBOS data from [34].
kaon multiplicities are very well reproduced, we can conclude that values of SW ≈ 0.4 − 0.5 GeV/fm 3 would be preferred for strange particle production.Also the proton spectrum, which was not part of the data calibration, matches quite well with the PHOBOS data, although slight overabundance of protons is evident from the figure.As the collision energy dependence of these smearing parameters lack an obvious physical motivation, it becomes necessary to impose more stringent theoretical constraints on the initial state for the future studies.

FIG. 1 :
FIG. 1: (Color online) The ratio of the Gaussian process prediction over the model result for positively charged pion yields, for the 6 test points at √ s N N = 19.6GeV.Top panel: Untrained emulator.Bottom panel: After emulator training.
2, 12, and 20.We visualize these 5-dimensional probability distributions using corner plots.The histograms on the diagonal panels represent the probability distributions of each input parameter integrated over the possible values of the other four parameters, while the off-diagonal panels show pairwise correlations between the parameters.Model results after the calibration are compared to experimental data in the following way: To illustrate how the uncertainty in parameter values reflects to the observables, 100 random parameter combinations were drawn from the posterior and the model outputs for these combinations were predicted by the Gaussian process emulator.These emulator predictions are drawn as box-and-whisker plots.The single red horizontal line within each box represents the median prediction; that is, half of the predictions are smaller and half are larger than this value.Boxes represent the values between first and third quartile Q 1 and Q 3 ; this means 25% of the predictions were smaller than the values in the box and 25% were larger.Finally, the whiskers represent either the complete range of predictions, or 1.5 times the interquartile range Q 3 − Q 1 , whichever is smaller.In the latter case, the predictions outside the whisker range ("outliers") are represented by individual points.Thus, the boxand-whisker representation illustrates the shape of the distribution of emulator predictions.To verify the quality of the calibration, full model simulations were run with the median values of the one-dimensional projections of each posterior probability distribution, represented as dashed red lines and red triangles in figures 2, 12, and 20.C.√ s N N = 19.6GeV

8 WFIG. 3 :
FIG. 2: (Color online) Representation of the posterior probability distribution of the input parameters at √ s N N = 19.6GeV.Panels on a diagonal show one-dimensional projections of the posterior for each parameter axis, while the off-diagonal panels display projections on the different 2-D planes of the input parameter space.Dashed red lines indicate median values of the one-dimensional distributions.

FIG. 4 :
FIG. 4: (Color online) Centrality dependence of average charged particle pseudorapidity η for η > 0.0 at √ s N N = 19.6GeV.Blue box-whisker plots are GP estimates for samples from the posterior distribution.Black open squares are emulator predictions for the model output at distribution median values (red triangles).Dashed grey boxes represent the range of values from the training data.Experimental values are estimated from PHOBOS data [24].

8 W
FIG. 12: (Color online) Representation of the posterior probability distribution of the input parameters at √ s N N = 39 GeV.See caption of Fig. 2 for details.

8 W
FIG. 20: (Color online) Representation of the posterior probability distribution of the input parameters at √ s N N = 62.4 GeV.See caption of Fig. 2 for details.

FIG. 32 :FIG. 33 :
FIG. 32: (Color online) Collision energy dependence of the smearing factor W trans .Error bars represent 90% confidence range around the median value (filled triangles).Open symbols indicate the peak position of the distribution.
29.In contrast to the results at lower collision energies, fully p T -integrated Ω multiplicity and mean transverse momentum were used in the data calibration, leading to a good agreement with measured p T spectra; the total multiplicity is slightly overestimated.As we observe an excess of strange particles both at 19.6 GeV and 62.4 GeV, while at √ s N N = 39 GeV the

TABLE I :
Calibration data at√s N N = 19.6GeV.