Parameterizing Smooth Viscous Fluid Dynamics With a Viscous Blast Wave

Blast waves are simple and effective tools that provide snapshots of a system that is close enough to local kinetic equilibrium so that macroscopic concepts like temperature, local collective velocity, and shear and bulk stress can be used. 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 system are known, fluid dynamics can evolve the system starting from a given initial state. In contrast, blast waves usually provide a static picture, typically the temperature field $T(x)$ and collective flow field $u^{\mu}(x)$ on a hypersurface $\Sigma$ consisting of points $r^{\mu}$. E.g., if the hypersurface of kinetic freeze-out in an expanding system is chosen, these fields can often be directly related to observable particles of the system and can be obtained by fitting to data.

In high energy nuclear collisions, blast waves have routinely been used to analyze the properties of the fireball of hadrons at the time of kinetic freeze-out [2, 3, 4, 5, 6, 7, 8, 9, 10]. After kinetic freeze-out hadrons free stream to the detectors and thus directly carry information about the freeze-out hypersurface. Blast wave parameters like the freeze-out temperature and average radial flow velocity can be determined by fits to transverse momentum spectra of hadrons. Additional information, like elliptic deformations of the fireball in coordinate and momentum space at the time of freeze-out can be extracted from fits of elliptic flow coefficients $v_{2}$. More recently, viscous corrections to blast waves have been considered [11, 12, 13, 1]. They can be used to extract the specific shear viscosity $\eta/s$ for hadronic matter at the freeze-out temperature [1] and parton matter at the pseudocritical temperature [14]. Viscous blast waves can also be used to extend the range of validity of ideal blast wave fits to larger transverse hadron momenta $P_{T}$.

The question arises to what extent blast waves, which use certain simplifying approximations, can faithfully capture important properties of the full dynamical system. The global picture is one of two successive approximations

where $\succ$ means ”approximated by”. Viscous fluid dynamic simulations are widely accepted to give accurate descriptions of the low transverse momentum region of high energy nuclear collisions, and analyses of experimental data with fluid dynamic simulations continue to be an active field of study [15, 16, 17]. There are several approximations that enter when describing the evolution of the system and its freeze-out with fluid dynamics. When dealing with the late hadronic phase the main issue is the instantaneous approximation made for the freeze-out. A sudden freeze-out corresponds to the mean free path of hadrons suddenly rising to infinity, for all hadron species at once. Simulations with hadronic transport models solving the Boltzmann equation show a more realistic picture of a gradual freeze-out process. A systematic comparison of fluid dynamics and Boltzmann transport in the hadronic phase could illuminate uncertainties arising from these approximations. However, the quantification of these types of uncertainties is outside the scope of this work.

Here, we focus on the second step in the approximation chain, and quantify uncertainties that arise when fluid dynamic systems are approximated by blast waves at freeze-out. One might ask whether the deployment of blast waves is still needed, given the proliferation of viscous fluid dynamic codes, and the cheap numerical cost, at least for 2+1D codes. The motivation lies only partly in the simplicity of blast wave fits. A second, important argument is their complementarity. Fluid dynamic calculations come with their own set of uncertainties, many of which are not shared by blast waves. As an important example, fluid dynamics computes flow fields using, among other input, initial conditions and an equation of state. The final flow field at freeze-out will depend on these inputs. On the other hand, blastwaves are independent of these specific inputs and rather find the final flow field by fits to data. Of course, blast waves suffer from other, complementary, uncertainties which will be discussed here. We refer the reader to our work [1] for an example that uses the complementarity of blast waves to extract properties of hadronic matter at freeze-out.

The outline of our work is as follows. We use smooth relativistic viscous fluid dynamics to create systems close to local equilibrium as they typically occur in the late stages of high energy nuclear collisions. Specifically, we utilize the viscous fluid dynamics code MUSIC [18, 15, 19] to generate simulation pseudodata. The setups of the calculations roughly reflect conditions in Au+Au collisions at the Relativistic Hadron Collider (RHIC) and Pb+Pb collisions at the Large Hadron Collider (LHC), as described below, although a precise description of data is not the point of this work. For direct comparisons of MUSIC simulations to data we refer the reader to [18, 19]. We subsequently use the viscous blastwave introduced by us in [1] to fit the transverse momentum spectra and elliptic flow computed in MUSIC for several species of identified hadrons. This blast wave is a generalization of the ideal blast wave by Retiere and Lisa (RL) [4]. The blast wave is defined by an ansatz for the temperature, flow field and space-time structure of the freeze-out hypersurface, as discussed in detail below. The ansatz contains several parameters with physical meaning, like the size of the fireball at freeze-out, and the shape of the collective flow field. We focus here on the extracted freeze-out temperature $T_{\mathrm{fo}}$ and the specific shear viscosity $\eta/s$, i.e. the ratio of shear viscosity $\eta$ and entropy density $s$, at freeze-out. The extracted values can be directly compared to their ”true” counterparts used in MUSIC. We also compare the flow field extracted by the blast wave fit to its counterpart in MUSIC.

Differences that are seen between ”true” fluid dynamic results and extracted values can be due to simplifications made in the blast wave ansatz, or due to limiations imposed on the range and error bars on pseudodata. We will discuss both of these below. For future applications of viscous blast wave fits it is is important to understand and quantify the uncertainties and systematic biases in the fit results. As an example, we introduce a map from the ”true” values of $T_{\mathrm{fo}}$ and $\eta/s$ set in fluid dynamic simulations to the corresponding values extracted from viscous blast wave fits of spectra and elliptic flow. The inverse map can be used to improve blast wave fits by systematic unfolding. Blast waves equipped with such procedures to remove systematic bias will have significantly improved precision. In this work we take a first step in this direction. Ref. [1] in which $\eta/s$ in the hadronic phase is extracted from experimental data serves as an example for the usefulness of such procedures.

The paper is organized as follows. In section II, we review the viscous Retiere-Lisa blast wave and discuss the approximations made. In section III, we describe the setup of the MUSIC hydrodynamic calculations and 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 quantify the deviations of blast wave fits. We conclude with a discussion of our results and possible improvements.

In this section we briefly review some basic concepts shared by both fluid dynamics freeze-out and blast waves. We will then discuss the particular ansatz for the viscous blast wave in [1], based on the work by Retiere and Lisa [4]. For a system close enough to local kinetic equilibrium one can assign a local temperature field $T(r)$ and a flow field $u^{\mu}(r)$ to describe the temperature and collective motion as a function of position 4-vector $r^{\mu}$. The particle distribution in the local rest frame of a fluid cell can then be written as

where $f_{0}$ is the equilibrium Bose or Fermi-distribution as a function of particle momentum $p^{\mu}=(E,\mathbf{p})$ for given chemical potential $\mu$ and local temperature $T$,

and $\delta f$ is a small correction that accounts for the out-of-equilibrium behavior. We neglect chemical potentials in this study and set $\mu=0$ here, but note that realistic chemical potentials for stable hadrons were used in [1]. The general form of the correction term is [20, 1]

Here the shear stress tensor $\pi^{\mu\nu}$ has been expressed by its Navier-Stokes approximation, $\pi^{\mu\nu}=2\eta\sigma^{\mu\nu}$, where $\sigma^{\mu\nu}$ is the traceless gradient tensor, defined as

We have used the notation $\nabla^{\mu}=\Delta^{\mu\nu}\partial_{\nu}$, with $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$, for the derivative perpendicular to the flow field vector $u^{\mu}$. In the following we will use the standard choice $\lambda=2$ for the residual momentum dependence of the correction, which is widely used in relativistic viscous hydrodynamics [21] including MUSIC. In order to guarantee the applicability of the equations in this section, the correction term $\delta f$ needs to be small. We have checked that $\delta f(r,p)$ is less than 20% of $f_{0}(r,p)$ for the majority of transverse momentum bins within the fit ranges, with very few bins receiving corrections up to 35%. These corrections are smaller than the necessary upper bound $\delta f/f_{0}\lesssim 1$ [22].

When solving the viscous fluid dynamics equations of motion, numerical stability requires second order gradient terms to be included, leading to equations of motion for the shear stress tensor $\pi^{\mu\nu}$ and bulk stress $\Pi$ [23, 24, 21, 25, 26]. At freeze-out, it is then convenient to compute $\delta f$ directly from the shear stress tensor $\pi^{\mu\nu}$. On the other hand, for the blast wave it is more practical to utilize the Navier-Stokes approximation $\pi^{\mu\nu}=2\eta\sigma^{\mu\nu}$ and to compute viscous corrections using Eq. (4). In that case $\delta f$ is calculated simply from the flow field, which can be independently constrained by fits to flow data and the specific shear viscosity $\eta/s$ of nuclear matter. The differences between the two approaches of calculating $\delta f$ at freeze-out, $\pi^{\mu\nu}$ vs Navier-Stokes, are parametrically small in situations of small gradients towards the end of the time evolution. However, they could still be noticeable at freeze-out in realistic systems 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 hypersurface $\Sigma$ in Minkowski space is given by the Cooper-Frye formula [27]

where $g$ is the degeneracy factor for a given particle and $d\Sigma^{\mu}$ is the forward normal vector on the freeze-out hypersurface. The momentum vector in the laboratory frame is written as usual, $p^{\mu}=(M_{T}\cosh Y,P_{T}\cos\psi,$ $P_{T}\sin\psi,M_{T}\sinh Y)$, in terms of the transverse momentum $P_{T}$, the longitudinal momentum rapidity $Y$ and the azimuthal angle $\psi$ in the transverse plane. $M_{T}^{2}=P_{T}^{2}+M^{2}$ 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_{\mathrm{fo}}$. In contrast, in fluid dynamics this isothermal hypersurface, as well as the flow field $u^{\mu}$ on it can be computed in the simulation itself. For the blast wave we have to choose ansätze for both.

Following Ref. [4] we assume that freeze-out from an isothermal hypersurface at temperature $T$ can be approximated by freeze-out from a hypersurface at constant proper longitudinal time $\tau$. We enforce boost invariance, which is a good approximation for nuclear collisions around midrapidity at top RHIC and LHC energies and is also often found in fluid dynamic calculations. To keep the blast wave simple we have to restrict ourselves to describing smooth fluid dynamics 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 the coordinate axes such that the impact parameter $b$ of the collision is measured along the $x$-axis. In the following we use the reduced radius $\rho=\sqrt{x^{2}/R_{x}^{2}+y^{2}/R_{y}^{2}}$ together with the azimuthal angle $\theta$, with $\tan\theta=(R_{x}y)/(R_{y}x)$, and space time rapidity $\eta_{s}=1/2\log[(t+z)/(t-z)]$ to carry out the integral over the hypersurface. Restricting ourselves to hadrons measured at midrapidity $Y=0$ and changing to convenient coordinates we obtain the final expression for the particles from the blastwave:

Next, we have to make an ansatz for the collective flow field. The general parameterization is

where $\eta_{T}$ is the transverse rapidity in the $x-y$-plane, and $\phi_{u}$ is the azimuthal angle of the flow vector in the transverse plane. Boost invariance fixes the longitudinal flow rapidity to be equal to the space time rapidity $\eta_{s}$. We follow Retiere and Lisa and choose to model the transverse flow velocity $v_{T}=\tanh\eta_{T}$ as [4]

which encodes a Hubble-like velocity ordering with an additional shape parameter $n$. $\alpha_{0}$ is the average velocity on the boundary $\rho=1$, and $\alpha_{2}$ parameterizes the elliptic deformation of the flow field built up from the initial elliptic spatial deformation of the system. Flow vectors tend to be tilted towards the smaller axis of the ellipse. In the RL approach they are chosen to be perpendicular to the elliptic surface at $\rho=1$, i.e. $\tan\phi_{u}=R_{x}^{2}/R_{y}^{2}\tan\phi$, where $\phi=\arctan y/x$ is the azimuthal angle of the position $r^{\mu}$.

With a parameterization of the flow field at hand, the next step is the calculation of the gradient tensor $\sigma^{\mu\nu}$. This has been carried out in [1] and we refer the reader to the details in that reference. We want to point out that temporal derivatives are calculated using ideal fluid dynamic equations of motion rather than the free-streaming approximation [11, 12]. This introduces the nuclear matter equation of state, specifically the speed of sound squared $c_{s}^{2}$ into our calculation of $\delta f$.

The blast wave ansatz has several parameters which allow us to adjust the flow field and the hypersurface, as well as the temperature, specific shear viscosity, and speed of sound squared at freeze-out. The full set of parameters is $\tilde{\mathcal{P}}=(\tau_{\mathrm{fo}},R_{x},R_{y},T_{\mathrm{fo}},n,\alpha_{0},\alpha_{2},\eta/s,c_{s}^{2})$. As in Ref. [1] we drop $c_{s}^{2}$ from this list and rather use guidance on the hadronic matter equation of state from existing literature which gives $c_{s}^{2}=0.16\sim 0.17$ for T = 110-140 MeV [28]. We also use a simple geometric argument for the time-dependence of the system size along the impact vector, $R_{x}=(R_{0}-b/2)+\tau c_{\tau}(\alpha_{0}+\alpha_{2})$, to determine $R_{x}$. Here $R_{0}$ is the radius of the colliding nucleus, $b$ is the impact parameter and $c_{\tau}$ relates the time-averaged surface velocity to the final velocities $\alpha_{0}$ and $\alpha_{2}$. The value of $c_{\tau}$ can be inferred from typical radial velocity-vs-time curves obtained in fluid dynamic simulations [21] and we set $c_{\tau}=0.65$ here. Note that the viscous blast wave depends on the parameters $\tau_{\mathrm{fo}}$, $R_{x}$ and $R_{y}$ separately, and not just on the total volume $\tau R_{x}R_{y}$ and the elliptic deformation $R_{y}/R_{x}$, as is the case for the ideal blast wave. With $R_{x}$ fixed in each instance we are left with the $R_{y}/R_{x}$ as a fit parameter.

Despite the large number of parameters it is clear that the blast wave has introduced significant simplifications compared to the freeze-out calculated in fluid dynamics. They come mostly from the simplified shape of the hypersurface (constant proper time) and the spatial structure of the flow field. Two more major approximations are made for sake of simplicity. First, resonance decays are usually neglected in blast wave calculations, and only hadrons stable under strong decays are taken into account. Secondly, correction terms to the particle distribution $f$ due to bulk stress have been ignored. They could in principle be added and we plan to do so in the future. We summarize the five major approximations compared to fluid dynamic freeze-out in the following list:

• •

  Navier-Stokes approximation used for $\delta 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 dynamic simulations compared to smooth fluid dynamics need to be considered separately. They are not included in the study below. The same is true for deviations of state-of-the-art 3+1D fluid dynamics from the boost-invariant 2+1D fluid dynamics used here. The effects of fluctuations and breaking of boost invariance have already been studied within fluid dynamics [29, 30] and can be added to the considerations in this work.

We use the viscous hydrodynamics code MUSIC to simulate averaged nuclear collisions at RHIC and LHC energies at various impact parameters. MUSIC is a relativistic second-order viscous hydrodynamics code for heavy ion collisions [18, 29, 19]. We choose boost-invariant (2+1)D mode consistent with the boost-invariant blast wave set-up. We use the built-in optical Glauber model to generate initial conditions with the appropriate nucleon-nucleon cross section and an overall normalization roughly consistent with pertinent multiplicity data for Au+Au collisions at RHIC at $\sqrt{s_{NN}}=200$ GeV and Pb+Pb collisions at the LHC at $\sqrt{s_{NN}}=2.76$ TeV. Other collision systems can be treated similarly. We use the equation of state (EOS) s95p-v1.2 in MUSIC, and the default MUSIC bulk viscosity. The shear viscosity over entropy ratio $\eta/s$ is chosen to be a constant which we vary as a parameter. We freeze out at pre-determined temperatures $T_{\mathrm{fo}}$ and compute the final spectra and elliptic flow for pion, kaons and protons, including resonance decays and including viscous corrections to freeze-out. The detailed MUSIC settings are documented in Appendix A.

Recall that we want to establish a map from the temperature $T_{\mathrm{fo}}^{(\text{extr})}$ and specific shear viscosity $(\eta/s)^{(\text{extr})}$ extracted from a blast wave fit of the pseudodata to the true values $T_{\mathrm{fo}}^{(\text{true})}$ and $(\eta/s)^{(\text{true})}$ used in the generation of these pseudodata. To focus on the relevant region tested in heavy ion collisions, we choose nine points in the $T_{\mathrm{fo}}^{(\text{true})}$-$(\eta/s)^{(\text{true})}$-plane for simulations at RHIC energy, such that the corresponding fitted values $(T_{\mathrm{fo}}^{(\text{extr})},(\eta/s)^{(\text{extr})})$ are roughly consistent with the values extracted from RHIC data in Ref. [1], for some impact parameter $b$. Thus, a single set of parameters to run MUSIC consists of an impact parameter $b$ and values $T_{\mathrm{fo}}^{(\text{true})}$, $(\eta/s)^{(\text{true})}$. We run MUSIC and perform a blast wave fit to the resulting hadron spectra and elliptic flow for all nine such sets at RHIC energies. Similarly, we choose eight points for Pb+Pb collisions at LHC. The 9+8 sets of parameters are shown in Tab. 1. We will refer to these points often as Set I.

The $P_{T}$-range of the pseudodata generated by MUSIC consists of the interval from 0 to 3 GeV/$c$. However, we have to restrict the range of data in which we fit the blast wave to pseudodata. At very low momenta hadron spectra tend to be dominated by resonance decays, which are not included in the blast wave. At momenta which are too high, hadron production receives viscous corrections larger than what can be reliably described by the Navier-Stokes approximation. Such restrictions of fit ranges for blast waves seem to be good practice in the literature. They are also used in Ref. [1]. The fit ranges used for blast wave fits of the pseudodata here are shown in Tab. 2. They are inspired by what was 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 we choose uncertainties in line with error bars in the pertinent available experimental data. We assign 5% uncertainty and 2% uncertainty to pseudodata spectra and $v_{2}$, respectively. We add a pedestal of 0.002 to the uncertainty of $v_{2}$ for realistic error bars at smaller $P_{T}$ where $v_{2}$ is very small. The choices of fit ranges and error bars introduce additional uncertainties in our analysis. The dependence of blast wave fits on the choice of fit ranges were studied in Ref. [1]. The additional uncertainty due to the choice of error bars for the pseudodata is studied later in this work by varying the size of the assumed error bars, see Tab. 7.

We utilize Bayesian inference to extract likelihoods for the relevant parameters. We use the statistical analysis package from the Models and Data Analysis Initiative (MADAI) project [31, 32]. The MADAI package includes a Gaussian process emulator and a Markov Chain Monte Carlo. We use $N=500$ training points for each Gaussian emulator. Closure tests find the errors in the Gaussian emulator to be negligible compared to the assumed uncertainties in the pseudodata.

We fit the following, reduced set of blast wave parameters in the Bayesian analysis: $\mathcal{P}=(T_{\mathrm{fo}},R_{y}/R_{x},\alpha_{0},\alpha_{2},\eta/s)$. Chemical potentials $\mu$ are set to zero in MUSIC and in the blast wave. We will rather fit a normalized $P_{T}$-spectrum and do not utilize the absolute yield of hadrons to reduce complexity. To further simplify the analysis we fix the radial shape parameters $n$ and the freeze-out time $\tau$ by choosing values close to those extracted from RHIC and LHC data in Ref. [1], see Tab. 3 and Tab. 4.

As an example, we take a look at the case of the input parameter set $T_{\mathrm{fo}}^{(\text{true})}=130$ MeV, $(\eta/s)^{(\text{true})}=2.51/(4\pi)$ in Au+Au collisions. Fig. 1 shows posterior distributions and correlations of the values $\mathcal{P}$ extracted from the Bayesian analysis. The likelihoods for all parameters exhibit well defined peaks. The preferred values (defined as the means) for the set $\mathcal{P}$ of parameters in this case are $T_{\mathrm{fo}}=116.4$ MeV, $\alpha_{0}$=0.74$c$, $R_{y}/R_{x}$=1.10, $\alpha_{2}$=0.045$c$, $\eta/s$=2.5/4$\pi$, with fixed parameters $\tau=8.4$ fm/$c$ and $n$=0.88.

We proceed analogously for the other $T_{\mathrm{fo}}^{(\text{true})}$, $(\eta/s)^{(\text{true})}$ Au+Au points, and for the Pb+Pb points from Set I. The results are summarized in Tabs. 3 and 4. Let us first discuss the quality of the fit results for the observables. Figs. 2 and 3 show the identified spectra and $v_{2}$ computed from the blastwave with the preferred values together with the MUSIC pseudodata for all impact parameters in Au+Au and Pb+Pb. The figures demonstrate that the fits to identified hadron spectra and elliptic flow are working quite well across the parameters chosen for Set I.

In this work, the declared measures for success are the accuracies of the extracted specific shear viscosity and temperature compared to the true values used in the simulations. For our example set $T_{\mathrm{fo}}^{(\text{true})}=130$ MeV, $(\eta/s)^{(\text{true})}=2.51/(4\pi)$ in Au+Au the extracted specific shear viscosity reflects the true value quite accurately — 2.47 vs 2.51 — while the extracted freeze-out temperature of 116 MeV is lower than the true value of 130 MeV. Looking across the entire set of fitted values we note that agreement between the extracted $\eta/s$ and the true value stays within 10%, and often much closer, in all cases except for the highest freeze-out temperatures. In contrast, there is a noteworthy trend for the extracted freeze-out temperatures. The agreement between true and extracted values is rather good for central collisions, in particular for Pb+Pb, and accuracy decreases toward peripheral collisions. For the latter, the extracted temperatures are underestimated by up to 20 MeV.

To understand this behavior further it is informative to compare some details of the two simulated fireballs. The left panel of Fig. 4 shows the freeze-out cells in fluid dynamics vs the corresponding freeze-out hypersurface in the blast wave in the $t-x$-plane in Minkowski space. The ”muffin” shape of the $T=T_{\mathrm{fo}}$ hypersurface in fluid dynamics, which emerges from the competition of cooling and radial flow, is well known. On the other hand, the choice made in blastwaves is to use an average constant time hypersurface. This comparison makes it clear that fields from fluid dynamics and blast waves should not be expected to line up point-by-point. Rather, what one can reasonably hope for is for a matching of global features and of the relevant numerical averages.

With that in mind we move to a comparison of components of the flow velocites $u^{\mu}(r)$ at freeze-out. $u^{x}(x)$ for $y=0$ and $u^{y}(y)$ for $x=0$ for our example set $T_{\mathrm{fo}}^{(\text{true})}=130$ MeV, $(\eta/s)^{(\text{true})}=2.51/(4\pi)$ in Au+Au are shown for fluid dynamic cells at freeze-out and for the viscous blast wave in the right panel of Fig. 4. We make note of two observations: The sizes $R_{x}$ and $R_{y}$ of the fireballs in $x$- and $y$-directions, respectively, indicated by the endpoints of the blast wave curves and the rightmost fluid cells, match up well. On the other hand, it looks like the blastwave overestimates the radial flow in both directions compared to fluid dynamics. However, the weights with which the different fluid cells contribute to the observed particle yields need to be taken into account. The sequences of MUSIC fluid cells fold in on themselves at large values of $x$ and $y$ because of the muffin shape of the hypersurface. Eq. (7) exhibits a Jacobian $\tau\rho$ which makes the outer part of the muffin edge the largest contributer to particle production. This is where the cells with the largest radial flow are located, closer to their counterpart in the blast wave. After accounting for this effect $\langle u^{x}\rangle\approx 0.65$ for the freeze-out cells along the $x$-axis in fluid dynamics, while the corresponding value for the blast wave is $\langle u^{x}\rangle\approx 0.76$. The full expression in Eq. (7) is used for the averaging here, with the restriction $\eta=0$, $\theta=0$. Thus a residual mismatch remains, with the blast wave overestimating the radial flow. This is consistent with the lower fitted temperatures in the blast wave. For momenta which are not too small the main effect of flow is an effective blue-shift of the temperature which is also reflected by the anti-correlation between $T$ and $\alpha_{0}$ in Fig. 1. This explains why the pseudodata for this set can be fit rather well in the chosen range, although the extracted temperature is too low.

The effect of lower temperatures being compensated by larger radial flow is largest in peripheral collisions and dissappears in more central collisions. This is confirmed by Fig. 5 which shows the average flow values $\langle u^{x}\rangle$ and $\langle u^{y}\rangle$ along the $x$- and $y$-axes, respectively, discussed earlier for the example set, for all sets in Au+Au as a function of the temperature used in the MUSIC simulation. The highest temperatures, corresponding to the most peripheral collisions, see the largest deviation between fluid dynamics and blastwave fits for both radial flow and the extracted temperature. It is not quite clear why the peripheral blastwave fits prefer to trade a lower temperature against a higher flow field. One can speculate that the simplistic $\tau$-dependence of the hypersurface in the blastwave obviously becomes a worse approximation in peripheral collisions. On a practical level, a better separation of flow and temperature effects on observables, and thus better fits for the temperature, could presumably be achieved by employing a wider fit range, in particular through the inclusion of smaller momenta. To this end, resonance decays would have to be taken into account for the viscous blast wave, which is possible but beyond the scope of the current model. Here, we will simply document the biases in the fits and prepare the tools to remove them.

Before we attempt to correct the extracted temperatures let us briefly look at the shear stress in both MUSIC and the blastwave. In Fig. 6 we compute the off-diagonal component $\pi^{xy}$ divided by the enthalpy $e+p$ as a function of the space-time angle $\phi$, averaged over the remaining coordinates of the hypersurface for both MUSIC fluid dynamics and the blastwave for the preferred parameters. We show the results for Au+Au events at $T_{\mathrm{fo}}^{(\text{true})}=115$ MeV and $140$ MeV, corresponding to rather central and peripheral collisions, respectively. Interestingly, the average shear stress in the fluid dynamic simulation is represented quite well by the viscous blast wave for both centralities. This comparison provides context for the good agreement we find for the extracted values of the specific shear viscosity.

We summarize the results of this section in Figs. 7 and 8 for Au+Au and Pb+Pb collisions, respectively. The plots show the correlation between true and extracted values of our two main observables, the freeze-out temperature (right panels) and specific shear viscosity (left panels). Diagonal lines for perfect correlation (”y=x”) are added to guide the eye. The uncertainties for each point are combined uncertainties from the posterior likelihoods from the Bayesian analysis, and systematic uncertainties, see Appendix B.

We now interpret the extracted values of freeze-out temperatures and specific shear viscosities as images of the original values $T_{\mathrm{fo}}^{(\text{true})}$ and $(\eta/s)^{(\text{true})}$ under a map that is determined by the approximations made in the blast wave approach. The left panels of Figs. 7 and 8 make it clear that — within reasonable accuracy — it is only necessary to consider the mapping of the true temperature onto the extracted temperature, $T_{\mathrm{fo}}^{(\text{true})}\mapsto T_{\mathrm{fo}}^{(\text{extr})}$ and that it is a sufficiently accurate assumption that $(\eta/s)^{(\text{true})}\approx(\eta/s)^{(\text{extr})}$. Moreover, Figs. 7 and 8 tell us that a linear map should have the desired level of accuracy.

By minimizing

for each input paramater point $i$, where $\sigma_{i}$ is the uncertainty of the extracted temperature at each point, we obtain an approximation $T_{\mathrm{fo}}^{(\text{extr,lin})}$ to the extracted temperature by linear regression. For Au+Au we obtain $T_{\mathrm{fo}}^{(\text{extr,lin})}=0.36T^{(\text{true})}_{\mathrm{fo}}+72.5\,\text{MeV}$. Similarly, for Pb+Pb we have $T_{\mathrm{fo}}^{(\text{extr,lin})}=0.47\,T_{\mathrm{fo}}^{(\text{true})}+62.4\,\text{MeV}$. These linear functions are shown in the right panels of Figs. 7 and 8 in addition to the true and extracted values. The linear maps provide a description of the extracted temperatures within the uncertainty bars. For all practical purposes we can thus identify $T_{\mathrm{fo}}^{(\text{extr,lin})}$ and $T_{\mathrm{fo}}^{(\text{extr})}$ for the rest of this section.

To further validate the above map, we introduce two more sets of points in the $T_{\mathrm{fo}}^{(\text{true})}$-$(\eta/s)^{(\text{true})}$-plane. We arrive at these sets by going through Set I and decreasing or increasing, respectively, $\eta/s$ by $0.5/(4\pi)$, keeping the impact parameters and freeze-out temperatures in Tab. 1. The resulting sets will be called Set II and III, respectively. Following the same process, we generate fluid dynamic pseudodata for these sets by running MUSIC, and extract $T_{\mathrm{fo}}$ and $\eta/s$ from subsequent blast wave fits of the pseudodata. The results are detailed in Appendix C, since we only need to discuss the final conclusion. We find that the extracted values of $\eta/s$ are again consistent with the true values within uncertainties. The general trends seen for the extracted freeze-out temperatures $T_{\mathrm{fo}}$ using Set I are the same for the two new sets. As a result, the maps we obtain for $T_{\mathrm{fo}}^{(\text{true})}\mapsto T_{\mathrm{fo}}^{(\text{extr})}$ for the new sets are very close to the ones obtained earlier, as shown in the right panels of Figs. 9 and 10 for the Au+Au case. The Pb+Pb case is analogous. We conclude that the map for set I already provides a reliable description of the extracted temperatures over a wide range of parameters.

We are now in a position to remove the bias on the extracted temperature that comes from using the approximations of the blast wave. By inverting the maps and applying them to the extracted temperatures, we can find ”corrected” values, $(T_{\mathrm{fo}}^{(\text{extr})},(\eta/s)^{(\text{extr})})\mapsto(T_{\mathrm{fo}}^{(\text{corr})},(\eta/s)^{(\text{corr})})$. The inverted maps are

By virtue of $T_{\mathrm{fo}}^{(\text{extr,lin})}\approx T_{\mathrm{fo}}^{(\text{extr})}$ we have $T_{\mathrm{fo}}^{(\text{corr})}\approx T_{\mathrm{fo}}^{(\text{true})}$. Thus we have an algorithm for predicting the correct temperature. Although our work has focused on freeze-out temperature and specific shear viscosity, the analysis in this work could be repeated for other quantities extracted from blast wave fits in a straight forward way.

As an application of the above procedure, we point the readers to our extraction of specific shear viscosity from data in [1]. For example, considering the extracted value from ALICE Pb+Pb at 2.76 TeV in the 50-60% centrality bin we extract the raw values $(T,\eta/s)^{\text{extr}}=(130.1,1.66/(4\pi))$ and subsequently obtain the corrected values $(T,\eta/s)^{\text{corr}}=(144.3,1.66/(4\pi))$. The same process can be applied to other centrality bins of ALICE Pb+Pb as well as RHIC Au+Au at 200 GeV. As a result of the correction $\eta/s$ drops more slowly with increasing temperature, and a value close to $4\pi\eta/s=1$ is only reached around $T\approx 150$ MeV. For further discussions of the implications of this particular result we refer the reader to reference [1].

In this paper we have discussed differences and complementarity between viscous blast wave descriptions and fluid dynamic simulations in the context of high energy nuclear collisions. We have carried out a systematic study of parameters at kinetic freeze-out set in fluid dynamic simulations which have subsequently been extracted by blast wave fits to hadron spectra and elliptic flow. We find that viscous blast wave fits correctly reproduce broad trends of the true values set in the simulations, and that they qualitatively agree with the numerical values set in fluid dynamics. To be more precise, extractions of the specific shear viscosity $\eta/s$ tend to be accurate within expected uncertainties. This is backed up a direct comparison of relevant average shear viscosity components between fluid dynamics and blast waves. On the other hand, we find the extracted temperatures and average radial flow velocities deviate from their true values. While still somewhat accurate in central collisions, freeze-out temperatures are underestimated, and average radial flow velocities overestimated in peripheral collisions. However, note that the overall fit quality of hadron spectra and elliptic flow remains excellent within the chosen ranges.

The quality of blast wave fits can be improved by understanding and quantifying these deviations. One can establish a map from the true temperature values to the ones extracted from the blast wave fits. The inverse of this map can be applied to arrive at corrected fit values for the temperatures. One can argue that these should be close to the true values and thus correct the biases in blast wave fits. Our study focuses on fitted shear viscosities and temperatures, but it can be readily extended to other physical parameters following the blueprint laid out here. The benefit of removing biases in blastwave fits has been demonstrated for the case of the extraction of the specific shear viscosity as a function of temperature from data in Ref. [1]. In the latter work, an added complexity is the application of the procedure to experimental data, which leaves the additional question of the quality of hydrodynamic modeling of nuclear collisions, as discussed in the original work.

The main point of this paper is the demonstration that average behavior in fluid dynamics is generally well described by viscous blastwaves, and that remaining deviations can be systematically removed. An obviouly interesting question emerges regarding the universality and systematic behavior of the specific corrections discussed. The maps for freeze-out temperature and shear viscosity discussed here might find direct application if the same blastwave and fit ranges are used for diffent fluid dynamic simulations. However, caution would dictate that interested readers should always check and if needed establish a custom map for their own setup. E.g., it is rather obvious that significant changes to the fit ranges, or additional observables, will change the map between true and extracted values.

The viscous blastwave itself could also be improved. The most straight forward changes could be the addition of bulk stress corrections and the inclusion of important hadronic resonances and their decays. The inclusion of resonances has been discussed in [9]. Indeed, the fits of ALICE data performed in that work lead to somewhat higher freeze-out temperatures than the raw values found in [1]. Qualitatively, this trend is compatiblen with the way the corrections found here move raw values. However, a closer inspection of the dependence on fit ranges, which are different in both cases, would be necessary for an qualitative comparison.

Finally, let us recall that simplicity is key to blast wave fits. While fluid dynamics has matured rapidly and is numerically feasible, blast waves are cheaper in labor and computational requirements. They remain a quick and useful tool for certain data analyses. The formalism laid out here can add significantly to their reliability.