Estimation of the shear viscosity at finite net-baryon density from A+A collision data at \sqrt{s_{\mathrm{NN}}}=7.7–200 GeV

# Estimation of the shear viscosity at finite net-baryon density from A+A collision data at √sNN=7.7–200 GeV

## Abstract

Hybrid approaches based on relativistic hydrodynamics and transport theory have been successfully applied for many years for the dynamical description of heavy ion collisions at ultrarelativistic energies. In this work a new viscous hybrid model employing the hadron transport approach UrQMD for the early and late non-equilibrium stages of the reaction, and 3+1 dimensional viscous hydrodynamics for the hot and dense quark-gluon plasma stage is introduced. This approach includes the equation of motion for finite baryon number, and employs an equation of state with finite net-baryon density to allow for calculations in a large range of beam energies. The parameter space of the model is explored, and constrained by comparison with the experimental data for bulk observables from SPS and the phase I beam energy scan at RHIC. The favored parameter values depend on energy, but allow to extract the effective value of the shear viscosity coefficient over entropy density ratio in the fluid phase for the whole energy region under investigation. The estimated value of increases with decreasing collision energy, which may indicate that of the quark-gluon plasma depends on baryochemical potential .

## I Introduction

Ultra-relativistic heavy ion collisions allow to investigate the properties of strongly interacting matter under extreme conditions. At high temperatures and/or high net-baryon densities a new state of matter, the so-called quark-gluon plasma, QGP, is formed. The two main goals of heavy ion research are the exploration of the phase diagram of quantum chromodynamics and the determination of transport coefficients of this new state of matter.

The studies of high energy heavy-ion collisions at the Large Hadron Collider (LHC) at CERN and the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory have revealed that the quark-gluon plasma behaves like an almost perfect fluid. In recent years, so-called hybrid approaches (1); (2); (3); (4); (5) based on (viscous) relativistic hydrodynamics for the hot and dense stage coupled to hadron transport approaches for the decoupling stage of the reaction have been applied with great success to extract average values of the shear viscosity over entropy ratio . The results are very close to the conjectured universal limit of , based on the anti–de Sitter/conformal field theory (AdS/CFT) correspondence (6). For example, the values extracted in Ref. (7) are for collisions at RHIC, and at the LHC.

One expects the formation of partonic matter in heavy ion collisions at ultra-relativistic energies (see, e.g., Ref. (8)). However, it is unknown at what collision energy the transition from hadronic to partonic matter sets in. In addition, as the collisions at lower energies probe the phase diagram at larger net-baryon densities, it may be possible to experimentally see signs of the theoretically predicted critical point (9) and the first-order phase transition beyond it. To investigate these questions the so-called beam energy scan (BES) programs at SPS (NA49, NA61 experiments) and at RHIC (STAR, PHENIX experiments) were started. One of the surprises of the stage I of the BES program at RHIC has been that the -differential elliptic flow, , of charged hadrons does not change significantly when the collision energy is reduced from to GeV (10). The large values of elliptic flow measured at GeV collisions were taken as a sign of very low shear viscosity of the matter formed in these collisions. Thus, the large measured in collisions at lower energy leads to the question how changes as function of net-baryon density and baryochemical potential  (11).

Unfortunately, many of the hydrodynamical and hybrid models used to model collisions at full RHIC and LHC energies are not directly applicable to collisions at RHIC BES and CERN SPS energies, nor to collisions at even lower energies in the future at Facility for Antiproton and Ion Research (FAIR), Nuclotron-based Ion Collider Facility (NICA) and the stage II of the BES program at RHIC. The simplifying approximations of boost invariance and zero net-baryon density are not valid, and different kinds of non-equilibrium effects play a larger role. To overcome these limitations, a novel hybrid approach has been developed. This approach is based on the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) transport (12) for the non-equilibrium early and late stages, and on a (3+1)-dimensional viscous hydrodynamical model (13) for the hot and dense stage of the reaction.

In this paper, this approach is applied to extract the shear viscosity coefficient over entropy density ratio of strongly interacting matter from the heavy-ion collision data at RHIC beam energy scan energies in the broad range GeV. The details of the model are explained in Section II, and the generic effects of finite shear viscosity on the hydrodynamical expansion are described in Section III. The sensitivity of particle yields and spectra to the parameters for the initial and final state transitions is explored in Section IV. Section V contains the main results of this work including the estimated values of the effective shear viscosity over entropy density ratio as a function of beam energy. Finally, the main conclusions are summarized and an outlook on future work is given in Section VI.

## Ii Model Description

Our hybrid approach combines the UrQMD transport model (12) for the early and late stages of the evolution with a dissipative hydrodynamical model, called vHLLE (13), for the hot and dense stage. The distinguishing features of the present model are that the fluid dynamical expansion is solved numerically in all three spatial dimensions without assuming boost invariance nor cylindrical symmetry, the equations of motion for finite net-baryon and charge densities are explicitly included and, in contrast to the standard UrQMD hybrid approach (UrQMD-3.4 at urqmd.org) (3); (14), dissipation in the form of shear viscosity is included in the hydrodynamical stage. Different to our previous studies (17); (18), event-by-event fluctuations are now included. The hadronic cascade operates with the full phase-space distribution of the final particles which allows for a proper comparison to experimental data.

### ii.1 Pre-thermal Phase

The UrQMD string/hadronic cascade is used to describe the primary collisions of the nucleons, and to create the initial state of the hydrodynamical evolution. The two nuclei are initialized according to Woods-Saxon distributions and the initial binary interactions proceed via string or resonance excitations, the former process being dominant in ultrarelativistic collisions (including the BES collision energies). All the strings are fragmented into hadrons before the transition to fluid phase (fluidization) takes place, although not all hadrons are yet fully formed at that time, i.e., they do not yet have their free-particle scattering cross sections, and thus do not yet interact at all. The hadrons before conversion to fluid should not be considered physical hadrons, but rather marker particles to describe the flow of energy, momentum and conserved charges during the pre-equilibrium evolution of the system. The use of UrQMD to initialise the system allows us to describe some of the pre-equilibrium dynamics and dynamically generates event-by-event fluctuating initial states for hydrodynamical evolution.

The interactions in the pre-equilibrium UrQMD evolution are allowed until a hypersurface of constant Bjorken proper time is reached, since the hydrodynamical code is constructed using the Milne coordinates , where  (13). The UrQMD evolution, however, proceeds in Cartesian coordinates , and thus evolving the particle distributions to constant means evolving the system until large enough time in such a way that the collisional processes and decays are only allowed in the domain . The resulting particles on surface are then propagated backwards in time to the surface along straight trajectories to obtain an initial state for the hydrodynamic evolution.

The lower limit for the starting time of the hydrodynamic evolution depends on the collision energy according to

 τ0=2R/√(√sNN/2mN)2−1, (1)

which corresponds to the average time, when two nuclei have passed through each other, i.e., all primary nucleon-nucleon collisions have happened. This is the earliest possible moment in time, where approximate local equilibrium can be assumed. The values calculated according to this formula are plotted in Fig. 1.

To perform event-by-event hydrodynamics using fluctuating initial conditions every individual UrQMD event is converted to an initial state profile. As mentioned, the hadron transport does not lead to an initial state in full local equilibrium, and the thermalization of the system at has to be artificially enforced. The energy and momentum of each UrQMD particle at is distributed to the hydrodynamic cells assuming Gaussian density profiles

 ΔPαijk =Pα⋅C⋅exp(−Δx2i+Δy2jR2⊥−Δη2kR2ηγ2ητ20) (2) ΔN0ijk =N0⋅C⋅exp(−Δx2i+Δy2jR2⊥−Δη2kR2ηγ2ητ20), (3)

where , , are the differences between particle’s position and the coordinates of the hydrodynamic cell , and is the longitudinal Lorentz factor of the particle as seen in a frame moving with the rapidity . The normalization constant is calculated from the condition that the discrete sum of the values of the Gaussian in all neighboring cells equals one. The resulting and are transformed into Milne coordinates and added to the energy, momentum and baryon number in each cell. This procedure ensures that in the initial transition from transport to hydrodynamics energy, momentum and baryon number are conserved.

For the present study energy and momentum of the initial particles are converted at into a perfectly equilibrated fluid, i.e., the initial values for the viscous terms in the energy-momentum tensor are set to zero: . In other words the components of the energy-momentum tensor stay the same, but the components change, when we switch from UrQMD to the fluid. Thus, we do not consider how much the energy-momentum tensor of UrQMD deviates from the ideal fluid energy-momentum tensor, but leave this topic for further studies.

One typical example of the initial energy density distributions in the transverse plane at midrapidity for one event is presented in Fig. 2. The parameters and regulate the granularity of the initial state. At the same time they influence the initial entropy of the hydrodynamic evolution, while the total initial energy and momentum are always fixed to be equal to the energy and momentum of the pre-equilibrium UrQMD event. The dependence of the final results on these two parameters is discussed later in Section IV.

### ii.2 Hydrodynamic Evolution

The (3+1)-dimensional viscous hydrodynamical code vHLLE is described in full detail in Ref. (13). We repeat here only its main features. The code solves the usual local energy-momentum conservation equations

 ∂;νTμν =0, (4) ∂;νNνB,Q =0, (5)

where and are net baryon and electric charge currents respectively, and the semicolon denotes the covariant derivative. The calculation1 is done in Milne coordinates , where and .

In the Israel-Stewart framework of causal dissipative hydrodynamics (19), the dissipative currents are independent variables. For the purpose of the present work we set the bulk viscosity to zero, . We work in the Landau frame, where the energy diffusion flow is zero, and neglect the baryon and charge diffusion currents, which is equivalent to zero heat conductivity. For the shear stress evolution we choose the relaxation time , the coefficient , and approximate all the other coefficients (20); (21) by zero. For the shear-stress tensor we obtain the evolution equation

 ⟨uγ∂;γπμν⟩=−πμν−πμνNSτπ−43πμν∂;γuγ, (6)

where the brackets denote the traceless and orthogonal to part of the tensor and is the Navier-Stokes value of the shear-stress tensor.

Another necessary ingredient for the hydrodynamic stage is the equation of state (EoS) of the medium. We use the chiral model EoS (22), which features correct asymptotic degrees of freedom, i.e., quarks and gluons in the high temperature and hadrons in the low-temperature limits, crossover-type transition between confined and deconfined matter for all values of and qualitatively agrees with lattice QCD data at .

The tests to confirm the accuracy of the code have been reported in Ref. (13). In particular the solutions have been compared to the ideal Gubser solution (23) and to a numerical solution of dissipative hydrodynamics calculated using the VISH2+1 hydro code (24).

### ii.3 Particlization and Hadronic Rescattering

It is well known that hydrodynamics loses its validity when the system becomes dilute. To deal with this problem we apply the conventional Cooper-Frye prescription (25) to particlize the system (convert the fluid to individual particles) at a hypersurface of constant local rest frame energy density, and use the UrQMD cascade to describe the further evolution of these particles. This switching hypersurface is evaluated during the hydrodynamic evolution using the Cornelius routine (26), and as a default value for the switching density we use  GeV/fm, which in the chiral model EoS corresponds to  MeV at . At this energy density the crossover transition is firmly on the hadronic side, but the density is still a little higher than the chemical freeze-out energy density suggested by the thermal models (27). Thus the hadronic transport can take care of both chemical and kinetic decoupling of hadrons. We discuss the sensitivity of the results to the value of the switching density in section IV.

As given by the Cooper-Frye prescription, the hadron distribution on each point of the hypersurface is

 p0d3Ni(x)d3p=dσμpμf(p⋅u(x),T(x),μi(x)). (7)

The phase space distribution function is usually assumed to be the one corresponding to a noninteracting hadron resonance gas in or close to the local thermal equilibrium. This is inconsistent with mean fields included in the chiral model EoS used during the evolution. To solve this inconsistency we evaluate the switching surface using the chiral model EoS, but use a free hadron resonance gas EoS to recalculate the energy density, pressure, flow velocity , temperature, and chemical potentials from the ideal part of the energy-momentum tensor and charge currents, and use these values to evaluate the particle distributions on the switching surface. For example the above mentioned temperature of MeV in chiral model EoS at zero baryon density and  GeV/fm, drops to MeV in the free hadron resonance gas EoS. This procedure ensures that the total energy of the produced particles is reasonably close to the overall energy flow through the particlization hypersurface (up to negative contributions to the Cooper-Frye formula), although a small error arises since we use a different energy density to evaluate the position of the surface, and the properties of the fluid on it2. We have checked that in a case of event-averaged initialization, this error is on the level of few percents. In addition, the conservation of energy and momentum in the 3+1 dimensional numerical solution of the fluid-dynamical equations using Milne coordinates is slightly violated as discussed in Refs. (13); (21).

To take into account the dissipative corrections to the distribution function , we use the well-known Grad’s 14-moment ansatz for a single component system, and assume that the correction is the same for all hadron species. We evaluate the particle distribution in the rest frame of the fluid at each surface element using the Cooper-Frye formula

 d3ΔNidp∗d(cosθ)dϕ=Δσ∗μp∗μp∗0Wresidualp∗2feq(p∗0;T,μi)isotropic×[1+(1∓feq)p∗μp∗νπ∗μν2T2(ϵ+p)]Wvisc. (8)

The distribution function in Eq. (8) is expressed in terms of temperature and chemical potential(s), which implies a grand canonical ensemble. This allows to do the particle sampling independently on each surface element. To create an ensemble for particles, we perform the following steps at each element :

• First, the average number of hadrons of every sort is calculated:

 ΔNi=Δσμuμni,th=Δσ∗0ni,th
• For a given , the number of particles to be created is generated according to a Poisson distribution with a mean value .

• For each created particle, the type is randomly chosen based on the probabilities .

• A momentum is assigned to the particle in two steps:

1. The modulus of the momentum is sampled in the local rest frame of the fluid, according to the isotropic part of Eq. (8), and the direction of momentum is picked randomly in solid angle.

2. The correction for or in Eq. (8) is applied via rejection sampling: A random number in the range is generated. If , the generated momentum is accepted, if not, the momentum generating procedure is repeated.

• The particle momentum is Lorentz boosted to the center of mass frame of the system.

• The particle position is taken to be equal to the coordinate of the centroid of the corresponding surface element, except for the spacetime rapidity of the particle, which is uniformly distributed within the longitudinal size of the volume element.

For the current study, no correction over the grand canonical procedure is made to effectively account for the exact conservation of the total baryon/electric charge, strangeness and total energy/momentum for every sampled event3. As a result, these quantities fluctuate from event to event.

The generated hadrons are then fed into the UrQMD cascade. Since the cascade accepts only a list of particles at an equal Cartesian time as an input, the created particles are propagated backwards in time to the time when the first particle was created. The particles are not allowed to interact in the cascade until their trajectories cross the particlization hypersurface.

## Iii Sensitivity to Shear Viscosity

The overall effects of shear viscosity on hydrodynamical expansion have been extensively discussed in the literature (24); (30); (31); (32). Here we show that the results from high energy collisions, e.g., entropy increase, enhancement of transverse and inhibition of longitudinal expansion, and suppression of anisotropies are also manifested in calculations at lower collision energies.

We have carried out event-by-event simulations for different collision energies, centralities, and two fixed values of shear viscosity: (ideal hydro evolution) and . For these simulations we use the values of the Gaussian radii for the particles’ energy/momentum deposition  fm (see Eqs. (2) and (3)). The initial time is chosen according to Eq. (1), however for the collisions at energies equal or higher than  GeV we set fm/.

To reduce the need for CPU time, we use so called oversampling technique, as in Ref. (33). For each collision energy, centrality, and parameter set we have created around 500 hydrodynamic events with randomly generated initial conditions. For each hydrodynamic event, or transition hypersurface, we generate final state events, which results in events used to calculate observables. We have checked that the oversampling procedure does not significantly affect the final observables by creating 1000 or 10000 hydrodynamic events, with respectively, for several parameter sets. In both cases the calculated observables agreed within statistical errors.

The available experimental data set for the basic bulk hadron observables at the BES energies is inhomogeneous. (Pseudo)rapidity spectra of all charged hadrons for Au-Au collisions are available from the PHOBOS analysis (34) for , 62.4 and 200 GeV energies only. The spectra are published for  GeV by the PHOBOS collaboration (35) and for  GeV by the PHENIX collaboration (36). To cover the lower collision energies we use and -spectra from the NA49 (29) collaboration for Pb-Pb collisions at and  A GeV, which correspond to and 17.6 GeV, and set up the simulations accordingly for Pb-Pb system. For the elliptic flow we compare to the STAR results at , 11.5, 19.6, 27, 39 GeV (10) and 200 GeV (37) collision energies. In the model we define the centrality classes as impact parameter intervals based on the optical Glauber model estimates (38); (39).

The transverse momentum distributions of identified particles at GeV ( A GeV) collision, and (pseudo)rapidity distributions of identified or charged particles at collision energies –200 GeV are shown in Figs. 3 and 4, respectively. As can be seen, the inclusion of shear viscosity in the hydrodynamic phase hardens the spectra, and increases (and similarly ) at midrapidity, squeezing the overall rapidity distribution. This effect can be attributed to the effect of shear viscosity on the strong longitudinal expansion of the system in the initial state for the hydrodynamic phase. Shear viscosity attempts to isotropize the expansion by decelerating it in the longitudinal direction and accelerating it in the transverse direction. The energy of the hydrodynamic system is always conserved, whereas additional entropy is produced during the viscous hydrodynamic evolution, which explains the increased total particle multiplicity. Comparing to the experimental data one observes that gives a good estimate of the rapidity and transverse momentum distributions at the lowest collision energy point GeV ( A GeV), while overestimating at midrapidity for the rest of collision energies except for the highest energy,  GeV, where we underestimate the PHOBOS results.

In Fig. 5 the -averaged elliptic and triangular flow coefficients and are shown as a function of collision energy. The flow coefficients are calculated using the event-plane method as described in Ref. (33), including the event plane resolution correction. As expected, the elliptic and triangular flow coefficients are suppressed by the shear viscosity. However, comparing the results for to the STAR experimental results at centrality we find that the suppression is too weak for  GeV and too strong otherwise. The latter is consistent with the fact that the optimal value of required to fit the elliptic flow data at  A GeV is assuming the initial energy density profile from Monte Carlo Glauber approach (41). Another particular feature of both curves is that, in the region –62 GeV the elliptic flow decreases as a function of . If we do not limit the initial time from below at energies GeV, but take it directly from Eq. (1), we do not see this decrease, but increases monotonously with increasing collision energy. Thus we expect that the reason for the nonmonotonous behavior is in our choice for the initial time of the hydrodynamic evolution.

The results from the standard UrQMD cascade (without intermediate hydrodynamic phase) are also shown for comparison on Figs. 3 and 4 with dotted lines. One may conclude that, whereas standard UrQMD does a good job for -spectra and rapidity distributions at the lowest energy, it clearly underestimates when the collision energy increases (which repeats the conclusion about the excitation function from Ref. (42), and later results from analysis in Ref. (14)). This is an indication of too large viscosity of the high-density medium and served historically as a motivation to introduce the intermediate hydrodynamic stage.

## Iv Investigation of Parameter space

After investigating the generic influence of a finite shear viscosity during the hydrodynamic evolution on basic bulk observables, it is clear that we cannot fit all the available experimental data using the same set of parameters4. Thus we have to adjust the model parameters according to the collision energy before drawing any conclusions about the physical properties of the system.

In this section we study systematically the sensitivities of the particle yield at midrapidity, which is a measure for the final entropy, the effective slope parameter that measures the strength of the transverse expansion, and the anisotropic flow to the main parameters of the model. Due to the limited space, and to emphasize the main features of the dependencies, we restrict ourselves to one collision energy,  GeV, in the middle of the investigated range. Since the influence of shear viscosity was discussed above, we now concentrate on the remaining parameters of the model: the two Gaussian radii and for the initial distribution of energy, momentum and charges, the starting time for the hydro phase , and the energy density when the switch to the hadronic cascade happens. The default case is  fm,  fm/c (calculated according to Eq. 1), (for simplicity) and  GeV/fm. The dependencies are presented in Figs. 6 and 7, where each curve corresponds to the variation of only one of the parameters, while keeping the default values for the others. All values are normalized to their default values to allow a direct comparison to each other. The effective temperatures of the hadron spectra in the lower panel of Fig. 6 are defined as the parameter of the exponential fit:

 dNmTdmTdy=Cexp(−mTTeff),

where the range is [0.2–1] GeV for pions and protons and [0.05–1] GeV for kaons5. In general we do observe only a very weak dependence on the parameters, that is less than 10% for a 10% change in parameters. The observed dependencies can be summarized as:

• Increased smoothens the initial energy density profile in the transverse plane, which leads to smaller gradients and less explosive transverse expansion. The latter leads to a decrease of the effective temperature (inverse slope) of the -spectra, see Fig. 6, lower panel. Larger also results in decreased ellipticity and triangularity of an initial energy density profile, which is hydrodynamically translated into smaller final elliptic (, see Fig. 7) and triangular () flow components.

• In a similar manner, the increase of leads to shallower longitudinal gradients and weaker longitudinal expansion. Thus more energy remain at midrapidity to form stronger transverse expansion, which increases and . On the other hand, larger also results in larger initial entropy of the system, which considerably increases the final particle multiplicity, see Fig. 6, upper panel.

• Increased has two effects:

1. It leads to a shorter lifetime of the hydrodynamic phase, as a result of longer pre-thermal phase.

2. At the same time enters the Gaussian energy/momentum smearing profile. Thus its increase acts opposite to the increase of .

From the response of the observables to the increase of we find that the second effect is stronger.

• Increased shortens the effective lifetime of the hydrodynamic phase. The shorter time to develop radial and elliptic flows is not fully compensated by the longer cascade phase, which results in the decrease of both final and final . Since the total entropy is conserved in the ideal hydrodynamic expansion, but increases in the cascade stage, the final particle multiplicity increases with the increase of .

The observed dependencies are schematically depicted in Table 1, where the signs of the responses of the observables to the increase of a particular model parameter are shown. As for the magnitudes of the response, one can see from the plots that by varying the parameters of the initialization procedure one has a nearly linear influence on the final , and . From Fig. 7 one can see that by choosing a larger value of it is possible to approach the experimental value of with zero shear viscosity in the hydrodynamic phase. However, such value is inconsistent with the spectra and charged particle multiplicity.

Investigating all the dependencies in detail allows us to choose parameter values which lead to a reasonable reproduction of the data. These values are shown in Table 2. For reasons of simplicity we keep  GeV/fm for all collision energies, since the other parameters provide enough freedom for adjustment. Note that since the model requires a lot of CPU time to obtain results for each particular collision energy and centrality, it is at the moment impractical to provide -optimized values of the model parameters and their errors. Thus the parameters are adjusted manually based on a visual correspondence to the data. A full fledged fit to the data is planned for the future using a model emulator, as suggested in Refs. (43); (44); (45).

## V Results for Bulk Observables

Finally, let us have a look at the results for bulk observables with the energy dependent parameters for the hydrodynamic description (see Table 2).

The (pseudo)rapidity spectra are presented in Fig. 8. One can see that whereas the parameters were adjusted to reproduce the total multiplicities, the resulting shapes of the pseudorapidity distributions are also in a reasonable agreement with the data. From the model results one can observe the change in shape from the single peak structure at  GeV to a doubly-peaked distribution (or from a Dromedary to a Bactrian camel shape) which starts to form at  GeV. At higher collision energies we observe a shallow dip around zero pseudorapidity.

The spectra of pions, kaons and protons in collisions at , 17.6, and 8.8 GeV energies are shown in Fig. 9. In general the spectra and especially the slopes are reproduced, which indicates that both the collective radial flow (generated in the hydrodynamic and cascade stages), and thermal motion are combined in the right proportion.

The elliptic and triangular flow coefficients for 20-30% central Au-Au collisions as a function of collision energy are presented in Fig. 10. As expected, the calculated values of the elliptic flow follow the data closely, since this quantity was used to fix the parameters. In contrast to that, triangular flow is calculated from the same simulated events, and thus can be considered as a prediction of the model. We expect that the non-monotonous behavior of is an artifact of our fitting procedure, and more careful adjustment of the model parameters would further smoothen the behavior of .

The 20-30% centrality class was chosen because the elliptic flow signal is strongest around this centrality class. Also, at this centrality nonflow contributions from minijets, which are not included in the model, are small. The centrality dependence of elliptic flow at  GeV is shown in Fig. 11. The parameters are the same at all centralities. In peripheral collisions the model significantly undershoots the data. This is due to the smoothening procedure used to convert individual particles to the fluid-dynamical initial state. With the present smearing parameters the eccentricity of the system is too small in peripheral collisions, where the size of the entire system is comparable to the smearing radius.

The most important conclusion from the adjustment procedure is that reproduction of the data requires an effective which decreases as a function of increasing collision energy, see Table 2 and Fig. 12. On Fig. 12 one can also see an estimated error band around the optimal values of . As mentioned, a proper determination of the error bars would require a fit. Currently the error band is estimated from the variations of two parameters of the model ( and ) which result in the same value of integrated elliptic flow and a 5% variation in the slope of proton spectrum, which is the most sensitive to a change in radial flow.

In the present calculations is taken to be constant during the evolution of the system, and its value changes only with the collision energy. However we expect that physical depends on both the temperature and baryon chemical potential, and that has a minimum around and zero (47); (48); (49); (50). The smaller the collision energy, the larger the average baryon chemical potential in the system. This indicates that the physical value of should increase with increasing .

In Ref. (51) it was argued that is not an appropriate measure of the fluidity of the system. However, the measure of fluidity proposed in that paper, , where is the total particle number density, enthalpy, and the speed of sound, is difficult to implement in the present fluid-dynamical calculation since is not well defined in our two-phase EoS. Instead, we use as an alternative measure of fluidity the combination , where are the charge densities (baryon, strange, electric) and the corresponding chemical potentials, and which approaches in the limit of small charge densities. We have performed an additional round of simulations, keeping and at all collision energies to see whether different measure of fluidity makes any difference. The resulting elliptic and triangular flow coefficients are shown in Fig. 13. One can see that at all considered collision energies there is no visible difference in the elliptic flow coefficient between the and cases. We have also checked that the two scenarios result in virtually same spectra and distributions. This indicates that the contribution from baryon/electric charge density to the entropy density does not induce strong enough baryon density dependence of the ratio to affect the hydrodynamic evolution.

## Vi Summary and Outlook

A hybrid model featuring a 3+1-dimensional viscous hydrodynamic phase with an explicit treatment of finite baryon and charge densities is introduced. The model employs a chiral model equation of state for the hydrodynamic stage. The initial and late non-equilibrium stages are modeled using the UrQMD hadron cascade on an event-by-event basis.

This hybrid model was applied to describe the dynamics of relativistic heavy ion collisions at energies ranging from the lowest RHIC beam energy scan energy to full RHIC energy,  GeV. After tuning the parameters, it was possible to reproduce the observed pseudorapidity and transverse momentum distributions of produced hadrons and their elliptic flow coefficients. The reproduction of the data requires a finite shear viscosity over entropy density ratio which depends on collision energy. This ratio was found to decrease from to as collision energy increases from to  GeV, and to stay at for  GeV. Since the average baryochemical potential at midrapidity decreases with increasing collision energy, the required collision energy dependence of the effective indicates that the physical -ratio may depend on baryochemical potential, and that increases with increasing . It was also found that a constant and collision energy independent and in hydrodynamic phase yield quantitatively similar results. This indicates that the term in entropy density does not induce the baryon density dependence of required to reproduce the data when is kept independent of collision energy.

In addition we have explored the parameter dependence of the model results and generally found a % variation of the results, when the individual parameters were varied by 10%. Of course, the proper evaluation of the effect of finite baryochemical potential on would require reproducing all the data using the same temperature and baryochemical potential dependent parametrization of at all energies and centralities. This will be addressed in future studies.

###### Acknowledgements.
The authors acknowledge the financial support by the Helmholtz International Center for FAIR and Hessian LOEWE initiative. The work of P.H. was supported by BMBF under contract no. 06FY9092. H.P. acknowledges funding by the Helmholtz Young Investigator Group VH-NG-822 from the Helmholtz Association and GSI. Computational resources have been provided by the Center for Scientific Computing (CSC) at the Goethe University Frankfurt.

### Footnotes

1. Typical grid spacing used in the calculations:  fm, and timestep  fm/c depending on the collision energy. A finer grid with  fm was taken to simulate peripheral collisions.
2. The exact procedure suggested in Ref. (28) requires a numerical solution of a cubic equation for each surface element, and is therefore too slow for event-by-event studies.
3. For a suggested procedure to impose the conservation laws, see Ref. (26).
4. The internal parameters of UrQMD, e.g., particle properties and cross sections, are fixed by experimental data as explained in Ref. (15). The effects of changes in resonance properties were studied in Ref. (16). It was found that if the changes stay within experimentally acceptable limits, the effects on final particle distributions are small.
5. Smaller range for pions and protons is taken since the lowest part of the spectrum has a different slope than the intermediate range.

### References

1. T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Lett. B 636, 299 (2006).
2. C. Nonaka and S. A. Bass, Nucl. Phys. A 774, 873 (2006).
3. H. Petersen, J. Steinheimer, G. Burau, M. Bleicher and H. Stocker, Phys. Rev. C 78, 044901 (2008).
4. K. Werner, I. Karpenko, T. Pierog, M. Bleicher and K. Mikhailov, Phys. Rev. C 82, 044904 (2010).
5. H. Song, S. A. Bass and U. Heinz, Phys. Rev. C 83, 024912 (2011).
6. G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 87, 081601 (2001); P. K. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005).
7. C. Gale, S. Jeon, B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 110, no. 1, 012302 (2013).
8. M. Gyulassy and L. McLerran, Nucl. Phys. A 750, 30 (2005).
9. M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999); Phys. Rev. Lett. 81, 4816 (1998).
10. L. Adamczyk et al. [STAR Collaboration], Phys. Rev. C 86, 054908 (2012).
11. G. S. Denicol, C. Gale, S. Jeon and J. Noronha, Phys. Rev. C 88, no. 6, 064901 (2013).
12. S. A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland and M. Hofmann et al., Prog. Part. Nucl. Phys. 41, 255 (1998); M. Bleicher et al., J. Phys. G G 25 1859 (1999).
13. I. Karpenko, P. Huovinen and M. Bleicher, Comput. Phys. Commun. 185, 3016 (2014).
14. J. Auvinen and H. Petersen, Phys. Rev. C 88, no. 6, 064908 (2013).
15. H. Petersen, M. Bleicher, S. A. Bass and H. Stocker, arXiv:0805.0567 [hep-ph].
16. J. Gerhard, B. Bauchle, V. Lindenstruth and M. Bleicher, Phys. Rev. C 85, 044912 (2012).
17. I. A. Karpenko, M. Bleicher, P. Huovinen and H. Petersen, J. Phys. Conf. Ser. 509, 012067 (2014).
18. I. A. Karpenko, M. Bleicher, P. Huovinen and H. Petersen, J. Phys. Conf. Ser. 503, 012040 (2014).
19. W. Israel, Annals Phys. 100, 310 (1976); W. Israel and J. M. Stewart, Annals Phys. 118, 341 (1979).
20. G. S. Denicol, H. Niemi, E. Molnar and D. H. Rischke, Phys. Rev. D 85, 114047 (2012).
21. E. Molnar, H. Holopainen, P. Huovinen and H. Niemi, Phys. Rev. C 90, no. 4, 044904 (2014).
22. J. Steinheimer, S. Schramm and H. Stocker, J. Phys. G 38, 035001 (2011).
23. S. S. Gubser, Phys. Rev. D 82, 085027 (2010).
24. H. Song and U. W. Heinz, Phys. Rev. C 77, 064901 (2008). Phys. Rev. C 78, 024902 (2008).
25. F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
26. P. Huovinen and H. Petersen, Eur. Phys. J. A 48, 171 (2012).
27. F. Becattini, J. Manninen and M. Gazdzicki, Phys. Rev. C 73, 044905 (2006).
28. Y. Cheng, L. P. Csernai, V. K. Magas, B. R. Schlei and D. Strottman, Phys. Rev. C 81, 064910 (2010).
29. S. V. Afanasiev et al. [NA49 Collaboration], Phys. Rev. C 66, 054902 (2002).
30. A. Muronga, Phys. Rev. Lett. 88, 062302 (2002) [Erratum-ibid. 89, 159901 (2002)].
31. R. Baier and P. Romatschke, Eur. Phys. J. C 51, 677 (2007).
32. D. Teaney, Phys. Rev. C 68, 034913 (2003).
33. H. Holopainen, H. Niemi and K. J. Eskola, Phys. Rev. C 83, 034901 (2011).
34. B. Alver et al. [PHOBOS Collaboration], Phys. Rev. C 83, 024913 (2011).
35. B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. C 75, 024910 (2007).
36. S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. C 69, 034909 (2004).
37. J. Adams et al. [STAR Collaboration], Phys. Rev. C 72, 014904 (2005).
38. K. J. Eskola, K. Kajantie and J. Lindfors, Nucl. Phys. B 323, 37 (1989).
39. D. Miskowiec. http://web-docs.gsi.de/â¼misko/overlap/
40. L. Adamczyk et al. [STAR Collaboration], Phys. Rev. C 88, no. 1, 014904 (2013).
41. H. Song, S. A. Bass, U. Heinz, T. Hirano and C. Shen, Phys. Rev. Lett. 106, 192301 (2011) [Erratum-ibid. 109, 139904 (2012)].
42. H. Petersen and M. Bleicher, Phys. Rev. C 79, 054904 (2009).
43. H. Petersen, C. Coleman-Smith, S. A. Bass and R. Wolpert, J. Phys. G 38, 045102 (2011).
44. J. Novak, K. Novak, S. Pratt, J. Vredevoogd, C. E. Coleman-Smith and R. L. Wolpert, Phys. Rev. C 89, 034917 (2014).
45. J. E. Bernhard, P. W. Marcy, C. E. Coleman-Smith, S. Huzurbazar, R. L. Wolpert and S. A. Bass, arXiv:1502.00339 [nucl-th].
46. T. Anticic et al. [NA49 Collaboration], Phys. Rev. C 83, 014901 (2011).
47. L. P. Csernai, J. I. Kapusta and L. D. McLerran, Phys. Rev. Lett. 97, 152303 (2006).
48. G. S. Denicol, T. Kodama and T. Koide, J. Phys. G 37, 094040 (2010).
49. H. Niemi, G. S. Denicol, P. Huovinen, E. Molnar and D. H. Rischke, Phys. Rev. Lett. 106, 212302 (2011).
50. H. Niemi, G. S. Denicol, P. Huovinen, E. Molnar and D. H. Rischke, Phys. Rev. C 86, 014909 (2012).
51. J. Liao and V. Koch, Phys. Rev. C 81, 014902 (2010).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters