Parameterizing Smooth Viscous Fluid Dynamics With a Viscous Blast Wave
PParameterizing Smooth Viscous Fluid Dynamics With a Viscous Blast Wave
Zhidong Yang ∗ and Rainer J. Fries † Cyclotron Institute and Department of Physics and Astronomy,Texas A&M University, College Station TX 77843, USA (Dated: August 17, 2020)Blast wave fits can capture essential features of global properties of systems near kinetic equilib-rium. They usually provide temperature fields and collective velocity fields on a given hypersurface.We investigate how faithful the viscous blast wave introduced in [1] can reproduce the given tem-perature and specific shear viscosity at freeze-out of a viscous fluid dynamic calculation, if the finalspectrum and elliptic flow of several particle species are fitted. We focus here on fluid dynamicsimulations appropriate for high energy nuclear collisions at current collider energies. We find thatspecific shear viscosities are reproduced to good accuracy by viscous blast wave fits while tempera-tures tend to be slightly underpredicted. We quantify the deviations of fitted from true quantitiesfor some examples. The maps we obtain can be used to improve raw results obtained from viscousblast wave fits. ∗ [email protected] † [email protected] a r X i v : . [ nu c l - t h ] A ug I. INTRODUCTION
Blast waves are simple and effective tools that provide snapshots of a system that is close enough to local kineticequilibrium so that macroscopic concepts like temperature, local collective velocity, and shear and bulk stress can beused. The full dynamics of such systems is usually very well captured by viscous fluid dynamic equations of motion.If the equation of state and a sufficient number of transport coefficients (e.g. shear and bulk viscosity) of the systemare known, fluid dynamics can evolve the system starting from a given initial state. In contrast, blast waves usuallyprovide a static picture, typically the temperature field T ( x ) and collective flow field u µ ( x ) on a hypersurface Σconsisting of points r µ . E.g., if the hypersurface of kinetic freeze-out in an expanding system is chosen, these fieldscan often be directly related to observables of the system and can be fitted to data.In high energy nuclear collisions, blast waves have routinely been used to analyze the properties of the fireball ofhadrons at the time of kinetic freeze-out [2–6]. After kinetic freeze-out hadrons free stream to the detectors and thusdirectly carry information about the freeze-out hypersurface. Blast wave parameters like the freeze-out temperatureand average radial flow velocity can be determined by fits to transverse momentum spectra of hadrons. Additionalinformation, like elliptic deformations of the fireball in coordinate and momentum space at the time of freeze-outcan be extracted from fits of elliptic flow coefficients v . More recently, viscous corrections to blast waves have beenconsidered [1, 7–9]. They can be used to extract the specific shear viscosity η/s of hadronic matter at the freeze-outtemperature [1]. Viscous blast waves can also be used to extend the range of validity of ideal blast wave fits to largertransverse hadron momenta P T .The question arises to what extent blast waves, which use certain simplifying approximations, can faithfully repro-duce important properties of the full dynamical system. The blast wave used here is defined by an ansatz for thetemperature, flow field and space-time structure of the hypersurface, as discussed below. The ansatz contains severalparameters with physical meaning, like the freeze-out temperature and specific shear viscosity which can be fitted todata. It is important to understand and quantify the uncertainties in the fit results that are due to the simplificationsmade in the blast wave ansatz. Once quantified, the systematic deviations could be corrected for and the extractionof important parameters from blast waves could turn into a tool with improved precision. In this work we take a firststep in this direction by comparing freeze-out conditions from smooth viscous fluid dynamics to blast wave fits. Theglobal picture is one of two successive approximationsreal collision system (cid:31) fluid dynamic simulation (cid:31) blast wave fitwhere (cid:31) means ”approximated by”. Viscous fluid dynamic simulations are widely accepted to give accurate descrip-tions of the low transverse momentum region of high energy nuclear collisions, although analyses of experimental datawith fluid dynamic simulations continue to be an active field of study [10–12]. There are several approximations thatenter when describing the evolution of the system and its freeze-out with fluid dynamics. Some relevant ones arebriefly discussed in [1], but their quantification is outside the scope of this work. Here we focus on the second step,and quantify uncertainties when approximating fluid dynamic systems by blast waves at freeze-out. We use smoothrelativistic viscous fluid dynamics to create systems close to local equilibrium as they typically occur in the late stagesof high energy nuclear collisions. Approximations used by deploying blast waves are discussed and then quantifiedin a specific case. We focus here on the accuracy of the extracted freeze-out temperature T fo and the specific shearviscosity η/s , i.e. the ratio of shear viscosity η and entropy density s , at freeze-out. The results from the analysishere can be used directly to improve the extraction of η/s from experimental data in [1].In our work we use the viscous generalization of the Retiere-Lisa (RL) blast wave [4] introduced in [1]. We utilizethe viscous fluid dynamics code MUSIC [10, 13, 14] to generate simulation pseudodata. The setup of the fluid dynamiccalculations roughly reflect conditions in Au+Au collisions at the Relativistic Hadron Collider (RHIC) and Pb+Pbcollisions at the Large Hadron Collider (LHC), as described below. For direct comparisons of MUSIC simulationsto data we refer the reader to [13, 14]. We will introduce a map from the ( T fo - η/s )-plane of ”true” values set influid dynamic simulations to the corresponding ( T fo - η/s )-plane extracted from viscous blast wave fits of fluid dynamicspectra and elliptic flow. We are providing a parameterization of this map and its inverse. The inverse map canbe used to improve blast wave fits by removing the bias introduced by the blast wave approximations. We will alsoestimate the uncertainties from this procedure.The paper is organized as follows. In section II, we review the viscous Retiere-Lisa blast wave and discuss theapproximations made in blast waves. In section III, we describe the setup of the MUSIC hydrodynamic calculationsand the pseudodata that is fitted. In section IV, we provide the results of the viscous blast wave fits. In section V,we discuss the relation between fluid dynamic parameters and blast wave fit parameters and utilize it to improve aprevious extraction of hadronic shear viscosity. II. FLUID DYNAMICS AND VISCOUS BLASTWAVE
In this section we briefly review some basic concepts shared by both fluid dynamics freeze-out and blast waves. Wewill then discuss the particular ansatz for the viscous blast wave in [1], based on the work by Retiere and Lisa [4]. Fora system close enough to local kinetic equilibrium one can assign a local temperature field T ( r ) and a flow field u µ ( r )to describe the temperature and collective motion as a function of position 4-vector r µ . The particle distribution inthe local rest frame of a fluid cell can then be written as f ( r, p ) = f ( r, p ) + δf ( r, p ) (1)where f is the equilibrium Bose or Fermi-distribution as a function of particle momentum p µ = ( E, p ) with the localtemperature T , f ( r, p ) = 1 e E/T ∓ , (2)and δf is a small correction that accounts for the out-of-equilibrium behavior. We neglect chemical potentials in thedistribution but their inclusion is straight forward [1]. The general form of the correction term is [1, 15, 16] δf ( r, p ) = ηs Γ(6)Γ(4 + λ ) (cid:18) ET (cid:19) λ − p µ p ν T σ µν f ( r, p ) . (3)Here the shear stress tensor π µν has been expressed by its Navier-Stokes approximation, π µν = 2 ησ µν , where σ µν isthe traceless gradient tensor, defined as σ µν = 12 ( ∇ µ u ν + ∇ ν u µ ) −
13 ∆ µν ∇ λ u λ . (4)We have used the notation ∇ µ = ∆ µν ∂ ν , with ∆ µν = g µν − u µ u ν , for the derivative perpendicular to the flow fieldvector u µ . In the following we will use the standard choice λ = 2 for the residual momentum dependence of thecorrection as in [1].When solving the viscous fluid dynamics equations of motion, numerical stability requires second order gradientterms to be included, leading to equations of motion for the shear stress tensor π µν and bulk stress Π [17–21]. Atfreeze-out, it is then convenient to compute δf directly from the shear stress tensor π µν . On the other hand, forthe blast wave it is more practical to utilize the Navier-Stokes approximation π µν = 2 ησ µν and to compute viscouscorrections using Eq. (3). In that case δf is calculated simply from the flow field, which can be independentlyconstrained by fits to flow data and the specific shear viscosity η/s of nuclear matter. The differences between thetwo approaches of calculating δf at freeze-out, π µν vs Navier-Stokes, are parametrically small in situations of smallgradients towards the end of the time evolution. However, they could still be noticeable at freeze-out in realisticsystems and are part of the uncertainties to be accounted for.In both blast wave and fluid dynamic freeze-out, the invariant particle momentum spectrum emitted from a hyper-surface Σ in Minkowski space is given by the Cooper-Frye formula [22] dNdyd P T = g (cid:90) p · d Σ(2 π ) f ( r, u · p ) (5)where g is the degeneracy factor for a given particle and d Σ µ is the forward normal vector on the freeze-outhypersurface. The momentum vector in the laboratory frame is written as usual, p µ = ( M T cosh y, P T cos ψ,P T sin ψ, M T sinh y ), in terms of the transverse momentum P T , the longitudinal momentum rapidity y and theazimuthal angle ψ in the transverse plane. M T = P T + M defines the transverse mass M T for a hadron of mass M .The final particle spectrum at freeze-out is usually calculated on a hypersurface at constant temperature T = T fo . Incontrast, in fluid dynamics this isothermal hypersurface, as well as the flow field u µ on it can be computed in thesimulation itself. For the blast wave we have to choose ans¨atze for both.Following Ref. [4] we assume that freeze-out from an isothermal hypersurface at temperature T can be approximatedby freeze-out from a hypersurface at constant proper longitudinal time τ . We enforce boost invariance, which is agood approximation for nuclear collisions around midrapidity at top RHIC and LHC energies and is also often foundin fluid dynamic calculations. To keep the blast wave simple we have to restrict ourselves to describing smooth fluiddynamics which corresponds to the averaging over many events. We can then assume that the hypersurface in the x − y -plane is approximately an ellipse with semi-axes R x and R y in x - and y -directions, respectively. We define thecoordinate axes such that the impact parameter b of the collision is measured along the x -axis. In the following weuse the reduced radius ρ = (cid:113) x /R x + y /R y together with the azimuthal angle θ , with tan θ = ( R x y ) / ( R y x ), andspace time rapidity η s = 1 / t + z ) / ( t − z )] to carry out the integral over the hypersurface. Restricting ourselvesto hadrons measured at midrapidity y = 0 and changing to convenient coordinates we obtain the final expression forthe particles from the blastwave: dNdyd P T = gτ R x R y M T (cid:90) dρ (cid:90) π dθ (cid:90) ∞−∞ dη s ρ cosh η s (2 π ) f ( ρ, θ, η s ; u · p ) (cid:20) ηs T p µ p ν σ µν (cid:21) . (6)Next we have to make an ansatz for the collective flow field. The general parameterization is u µ = (cosh η s cosh η T , sinh η T cos φ u , sinh η T sin φ u , sinh η s cosh η T ) (7)where η T is the transverse rapidity in the x − y -plane, and φ u is the azimuthal angle of the flow vector in the transverseplane. Boost invariance fixes the longitudinal flow rapidity to be equal to the space time rapidity η s . We follow Retiereand Lisa and choose to model the transverse flow velocity v T = tanh η T as [4] v T = ρ n ( α + α cos(2 φ u )) (8)which encodes a Hubble-like velocity ordering with an additional shape parameter n . α is the average velocity on theboundary ρ = 1, and α parameterizes the elliptic deformation of the flow field built up from the initial elliptic spatialdeformation of the system. Flow vectors tend to be tilted towards the smaller axis of the ellipse. In the RL approachthey are chosen to be perpendicular to the elliptic surface at ρ = 1, i.e. tan φ u = R x /R y tan φ , where φ = arctan y/x is the azimuthal angle of the position r µ .With a parameterization of the flow field at hand, the next step is the calculation of the gradient tensor σ µν . Thishas been carried out in [1] and we refer the reader to the details in that reference. We want to point out that temporalderivatives are calculated using ideal fluid dynamic equations of motion rather than the free-streaming approximation[7, 8]. This introduces the nuclear matter equation of state, specifically the speed of sound squared c s into ourcalculation of δf .The blast wave ansatz has several parameters which allow us to adjust the flow field and the hypersurface, as wellas the temperature, specific shear viscosity, and speed of sound squared at freeze-out. The full set of parameters is˜ P = ( τ fo , R x , R y , T fo , n, α , α , η/s, c s ). We will drop c s from this list and rather use guidance on the hadronic matterequation of state from existing literature [27]. Note that the viscous blast wave depends on the parameters τ fo , R x and R y separately, and not just on the total volume τ R x R y and the ellipticity R y /R x [1], as is the case for the idealblast wave.Despite the large number of parameters it is clear that the blast wave has introduced significant simplificationscompared to a realistic freeze-out calculated in fluid dynamics. They come mostly from the simplified shape of thehypersurface (constant proper time) and the spatial structure of the flow field. Two more major approximations aremade for sake of simplicity. First, resonance decays are usually neglected in blast wave calculations, and only hadronsstable under strong decays are taken into account (see [28] for recent work on implementing decays). Secondly,correction terms to the particle distribution f due to bulk stress have been ignored. They could in principle beadded and we plan to do so in the future. We summarize these five major approximations compared to fluid dynamicfreeze-out in the following list: • Navier-Stokes approximation used for δf . • Certain aspects of the shape of the hypersurface are fixed. • General shape of the flow field is fixed. • Lack of resonance production and decay. • Missing bulk stress effects on particle distributions.Since we have eliminated event-by-event fluctuations from the comparison, the effects of event-by-event fluid dy-namic simulations compared to smooth fluid dynamics need to be considered separately. They are not included in thestudy below. The same is true for deviations of state-of-the-art 3+1D fluid dynamics from the boost-invariant 2+1Dfluid dynamics used here. The effects of fluctuations and breaking of boost invariance have already been studiedwithin fluid dynamics [23, 24] and can be added to the considerations in this work. b (fm) (Au+Au) 5 5 6 6.5 7 8 9 10.5 10.5 b (fm) (Pb+Pb) 5.3 6.3 6.9 7.4 8.5 9.6 11.1 11.1 T (true)fo (MeV) 105 110 115 120 125 130 135 140 1454 π ( η/s ) (true) T (true)fo , ( η/s ) (true) chosen for MUSIC simulations of Au+Au and Pb+Pb collisions, together withthe impact parameter b . T (true)fo (MeV) P T -range spectra (GeV/c) P T -range v (GeV/c)pion kaon proton pion kaon proton105 0.34-1.95 0.34-2.23 0.76-2.52 0.53-3.0 0.34-3.0 0.34-3.0110 0.34-2.37 0.34-2.68 0.34-3.0 0.34-3.0 0.34-3.0 0.34-3.0115 0.34-1.95 0.34-2.23 0.34-2.52 0.34-3.0 0.34-3.0 0.34-3.0120 0.40-1.95 0.40-2.09 0.34-2.37 0.34-2.84 0.34-3.0 0.34-3.0125 0.40-1.95 0.40-2.09 0.34-2.37 0.34-2.68 0.34-2.84 0.34-3.0130 0.46-1.95 0.40-2.09 0.29-2.23 0.34-2.68 0.34-2.68 0.34-2.84135 0.34-1.82 0.29-1.95 0.24-2.09 0.34-2.52 0.34-2.68 0.34-2.68140 0.24-1.57 0.20-1.69 0.20-1.82 0.24-2.23 0.53-2.23 0.20-2.37145 0.24-1.57 0.20-1.69 0.20-1.82 0.24-2.23 0.53-2.23 0.20-2.37TABLE II. Preferred fit ranges for the pseudodata from Set I for Au+Au collisions. Similar fit ranges are used for Pb+Pbpseudodata. III. GENERATION OF MUSIC PSEUDODATA
We use the viscous hydrodynamics code MUSIC to simulate averaged nuclear collisions at RHIC and LHC energiesat various impact parameters. MUSIC is a relativistic second-order viscous hydrodynamics code for heavy ion collisions[13, 14, 23]. We choose boost-invariant (2+1)D mode consistent with the boost-invariant blast wave set-up. We use thebuilt-in optical Glauber model to generate initial conditions with the appropriate nucleon-nucleon cross section and anoverall normalization roughly consistent with pertinent multiplicity data from RHIC and LHC. We use the equationof state (EOS) s95p-v1.2 in MUSIC. We set a constant shear viscosity η/s and the default MUSIC bulk viscosity.We freeze out at pre-determined temperatures T fo and compute the final hadron spectra and elliptic flow, includingresonance decays and including viscous corrections to freeze-out. The detailed MUSIC settings are documented inAppendix A. With later applications in mind, we focus on Au+Au collisions at the top RHIC energy of √ s NN = 200GeV, as well as Pb+Pb collisions at the LHC energy of √ s NN = 5 .
02 GeV.As discussed above, our main goal here are the relations of the temperature T (extr)fo and specific shear viscosity( η/s ) (extr) extracted from a blast wave fit of the pseudodata to the true values T (true)fo and ( η/s ) (true) used in thegeneration of these pseudodata sets. Generally, we expect η/s to decrease with temperature in the hadronic phase,and T fo to increase with impact parameter b [1]. To focus on the physical region we investigate a band in the T (extr)fo -( η/s ) (extr) -plane which covers the values extracted from RHIC and LHC data in Ref. [1]. For the case of simulationsat RHIC energy we choose 27 points, in the T (true)fo -( η/s ) (true) -plane, organized in 3 sets of nine points each. ForLHC we choose a similar set of points. For each point we run MUSIC and perform a blast wave fit to the resultinghadron spectra and elliptic flow. The value of the impact parameter b for the fluid dynamic simulation for each T (true)fo -( η/s ) (true) -point is chosen using guidance from the fits to data carried out in [1]. The nine points (eight pointsfor Pb+Pb) of the first set (Set I) used in the T (true)fo - ( η/s ) (true) -plane are shown in Tab. I together with the impactparameter for each simulation.We briefly discuss the details of this process for Set I for Au+Au collisions as an example. For a given parametersvalues MUSIC runs and computes the transverse momentum spectra and elliptic flow v at rapidity y = 0 for pions,kaons and protons. The P T -range of the pseudodata generated is from 0 to 3 GeV/ c . The fit ranges used for blastwave fits of the pseudodata in Set I are shown in Tab. II. The rationale behind the choice of fit ranges was discussedin detail in Ref. [1]. To summarize briefly, very low momenta tend to be not directly described by blast wave fits asthey are dominated by resonance decays, while higher momenta receive viscous corrections larger than what can bereliably described by the Navier-Stokes approximation. Overall, the fit ranges for fitting the simulation pseudodataare very similar to the ranges used for good quality fits of experimental data in [1].The fluid dynamic simulations do not provide useful uncertainty estimates on the pseudodata. For this study wechoose uncertainties in line with uncertainties in the pertinent available experimental data. We assign 5% uncertaintyand 2% uncertainty to pseudodata spectra and v , respectively. We add a pedestal of 0.002 to the uncertainty of v for realistic error bars at smaller P T where v is very small. To quantify the uncertainty from this choice we vary theuncertainties, similar to the procedure in [1].For fits we use the statistical analysis package from the Models and Data Analysis Initiative (MADAI) project[25, 26]. The MADAI package includes a Gaussian process emulator and a Markov Chain Monte Carlo for Bayesiananalyses. One example of such a Bayesian analysis can be found in Fig. 1. We use N = 500 training points for eachGaussian emulator. We check the Gaussian emulators and find their errors to be acceptably low. IV. BLAST WAVE FITS
We fit the following blast wave parameters in the Bayesian analysis: P = ( T fo , R y /R x , α , α , η/s ). Chemicalpotentials µ are set to zero in MUSIC. To simplify the analysis we fix the radial shape parameters n and the freeze-out time τ by choosing values close to those extracted from RHIC and LHC data in Ref. [1], see Tab. III and Tab.IV. As in [1] we also use geometric arguments to determine R x and R y independently from the fitted ratio. A l pha0 R y O v e r R x A l pha2 E t a O v e r S i Temp
Alpha0
RyOverRx
Alpha2 L i k e li hood L i k e li hood EtaOverS4pi
FIG. 1. Likelihood analysis for the MUSIC run for Au+Au with T (true)fo = 130 MeV, ( η/s ) (true) = 2 . / (4 π ). The axesshow (from left to right, and top to bottom): T fo (MeV), α , R y /R x , α , η/s .The diagonal plots show the posterior likelihooddistributions. The off-diagonal plots show correlations between parameters. As an example, let us take a brief look at the case of parameters T (true)fo = 130 MeV, ( η/s ) (true) = 2 . / (4 π )in Au+Au Set I. Fig. 1 shows the result from the statistical analysis. The likelihood plots exhibit well definedpeaks. The extracted values for the full set of parameters in this case are T fo = 117 . α =0.753 c , R y /R x =1.10, α =0.048 c , η/s =2.6/4 π , with τ = 8 . c and n =0.88 fixed. We proceed analogously for the other T (true)fo , ( η/s ) (true) points of Au+Au Set I in Tab. I, and for the Pb+Pb points. The results are summarized in Tab. III and Tab. IV. Hydro Au+Au Blast Wave T (true)fo π ( η/s ) (true) T fo (MeV) α /c R y /R x α /c πη/s τ (fm/ c ) n
105 6.03 111.2 0.824 0.99 0.021 5.83 12.2 0.86110 5.28 114.0 0.822 1.01 0.021 5.43 11.4 0.87115 4.52 112.7 0.833 1.04 0.025 4.67 10.6 0.81120 3.77 113.9 0.820 1.06 0.028 3.75 9.8 0.84125 3.02 117.7 0.786 1.08 0.037 3.01 9.1 0.88130 2.51 117.2 0.753 1.10 0.048 2.59 8.4 0.88135 2.01 120.0 0.715 1.15 0.059 2.07 7.8 0.92140 1.51 123.0 0.654 1.27 0.069 1.55 7.2 0.96145 1.01 126.3 0.604 1.35 0.080 1.23 6.8 1.00TABLE III. Extracted parameter values P for points in Set I for Au+Au collisions, together with values set for τ and n .Hydro Pb+Pb Blast Wave T (true)fo π ( η/s ) (true) T fo (MeV) α /c R y /R x α /c πη/s τ (fm/ c ) n
110 5.28 111.4 0.822 0.99 0.020 5.74 13.2 0.84115 4.52 117.5 0.827 1.00 0.023 4.73 12.6 0.87120 3.77 120.0 0.818 1.03 0.026 3.98 12 0.85125 3.02 123.3 0.822 1.07 0.032 3.34 11.6 0.88130 2.51 123.7 0.787 1.09 0.043 2.62 10.8 0.90135 2.01 125.6 0.750 1.13 0.054 2.01 10.0 0.94140 1.51 127.2 0.689 1.19 0.063 1.48 9.2 0.98145 1.01 130.3 0.642 1.24 0.075 1.18 8.6 1.00TABLE IV. The same as Tab. III for points in Set I for Pb+Pb collision.
Interestingly, the scaling between true and extracted events does not seem to depend very much on the initial energyof the fireball. For example, for (
T, η/s ) (true) = (125 MeV, 3.0/4 π ), we extract ( T, η/s ) (extr) = (118 MeV, 3.0/4 π ) forAu+Au collisions and ( T, η/s ) = (122 MeV, 3.3/4 π ) for Pb+Pb collisions. Spectra and v calculated using the blastwave with the preferred parameter values from Tab. III are compared to the MUSIC Au+Au pseudodata in Figs. 2and 3, respectively. Comparisons with the Pb+Pb pseudodata lead to similar conclusions but are not shown here.The agreement with the pseudodata is generally good.From here on we focus on the extracted freeze-out temperatures and specific shear viscosities in Tabs. III and IV,i.e. T (extr)fo and ( η/s ) (extr) . We interpret them as the image of the original values T (true)fo and ( η/s ) (true) under a map M . Before we discuss this map in detail in the next section let us discuss uncertainties for the extracted values. At theend of Sec. II we had summarized important approximations made by the blast wave approach. These approximationslead to differences between the true and extracted values whose quantification is the goal of this and the followingsection. The extracted values themselves already come with a set of uncertainties. We have quantified three of them:(i) Uncertainties due to using a fixed value of the shape parameter n in our analysis. (ii) Uncertainties due to theerrors assigned to MUSIC pseudodata; while our choice of errors for the pseudodata is motivated by experiment, thesize of the error is expected to influence fit results. Finally, (iii) uncertainties of the Bayesian analysis itself, expressedin the posterior distributions. The latter ones can be extracted directly from the analysis and can be seen in Fig.1. We have studied the uncertainties from sources (i) and (ii) by systematically varying n and the error given topseudodata. We combine all three sources into combined uncertainties in temperature and specific shear viscosity foreach point. An example for this analysis is given in appendix B.The final extracted freeze-out temperature and specific viscosities for the Au+Au points from Tab. I are shownin Fig. 4 (left panel) with their combined uncertainties. We find that the extracted specific shear viscosities aremostly consistent with true values within uncertainties. However the extracted temperatures underestimate the truetemperatures significantly at higher values of T fo . The extracted T fo is about 15 MeV smaller for hydro events at hightemperature (or peripheral collisions). For events at low freeze-out temperature (or central collisions), the extracted T fo are within 5 MeV of the true temperature. The deviations in freeze-out temperature for peripheral collisions aresignificant with respect to the estimated uncertainties also shown in Fig. 4. A very similar picture emerges for theresults for Pb+Pb points also shown in Fig. 4 (right panel). - - - d N / d P T d y ( c / G e V ) P + P (× ) K ± (× ) π ± MUSIC T = - - - P + P (× ) K ± (× ) π ± MUSIC T = - - - P + P (× ) K ± (× ) π ± MUSIC T = - - - d N / d P T d y ( c / G e V ) P + P (× ) K ± (× ) π ± MUSIC T = - - - P + P (× ) K ± (× ) π ± MUSIC T = - - - P + P (× ) K ± (× ) π ± MUSIC T = - - - - P T ( GeV / c ) d N / d P T d y ( c / G e V ) P + P (× ) K ± (× ) π ± MUSIC T = - - - - P T ( GeV / c ) P + P (× ) K ± π ± (× ) MUSIC T = - - - - P T ( GeV / c ) P + P (× ) K ± π ± (× ) MUSIC T = FIG. 2. Transverse momentum spectra (symbols) for protons, pions and kaons in Au+Au collisions calculated in MUSICtogether with blast wave calculations (solid lines) using the extracted parameter values in Tab. III. The nine plots correspondto the nine points of freeze-out temperature and specific shear viscosity in Set I.
V. QUANTIFYING THE ACCURACY OF BLAST WAVE FITS
We are now ready to quantify the deviations between true and fitted values. We define a map M as (cid:18) T (extr) ( η/s ) (extr) (cid:19) = M (cid:18) T (true) ( η/s ) (true) (cid:19) . (9)We first focus on Au+Au collisions. We find that a linear map M in this 2-dimensional parameter space has sufficientaccuracy given the uncertainties established for the fitted values. For a given set of points we apply the method ofnormal equations which gives the least square approximate solution to M . For Au+Au Set I we find, see AppendixC, M = (cid:20) . . . . (cid:21) . ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● V P + PK ± π ± MUSIC T = ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● P + PK ± π ± MUSIC T = ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● P + PK ± π ± MUSIC T = ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● V P + P K ± π ± MUSIC T = ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● P + P K ± π ± MUSIC T = ●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ● P + P K ± π ± MUSIC T = ●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ● ● P T ( Gev / c ) V P + PK ± π ± ↖ MUSIC T = ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● P T ( Gev / c ) P + P K ± π ± MUSIC T = ●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● P T ( Gev / c ) P + P K ± π ± MUSICT = FIG. 3. Same as Fig. 2 for the elliptic flow v . This M is only an approximate solution to Eq. (9). We define the predictions from the map to be (cid:18) T ( M ) ( η/s ) ( M ) (cid:19) = M (cid:18) T (true) ( η/s ) (true) (cid:19) . (10)and show them together with the true and extracted values in Fig. 4 (left panel). Overall, the map M provides agood approximation of extracted blast wave values. For Pb+Pb collisions, we obtain a similar map through the sametechnique, M Pb+Pb = (cid:20) . . − . . (cid:21) , which is shown in Fig. 4 (right panel) together with the true and extracted values in that case.So far we have only used Set I of points. Appendix C contains some details on the other sets, as well as variationsof the map M which are needed to capture the uncertainties of the extracted values. The results for Sets II (Tab. IX)and III (Tab. X) in the Au+Au case are summarized in Fig. 5.We are now in a position to remove biases that come with the extraction of properties of nuclear collisions usingblast wave fits. As an example, we will attempt to correct the bias introduced by the blast wave fit from the analysis0 ● ●●● ●● ● ● ● ▲▲▲▲▲ ▲ ▲ ▲ ▲◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆
100 110 120 130 1402468 T ( MeV ) π η / s ◆◆◆ MUSIC Au + Au Blastwave ▲▲▲
Map using M ● ● ● ●● ● ● ● ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆
110 120 130 1402468 T ( MeV ) π η / s ◆◆◆ MUSIC Pb + Pb Blastwave ▲▲▲
Map using M Pb + Pb FIG. 4. Left panel: Comparison of the true freeze-out temperature and specific shear viscosity used to create Au+Au SetI pseudodata (red diamonds), their counterparts extracted from blastwave fits (green circles), and the values using the linearmap M discussed in the text (gray triangles). The lines are drawn to guide the eye. Right panel: The same comparison forPb+Pb collisions. ● ●●● ●● ● ● ● ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲▲▲▲ ▲ ▲ ▲ ▲ ▲
100 110 120 130 1402468 T ( MeV ) π η / s ◆◆◆ MUSIC smaller η / s Blastwave ▲▲▲
Map using M ● ●●● ●● ● ●● ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆▲▲▲▲▲▲ ▲ ▲ ▲
100 110 120 130 1402468 T ( MeV ) π η / s ◆◆◆ MUSIC larger η / s Blastwave ▲▲▲
Map using M FIG. 5. The same comparison as shown in Fig. 4 for Set II (left panel) and Set III (right panel) of Au+Au collisions. Thedashed line indicates the location of points in Set I. of the specific shear viscosity in [1]. Although our work has focused on those two fit parameters, the analysis inthis work could be repeated for other quantities extracted from blast wave fits in a straight forward way. We definecorrected blast wave fit values as (cid:18) T (corr) ( η/s ) (corr) (cid:19) = M − (cid:18) T (extr) ( η/s ) (extr) (cid:19) . (11)Since Eq. (9) assures us that (cid:18) T (true) ( η/s ) (true) (cid:19) ≈ M − (cid:18) T (extr) ( η/s ) (extr) (cid:19) . (12)1we postulate that this correction step takes the blast wave fitted values closer to physical reality. We will revisit thisargument in the discussion in the next section.As an example, for Au+Au collisions at RHIC the inverse matrix is M − = (cid:20) . − . − . . (cid:21) . In Ref. [1] T and η/s have been extracted from data with carefully quantified Gaussian uncertainties. These un-certainties are transported by M − to Gaussian uncertainties on T (corr) , ( η/s ) (corr) , see left panel of Fig. 7 for anexample. We add the uncertainty on the map M − itself which we do by applying additional Gaussian smearing usingthe matrices M ( u ) − , M ( d ) − , M ( l ) − , M ( r ) − , and M − , M − introduced in Appendix C. The right panel of Fig. 7gives an impression of the size of the effect by plotting the variation of the center value under these inverse maps.The final results of this correction, including all uncertainties combined, is shown in Fig. 6. We notice that aftercorrection the temperature dependence is less steep, and values close to 4 πη/s = 1 are only reached around T ≈ ● ● ● ● ●◆ ◆◆ ◆ ◆
120 140 160 T ( MeV ) π η / s ALICE unfold ◆◆◆
PHENIX unfold ● ● ● ● ● ◆◆◆ ◆◆
100 120 140 1602468 T ( MeV ) π η / s ALICE fit ◆◆◆
PHENIX fit
FIG. 6. Specific shear viscosity η/s at corresponding kinetic freeze-out temperature T extracted from ALICE and PHENIXbefore removing the blast wave bias (left panel) and final values after correction with combined uncertainties (right panel). VI. SUMMARY
We have discussed the shortcomings of blast wave fits compared to fluid dynamic simulations applied to high energynuclear collisions. We have carried out a systematic study of parameters at kinetic freeze-out set in fluid dynamicsimulations which have subsequently been extracted from blast wave fits to hadron spectra and elliptic flow. We findthat blast wave fits correctly reproduce trends of the true values set in the simulations, and they qualitatively agreewith the numerical values. However, some deviations between extracted and true values exist are significant. Thequality of blast wave fits can be improved by understanding and quantifying these deviations. We have done so forboth the freeze-out temperature and specific shear viscosity at freeze-out and have presented the results in this paper.We have established a map from true parameter values to blast wave fitted parameters. This map can be invertedand can be applied to correct the biases in blast wave fits. This has been demonstrated for the case of the extraction2of the specific shear viscosity as a function of temperature from data. We have also estimated uncertainties for thecorrection which have been added to the uncertainties of the final corrected results. There are good reasons to expectthat the corrected results are closer to the true physical values than the raw blast wave fit values. A rigorous argumentcan only be made for the correction to improve fits to hydrodynamic pseudodata. This leaves the questions how closehydrodynamic modeling is to the true physical situation in high energy nuclear collisions. We can only approach ananswer to that by continued studies of nuclear collisions with fluid dynamic models, which are outside of the scope ofthis work.One interesting question regards the universality and systematic behavior of the corrections discussed here. Wefound a large dependence on the temperature, which might be a reflection of changing system size. On the other handwe found very little variation going from RHIC to LHC energies. The systematic study of the deviations as a functionof various changes that can be made to the fluid dynamic calculations would be a worthwhile goal for a follow upstudy. Future improvements to the blast wave itself can be made by adding bulk stress corrections, and to includesome resonances and their decays. However, simplicity is key to blast wave fits and perhaps a simpler model withknown corrections and a study of uncertainties is more useful. Blast waves with corrections applied are still orders ofmagnitude cheaper in labor and computational cost compared to fluid dynamics. They remain a quick tool for dataanalysis. The mechanism to reduce bias laid out here should inspire confidence in their ability, as long as realisticuncertainty estimates are added to the analysis.
ACKNOWLEDGMENTS
Z. Y. would like to thank Michael Kordell, Mayank Singh and Jean-Francois Paquet for help and support regardingMUSIC. This work was supported by the US National Science Foundation under award nos. 1550221 and 1812431.
Appendix A: MUSIC Settings
MUSIC settings used are documented in Tab. V for RHIC energies and for LHC (in parentheses where different).
Parameter SetTarget+Projectile Au+Au (Pb+Pb)Maximum energy density 54.0 (96.0) a SigmaNN 42.1 (70.0)Initial profile Optical Glauber modelEOS to use lattice EOS s95p-v1.2reconst type solve flow velocity for hydro eqnsboost invariant 1Viscosity Flag 1Include Shear Visc 1T dependent Shear to S ratio 0Include Bulk Visc 1Include second order terms 0Do FreezeOut 1freeze out method complex methoduse eps for freeze out use temperaturept steps 36min pt 0.01max pt 3.0phi steps 40Include deltaf in Cooper-Frye formula 1Include deltaf bulk 1 a e max = 75.0 was used for T fo = 110 MeV TABLE V. Parameter set for MUSIC runs generating the pseudodata. 1 and 0 are flags corresponding to YES and NO. Appendix B: Blast wave extraction error estimates
If ones assume the blast wave model as given, there are uncertainties to consider when it is deployed to describeMUSIC pseudodata. As described in the main text they mainly come from three sources: (i) Uncertainties due tousing a fixed value of the radial shape parameter n in our analysis. (ii) Uncertainties due to the errors assigned toMUSIC pseudodata. (iii) Uncertainties of the fit itself as encoded in the posterior distributions. We have studied theuncertainties from source (i) by systematically varying n . An example for one point in Au+Au Set I is given in Tab.VI. Similarly, we vary the error given to the pseudodata before the fit. Tab. VII gives an example for the same pointin Set I when the error on the hadron spectra is varied. n T fo (MeV) 4 πη/s small 0.83 116.4 3.26regular 0.88 117.2 2.59large 0.93 117.8 1.81TABLE VI. The extracted values of T fo and η/s for different values set for n . This example uses MUSIC Au+Au pseudodatafor T fo =130 MeV. spectra uncertainty T fo (MeV) 4 πη/s small 4% 116.2 2.43regular 5% 117.2 2.59large 6% 117.8 2.62TABLE VII. The extracted values of T fo and η/s for different uncertainties assigned to MUSIC pseudodata hadron spectra forAu+Au with T fo =130 MeV. The uncertainty for v is fixed at 2% with a pedestal 0.002. The uncertainty for spectra is variedas shown in the table. In general, we find that T fo varies only within a few MeV while η/s shows larger changes when varying n . Variationsin the assigned pseudodata error lead to small changes of both the extracted temperature and η/s . Uncertaintiesfrom the fit can be taken directly from the MADAI code. We treat the sources of uncertainties as independent andadd them quadratically to arrive at estimates of the total uncertainty assigned to the extracted values of T fo and η/s .For the point selected for this example these errors are summarized in Tab. VIII. Origin of uncertainty n error assigned stat. analysis total σσ T (MeV) 0.57 0.65 1.47 1.70 σ η (4 π ) 0.59 0.08 0.22 0.64TABLE VIII. A summary of uncertainties for temperature and specific shear viscosity extracted from MUSIC pseudodata forAu+Au with T fo =130 MeV. Appendix C: Establishing The Map M
To find an approximate solution of Eq. (9) we write the true and extracted values of a given set (we use Au+AuSet I as an example here, see Tab. I and Tab. III) in numerical matrix forms A T = (cid:20) . . . . . . . . . .
03 5 .
28 4 .
52 3 .
77 3 .
02 2 .
51 2 .
01 1 .
51 1 . (cid:21) (C1) B T = (cid:20) . . . . . . . . . .
83 5 .
43 4 .
67 3 .
75 3 .
01 2 .
59 2 .
07 1 .
55 1 . (cid:21) (C2)respectively. We are looking for an approximate solution of the overdetermined equation B T ≈ M A T . (C3)Here T means transposition and the first columns have units of MeV and the second columns units of 1 / π .We proceed by solving the associated normal equations B T A = M A T A . (C4)4yielding the solution M = (cid:20) . . . . (cid:21) . Using M , we can calculate the values B ( M ) = ( M A T ) T which are now predictions for the blast wave fit results usingthe original values,To capture the uncertainties in the extracted values we parameterize auxiliary maps which approximate the verticesand co-vertices of the ellipses for each point (called u , d , l , r for up, down, left, right (co)-vertices). They define anerror corridor around the map M . These auxiliary maps for Au+Au Set I are M ( u ) = (cid:20) . . . . (cid:21) M ( d ) = (cid:20) . . . . (cid:21) M ( l ) = (cid:20) . . . . (cid:21) M ( r ) = (cid:20) . . . . (cid:21) . We proceed to process Sets II and III of parameter points in the same way as Set I. The true parameter valuesand the MUSIC pseudodata and extracted blast wave fit parameters are listed in Tabs. IX and X. Sets II and III aregenerated by varying η/s compared to the values of Set I. The parameterizations of the maps from true to extracted
Hydro (true) T fo (MeV) 105 110 115 120 125 130 135 140 1454 πη/s T fo (MeV) 110.7 113.8 111.9 113.4 118.5 117.5 121.4 124.6 128.44 πη/s T fo and η/s values extracted from the blastwave are also shown.Hydro (true) T fo (MeV) 105 110 115 120 125 130 135 140 1454 πη/s T fo (MeV) 111.9 114.4 113.0 113.8 118.2 117.6 120.7 123.7 124.44 πη/s T fo and η/s values extracted from the blastwave are also shown. values for these two sets are M = (cid:20) . . . . (cid:21) M = (cid:20) . . − . . (cid:21) . Using M and M , we can calculate the values predicted for the blastwave fits and compare to the actual numbers,shown in Fig. 5. The maps M and M are close to M but differences are present and large enough to warrant theintroduction of an additional uncertainty for the map M . Fig. 7 shows an example how uncertainties in raw blastwave values are transported by M − , and how the center value varies due to different maps M ( u ) − , M ( d ) − , M ( l ) − , M ( r ) − , M − and M − . [1] Z. Yang and R. J. Fries, “Extraction of the Specific Shear Viscosity of Hot Hadron Gas,” arXiv:1807.03410 [nucl-th].[2] P. J. Siemens and J. O. Rasmussen, “Evidence for a blast wave from compress nuclear matter,” Phys. Rev. Lett. , 880(1997).[3] E. Schnedermann, J. Sollfrank and U. W. Heinz, “Thermal phenomenology of hadrons from 200-A/GeV S+S collisions”,Phys. Rev. C , 2462 (1993).[4] F. Retiere and M. A. Lisa, “Observable implications of geometrical and dynamical aspects of freeze out in heavy ioncollisions”, Phys. Rev. C , 044907 (2004).[5] W. Florkowski and W. Broniowski, “Hydro-inspired parameterizations of freeze-out in relativistic heavy-ion collisions,”Acta Phys. Polon. B , 2895 (2004).[6] A. Mazeliauskas and V. Vislavicius, Phys. Rev. C , 014910 (2020). ●● ◆◆
120 130 1401234 T ( MeV ) π η / s Phenix fit ◆ Unfold using M - ●● ▲▲ ◆◆▲▲
120 130 140 150 T ( MeV ) Unfold using M ( u )- etc. ▲▲ Unfold using M - ( M - ) FIG. 7. Example for the bias correction for one blast wave fit result from Ref. [1] (PHENIX 20-30% fit result, T fo = 123.5 ± πη/s = 1.98 ± T (extr) , ( η/s ) (extr) extracted from the blast wave fit of PHENIX data and itsuncertainty (green circle) and the corrected value with its uncertainty using M − (purple diamond and ellipse). Right panel:Variance of the corrected center value due to the six maps M ( u ) − , etc., described in the text.[7] D. Teaney, “The Effects of viscosity on spectra, elliptic flow, and HBT radii”, Phys. Rev. C , 034913 (2003).[8] A. Jaiswal and V. Koch, “A viscous blast-wave model for relativistic heavy-ion collisions”, arXiv:1508.05878 [nucl-th].[9] Z. Yang and R. J. Fries, “A Blast Wave Model With Viscous Corrections,” J. Phys. Conf. Ser. , 012056 (2017).[10] C. Gale, S. Jeon and B. Schenke, “Hydrodynamic Modeling of Heavy-Ion Collisions,” Int. J. Mod. Phys. A , 1340011(2013).[11] J. E. Bernhard, J. S. Moreland and S. A. Bass, “Bayesian estimation of the specific shear and bulk viscosity of quark gluonplasma,” Nature Phys. , 1113 (2019).[12] J. F. Paquet et al. [JETSCAPE], “Revisiting Bayesian constraints on the transport coefficients of QCD,” arXiv:2002.05337[nucl-th].[13] B. Schenke, S. Jeon and C. Gale, “(3+1)D hydrodynamic simulation of relativistic heavy-ion collisions”, Phys. Rev. C ,014903 (2010).[14] S. Ryu, J.-F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon and C. Gale, “Importance of the Bulk Viscosity of QCDin Ultra-relativistic Heavy-Ion Collisions”, Phys. Rev. Lett. , 132301 (2015).[15] M. Damodaran, D. Molnar, G. G. Barnafldi, D. Bernyi and M. Ferenc Nagy-Egri, “Testing and improving shear viscousphase space correction models”, arXiv:1707.00793 [nucl-th].[16] M. Damodaran, D. Molnar, G. G. Barnafldi, D. Bernyi and M. F. Nagy-Egri, “Improved single-particle phase-spacedistributions for viscous fluid dynamic models of relativistic heavy ion collisions,” Phys. Rev. C , 014907 (2020).[17] W. Israel and J. M. Stewart, “Transient relativistic thermodynamics and kinetic theory,” Annals Phys. , 341 (1979).[18] R. Baier, P. Romatschke and U. A. Wiedemann, “Dissipative hydrodynamics and heavy ion collisions,” Phys. Rev. C ,064903 (2006).[19] H. Song and U. W. Heinz, “Causal viscous hydrodynamics in 2+1 dimensions for relativistic heavy-ion collisions,” Phys.Rev. C , 064901 (2008).[20] I. Karpenko, P. Huovinen and M. Bleicher, “A 3+1 dimensional viscous hydrodynamic code for relativistic heavy ioncollisions,” Comput. Phys. Commun. , 3016 (2014).[21] P. Romatschke and U. Romatschke, “Relativistic Fluid Dynamics In and Out of Equilibrium – Ten Years of Progress inTheory and Numerical Simulations of Nuclear Collisions”, arXiv:1712.05815 [nucl-th].[22] F. Cooper and G. Frye, “Comment on the Single Particle Distribution in the Hydrodynamic and Statistical ThermodynamicModels of Multiparticle Production”, Phys. Rev. D , 186 (1974).[23] B. Schenke, S. Jeon and C. Gale, “Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics,” Phys.Rev. Lett. , 042301 (2011).[24] Z. Qiu and U. W. Heinz, “Event-by-event shape and flow fluctuations of relativistic heavy-ion collision fireballs,” Phys.Rev. C , 024911 (2011). [25] MADAI collaboration, https://madai-public.cs.unc.edu/ , accessed June 20, 2016.[26] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu and U. Heinz, “Applying Bayesian parameter estimation to relativisticheavy-ion collisions”, Phys. Rev. C , 024907 (2016).[27] D. Teaney, “Chemical freezeout in heavy ion collisions”, arXiv:0204023[nucl-th].[28] A. Mazeliauskas, S. Floerchinger, E. Grossi and D. Teaney, Eur. Phys. J. C79