Viscous Hydrodynamic Model for Relativistic Heavy Ion Collisions

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 themodeling, 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.


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][2][3][4][5][6], with the crossover temperature (or pseudocritical temperature)   ≈ 150-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 √  = 200 GeV Au+Au collisions at Relativistic Heavy Ion Collider (RHIC) [7][8][9][10] and √  = 2. 76 TeV Pb+Pb collisions at Large Hadron Collider (LHC) [11][12][13][14][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 Advances in High Energy Physics 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 freezeout.It is called kinetic freeze-out.Hadrons from the kinetic freeze out surface will be detected in the detector.
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 √  = 200 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 / ≥ 1/4 [18][19][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   ≈ 1.5 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, / ≥ 1/4, nevertheless it is nonzero.Indeed, much earlier, Danielewicz and Gyulassy [23], from the uncertainty principle, estimated the lower bound of viscosity to entropy ratio, / ≥ 1/12, 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].

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  2 = 1).One also defines a projector, Δ ] =  ] −    ] , orthogonal to the 4-velocity (Δ ]  ] = 0).In equilibrium, an unique 4-velocity (  ) exists such that the particle density (), energy density (), and the entropy density () can be obtained from eq =   . ( 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, (  ) =   eq  −   eq   , following relations can be established on the equilibrium hypersurface Σ eq (,   ),   eq = −  eq +     eq .
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 Φ  .Π = −(1/3)Δ ]  ] −  is the bulk viscous pressure, and is the shear stress tensor.Hydrodynamic 4velocity 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,   =   −   eq , and  ] =  ] −  ] eq .Detail form of   is constrained by the second law     ≥ 0. 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,   = 0, and entropy current contains terms up to first order in deviations,   and  ] .Entropy production rate can be written as where  = −∇ ⋅ ,   = ∇  / −  ]  ]   , and The second law,     ≥ 0, 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.

Advances in High Energy Physics
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,   ̸ = 0.The most general   containing terms up to second order in deviations can be written as where  0 ,  1 , and  2 are thermodynamic coefficients for the bulk stress (Π), heat flow (  ), and shear viscous stress ( ] ), respectively. 0 and  1 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][32][33][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][38][39][40][41][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.

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  ≥ 150 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 confinementdeconfinement transition is a crossover, and the crossover or the pseudocritical temperature for the transition is   ≈ 170 MeV [3][4][5][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  =   ≈ 170 MeV, with the entropy density of the hadronic resonance gas result.The energy density and pressure are then obtained from the thermodynamical relation:

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 u(, , ) = v(, , ) 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 ( = √  2 −  2 , , ,   = (1/2) ln( + )/( − )).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][47][48].I briefly discuss the Glauber and CGC initial conditions.In Glauber model, in A+B collisions at impact parameter b, at the initial time   , the initial entropy density ((, )) or energy density ((, )) is parameterised as or where  part (, ) and  coll (, ) 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][50][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][51][52][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. 1,2 = (  /√) ± 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.2-0.3 is motivated from DIS experiments. ()  part 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   ∼ 1/  , 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 (V 2 /) and the multiplicity density (1/)( ch /) is approximately universal function, only depending on the value of QGP viscosity over entropy ratio, (/) QGP , but not on any details of the model from which  and  =  √ ⟨⟩ 2 ⟨⟩ 2 is computed.Results of their simulations are shown in Figure 3 (for details see [57]).
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, Π(, , ) =  ] (, , ) = 0 [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 √  = 200 GeV Au+Au and √  = 2.76 TeV Pb+Pb collisions is that hydrodynamical simulations require small thermalisation time   ≈ 0.5-1.0fm.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  ∼ √/ hard , with  being the density and  hard 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 ( = / → ∞),  Π =  0 ≈ (6/5)( 2 / 2 )(/) and   = 2 2 ≈ /.In the relativistic limit ( → 0),  Π =  0 = (216/)(/) 4 and   = 2 2 = 3/2.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.

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:  [62,63], lattice QCD calculation [64].AdS/CFT bound / = 1/4 is also shown.
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 / ∼ 1, 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 crosssection.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 AdS 5 ×  5 is equivalent to quantum field theory without gravity defined on the conformal boundary of AdS 5 ; here AdS 5 is the 5dimensional anti-deSitter space, and  5 is a five-dimensional sphere.For more details, see [69].In this approach, energymomentum 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  = √ 8 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  = /4, one can easily obtain the AdS/CFT limit for shear viscosity over entropy ratio [18][19][20]: Bulk viscosity is proportional to interaction measure  − 3, 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 0.06 ≤   ≤ 0.3 / = 0.020 2  [71], and is consistent with the old relation, , with  2  being the speed of sound in the medium.AdS/CFT calculations that incorporate deviations from conformality at strong coupling find  ≈ 2(1/3 −  2  ) [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][76][77][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.

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   = 110-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  = 1/ 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  eq (, ).In a highly nonequilibrium system, distribution function (, ) is unknown.If the system is slightly offequilibrium, 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, −1, and 0 for Fermi, Bose, and Boltzmann gas, respectively.(, ) ≪ 1 is the deviation from equilibrium distribution function  eq .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  0 ,  0 , and  0 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   > 1 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][84][85][86][87].In the hybrid models, below a switching temperature  sw <   , motion of hadrons is governed by transport equations, for example, UrQMD model.At the switching temperature  sw 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.

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][89][90].
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 Ψ RP (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 V  are easily obtained, is called (integrated) hexadecapole flow, and so forth.Similar to (35), one can Fourier expand the invariant distribution The coefficients V  (  ) 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 V 2 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.
The second harmonic coefficient or the elliptic flow (V 2 ) has been studied extensively in √  = 200 GeV Au+Au collisions at RHIC [9,10].Recently, ALICE collaboration measured elliptic flow in √  = 2.76 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][98][99][100][101][102]].It appears that QGP viscosity over entropy ratio is close to / ≈ 1/4.

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  +1/2 .In the second step, using the new velocity and the new source term in the th step, one calculates the quantities at the +1th 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 +1th step to get the final result at the  + 1th 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 √  = 200 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, / = 0 (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, 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 √  = 200 GeV Au+Au collision data, the viscosity over entropy ratio was estimated as / = 0.1 ± 0.1 (theory) ±0.08 (experiment) [47].
In Figure 9, viscous hydrodynamic simulation results [48] for transverse momentum spectra in √  = 200 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.08-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, / ≈ 0.12, and still higher value, / ≈ 0.16, 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.08-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.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][106][107][108].In [105], four different temperatures dependence for / were considered: (i) LH-LQ: / = 0.08 for all temperatures, (ii) LH-HQ: / 0.08 in the hadron gas and above  = 180 MeV / increases according to lattice QCD data [68], (iii) HH-LQ: below  = 180 MeV, / is that of a hadron gas [76], and / = 0.08, 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 √  = 200 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 √  = 2.76 TeV.
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 √ = 200 GeV Au+Au and √  = 2.76 TeV Pb+Pb collisions were confronted with three slope parameters: (i)  = 0 (ii)  = 0.2 and (iii)  = 0.4.Slope parameter  = 0-0.2appears to be consistent with experimental data in √  = 2.76 TeV Pb+Pb collisions.LHC data clearly disfavor very strong temperature dependence.At RHIC energy, however, strong temperature dependence  = 0.4 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.

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][111][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 V 3 develops in the collisions [115][116][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].
For example, the azimuthal correlation function can be Fourier decomposed as where the first component  1Δ is understood to be due to momentum conservation and directed flow, the second component  2Δ is dominated by the contribution from elliptic flow, and the third component  3Δ is dominated by the triangular flow.In Figure 14, PHOBOS and STAR measurements of long range azimuthal correlation in 200 GeV Au+Au collisions are shown.In the top panel of Figure 14, the first three Fourier components of the azimuthal correlations are shown.The bottom panels show the residual after these components are taken out.Evidently, experimental data on two particle correlation are very well described by the three Fourier components.The analysis indicates that the two particle correlation in Δ − Δ plane is consistent with hydrodynamic models, if triangular flow develops during the evolution, which is possible only when initial condition fluctuates.Recently, ALICE collaboration has measured triangular flow in √  = 2.76 TeV Pb+Pb collisions [15].In most central collisions, the elliptic flow (V 2 ) and triangular flow (V 3 ) are of similar magnitude.In peripheral collisions, however, elliptic flow dominates.More recently, PHENIX collaboration [118][119][120] measured triangular flow in √  = 200 GeV Au+Au collisions.
In event-by-event hydrodynamics, one generally uses Monte-Carlo Glauber model or Monte-Carlo CGC (the Kharzeev-Levin-Nardi version) to obtain the initial conditions, event-by-event.Recently, Monte-Carlo CGC model is improved by combining the impact parameter-dependent saturation model [133,134] of high energy nucleon or nuclear wave function with the classical Yang-Mills description of Glasma fields [135].The model is called IP-Glasma model [122] for initial condition.In addition to fluctuations of nucleon positions, IP-Glasma description includes quantum fluctuations of color charges on the length scale determined by the inverse nuclear saturation scale   , missed in MC-KLN models [136].In Figure 15, initial energy density distribution in three types of model, MC-Glauber, MC-KLN, and IP-Glasma is shown.Color charge fluctuations in the length scale  −1  introduces addition fluctuations, and in IP-Glasma model, initial density has finer structures than in MC-KLN model or in MC-Glauber model.
Fluid dynamical models require continuous density distribution; however, in event-by-event hydrodynamics, initial conditions are generated using Monte-Carlo algorithm, and in general the distribution is not continuous.For example, in MC-Glauber model, in an event, participants nucleons positions can be obtained.If a particular MC event has  part participants, participants positions in the transverse plane can be labeled as ( 1 ,  1 ), ( 2 ,  2 ) ⋅ ⋅ ⋅ ( part ,  part ).The energy density in the transverse plane can be approximated as However, the discrete distribution as in ( 42) cannot be evolved.To use in a hydrodynamic model, the discrete density distribution has to be converted into a smooth energy-density distribution.This can be done by smearing the discrete participant positions by some smoothing function, ( −   ,  −   ) → ( −   ,  −   ,  1 ,  2 . ..), with   being parameters of the smoothing function .Effect of smoothing of was studied in [137].In Monte-Carlo Glauber model initial participant positions were smoothed with a Gaussian function of various width and Woods-Saxon function of various diffuseness.It was shown that effect of smoothing of participant positions, on elliptic and triangular flow, is minimum.
One important aspect of event-by-event hydrodynamics is the characterisation of the asymmetry of the initial collision zone as well as the azimuthal angle of the participant plane.Each flow harmonics can have their own participant plane.One can generalise the definition of eccentricity to give a simple ansatz to characterise the asymmetry of the initial collision zone [129], where  =  cos  and  =  sin .Equation ( 43) also determines the participant plane angle    .Asymmetry measures,  2 and  3 , are called eccentricity and triangularity. 4 and  5 Advances in High Energy Physics essentially measure the squareness and five sidedness of the initial distribution.They may be called rectangularity and pentangularity, respectively.Fourth flow coefficient V 4 is generally referred as hexadecapole flow, and rectangular flow may be more appropriate.Similarly V 5 may be referred to as the pentangular flow.
Viscous effects on elliptic and triangular flow, in eventby-event hydrodyanmics, were studied in [92].The results are shown in Figures 16 and 17.In [92], a 3 + 1 dimensional hydrodynamical model was used.As it was in smooth hydrodynamics, effect of viscosity is to reduce the elliptic flow.More viscous is the fluid, and less is the elliptic flow.Except in p T (GeV) Figure 16: Charged hadron elliptic flow (V 2 ) for different centralities as a function of transverse momentum for averaged initial conditions (avg) and event-by-event simulations (e-b-e) using different viscosity to entropy density ratios compared to STAR [33] and PHENIX [34] data.The figures are from [92].
very central collision, where event-by-event hydrodynamics produces less flow than that in smooth hydrodynamics (the curve labeled as "avg"), the increase of V 2 in central 0-5% collision is understood.In event-by-event hydrodynamics, V 2 , is determined in each event.Single events have a larger anisotropy with respect to the event plane than the average with respect to the reaction plane, and V 2 is increased.In more peripheral collisions, the effect reduces, and event-by-event V 2 becomes smaller than that for smooth initial condition.Triangular flow is completely fluctuation driven and depend less strongly on collision centrality than the elliptic flow.Triangular flow also reduces with viscosity.One also understands that triangular flow is more sensitive to viscosity than the elliptic flow.For example, in 30-40% collision, at   ≈ 1.5 GeV, elliptic flow reduces by ∼15% and 30% when / is increased from 0 to 0.08 and 0.12, respectively.The reduction in triangular flow is much larger, respectively, ∼30% and 45%.One expects that experimental data on triangular flow can constrain viscosity over entropy ratio much better than the elliptic flow data.
In Figure 18, event-by-event hydrodynamic simulations, with IP-Glasma initial condition, are compared with ALICE measurements for   spectra of pion, kaon, and proton [92].Viscosity over entropy ratio was fixed at / = 0.2.One find that simulations with IP-Glasma initial condition explains the   spectra of charged particles rather nicely.IP-Glasma initial condition also explain the measured flow coefficients rather nicely.In Figure 19, simulations results are compared with ATLAS measurements.Charged particles flow coefficients V 2 , V 3 , V 4 , and V 5 appear to agree with the experimental measurements.However, whether or not similar fit can be obtained with other values of viscosity over entropy ratio was not explored in [92].
In Figures 16 and 17, only event averaged flows are shown.In event-by-event hydrodynamics, participant plane fluctuates.It is natural then the flow coefficient also fluctuates.Unless the fluctuations are within some reasonable limit, sensitivity of the flow to viscosity over entropy ratio will reduce greatly.In [138] fluctuations in elliptic and triangular flow in event-by-event hydrodynamic were studied.30-40% Pb+Pb collisions at √  = 2.76 TeV was simulated in event-byevent hydrodynamics.Four values of viscosity over entropy ratio were considered, / = 0, 0.08, 0.12, and 0.16.Simulation results are shown in Figures 20 and 21.In Figures 20 and 21, the symbols represent the event averaged values and the bars the variances.As it was shown earlier, event-averaged flow coefficients smoothly decrease with increasing viscosity.However, both elliptic and triangular flow fluctuates strongly.For example, in the   range 1-3 GeV, fluctuations in V 2 (  ) are ∼15-20%.Fluctuations in V 3 (  ) are even more ∼70-80%.The fluctuations in anisotropic flow greatly reduce their efficacy as a diagnostic tool.For example, within the fluctuations, differential elliptic flow does not distinguish between ideal fluid and fluid with viscosity to entropy ratio 0.08.Triangular flow is even more insensitive.Fluid viscosity over entropy ratio / varying between 0 and 0.16 is not distinguished.It was also shown in [138] that the fluctuations are not statistical and cannot be eliminated by increasing the number of events.Reduced sensitivity of elliptic and triangular flow towards viscosity belies the possibility of constraining shear viscosity from a simultaneous fit to elliptic and triangular flow.In a hydrodynamic model, collective flow is a response of the spatial asymmetry of the initial state.For example, elliptic flow is the response of ellipticity of the initial medium.If ellipticity in the initial medium is characterized by spatial eccentricity,  2 = ⟨ 2 −  2 ⟩/⟨ 2 −  2 ⟩, more eccentric is the initial medium, and more flow will be generated, V 2 ∝  2 .Similar correlation is expected between the triangular flow and initial triangularity.Recently in [131], the correlation between asymmetry measures   and flow coefficients V  was studied in event-by-event hydrodynamics.As it was in smooth hydrodynamics, in event-by-event hydrodynamics also, elliptic flow is strongly correlated with initial eccentricity.Triangular flow was also found to be correlated with initial triangularity.Higher flow harmonics, however, are only weakly correlated with initial eccentricity measured.It was shown that to correctly predict V 4 and V 5 , one must take into account nonlinear terms proportional  2 2 and  2  3 , respectively.Correlation between initial (spatial) asymmetry measures and flow coefficients was also studied in [140].Similar results were obtained.Higher flow coefficients are only weakly correlated with initial asymmetry measures, and the correlations are weakened in more peripheral collisions.However, in [131,140] only ideal fluid was considered, dependence on viscosity was not investigated.Since viscosity introduces additional length scales, qualitatively, one can argue that correlation between momentum anisotropy of the produced particles and asymmetry of initial density distribution will reduce in presence of viscosity.Effect of viscosity on the correlation between elliptic flow and initial eccentricity and between triangular flow and initial triangularity was studied in [139].In Figure 22, event-by-event hydrodynamic simulations for elliptic flow in 30-40% Pb+Pb collisions are plotted against the initial eccentricity.For perfect correlation V 2 ∝  2 , all the points should lie on a straight line.In ideal hydrodynamics, elliptic flow is strongly correlated with initial eccentricity.Scattering of the points around a common straight line is very small.Scattering increases in viscous evolution indicating that the correlation is weakened.Even then scattering of the points is not large even for large viscosity / = 0.16, and one understands that the correlation between initial eccentricity and elliptic flow remains strong in ideal  and viscous hydrodynamics.Simulation result for triangular flow is shown in Figure 23.Comparatively large scattering of the points indicates that the correlation between triangular flow and initial triangularity is comparatively weak.Viscosity further weakens the correlation.In [139], a viscous effect on higher flow correlations was not studied.It will be interesting to see the effect of viscosity on higher flow coefficients and in other collision centralities.
Comparatively weak correlation in higher flow harmonics with the eccentricity measures can be understood, qualitatively.In Figure 24, primary shapes of the collision zone, responsible elliptic, triangular, rectangular, and pentangular flow are shown.Thick arrows attached to the shapes indicate the primary direction of particle emission.One notices that elliptic shape for elliptic flow can contribute both to elliptic and to rectangular flow (indicated by the thin arrows) and vice versa.However, since V 4 ≪ V 2 , elliptic flow will be less affected by the contamination than the rectangular flow.Then while elliptic flow and initial eccentricity will be correlated strongly, correlation between rectangular flow and initial rectangularity will be comparatively weak.Similarly, triangular shape can contribute both to triangular flow and pentangular flow.Here again, since V 5 ≪ V 3 , triangular flow will be less affected than the pentangular flow.

Summary
In the preceding sections, I have briefly discussed some aspects of viscous hydrodynamics.Viscous hydrodynamics has been greatly successful in explaining many experimental Figure 22: Correlation between elliptic flow (V 2 ) and initial eccentricity ( 2 ) is shown for four values of fluid viscosity to entropy ratio / = 0, 0.08, 0.12, and 0.16.The figure is from [139].data in high energy collisions, in particular, data related to the bulk production.However, there are several issues that need to be investigated.For example, hydrodynamics modeling of nuclear collisions require a number of parameters to be fixed initially, for example, initial time, initial energy density/velocity/shear (bulk) stress tensor profile, freeze-out, and so forth.Intercorrelation between these parameters is not fully investigated.How the intercorrelations affect sensitive observables like elliptic or triangular flow is not known.Several attempts have been made to extract QGP viscosity over entropy ratio from hydrodynamical analysis of experimental data.However, how the extracted viscosity over entropy ratio depends on the intercorrelation between different initial conditions was not investigated.In [98], an attempt was made to obtain systematic uncertainty in extracted / due to uncertain initial conditions, but the investigation was not complete.The small thermalisation time used in most of the simulations is also an issue.The mechanism of fast thermalisation is not understood.It was mentioned previously that chromo-Weibel instability was found to be inadequate to thermalise the medium in the time scale ∼1 fm.If indeed, thermalisation takes longer time, for example, ∼5 fm in LHC energy and ∼10 fm in RHIC energy [61], then many of our understanding, in particular regarding elliptic flow, need to be readjusted.Freeze-out is also an issue, in particular in event-by-event hydrodynamics.Unlike in smooth initial condition, in event-by-event hydrodynamic, the initial condition is granular.Some interior parts of the fluid can freeze-out early, but the hadrons emitted from that part may not come out but absorbed in the medium.The medium can be reheated.The issue needs to be investigated.Freeze out issue may be less important in hybrid models, where the late stage of evolution is governed by transport equations.But hybrid model requires an arbitrary parameter, the switching temperature, below which hydrodynamic description give way to transport description.The issue of temperature dependence of viscosity over entropy ratio is also not fixed.Is the temperature dependence linear or more complex?Whether or not the experimental data on bulk production can fix empirically the temperature dependence of QGP viscosity over entropy ratio is uncertain.In the present review, bulk viscosity was not discussed in details.Bulk viscosity may also become an important issue in understanding experimental data in terms of hydrodynamic models.As it was a noted earlier, correct form for the nonequilibrium correction due to bulk viscosity is still an issue and needed to be investigated in detail.

Figure 1 :
Figure 1: Different stages of nuclear collisions at ultrarelativistic energy.

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

Figure 3 :
Figure3: Hybrid model simulation for the eccentricity-scaled elliptic flow V 2 / as a function of multiplicity density (1/)( ch /), for different values of (/) QGP .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.

Figure 5 :
Figure 5: Schematic picture of a nonzero impact parameter collision.

Figure 6 :
Figure 6: Ideal and viscous hydrodynamic simulations for evolution of spatial eccentricity (  ) and momentum anisotropy (  ) in  = 7.4 fm Pb+Pb collisions.

16 Figure 7 :
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 :
Figure 8: Centrality dependence of charged particles (integrated) elliptic flow compared with experimental data.The figure is from [47].

Figure 9 :Figure 10 :Figure 11 :
Figure 9: PHENIX measurements of charged particles   spectra in √  = 200 GeV Au+Au collisions are compared with viscous hydrodynamic simulations.The figures are from [48].

Figure 12 :
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 ( PP ) is tilted with respect to the reaction plane ( RP ).

Figure 13 :Figure 14 :Figure 15 :
Figure 13: Two particle correlations in Δ − Δ in d+Au and Au+Au collisions are shown.The results are from STAR experiment.

Figure 17 :
Figure17: Same as in Figure16but for triangular flow.The figures are from[92].

Figure 18 :Figure 19 :
Figure18: ALICE measurements for the transverse momentum spectra of pion, kaon, and protons in 0-5% Pb+Pb collisions at LHC energy are compared with IP-Glasma model predictions with viscosity over entropy ratio / = 0.08.The figures are from[123].

Figure 21 :
Figure 21: Same as in Figure 20 but for triangular flow.The figures are from [138].

Table 1 :
[95]ous viscous hydrodynamics models are listed.Glauber model, CGC: color glass condensate, MC-G: Monte-Carlo Glauber, MC-CGC: Monte Carlo color glass condensate, and CD: central difference.The table is from[95].I have added the last reference, which uses smooth particle hydrodynamics (SPH) for solving hydrodynamic equations.