Abstract

Viscous hydrodynamical modeling of relativistic heavy ion collisions has been highly successful in explaining bulk of the experimental data in RHIC and LHC energy collisions. We briefly review viscous hydrodynamics modeling of high energy nuclear collisions. Basic ingredients of the modeling, the hydrodynamic equations, relaxation equations for dissipative forces, are discussed. Hydrodynamical modeling being a boundary value problem, we discuss the initial conditions, freeze-out process. We also show representative simulation results in comparison with experimental data. We also discuss the recent developments in event-by-event hydrodynamics.

1. Introduction

Lattice QCD simulations of strongly interacting nuclear matter suggests that, under appropriate conditions (e.g., high temperature, high density), the normally confining nuclear matter can exist in a deconfined state, commonly called Quark-Gluon Plasma (QGP). Lattice simulations also predict that the confinement-deconfinement transition is not a phase transition rather a crossover transition [1–6], with the crossover temperature (or pseudocritical temperature) –170 MeV. It also appears that the deconfined matter is not entirely interaction free even at very high temperature. Thermodynamic quantities, for example, energy density, pressure, and entropy density, are significantly below the interaction free values. It is thus appropriate to call the deconfined matter strongly interacting Quark-Gluon Plasma (sQGP). Recent experiments in  GeV Au+Au collisions at Relativistic Heavy Ion Collider (RHIC) [7–10] and  TeV Pb+Pb collisions at Large Hadron Collider (LHC) [11–15] produced convincing evidences that in RHIC and LHC collisions a deconfined medium or Quark-Gluon Plasma is produced.

A nucleus-nucleus collision at relativistic energy passes through different stages. Schematic picture of different stages of the collision are shown in Figure 1. One can broadly classify the following stages.(i)Preequilibrium stage: at relativistic energy, initial collisions are expected to be in the partonic level. Initial partonic collisions produce a fireball in a highly excited state. In all possibility, the fireball is not in equilibrium. Constituents of the fireball collide frequently to establish a “local” equilibrium’ state. The time taken to establish the local equilibrium is called the “thermalisation” or “equilibration” time.(ii)Expansion stage: in the equilibrium or the thermalised state, the constituents of the fireball, the partons or more specifically, the quarks and gluons are in the deconfined state. The system has thermal pressure which acts against the surrounding vacuum. The fireball then undergoes collective (hydrodynamic) expansion. As it expands, density (energy density) decreases, and the system cools. If the system supports phase transition, then below a critical temperature, the deconfined quarks and gluons will hadronise. In the hadronisation stage, over a small temperature interval, entropy density will decrease very fast. Since total entropy cannot decrease, it implies that the fire ball will expand rapidly, while temperature remains approximately constant.(iii)Freeze-out: even after the hadronisation, the matter can also be in thermal equilibrium. Constituent hadrons will collide to maintain local equilibrium. The system will expand and cool. A stage will come when inelastic collisions, in which hadron changes identity, become too small to keep up with expansion. The stage is called chemical freeze-out. Hadron abundances will remain fixed after the chemical freeze-out. However, due to elastic collisions, local equilibrium can still be maintained, and the system will further expands and cool, with fixed hadron abundances. Eventually a stage will come when average distance between the constituents will be larger than strong interaction range. Collisions between the constituents will be so infrequent that “local” thermal equilibrium cannot be maintained. The hydrodynamic description will break down. The hadrons decouple or freeze-out. It is called kinetic freeze-out. Hadrons from the kinetic freeze out surface will be detected in the detector.

Figure 1 

Different stages of nuclear collisions at ultrarelativistic energy.

From the expansion stage to freeze-out, the collision process can be modeled by relativistic hydrodynamic equations. If the macroscopic properties of the fluid, for example, local energy density, pressure, fluid velocity, and so forth, are known at the initial time , hydrodynamic equations can be solved to obtain the space-time evolution of the fireball till the freeze-out. At the freeze-out suitable algorithm (e.g., Cooper-Frye [16]) can be used to convert the fluid information at the freeze-out into particle’s invariant yield and compare with experimental data. Hydrodynamic models are unique that it is imperative to use an equation of state (EoS), a thermodynamic relation between the energy density, pressure, and number density of the fluid. By explicitly incorporating phase transition in the equation of state, one can study, dynamically, the effect of phase transition on the fluid evolution and associated particle production. Ideal hydrodynamic models have been largely successful in explaining a variety of experimental data, for example, transverse momentum spectra, elliptic flow of identified particles in  GeV Au+Au collisions [17]. Indeed, success of ideal fluid dynamics in explaining several experimental data, together with the string theory motivated lower limit of shear viscosity over entropy ratio [18–20], leads to a paradigm that in Au+Au collisions, a nearly “perfect” fluid, is created.

However, the paradigm of “perfect fluid” produced in Au+Au collisions at RHIC needs to be clarified. It so happens that the ideal fluid dynamic models do have their limitations [21]. For example, experimentally, elliptic flow tends to saturate at large transverse momentum. The ideal fluid dynamics on the other hand predicts a continually increasing elliptic flow. The transverse momentum spectra of identified particles also start to deviate from ideal fluid dynamics prediction beyond  GeV. Experimentally determined HBT radii are not reproduced in the ideal fluid dynamic models, the famous “HBT puzzle” [22]. Ideal fluid dynamics also works best in central collisions and gets poorer in more peripheral collisions. The shortcomings of ideal fluid dynamics possibly indicate greater importance of dissipative effects in the ranges greater than 1.5 GeV or in more peripheral collisions. Indeed, ideal fluid is a concept, which is never realized in nature. As suggested in string theory motivated models, QGP viscosity could be small, , nevertheless it is nonzero. Indeed, much earlier, Danielewicz and Gyulassy [23], from the uncertainty principle, estimated the lower bound of viscosity to entropy ratio, , a value very close to the AdS/CFT estimate. It is thus important to study the effect of viscosity, even if small, on space-time evolution of QGP fluid and quantify its effect. Furthermore, QGP fluid has to be characterized by its transport coefficients, for example, heat conductivity, bulk and shear viscosity. Theoretically, it is possible to obtain those transport coefficients in a kinetic theory model. However, in the present status of theory, the goal cannot be achieved immediately, even more so for a strongly interacting QGP (sQGP). Alternatively, one can compare viscous hydrodynamic simulations to experimental data and obtain a “phenomenological” limit to the transport coefficients of sQGP. There is another incentive to study dissipative hydrodynamics. Ideal hydrodynamics depends on the assumption of local equilibrium. In dissipative hydrodynamics, the strict assumption of local thermal equilibrium is relaxed to the assumption of “near” local thermal equilibrium, extending the range of validity of hydrodynamic description. Indeed, one can explore early times of fluid evolution better in dissipative hydrodynamics. It may be mentioned here that success of hydrodynamical modeling of relativistic heavy ion collisions does not necessarily imply the realisation of local thermodynamic equilibrium, although the inverse is true. Whether or not local equilibrium is achieved in high energy nuclear collisions is still in debate.

Theory of dissipative relativistic fluid has been formulated quite early. The original dissipative relativistic fluid equations were given by Eckart [24] and Landau and Lifshitz [25]. They are called first order theories. Formally, relativistic dissipative hydrodynamic equations are obtained from an expansion of entropy 4-current, in terms of the dissipative fluxes. In first order theories, entropy 4-current contains terms linear in dissipative quantities. First order theory of dissipative hydrodynamics suffers from the problem of causality violation and instabilities [26, 27]. Causality violation is unwarranted in any theory, even more in a relativistic theory. The problem of causality violation is removed in the Israel-Stewart’s second order theory of dissipative fluid [28, 29]. In second order theory, expansion of entropy 4-current contains terms second order in dissipative fluxes. However, these leads to complications that dissipative fluxes are no longer function of the state variables only. They become dynamic. The space of thermodynamic variables has to be extended to include the dissipative fluxes (e.g., heat conductivity, bulk, and shear viscosity). In the following section Israel-Stewarts phenomenological theory of dissipative hydrodynamics is briefly discussed. More detailed exposition can be found in [29, 30].

2. Dissipative Fluid Dynamics

A simple fluid, in an arbitrary state, is fully specified by primary variables: particle current , energy-momentum tensor , entropy current , and a number of additional (unknown) variables. Primary variables satisfy the conservation laws: and the second law of thermodynamics,

In relativistic fluid dynamics, one defines a time-like hydrodynamic 4-velocity, (normalized as ). One also defines a projector, , orthogonal to the 4-velocity . In equilibrium, an unique 4-velocity exists such that the particle density , energy density , and the entropy density can be obtained from

An equilibrium state is assumed to be fully specified by 5 parameters or equivalently by the thermal potential ( being the chemical potential) and inverse 4-temperature, . Given an equation of state, , pressure can be obtained from the generalized thermodynamic relation:

Using the Gibbs-Duhem relation, , following relations can be established on the equilibrium hypersurface ,

In a nonequilibrium system, no 4-velocity can be found such that (3) remains valid. Tensor decomposition leads to additional terms:

The new terms describe a net flow of charge , energy flow (where is the heat flow vector), and entropy flow . is the bulk viscous pressure, and is the shear stress tensor. Hydrodynamic 4-velocity can be chosen to eliminate either (the Eckart frame, is parallel to particle flow) or the energy flow (the Landau frame, is eigenvector of energy-momentum tensor). In relativistic heavy ion collisions, central rapidity region is nearly baryon free, and Landau’s frame is more appropriate than the Eckart’s frame. Dissipative flows are transverse to , and additionally, shear stress tensor is traceless. Thus a nonequilibrium state requires 1  +  3  +  5  =  9 additional quantities, the dissipative flows , (or ), and . In kinetic theory, and are the first and second moments of the distribution function. Unless the function is known a-priori, two moments do not furnish enough information to enumerate the microscopic states required to determine , and in an arbitrary nonequilibrium state, no relation exists between , , and . Only in a state, close to an equilibrium one, such a relation can be established. Assuming that the equilibrium relation (5) remains valid in a “near equilibrium state,” also the entropy current can be generalized as where is an undetermined quantity in second order in deviations, , and . Detail form of is constrained by the second law . With the help of conservation laws and Gibbs-Duhem relation, entropy production rate can be written as

Choice of leads to first order or second order theories of dissipative hydrodynamics. In first order theories the simplest choice is made, , and entropy current contains terms up to first order in deviations, and . Entropy production rate can be written as where , , and .

The second law, , can be satisfied by postulating a linear relation between the dissipative flows and thermodynamic forces: where , , and are the positive transport coefficients, bulk viscosity, heat conductivity, and shear viscosity, respectively.

In first order theories, causality is violated [26, 27]. Causality violation is corrected in second order theories [29]. In second order theories, entropy current contain terms up to second order in the deviations, . The most general containing terms up to second order in deviations can be written as where , , and are thermodynamic coefficients for the bulk stress , heat flow , and shear viscous stress , respectively. and are thermodynamic coefficients for the coupling between the heat flow and bulk and shear stress, respectively. As before, one can cast the entropy production rate in the form of (9). Neglecting the terms involving dissipative flows with gradients of equilibrium thermodynamic quantities (both are assumed to be small) and demanding that a linear relation exists between the dissipative flows and thermodynamic forces, the following relaxation equations for the dissipative flows can be obtained:

The relaxation times are

s are coupling coefficients between different dissipative flows, specifically,

Unlike in the first order theories, in second order theories, dynamical equations control the dissipative flows. Even if thermodynamic forces vanish, dissipative flows do not vanish instantly. It is important to mention that the parameters, and , are not connected to the actual state . The pressure in (7) is also not the “actual” thermodynamics pressure, that is, not the work done in an isentropic expansion. Chemical potential and 4-inverse temperature have meaning only for the equilibrium state. Their meaning need not to be extended to nonequilibrium states also. However, it is possible to fit a fictitious “local equilibrium” state, point by point, such that pressure in (7) can be identified with the thermodynamic pressure, at least up to first order. The conditions of fit fix the underlying nonequilibrium phase-space distribution.

It may be mentioned here that relaxation equations for the dissipative fluxes can also be derived in kinetic theory [29, 31–34]. In kinetic theory, relaxation equations as well as explicit expressions for transport coefficients, relaxation times, and so forth can be obtained. There is no unique method to obtain the relaxation equations from kinetic theory. For example, Israel-Stewart obtained the relaxation equations from the second moment of kinetic equation. Betz et al. [31] on the other hand obtained them from the kinetic equation itself. Both the methods gave identical relaxation equations. However, while the relaxation equations remain unchanged, the relaxation time and coupling coefficients do depend on whether the zeroth moment or the second moment of the kinetic equation is used to derive the relaxation equation.

Among the three dissipative coefficients, shear viscosity appears to be most important in heavy ion collisions. In a collision, initial momentum is predominantly longitudinal. Some shear force must act on it to isotropise the momentum distribution (as required by the assumption of thermal equilibrium). Most of the dissipative hydrodynamic studies for heavy ion collisions thus concern with the effect of shear viscosity on fluid evolution and subsequent particle emission. Bulk viscosity, in general is order of magnitude less than shear viscosity. However, there are indications that, in QGP, near the transition point, bulk viscosity can be large [35, 36]. Recently, effect of bulk viscosity on particle production is being investigated [37–42]. There is some uncertainty about the correct form of the nonequilibrium correction to the equilibrium distribution function in presence of bulk viscosity. Grad’s 14-moment method for nonequilibrium correction appears to give very large correction even for small bulk viscosity [42]. In [41], a different conclusion was reached. Using Grad’s 14-moment method in orthogonal basis form [32] along with bulk viscosity relaxation time of [34] does not lead to such problems. Conductivity of QGP fluid is the least studied dissipative coefficient. The central rapidity region in ultrarelativistic heavy ion collisions is essentially baryon free, and effect of conductivity is minimum. However, in future FAIR energy collisions, effect of conductivity may be important.

3. Equation of State

Hydrodynamic equations are closed with an equation of state (EOS) . As it was mentioned earlier, hydrodynamic models are unique that the effect of phase transition can be studied, dynamically, by including the phase transition in the equation of state. Earlier, hydrodynamical simulations of relativistic heavy ion collisions used EOS with first-order confinement-deconfinement phase transition. The simple Bag model equation of state was used extensively to model EOS of the QGP phase. The confined phase is generally modeled as interaction free hadron resonance gas. At sufficiently low temperature, thermodynamics of a strongly interacting matter is dominated by pions. As the temperature increases, larger and larger fraction of available energy goes into excitation of more and more heavier resonances. For temperature  MeV, heavy states dominate the energy density. However, densities of heavy particles are still small, . There mutual interaction, being proportional to , is suppressed. One can use Virial expansion to obtain an effective interaction. Virial expansion together with experimental phase shifts was used by Venugopalan and Prakash to study thermodynamics of low temperature hadronic matter [43]. It was shown that interplay of attractive interactions (characterised by positive phase shifts) and repulsive interactions (characterised by negative phase shifts) is such that, effectively, theory is interaction free. One can then consider that interaction-free resonances constitute the hadronic matter at low temperature.

In recent years, there is much progress in lattice simulation of QCD. In particular, thermodynamic properties of baryon free QCD have been studied accurately and extensively. Currently, there is consensus that the confinement-deconfinement transition is a crossover, and the crossover or the pseudocritical temperature for the transition is  MeV [3–6]. Accurate lattice simulations also show that at low temperature, interaction-free hadron resonance gas correctly reproduces lattice simulation results. Recent hydrodynamical simulations generally uses EOS where the lattice simulations results for the deconfined phase are smoothly joined at MeV, with hadronic resonance gas EOS [44, 45]. In Figure 2, Wuppertal-Budapest simulations for entropy density and pressure are shown. The lines in Figure 2 are for a lattice-based EOS [45], where the lattice simulations for entropy density are smoothly joined at MeV, with the entropy density of the hadronic resonance gas result. The energy density and pressure are then obtained from the thermodynamical relation:

Figure 2 

Lattice simulations [6] for energy density and entropy density are shown. The lines are parameterised EOS. For details of the parametric EOS, see [45].

4. Initial Condition

One understands that hydrodynamics is an initial value problem. To solve the conservation one has to initialise the (baryon) number density , energy density , and velocity distributions at the initial time , beyond which hydrodynamics is applicable. In viscous hydrodynamics, the viscous stresses have also to be initialised. In relativistic heavy ion collisions, rather than , more appropriate coordinates are . In coordinate system, for the initial energy density, a common practice is to assume a factorised form: where is the initial energy density in the transverse plane and is in the direction of (spatial) rapidity . One can use a Gaussian distribution for . Transverse energy distribution can be conveniently parameterised in a Glauber model [17] or in color glass condensate (CGC) model [46–48]. I briefly discuss the Glauber and CGC initial conditions. In Glauber model, in A+B collisions at impact parameter , at the initial time , the initial entropy density () or energy density () is parameterised as or where and are the transverse profile of participant nucleons and binary collision number and can easily be computed in Glauber model: where the thickness function is , with being the density distribution. in (20), (21) is the hard scattering fraction (to be fitted to experimental data).

CGC model initial condition relies on the physical argument that gluon production saturates at high energies [49–51], and prior to QGP, a new form of matter, “Color Glass Condensate (CGC)” is formed. The new form of matter then evolves into QGP. Evolution of CGC is governed by complicated nonlinear JIMWLK equation [50–53]. For hydrodynamical application, one generally uses the Kharzeev-Levin-Nardi (KLN) approach (which captures the essential features of the gluon saturation) [54, 55] and computes the transverse profile of the gluons in A+A collisions. Transverse profile of gluons in A+B collision can be computed as where and are the transverse momentum and rapidity of the produced gluon. is the momentum fraction of colliding gluons ladders at c.m. energy is the strong coupling constant at the momentum scale are unintegrated gluon distribution functions in nuclei and . In KLN approach [56], the unintegrated gluon distribution function is taken as [55, 56]: where is saturation momentum at the given momentum fraction and at the transverse position . In the KLN approach, the saturation scale in AB collision is parameterised as [55, 56]:

The form with –0.3 is motivated from DIS experiments. in the above equation is the transverse density of participant nucleons, which can be calculated in a Glauber model.

In the CGC model, the transverse energy density should follow the gluon density distribution (23). However, (23) is valid in the time scale , when the medium may not be in thermal equilibrium. One assumes that the medium undergoes one-dimensional Bjorken (longitudinal, isentropic) expansion during the period to . The density at the time , when hydrodynamic becomes applicable, is easily obtained as . The transverse energy density profile at the initial time is

It may be mentioned here that, in general, in a given impact parameter, hydrodynamical simulation with CGC initial condition generates more elliptic flow than the Glauber model initial condition. The reason is also understood. Elliptic flow is proportional to initial eccentricity of the medium. In CGC model initial condition, initial eccentricity is comparatively larger than that in the Glauber model. Recently, in an important publication, Song et al. [57] resolved the uncertainty in the initial conditions in viscous hydrodynamics. They used a hybrid model, where viscous fluid dynamics for the QGP was coupled with a microscopic transport model (UrQMD) for hadronic freeze-out. Hybrid models, without additional parameters, govern the complex hadron kinetic freeze-out and allow quantitative exploration of the transport properties of the earlier QGP phase using measured final hadron spectra. It was shown that at least in a fixed energy collision, the relation between the eccentricity scaled elliptic flow and the multiplicity density is approximately universal function, only depending on the value of QGP viscosity over entropy ratio, , but not on any details of the model from which and is computed. Results of their simulations are shown in Figure 3 (for details see [57]).

Figure 3 

Hybrid model simulation for the eccentricity-scaled elliptic flow as a function of multiplicity density , for different values of . Eccentricity-scaled elliptic flow does not depend on the details of the model for initial density distribution but only on the QGP viscosity. See [57] for details.

The number density distribution at the initial time can be similarly parameterised. In general, one assumes zero initial fluid velocity at the initial time, though it is possible that fluid has nonzero velocity, especially near the surface of the fluid. The reasoning is simple. Fluid constituents have random velocity. In the interior of the fluid, the random velocities will balance to produce net zero velocity, but near the surface random velocities will not be balanced. In viscous hydrodynamics, one also needs to initialise the viscous stresses. One common practice is to assume that at the initial time , viscous stresses are zero, [58]. In hydrodynamic models with the assumption of boostinvariance, boost-invariant values can also be used [59].

It may be mentioned here that the initial time is also a parameter of the model. It is essentially the thermalisation time. One important finding of hydrodynamical analysis of  GeV Au+Au and  TeV Pb+Pb collisions is that hydrodynamical simulations require small thermalisation time –1.0 fm. Origin of small thermalisation time is not properly understood. In QGP, nonabelian version of Weibel instabilities can grow very fast, isotropizing the medium [60]. The maximum growth rate is , with being the density and the characteristic momentum scale dominating the excitations in nonequilibrium QGP [60]. In high energy collisions, a large number of gluons will be freed in the first moments after the collisions, is large, and instabilities can grow fast. However, it must be mentioned here that Chromo-Weibel instabilities have never been shown to be able to isotropize the system in the time scale 0.5–1 fm/c. In a recent publication it was shown that in nucleus-nucleus collisions at LHC energy, Chromo-Wiebel instabilities can isotropise the system in the time as large as ~5 fm [61]. It is even more ~10 fm in RHIC energy scale.

Another important input to viscous hydrodynamics simulation is the relaxation time for bulk stress and shear stress. In principle, relaxation times and could be calculated from the underlying kinetic theory, which for strongly coupled QCD plasma, is a complex problem. Relaxation times and were calculated in [28, 29], for simple relativistic Boltzmann, Bose, and Fermi gases with mass using Grad 14 moment approximation in relativistic kinetic theory. For a Boltzmann gas, in the nonrelativistic limit , and . In the relativistic limit , and . Note that in the relativistic limit, the mass term appears in the denominator with a quartic power. Two phases in heavy ion collision, that is, QGP and hadronic phases, consist of quasiparticles of different masses; the dependence of this mass with temperature as well as the ambiguity around the crossover makes it difficult to use as given here. In most of the simulation, is either taken as a constant or kept same as . Recently, Jaiswal et al. [34] solved the long-standing ambiguity in the relaxation time for bulk viscosity. They derived causal dissipative hydrodynamic equations by invoking second law of thermodynamics from entropy four-current expressed in terms of distribution function. The bulk relaxation time thus obtained shows critical slowing down near transition temperature and does not lead to cavitation.

5. Shear and Bulk Viscosity Coefficients

In viscous hydrodynamic simulations, one also needs to specify the dissipative coefficients. The shear or bulk viscosity for a strongly interacting system can be calculated by using Kubo formula [65, 66] as given in the following:

However, first principle calculations of viscosity coefficients are very difficult for sQGP. See [67] for a detailed discussion of viscous properties of QCD matter. In Figure 4, I have shown some results from first principle attempts to obtain (shear) viscosity over entropy ratio of sQGP. Applicability of perturbative QCD in sQGP is questionable. Still, in [62, 63], was obtained in pQCD approach. pQCD leads to very large (shear) viscosity over entropy ratio , too large to explain experimental elliptic flow. Lattice QCD calculations [64, 68] require analytic continuation from imaginary to real time and results can be obtained only with large uncertainty. As shown in Figure 4, QGP viscosity over entropy ratio appears to increase with temperature, though very large error bars preclude any definite conclusion. In AdS/CFT approach, one relates the correlator to Graviton absorption cross-section. AdS/CFT correspondence is based on Malcedena’s gauge/gravity duality conjecture. It is conjectured that string theory with gravity defined on a product space is equivalent to quantum field theory without gravity defined on the conformal boundary of ; here is the 5-dimensional anti-deSitter space, and is a five-dimensional sphere. For more details, see [69]. In this approach, energy-momentum tensor couples to gravity at the boundary, and one can calculate the absorption cross-section of a graviton. Absorption cross-section of a graviton of frequency , polarised in direction, is the imaginary part of the retarded Green function: where appears due to normalisation of the graviton action. One can show that absorption cross-section of graviton is equal to that of a minimally coupled scalar, in the low energy limit which can be identified with the area of the horizon [70]. Considering that black hole entropy is , one can easily obtain the AdS/CFT limit for shear viscosity over entropy ratio [18–20]:

Figure 4 

Viscosity over entropy ratio as a function of temperature from perturbative calculation [62, 63], lattice QCD calculation [64]. AdS/CFT bound is also shown.

Bulk viscosity is proportional to interaction measure , and in a purely conformal system it is identical to zero. In QCD conformal symmetry is broken by the running coupling constant, , and can be nonzero. Leading log perturbative limit for bulk viscosity is very small, to within 10% for [71], and is consistent with the old relation, [72], with being the speed of sound in the medium. AdS/CFT calculations that incorporate deviations from conformality at strong coupling find [73, 74]. However, there are indications that in QGP, near the transition point, bulk viscosity can be large [35, 36].

Bulk and shear viscosity of the hadronic matter has been also been calculated, either in effective hadronic interaction model or in the transport model [75–78]. QGP viscosity over entropy ratio of the hadronic gas decreases with increasing temperature, bulk viscosity over entropy ratio and on the other hand shows opposite behavior increase with temperature. Present paradigm is that while shear viscosity over entropy ratio has a minimum at [79, 80], bulk viscosity over entropy ratio has a maxima at .

Considering the limited theoretical knowledge of the viscosity over entropy ratio, an alternative approach is to extract the ratio by comparing viscous fluid dynamical simulations with experimental data. This is the approach taken in most of the hydrodynamical simulations.

6. Freeze-Out

In addition to the initialisation of the fluid, hydrodynamic models also require a freeze-out condition. As it was mentioned earlier, there are two types of freeze-out, in hydrodynamical evolution, the chemical freeze-out and the kinetic freeze-out. Hadron abundances remain fixed after the chemical freeze-out. Statistical models, based on Grand Canonical ensemble, surprisingly give correct values for hadron ratios [81, 82]. Chemical freeze-out temperature and baryonic chemical potentials have been parameterised as a function of collision energy. The parameterised forms can be used to obtain the chemical freeze-out temperature for application in hydrodynamical modeling. The kinetic freeze-out temperature, generally, is treated as a parameter, and its value is obtained from fitting experimental data. Freeze-out temperature –150 MeV has been used in hydrodynamic simulations.

At kinetic freeze-out, Cooper-Frye prescription [16] is used to obtained the particle (invariant) distribution: where is the freeze-out hypersurface and is the one-body distribution function. Now in ideal dynamics, the fluid is in local equilibrium, and the one-body distribution function is well approximated by the equilibrium distribution function: with inverse temperature and chemical potential . is the degeneracy factor. In viscous dynamics on the other hand, the fluid is not in equilibrium, and cannot be approximated by the equilibrium distribution function . In a highly nonequilibrium system, distribution function is unknown. If the system is slightly off-equilibrium, then it is possible to calculate correction to equilibrium distribution function due to (small) nonequilibrium effects. Slightly off-equilibrium distribution function can be approximated as where , −1, and 0 for Fermi, Bose, and Boltzmann gas, respectively. is the deviation from equilibrium distribution function . With shear viscosity as the only dissipative forces, can be locally approximated by a quadratic function of 4-momentum: completely specifying the nonequilibrium distribution functions. As expected, correction factor increases with increasing viscosity. We also note that nonequilibrium correction depends quadratically on particle momentum. The effect of dissipation is more on large momentum particles. For bulk viscosity, the situation is more complicated. The nonequilibrium correction can be written as where , , and are space-time-dependent parameters, explicit expression for which can be found in [38]. It appears that even for small bulk viscosity, the correction factor can be very large, invalidating viscous hydrodynamics at momenta  GeV or so [42].

In the late stage of evolution, applicability of hydrodynamics may be poor. In more realistic models, hydrodynamic models are coupled with hadron-based transport models [83–87]. In the hybrid models, below a switching temperature , motion of hadrons is governed by transport equations, for example, UrQMD model. At the switching temperature hadron abundances can be computed using the Cooper-Frye formalism. The initial state of the transport models is produced through the Monte-Carlo sampling. Agreement with data appears to improve in hybrid models.

7. Collective Flow

In relativistic heavy ion collisions, one of the important observables is the azimuthal distribution of produced particles. In Figure 5, geometry of a collision at nonzero impact parameter is shown. The overlap region of the two nuclei is the participant region, where most of the collisions occur. The target and projectile remnants on the periphery act as spectator. It is obvious from Figure 5 that, in nonzero impact parameter collisions, the participant or the reaction zone in coordinate space does not possess azimuthal symmetry. Multiple collisions among the constituent particles translate this spatial anisotropy into momentum anisotropy of the produced particles. The observed momentum anisotropy is called collective flow and has a natural explanation in a hydrodynamic model [17]. In the following we briefly discuss collective flow phenomena. More detailed expositions can be found in [88–90].

Figure 5 

Schematic picture of a nonzero impact parameter collision.

Momentum anisotropy is best studied by decomposing the invariant distribution in a Fourier series. For example, the momentum integrated invariant distribution of a particle can be expanded as where is the azimuthal angle of the detected particle and is the plane of the symmetry of initial collision zone. For smooth initial matter distribution (obtained from geometric overlap of the colliding nuclei, e.g., Glauber model), plane of symmetry of the collision zone coincides with the reaction plane (the plane containing the impact parameter and the beam axis). The sine terms of the Fourier expansion does not contribute due to symmetry with respect to the reaction plane.

Flow coefficients are easily obtained,

is called (integrated) directed flow, is called (integrated) elliptic flow, is called (integrated) triangular flow, is called (integrated) hexadecapole flow, and so forth.

Similar to (35), one can Fourier expand the invariant distribution

The coefficients are called differential flow coefficient. Second flow coefficient or elliptic flow has been studied extensively in RHIC and LHC energy collisions. Finite-nonzero value of is thought to be direct signature of production of thermalised medium. Elliptic flow is best understood in a hydrodynamic model [17]. Elliptic flow measures the momentum anisotropy. In nonzero impact parameter collisions, the reaction zone is spatially asymmetric (see Figure 5). Spatial asymmetry of the initial reaction zone can be quantified in terms of eccentricity, defined as where indicates energy/entropy density weighted averaging. In nonzero impact parameter collision, initial eccentricity is positive. If a thermalised medium is produced in the reaction zone, due to thermodynamic pressure, the medium will expand against the outside vacuum. One can immediately see that pressure gradient will be more along the minor axis than along the major axis. Due to differential pressure gradient, as the system evolves with time, eccentricity will reduce. Momentum distribution of particles is isotropic initially. If momentum anisotropy is measured as initially will be zero. However, as the fluid evolves, rescattering of particles will introduce asymmetry, and will grow. Beyond certain time, when reaction zone attains azimuthal symmetry, it is expected to saturate. In a sense, elliptic flow is self-quenching phenomena, and driving force of the flow (the reaction zone asymmetry) continuously reduces to quench the flow. In Figure 6, hydrodynamic model simulations for temporal evolution of spatial eccentricity and momentum anisotropy, in ideal and viscous fluid, are shown [96]. As expected, in simulations also, while the eccentricity decreases with time, the momentum anisotropy increases and tends to saturate at large time. Effect of viscosity is not large on the evolution of spatial eccentricity. In viscous fluid, eccentricity decrease marginally faster. The effect is, however, more pronounced on momentum anisotropy. Viscosity hinders the growth of the momentum anisotropy. Saturation value of momentum anisotropy in viscous fluid is lower than that in ideal fluid.

Figure 6 

Ideal and viscous hydrodynamic simulations for evolution of spatial eccentricity and momentum anisotropy in  fm Pb+Pb collisions.

The second harmonic coefficient or the elliptic flow has been studied extensively in  GeV Au+Au collisions at RHIC [9, 10]. Recently, ALICE collaboration measured elliptic flow in  TeV Pb+Pb collisions at LHC [14, 15]. Large elliptic flow has provided compelling evidence that at RHIC and LHC, nearly perfect fluid is produced. Deviation from the ideal fluid behavior is controlled by shear viscosity to entropy ratio . Effect of shear viscosity is to dampen the flow coefficients. Elliptic flow has sensitive dependence on . In smooth hydrodynamics, sensitivity of elliptic flow has been utilised to obtain phenomenological estimates of [47, 85, 97–102]. It appears that QGP viscosity over entropy ratio is close to .

8. Hydrodynamic Simulation of Heavy Ion Collisions and Comparison with Experiment

Even though Israel-Stewart formulated second order theory for dissipative fluid some 30 years back, significant progress towards its numerical implementation has only been made very recently. In Table 1, I have listed the various viscous hydrodynamic models, developed across the world. Several algorithms can be used to solve the energy momentum conservation equations and the relaxation equations. One widely used algorithm is the Smooth and Sharp Transport Algorithm (SHASTA) followed by Flux Corrected Transport (FCT) [103]. SHATSA-FCT algorithm is known to be an accurate and effective algorithm to solve nonlinear generalized continuity equation of the type which occurs in fluid dynamics [103, 104]. It is essentially a three-step process. In the first step, called transport, the velocities of the fluid and the source term are first calculated at half time steps . In the second step, using the new velocity and the new source term in the th step, one calculates the quantities at the th step. The third step is called the antidiffusive step, designed to remove the numerical diffusion inherent to the transport scheme. This is done by calculating an antidiffusive flux, which is subtracted from the time-advanced quantities at the th step to get the final result at the th step. The calculation of the antidiffusion is carried out by a method called “flux correction.” For details of the method, see [103, 104].

In the following, some representative viscous hydrodynamics simulation results are discussed in brief. In Figure 7 viscous hydrodynamic simulation results [47] for the centrality dependence of pions, kaons, and protons multiplicity, in  GeV Au+Au collisions, are compared with experimental data. In the simulations, it was assumed that viscosity over entropy ratio remains a constant throughout the evolution. Three values, (ideal fluid), 0.08 (AdS/CFT limit), and 0.16, were used in the simulation. In [47], both Glauber model initial condition and CGC initial conditions were used to simulate Au+Au collisions. Figure 7 was obtained with Glauber model initial condition. Simulation results indicate that particle’s multiplicity is not sensitive to viscosity when fluid is appropriately initialised (see [47] for details), viscosity over entropy ratio , 0.08 and 0.16 describes the data equally well. In Figure 8, viscous hydrodynamic simulations for charged particle’s integrated elliptic flow are compared with the experimental data. Simulation results for integrated elliptic flow, however, show sensitive dependence on viscosity over entropy ratio, decreasing with increasing viscosity. Similar results are also obtained with CGC initial condition, and particle multiplicity is not sensitive to , but elliptic flow is. We will not discuss in detail but from the analysis of  GeV Au+Au collision data, the viscosity over entropy ratio was estimated as (theory) (experiment) [47].

Figure 7 

Centrality dependence of particle’s multiplicity in a 2  +  1D hydrodynamical simulation is compared with experimental data. The figure is from [47].

Figure 8 

Centrality dependence of charged particles (integrated) elliptic flow compared with experimental data. The figure is from [47].

In Figure 9, viscous hydrodynamic simulation results [48] for transverse momentum spectra in  GeV Au+Au collisions are compared with the PHENIX experiment. Glauber model initial condition was used for these simulations. Here also, viscosity over entropy ratio assumed to be temperature-independent. Effect of viscosity is manifest, and spectra is hardened. More viscous is the fluid, and more is the hardening. From the figures, one conclude that charged particles spectra in 0–50% collision are best described with viscosity to entropy ratio –0.12. In Figure 10, hydrodynamic simulations for elliptic flow are compared with PHENIX data. Elliptic flow is reduced in viscous fluid. In more viscous fluid elliptic flow is reduced even more. One observes that the experimental data on elliptic flow prefers higher values of in more peripheral collisions than in central collisions. For example, elliptic flow in 0–10% collisions is consistent with ideal fluid simulation. Elliptic flow in 20–30% collisions however prefer, , and still higher value, , in 40–50% collision. In [48], PHENIX data were also analysed with CGC initial condition. Conclusions were similar. While charged particles spectra prefer fluid viscosity –0.12, elliptic flow data demand more viscous fluid in more peripheral collision. That the elliptic flow data demand-ependent QGP viscosity was obtained in an earlier analysis [99] also.

(a) 0–10% 

(b) 10–20%

(c) 20–30%

(d) 30–40%

(e) 40–50%

(a) 0–10% 
(b) 10–20%
(c) 20–30%
(d) 30–40%
(e) 40–50%

Figure 9 

PHENIX measurements of charged particles spectra in  GeV Au+Au collisions are compared with viscous hydrodynamic simulations. The figures are from [48].

(a)  0–10%

(b) 10–20%

(c) 20–30%

(d) 30–40%

(e) 40–50%

(a)  0–10%
(b) 10–20%
(c) 20–30%
(d) 30–40%
(e) 40–50%

Figure 10 

Same as in Figure 9 but for elliptic flow. The figures are from [48].

Centrality dependence of is essentially manifest of temperature dependence of viscosity over entropy ratio. As mentioned earlier, simulation results depicted in Figures 7–10 were obtained with temperature-independent . However, viscosity over entropy ratio does have temperature dependence [79, 80]. Viscosity over entropy ratio of hadronic matter increases with lowering of temperature [76]. Lattice QCD simulations do indicate that for QGP increases with temperature [68]. Present paradigm is that the viscosity over entropy ratio for strongly interacting matter has a minimum, possibly with a cusp, around the critical temperature . Demand of higher in peripheral collisions in simulations with constant is not at variance with the prevailing paradigm. Rather it indicates the increasingly important role of hardronic matter in the development of elliptic flow in peripheral collisions. of hadronic matter increase with lowering of temperature.

In recent years, temperature dependence of QGP viscosity over entropy ratio was investigated by several authors [85, 105–108]. In [105], four different temperatures dependence for were considered: (i)   LH-LQ: for all temperatures, (ii) LH-HQ: 0.08 in the hadron gas and above  MeV increases according to lattice QCD data [68], (iii) HH-LQ: below  MeV, is that of a hadron gas [76], and , and (iv) HH-HQ: a realistic parametrization for both the hadron gas and the QGP viscosity [68, 76]. Simulation results for transverse momentum spectra and elliptic flow for pions are shown in Figure 11. One finds that the elliptic flow in  GeV Au+Au collisions at RHIC is dominated by the viscosity in the hadronic phase. It is largely insensitive to QGP viscosity. QGP viscosity dominates only at LHC energy  TeV.

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

Figure 11 

Transverse momentum spectra of positive pions in the 0–5% most central collisions and elliptic flow coefficients in the 20–30% centrality class at RHIC and LHC. Different curves correspond to the different parametrizations of the temperature dependence of described in the text. For details, see [105].

In [108], attempt was made to obtain the temperature dependence of empirically. A linear temperature dependence of for the QGP phase was assumed: and it was smoothly joined with hadronic viscosity [68] at . Experimental data in  GeV Au+Au and  TeV Pb+Pb collisions were confronted with three slope parameters: (i) (ii) and (iii) . Slope parameter –0.2 appears to be consistent with experimental data in  TeV Pb+Pb collisions. LHC data clearly disfavor very strong temperature dependence. At RHIC energy, however, strong temperature dependence was indicated. It was concluded that a unique, linear temperature-dependent QGP viscosity over entropy ratio fails to explain, simultaneously, the experimental elliptic flow data at RHIC and LHC energies. In [85, 107] a similar conclusion that was reached the temperature dependence of QGP viscosity over entropy ratio cannot be constrained by fitting spectra and elliptic flow data alone.

9. Event-by-Event Viscous Hydrodynamics

In the preceding section, hydrodynamical simulations discussed used smooth initial condition, either Glauber model or CGC model for initial condition. In recent years, there is much interest in event-by-event hydrodynamics. Unlike in smooth hydrodynamics, in event-by-event hydrodynamics, initial conditions fluctuate event-by-event. Event-by-event hydrodynamics takes into account the possibility that participant positions can fluctuate from event to event. It was also realised in recent years that the participating nucleons, rather than the reaction plane, determine the symmetry plane of the initial collision zone [109]. In Figure 12, a schematic representation of participating nucleons in the transverse plane in a Monte-Carlo event is shown. Geometric overlap region does not coincide with the participating nucleons. It is obvious that the symmetry plane of the participating nucleons is tilted with respect to the reaction plane. In such a situation, azimuthal angle should be measured with respect to the participant plane rather than the reaction plane. Participant plane can fluctuate, event-by-event, and give rise to the novel phenomena like triangular flow, the third flow harmonic in the Fourier expansion of azimuthal distribution, which will be absent in smooth initial condition. Triangular flow is instrumental in explaining the peculiar structures (known as Ridge) seen in the two particle correlation in plane [110–112]. In Figure 13 two particle correlations in plane, in d+Au and Au+Au collisions are shown. One notices the marked differences in the correlation. Unlike in d+Au collisions, in Au+Au collisions, two particles are correlated over a many units of pseudorapidity . Correlation in the azimuth, however, is narrow. The peculiar structure in two particle correlations known as “ridge” observed both in STAR and PHENIX experiments. The ridge like structure that is also observed in collisions [113, 114]. The ridge structure has most compelling explanation provided the third flow harmonic, and the triangular flow develops in the collisions [115–117]. Specifically, if initial condition is parameterized with quadrapole and triangular moments, response of the medium to these anisotropies is reflected in the two body correlation as ridge [115, 116].

Figure 12 

A schematic representation for distribution of participant nucleons on the transverse plane in Au+Au collisions in a Monte-Carlo event. The participant plane angle is tilted with respect to the reaction plane .