Relativistic hydrodynamics in heavy-ion collisions: general aspects and recent developments

Relativistic hydrodynamics has been quite successful in explaining the collective behaviour of the QCD matter produced in high energy heavy-ion collisions at RHIC and LHC. We briefly review the latest developments in the hydrodynamical modeling of relativistic heavy-ion collisions. Essential ingredients of the model such as the hydrodynamic evolution equations, dissipation, initial conditions, equation of state, and freeze-out process are reviewed. We discuss observable quantities such as particle spectra and anisotropic flow and effect of viscosity on these observables. Recent developments such as event-by-event fluctuations, flow in small systems (proton-proton and proton-nucleus collisions), flow in ultra central collisions, longitudinal fluctuations and correlations and flow in intense magnetic field are also discussed.


I. INTRODUCTION
The existence of both confinement and asymptotic freedom in Quantum Chromodynamics (QCD) has led to many speculations about its thermodynamic and transport properties. Due to confinement, the nuclear matter must be made of hadrons at low energies, hence it is expected to behave as a weakly interacting gas of hadrons. On the other hand, at very high energies asymptotic freedom implies that quarks and gluons interact only weakly and the nuclear matter is expected to behave as a weakly coupled gas of quarks and gluons. In between these two configurations there must be a phase transition where the hadronic degrees of freedom disappear and a new state of matter, in which the quark and gluon degrees of freedom manifest directly over a certain volume, is formed. This new phase of matter, referred to as Quark-Gluon Plasma (QGP), is expected to be created when sufficiently high temperatures or densities are reached [1][2][3].
The QGP is believed to have existed in the very early universe (a few microseconds after the Big Bang), or some variant of which possibly still exists in the inner core of a neutron star where it is estimated that the density can reach values ten times higher than those of ordinary nuclei. It was conjectured theoretically that such extreme conditions can also be realized on earth, in a controlled experimental environment, by colliding two heavy nuclei with ultra-relativistic energies [4]. This may transform a fraction of the kinetic energies of the two colliding nuclei into heating the QCD vacuum within an extremely small volume where temperatures million times hotter than the core of the sun may be achieved.
With the advent of modern accelerator facilities, ultrarelativistic heavy-ion collisions have provided an opportunity to systematically create and study different phases of the bulk nuclear matter. It is widely believed that the QGP phase is formed in heavy-ion collision experi-ments at Relativistic Heavy-Ion Collider (RHIC) located at Brookhaven National Laboratory, USA and Large Hadron Collider (LHC) at European Organization for Nuclear Research (CERN), Geneva. A number of indirect evidences found at the Super Proton Synchrotron (SPS) at CERN, strongly suggested the formation of a "new state of matter" [5], but quantitative and clear results were only obtained at RHIC energies [6][7][8][9][10][11][12][13][14], and recently at LHC energies [15][16][17][18]. The regime with relatively large baryon chemical potential will be probed by the upcoming experimental facilities like Facility for Anti-proton and Ion Research (FAIR) at GSI, Darmstadt. An illustration of the QCD phase diagram and the regions probed by these experimental facilities is shown in Fig. 1 [19].
It is possible to create hot and dense nuclear matter with very high energy densities in relatively large volumes by colliding ultra-relativistic heavy ions. In these conditions, the nuclear matter created may be close to (local) thermodynamic equilibrium, providing the opportunity to investigate the various phases and the thermodynamic and transport properties of QCD. It is important to note that, even though it appears that a deconfined state of matter is formed in these colliders, investigating and extracting the transport properties of QGP from heavy-ion collisions is not an easy task since it cannot be observed directly. Experimentally, it is only feasible to measure energy and momenta of the particles produced in the final stages of the collision, when nuclear matter is already relatively cold and non-interacting. Hence, in order to study the thermodynamic and transport properties of the QGP, the whole heavy-ion collision process from the very beginning till the end has to be modelled: starting from the stage where two highly Lorentz contracted heavy nuclei collide with each other, the formation and thermalization of the QGP or de-confined phase in the initial stages of the collision, its subsequent spacetime evolution, the phase transition to the hadronic or arXiv:1605.08694v3 [nucl-th] 27 Oct 2016 confined phase of matter, and eventually, the dynamics of the cold hadronic matter formed in the final stages of the collision. The different stages of ultra-relativistic heavy-ion collisions are schematically illustrated in Fig. 2 [20].
Assuming that thermalization is achieved in the early stages of heavy-ion collisions and that the interaction between the quarks is strong enough to maintain local thermodynamic equilibrium during the subsequent expansion, the time evolution of the QGP and hadronic matter can be described by the laws of fluid dynamics [21][22][23][24]. Fluid dynamics, also loosely referred to as hydrodynamics, is an effective approach through which a system can be described by macroscopic variables, such as local energy density, pressure, temperature and flow velocity. The most appealing feature of fluid dynamics is the fact that it is simple and general. It is simple in the sense that all the information of the system is contained in its thermodynamic and transport properties, i.e., its equation of state and transport coefficients. Fluid dynamics is also general because it relies on only one assumption: the system remains close to local thermodynamic equilibrium throughout its evolution. Although the hypothesis of proximity to local equilibrium is quite strong, it saves us from making any further assumption regarding the description of the particles and fields, their interactions, the classical or quantum nature of the phenomena involved etc. Hydrodynamic analysis of the spectra and azimuthal anisotropy of particles produced in heavy-ion collisions at RHIC [25,26] and recently at LHC [27,28] suggests that the matter formed in these collisions is strongly-coupled quark-gluon plasma (sQGP).
Application of viscous hydrodynamics to high-energy heavy-ion collisions has evoked widespread interest ever since a surprisingly small value of η/s was estimated from the analysis of the elliptic flow data [25]. It is interesting to note that in the strong coupling limit of a large class of holographic theories, the value of the shear viscosity to entropy density ratio η/s = 1/4π. Kovtun, Son and Starinets (KSS) conjectured this strong coupling result to be the absolute lower bound for all fluids, i.e., η/s ≥ 1/4π [29,30]. This specific combination of hydrodynamic quantities, η/s, accounts for the difference between momentum and charge diffusion such that even though the diffusion constant goes to zero in the strong coupling limit, the ratio η/s remains finite [31]. Similar result for a lower bound also follows from the quantum mechanical uncertainty principle. The kinetic theory prediction for viscosity, η = 1 3 nl mfpp , suggests that low viscosity corresponds to short mean free path. On the other hand, the uncertainty relation implies that the product of the mean free path and the average momentum cannot be arbitrarily small, i.e., l mfpp 1. For a weakly interacting relativistic Bose gas, the entropy per particle is given by s/n = 3.6. This leads to η/s 0.09 which is very close to the lower KSS bound.
Indeed the estimated η/s of QGP was so close to the KSS bound that it led to the claim that the matter formed at RHIC was the most perfect fluid ever observed. A precise estimate of η/s is vital to the understanding of the properties of the QCD matter and is presently a topic of intense investigation, see [32] and references therein for more details. In this review, we shall discuss the general aspects of the formulation of relativistic fluid dynamics and its application to the physics of high-energy heavy-ion collisions. Along with these general concepts we shall also discuss here some of the recent developments in the field. Among the recent developments the most striking feature is the experimental observation of flow like pattern in the particle azimuthal distribution of high multiplicity proton-proton (p-p) and proton-lead (p-Pb)collisions. We will discuss in this review the success of hydrodynamics model in describing these recent experimental measurements by assuming hydrodynamics flow of small systems.
The review is organised as follows, In section II we discuss the general formalism of a causal theory of relativistic dissipative fluid dynamics. Section III deals with the initial conditions necessary for the modeling of relativistic heavy-ion, p-p and p-Pb collisions. In section IV, we discuss some models of pre-equilibrium dynamics employed before hydrodynamic evolution. Section V briefly covers various equations of state, necessary to close the hydrodynamic equations. In section VI, particle production mechanism via Cooper-Frye freeze-out and anisotropic flow generation is discussed. In Section VII, we discuss models of hadronic rescattering after freezeout and contribution to particle spectra and flow from resonance decays. Section VIII deals with the extraction of transport coefficients from hydrodynamic analysis of flow data. Finally, in section IX, we discuss recent developments in the hydrodynamic modeling of relativistic collisions.
In this review, unless stated otherwise, all physical quantities are expressed in terms of natural units, where, = c = k B = 1, with = h/2π, where h is the Planck constant, c the velocity of light in vacuum, and k B the Boltzmann constant. Unless stated otherwise, the spacetime is always taken to be Minkowskian where the metric tensor is given by g µν = diag(+1, −1, −1, −1). Apart from Minkowskian coordinates x µ = (t, x, y, z), we will also regularly employ Milne coordinate system x µ = (τ, x, y, η s ) or x µ = (τ, r, ϕ, η s ), with proper time τ = √ t 2 − z 2 , the radial coordinate r = x 2 + y 2 , the azimuthal angle ϕ = tan −1 (y/x), and spacetime rapidity η s = tanh −1 (z/t). Hence, t = τ cosh η s , x = r cos ϕ, y = r sin ϕ, and z = τ sinh η s . For the coordinate system x µ = (τ, x, y, η s ), the metric becomes g µν = diag(1, −1, −1, −τ 2 ), whereas for x µ = (τ, r, ϕ, η s ), the metric is g µν = diag(1, −1, −r 2 , −τ 2 ). Roman letters are used to indicate indices that vary from 1-3 and Greek letters for indices that vary from 0-3. Covariant and contravariant four-vectors are denoted as p µ and p µ , respectively. The notation p · q ≡ p µ q µ represents scalar product of a covariant and a contravariant four-vector. Tensors without indices shall always correspond to Lorentz scalars. We follow Einstein summation convention, which states that repeated indices in a single term are implicitly summed over all the values of that index.
We denote the fluid four-velocity by u µ and the Lorentz contraction factor by γ. The projector onto the space orthogonal to u µ is defined as: ∆ µν ≡ g µν −u µ u ν . Hence, ∆ µν satisfies the conditions ∆ µν u µ = ∆ µν u ν = 0 with trace ∆ µ µ = 3. The partial derivative ∂ µ can then be decomposed as: In the fluid rest frame, D reduces to the time derivative and ∇ µ reduces to the spacial gradient. Hence, the notationḟ ≡ Df is also commonly used. We also frequently use the symmetric, anti-symmetric and angular brackets notations defined as is the traceless symmetric projection operator orthogonal to u µ satisfying the conditions ∆ αβ µν ∆ αβ = ∆ αβ µν ∆ µν = 0.

II. RELATIVISTIC FLUID DYNAMICS
The physical characterization of a system consisting of many degrees of freedom is in general a non-trivial task. For instance, the mathematical formulation of a theory describing the microscopic dynamics of a system containing a large number of interacting particles is one of the most challenging problems of theoretical physics. However, it is possible to provide an effective macroscopic description, over large distance and time scales, by taking into account only the degrees of freedom that are relevant at these scales. This is a consequence of the fact that on macroscopic distance and time scales the actual degrees of freedom of the microscopic theory are imperceptible. Most of the microscopic variables fluctuate rapidly in space and time, hence only average quantities resulting from the interactions at the microscopic level can be observed on macroscopic scales. These rapid fluctuations lead to very small changes of the average values, and hence are not expected to contribute to the macroscopic dynamics. On the other hand, variables that do vary slowly, such as the conserved quantities, are expected to play an important role in the effective description of the system. Fluid dynamics is one of the most common examples of such a situation. It is an effective theory describing the long-wavelength, low frequency limit of the underlying microscopic dynamics of a system.
A fluid is defined as a continuous system in which every infinitesimal volume element is assumed to be close to thermodynamic equilibrium and to remain near equilibrium throughout its evolution. Hence, in other words, in the neighbourhood of each point in space, an infinitesimal volume called fluid element is defined in which the matter is assumed to be homogeneous, i.e., any spatial gradients can be ignored, and is described by a finite number of thermodynamic variables. This implies that each fluid element must be large enough, compared to the microscopic distance scales, to guarantee the proximity to thermodynamic equilibrium, and, at the same time, must be small enough, relative to the macroscopic distance scales, to ensure the continuum limit. The coexistence of both continuous (zero volume) and thermodynamic (infinite volume) limits within a fluid volume might seem paradoxical at first glance. However, if the microscopic and the macroscopic length scales of the system are sufficiently far apart, it is always possible to establish the existence of a volume that is small enough compared to the macroscopic scales, and at the same time, large enough compared to the microscopic ones. Here, we will assume the existence of a clear separation between microscopic and macroscopic scales to guarantee the proximity to local thermodynamic equilibrium.
Relativistic fluid dynamics has been quite successful in explaining the various collective phenomena observed in astrophysics, cosmology and the physics of high-energy heavy-ion collisions. In cosmology and certain areas of astrophysics, one needs a fluid dynamics formulation consistent with the General Theory of Relativity [33]. On the other hand, a formulation based on the Special Theory of Relativity is quite adequate to treat the evolution of the strongly interacting matter formed in high-energy heavy-ion collisions when it is close to a local thermodynamic equilibrium. In fluid dynamical approach, although no detailed knowledge of the microscopic dynamics is needed, however, knowledge of the equation of state relating pressure, energy density and baryon density is required. The collective behaviour of the hot and dense quark-gluon plasma created in ultra-relativistic heavyion collisions has been studied quite extensively within the framework of relativistic fluid dynamics. In application of fluid dynamics, it is natural to first employ the simplest version which is ideal hydrodynamics which neglects the viscous effects and assumes that local equilibrium is always perfectly maintained during the fireball expansion. Microscopically, this requires that the microscopic scattering time be much shorter than the macroscopic expansion (evolution) time. In other words, ideal hydrodynamics assumes that the mean free path of the constituent particles is much smaller than the system size. However, as all fluids are dissipative in nature due to the quantum mechanical uncertainty principle [34], the ideal fluid results serve only as a benchmark when dissipative effects become important.
When discussing the application of relativistic dissipative fluid dynamics to heavy-ion collision, one is faced with yet another predicament: the theory of relativistic dissipative fluid dynamics is not yet conclusively established. In fact, introducing dissipation in relativistic fluids is not at all a trivial task and still remains one of the important topics of research in high-energy physics. Ideal hydrodynamics assumes that local thermodynamic equilibrium is perfectly maintained and each fluid element is homogeneous, i.e., spatial gradients are absent (zeroth order in gradient expansion). If this is not satisfied, dissipative effects come into play. The earliest theoretical formulations of relativistic dissipative hydrodynamics also known as first-order theories, are due to Eckart [35] and Landau-Lifshitz [36]. However, these formulations, collectively called relativistic Navier-Stokes (NS) theory, suffer from acausality and numerical instability. The reason for the acausality is that in the gradient expansion the dissipative currents are linearly proportional to gradients of temperature, chemical potential, and velocity, resulting in parabolic equations. Thus, in Navier-Stokes theory the gradients have an instantaneous influence on the dissipative currents. Such instantaneous effects tend to violate causality and cannot be allowed in a covariant setup, leading to the instabilities investigated in Refs. [37][38][39].
The second-order Israel-Stewart (IS) theory [40], restores causality but may not guarantee stability [41]. The acausality problems were solved by introducing a time delay in the creation of the dissipative currents from gradients of the fluid-dynamical variables. In this case, the dissipative quantities become independent dynamical variables obeying equations of motion that describe their relaxation towards their respective Navier-Stokes values. The resulting equations are hyperbolic in nature which preserves causality. Israel-Stewart theory has been widely applied to ultra-relativistic heavy-ion collisions in order to describe the time evolution of the QGP and the subsequent freeze-out process of the hadron resonance gas. Although IS hydrodynamics has been quite successful in modelling relativistic heavy-ion collisions, there are several inconsistencies and approximations in its formulation which prevent proper understanding of the thermodynamic and transport properties of the QGP. Moreover, the second-order IS theory can be derived in several ways, each leading to a different set of transport coefficients. Therefore, in order to quantify the transport properties of the QGP from experiment and confirm the claim that it is indeed the most perfect fluid ever created, the theoretical foundations of relativistic dissipative fluid dynamics must be first addressed and clearly understood. In this section, we review the basic aspects of thermodynamics and discuss the formulation of relativistic fluid dynamics from a phenomenological perspective. The salient features of kinetic theory in the context of fluid dynamics will also be discussed.

A. Thermodynamics
Thermodynamics is an empirical description of the macroscopic or large-scale properties of matter and it makes no hypotheses about the small-scale or microscopic structure. It is concerned only with the average behaviour of a very large number of microscopic constituents, and its laws can be derived from statistical mechanics. A thermodynamic system can be described in terms of a small set of extensive variables, such as volume (V ), the total energy (E), entropy (S), and number of particles (N ), of the system. Thermodynamics is based on four phenomenological laws that explain how these quantities are related and how they change with time [42][43][44].
• Zeroth Law: If two systems are both in thermal equilibrium with a third system then they are in thermal equilibrium with each other. This law helps define the notion of temperature.
• First Law: All the energy transfers must be accounted for to ensure the conservation of the total energy of a thermodynamic system and its surroundings. This law is the principle of conservation of energy.
• Second Law: An isolated physical system spontaneously evolves towards its own internal state of thermodynamic equilibrium. Employing the notion of entropy, this law states that the change in entropy of a closed thermodynamic system is always positive or zero.
• Third Law: Also known an Nernst's heat theorem, states that the difference in entropy between systems connected by a reversible process is zero in the limit of vanishing temperature. In other words, it is impossible to reduce the temperature of a system to absolute zero in a finite number of operations.
The first law of thermodynamics postulates that the changes in the total energy of a thermodynamic system must result from: (1) heat exchange, (2) the mechanical work done by an external force, and (3) from particle exchange with an external medium. Hence the conservation law relating the small changes in state variables, E, V , and N is where P and µ are the pressure and chemical potential, respectively, and δQ is the amount of heat exchange. The heat exchange takes into account the energy variations due to changes of internal degrees of freedom that are not described by the state variables. The heat itself is not a state variable since it can depend on the past evolution of the system and may take several values for the same thermodynamic state. However, when dealing with reversible processes (in time), it becomes possible to assign a state variable related to heat. This variable is the entropy, S , and is defined in terms of the heat exchange as δQ = T δS, with the temperature T being the proportionality constant. Then, when considering variations between equilibrium states that are infinitesimally close to each other, it is possible to write the first law of thermodynamics in terms of differentials of the state variables, Hence, using Eq. (2), the intensive quantities, T , µ and P , can be obtained in terms of partial derivatives of the entropy as (3) The entropy is mathematically defined as an extensive and additive function of the state variables, which means that Differentiating both sides with respect to λ, we obtain which holds for any arbitrary value of λ. Setting λ = 1 and using Eq. (3), we obtain the so-called Euler's relation Using Euler's relation, Eq. (6), along with the first law of thermodynamics, Eq. (2), we arrive at the Gibbs-Duhem relation In terms of energy, entropy and number densities defined as ≡ E/V , s ≡ S/V , and n ≡ N/V respectively, the Euler's relation, Eq. (6) and Gibbs-Duhem relation, Eq. (7), reduce to = −P + T s + µ n (8) dP = s dT + n dµ.
Differentiating Eq. (8) and using Eq. (9), we obtain the relation analogous to first law of thermodynamics It is important to note that all the densities defined above ( , s, n) are intensive quantities. The equilibrium state of a system is defined as a stationary state where the extensive and intensive variables of the system do not change. We know from the second law of thermodynamics that the entropy of an isolated thermodynamic system must either increase or remain constant. Hence, if a thermodynamic system is in equilibrium, the entropy of the system being an extensive variable, must remain constant. On the other hand, for a system that is out of equilibrium, the entropy must always increase. This is an extremely powerful concept that will be extensively used in this section to constrain and derive the equations of motion of a dissipative fluid. This concludes a brief outline of the basics of thermodynamics; for a more detailed review, see Ref. [44]. In the next section, we introduce and derive the equations of relativistic ideal fluid dynamics.

B. Relativistic ideal fluid dynamics
An ideal fluid is defined by the assumption of local thermal equilibrium, i.e., all fluid elements must be exactly in thermodynamic equilibrium [36,45]. This means that at each space-time coordinate of the fluid x ≡ x µ , there can be assigned a temperature T (x), a chemical potential µ(x), and a collective four-velocity field, The proper time increment dτ is given by the line element where v ≡ d x/dt. This implies that where γ( v) = 1/ √ 1 − v 2 . In the non-relativistic limit, we obtain u µ (x) = (1, v). It is important to note that the four-vector u µ (x) only contains three independent components since it obeys the relation The quantities T , µ and u µ are often referred to as the primary fluid-dynamical variables. The state of a fluid can be completely specified by the densities and currents associated with conserved quantities, i.e., energy, momentum, and (net) particle number. For a relativistic fluid, the state variables are the energy-momentum tensor, T µν , and the (net) particle four-current, N µ . To obtain the general form of these currents for an ideal fluid, we first define the local rest frame (LRF) of the fluid. In this frame, v = 0, and the energy-momentum tensor, T µν LRF , the (net) particle four-current, N µ LRF , and the entropy four-current, S µ LRF , should have the characteristic form of a system in static equilibrium. In other words, in local rest frame, there is no flow of energy (T i0 LRF = 0), the force per unit surface element is isotropic (T ij LRF = δ ij P ) and there is no particle and entropy flow ( N = 0 and S = 0). Consequently, the energy-momentum tensor, particle and entropy fourcurrents in this frame take the following simple forms For an ideal relativistic fluid, the general form of the energy-momentum tensor, T µν (0) , (net) particle fourcurrent, N µ (0) , and the entropy four-current, S µ (0) , has to be built out of the hydrodynamic tensor degrees of freedom, namely the vector, u µ , and the metric tensor, g µν . Since T µν (0) should be symmetric and transform as a tensor, and, N µ (0) and S µ (0) should transform as a vector, under Lorentz transformations, the most general form allowed is therefore In the local rest frame, v = 0 ⇒ u µ = (1, 0). Hence in this frame, Eq. (16) takes the form By comparing the above equation with the corresponding general expressions in the local rest frame, Eq. (15), one obtains the following expressions for the coefficients The conserved currents of an ideal fluid can then be expressed as where ∆ µν = g µν − u µ u ν is the projection operator onto the three-space orthogonal to u µ , and satisfies the following properties of an orthogonal projector, The dynamical description of an ideal fluid is obtained using the conservation laws of energy, momentum and (net) particle number. These conservation laws can be mathematically expressed using the four-divergences of energy-momentum tensor and particle four-current which leads to the following equations, where the partial derivative ∂ µ ≡ ∂/∂x µ transforms as a covariant vector under Lorentz transformations. Using the four-velocity, u µ , and the projection operator, ∆ µν , the derivative, ∂ µ , can be projected along and orthogonal to u µ Projection of energy-momentum conservation equation along and orthogonal to u µ together with the conservation law for particle number, leads to the equations of motion of ideal fluid dynamics, where θ ≡ ∂ µ u µ . It is important to note that an ideal fluid is described by four fields, , P , n, and u µ , corresponding to six independent degrees of freedom. The conservation laws, on the other hand, provide only five equations of motion. The equation of state of the fluid, P = P (n, ), relating the pressure to other thermodynamic variables has to be specified to close this system of equations. The existence of equation of state is guaranteed by the assumption of local thermal equilibrium and hence the equations of ideal fluid dynamics are always closed.

C. Covariant thermodynamics
In the following, we re-write the equilibrium thermodynamic relations derived in Sec. 2.1, Eqs. (8), (9), and (10), in a covariant form [40,46]. For this purpose, it is convenient to introduce the following notations In these notations, the covariant version of the Euler's relation, Eq. (8), and the Gibbs-Duhem relation, Eq. (9), can be postulated as, respectively. The above equations can then be used to derive a covariant form of the first law of thermodynamics, Eq. (10), The covariant thermodynamic relations were constructed in such a way that when Eqs. (27), (28) and (29) are contracted with u µ , we obtain the usual thermodynamic relations, Eqs. (8), (9), and (10). Here we have used the property of the fluid four-velocity, u µ u µ = 1 ⇒ u µ du µ = 0. The projection of Eqs. (27), (28) and (29) onto the three-space orthogonal to u µ just leads to trivial identities, From the above equations we conclude that the covariant thermodynamic relations do not contain more information than the usual thermodynamic relations. The first law of thermodynamics, Eq. (29), leads to the following expression for the entropy four-current divergence, After employing the conservation of energy-momentum and net particle number, Eq. (21), the above equation leads to the conservation of entropy, ∂ µ S µ (0) = 0. It is important to note that within equilibrium thermodynamics, the entropy conservation is a natural consequence of energy-momentum and particle number conservation, and the first law of thermodynamics. The equation of motion for the entropy density is then obtained as We observe that the rate equation of the entropy density in the above equation is identical to that of the net particle number, Eq. (25). Therefore, we conclude that for ideal hydrodynamics, the ratio of entropy density to number density (s/n) is a constant of motion.

D. Relativistic dissipative fluid dynamics
The derivation of relativistic ideal fluid dynamics proceeds by employing the properties of the Lorentz transformation, the conservation laws, and most importantly, by imposing local thermodynamic equilibrium. It is important to note that while the properties of Lorentz transformation and the conservation laws are robust, the assumption of local thermodynamic equilibrium is a strong restriction. The deviation from local thermodynamic equilibrium results in dissipative effects, and, as all fluids are dissipative in nature due to the uncertainty principle [34], the assumption of local thermodynamic equilibrium is never strictly realized in practice. In the following, we consider a more general theory of fluid dynamics that attempts to take into account the dissipative processes that must happen, because a fluid can never maintain exact local thermodynamic equilibrium throughout its dynamical evolution.
Dissipative effects in a fluid originate from irreversible thermodynamic processes that occur during the motion of the fluid. In general, each fluid element may not be in equilibrium with the whole fluid, and, in order to approach equilibrium, it exchanges heat with its surroundings. Moreover, the fluid elements are in relative motion and can also dissipate energy by friction. All these processes must be included in order to obtain a reasonable description of a relativistic fluid.
The earliest covariant formulation of dissipative fluid dynamics were due to Eckart [35], in 1940, and, later, by Landau and Lifshitz [36], in 1959. The formulation of these theories, collectively known as first-order theories (order of gradients), was based on a covariant generalization of the Navier-Stokes theory. The Navier-Stokes theory, at that time, had already become a successful theory of dissipative fluid dynamics. It was employed efficiently to describe a wide variety of non-relativistic fluids, from weakly coupled gases such as air, to strongly coupled fluids such as water. Hence, a relativistic generalisation of Navier-Stokes theory was considered to be the most effective and promising way to describe relativistic dissipative fluids.
The formulation of relativistic dissipative hydrodynamics turned out to be more subtle since the relativistic generalisation of Navier-Stokes theory is intrinsically unstable [37][38][39]. The source of such instability is attributed to the inherent acausal behaviour of this theory [47,48]. A straightforward relativistic generalisation of Navier-Stokes theory allows signals to propagate with infinite speed in a medium. While in non-relativistic theories, this does not give rise to an intrinsic problem and can be ignored, in relativistic systems where causality is a physical property that is naturally preserved, this feature leads to intrinsically unstable equations of motion. Nevertheless, it is instructive to review the first-order theories as they are an important initial step to illustrate the basic features of relativistic dissipative fluid-dynamics.
As in the case of ideal fluids, the basic equations governing the motion of dissipative fluids are also obtained from the conservation laws of energy-momentum and (net) particle number, However, for dissipative fluids, the energy-momentum tensor is no longer diagonal and isotropic in the local rest frame. Moreover, due to diffusion, the particle flow is expected to appear in the local rest frame of the fluid element. To account for these effects, dissipative currents τ µν and n µ are added to the previously derived ideal currents, T µν (0) and N µ (0) , where, τ µν is required to be symmetric (τ µν = τ νµ ) in order to satisfy angular momentum conservation. The main objective then becomes to find the dynamical or constitutive equations satisfied by these dissipative currents.

Matching conditions
The introduction of the dissipative currents causes the equilibrium variables to be ill-defined, since the fluid can no longer be considered to be in local thermodynamic equilibrium. Hence, in a dissipative fluid, the thermodynamic variables can only be defined in terms of an artificial equilibrium state, constructed such that the thermodynamic relations are valid as if the fluid were in local thermodynamic equilibrium. The first step to construct such an equilibrium state is to define and n as the total energy and particle density in the local rest frame of the fluid, respectively. This is guaranteed by imposing the so-called matching or fitting conditions [40], These matching conditions enforces the following constraints on the dissipative currents Subsequently, using n and , an artificial equilibrium state can be constructed with the help of the equation of state. It is however important to note that while the energy and particle densities are physically defined, all the other thermodynamic quantities (s, P, T, µ, · · · ) are defined only in terms of an artificial equilibrium state and do not necessarily retain their usual physical meaning.

Tensor decompositions of dissipative quantities
To proceed further, it is convenient to decompose τ µν in terms of its irreducible components, i.e., a scalar, a four-vector, and a traceless and symmetric second-rank tensor. Moreover, this tensor decomposition must be consistent with the matching or orthogonality condition, Eq. (42), satisfied by τ µν . To this end, we introduce another projection operator, the double symmetric, traceless projector orthogonal to u µ , with the following properties, The parentheses in the above equation denote symmetrization of the Lorentz indices, i.e., A (µν) ≡ (A µν + A νµ )/2. The dissipative current τ µν then can be tensor decomposed in its irreducible form by using u µ , ∆ µν and ∆ µν αβ as where we have defined The scalar Π is the bulk viscous pressure, the vector h µ is the energy-diffusion four-current, and the second-rank tensor π µν is the shear-stress tensor. The properties of the projection operators ∆ µ α and ∆ µν αβ imply that both h µ and π µν are orthogonal to u µ and, additionally, π µν is traceless. Armed with these definitions, all the irreducible hydrodynamic fields are expressed in terms of N µ and T µν as where the angular bracket notations are defined as, We observe that T µν is a symmetric second-rank tensor with ten independent components and N µ is a fourvector; overall they have fourteen independent components. Next we count the number of independent components in the tensor decompositions of T µν and N µ . Since n µ and h µ are orthogonal to u µ , they can have only three independent components each. The shear-stress tensor π µν is symmetric, traceless and orthogonal to u µ , and hence, can have only five independent components. Together with u µ , , n and Π, which have in total six independent components (P is related to via equation of state), we count a total of seventeen independent components, three more than expected. The reason being that so far, the velocity field u µ was introduced as a general normalized four-vector and was not specified. Hence u µ has to be defined to reduce the number of independent components to the correct value.

Definition of the velocity field
In the process of formulating the theory of dissipative fluid dynamics, the next important step is to fix u µ . In the case of ideal fluids, the local rest frame was defined as the frame in which there is, simultaneously, no net energy and particle flow. While the definition of local rest frame was unambiguous for ideal fluids, this definition is no longer possible in the case of dissipative fluids due to the presence of both energy and particle diffusion. From a mathematical perspective, the fluid velocity can be defined in numerous ways. However, from the physical perspective, there are two natural choices. The Eckart definition [35], in which the velocity is defined by the flow of particles and the Landau definition [36], in which the velocity is specified by the flow of the total energy We note that the above two definitions of u µ impose different constraints on the dissipative currents. In the Eckart definition the particle diffusion is always set to zero, while in the Landau definition, the energy diffusion is zero. In other words, the Eckart definition of the velocity field eliminates any diffusion of particles whereas the Landau definition eliminates any diffusion of energy. In this review, we shall always use the Landau definition, Eq. (50). The conserved currents in this frame take the following form As done for ideal fluids, the energy-momentum conservation equation in Eq. (38) is decomposed parallel and orthogonal to u µ . Using Eq. (51) together with the conservation law for particle number in Eq. (38), leads to the equations of motion for dissipative fluids. For u µ ∂ ν T µν = 0, ∆ α µ ∂ ν T µν = 0 and ∂ µ N µ = 0, one ob-tains˙ respectively. HereȦ ≡ DA = u µ ∂ µ A, and the shear tensor σ µν ≡ ∇ µ u ν = ∆ µν αβ ∇ α u β . We observe that while there are fourteen total independent components of T µν and N µ , Eqs. (52)-(54) constitute only five equations. Therefore, in order to derive the complete set of equations for dissipative fluid dynamics, one still has to obtain the additional nine equations of motion that will close Eqs. (52)- (54). Eventually, this corresponds to finding the closed dynamical or constitutive relations satisfied by the dissipative tensors Π, n µ and π µν .

Relativistic Navier-Stokes theory
In the presence of dissipative currents, the entropy is no longer a conserved quantity, i.e., ∂ µ S µ = 0. Since the form of the entropy four-current for a dissipative fluid is not known a priori, it is not trivial to obtain its equation. We proceed by recalling the form of the entropy fourcurrent for ideal fluids, Eq. (27), and extending it for dissipative fluids, The above extension remains valid because an artificial equilibrium state was constructed using the matching conditions to satisfy the thermodynamic relations as if in equilibrium. This was the key step proposed by Eckart, Landau and Lifshitz in order to derive the relativistic Navier-Stokes theory [35,36]. The next step is to calculate the entropy generation, ∂ µ S µ , in dissipative fluids.
To this end, we substitute the form of T µν and N µ for dissipative fluids from Eq. (51) in Eq. (55). Taking the divergence and using Eqs. (52)-(54), we obtain The relativistic Navier-Stokes theory can then be obtained by applying the second law of thermodynamics to each fluid element, i.e., by requiring that the entropy production ∂ µ S µ must always be positive, The above inequality can be satisfied for all possible fluid configurations if one assumes that the bulk viscous pressure Π, the particle-diffusion four-current n µ , and the shear-stress tensor π µν are linearly proportional to θ, ∇ µ α, and σ µν , respectively. This leads to where the proportionality coefficients ζ, κ and η refer to the bulk viscosity, the particle diffusion, and the shear viscosity, respectively. Substituting the above equation in Eq. (56), we observe that the source term for entropy production becomes a quadratic function of the dissipative currents In the above equation, since n µ is orthogonal to the timelike four-vector u µ , it is spacelike and hence n µ n µ < 0. Moreover, π µν is symmetric in its Lorentz indices, and in the local rest frame π 0µ = π µ0 = 0. Since the trace of the square of a symmetric matrix is always positive, therefore π µν π µν > 0. Hence, as long as ζ, κ, η ≥ 0, the entropy production is always positive. Constitutive relations for the dissipative quantities, Eq. (58), along with Eqs. (52)- (54) are known as the relativistic Navier-Stokes equations.
The relativistic Navier-Stokes theory in this form was obtained originally by Landau and Lifshitz [36]. A similar theory was derived independently by Eckart, using a different definition of the fluid four-velocity [35]. However, as already mentioned, the Navier-Stokes theory is acausal and, consequently, unstable. The source of the acausality can be understood from the constitutive relations satisfied by the dissipative currents, Eq. (58). The linear relations between dissipative currents and gradients of the primary fluid-dynamical variables imply that any inhomogeneity of α and u µ , immediately results in dissipative currents. This instantaneous effect is not allowed in a relativistic theory which eventually causes the theory to be unstable. Several theories have been developed to incorporate dissipative effects in fluid dynamics without violating causality: Grad-Israel-Stewart theory [40,46,49], the divergence-type theory [50,51], extended irreversible thermodynamics [52], Carter's theory [53],Öttinger-Grmela theory [54], among others. Israel and Stewart's formulation of causal relativistic dissipative fluid dynamics is the most popular and widely used; in the following we briefly review their approach.

Causal fluid dynamics: Israel-Stewart theory
The main idea behind the Israel-Stewart formulation was to apply the second law of thermodynamics to a more general expression of the non-equilibrium entropy fourcurrent [40,46,49]. In equilibrium, the entropy fourcurrent was expressed exactly in terms of the primary fluid-dynamical variables, Eq. (27). Strictly speaking, the nonequilibrium entropy four-current should depend on a larger number of independent dynamical variables, and, a direct extension of Eq. (27) to Eq. (55) is, in fact, incomplete. A more realistic description of the entropy four-current can be obtained by considering it to be a function not only of the primary fluid-dynamical variables, but also of the dissipative currents. The most general off-equilibrium entropy four-current is then given by where Q µ is a function of deviations from local equilib- . Using Eq. (51) and Taylor-expanding Q µ to second order in dissipative fluxes, we obtain where O(δ 3 ) denotes third order terms in the dissipative currents and β 0 , β 1 , β 2 , α 0 , α 1 are the thermodynamic coefficients of the Taylor expansion and are complicated functions of the temperature and chemical potential. We observe that the existence of second-order contributions to the entropy four-current in Eq. (61) should lead to constitutive relations for the dissipative quantities which are different from relativistic Navier-Stokes theory obtained previously by employing the second law of thermodynamics. The relativistic Navier-Stokes theory can then be understood to be valid only up to first order in the dissipative currents (hence also called first-order theory). Next, we re-calculate the entropy production, ∂ µ S µ , using the more general entropy four-current given in Eq. (61), As argued before, the only way to explicitly satisfy the second law of thermodynamics is to ensure that the entropy production is a positive definite quadratic function of the dissipative currents. The second law of thermodynamics, ∂ µ S µ ≥ 0, is guaranteed to be satisfied if we impose linear relationships between thermodynamical fluxes and extended thermodynamic forces, leading to the following evolution equations for bulk pressure, particle-diffusion four-current and shear stress tensor, where λ ≡ κ/T . This implies that the dissipative currents must satisfy the dynamical equations, The above equations for the dissipative quantities are relaxation-type equations with the relaxation times defined as Since the relaxation times must be positive, the Taylor expansion coefficients β 0 , β 1 and β 2 must all be larger than zero.
The most important feature of the Israel-Stewart theory is the presence of relaxation times corresponding to the dissipative currents. These relaxation times indicate the time scales within which the dissipative currents react to hydrodynamic gradients, in contrast to the relativistic Navier-Stokes theory where this process occurs instantaneously. The introduction of such relaxation processes restores causality and transforms the dissipative currents into independent dynamical variables that satisfy partial differential equations instead of constitutive relations. However, it is important to note that this welcome feature comes with a price: five new parameters, β 0 , β 1 , β 2 , α 0 and α 1 , are introduced in the theory. These coefficients cannot be determined within the present framework, i.e., within the framework of thermodynamics alone, and as a consequence the evolution equations remain incomplete. Microscopic theories, such as kinetic theory, have to be invoked in order to determine these coefficients. In the next section, we review the basics of relativistic kinetic theory and Boltzmann transport equation, and discuss the details of the coarse graining procedure to obtain dissipative hydrodynamic equations.

E. Relativistic kinetic theory
Macroscopic properties of a many-body system are governed by the interactions among its constituent particles and the external constraints on the system. Kinetic theory presents a statistical framework in which the macroscopic quantities are expressed in terms of singleparticle phase-space distribution function. The various formulations of relativistic dissipative hydrodynamics, presented in this review, are obtained within the framework of relativistic kinetic theory. In the following, we briefly outline the salient features of relativistic kinetic theory and dissipative hydrodynamics which have been employed in the subsequent calculations [55].
Let us consider a system of relativistic particles, each having rest mass m, momentum p and energy p 0 . Therefore from relativity, we have, p 0 = ( p) 2 + m 2 . For a large number of particles, we introduce a single-particle distribution function f (x, p) which gives the distribution of the four-momentum p = p µ = (p 0 , p) at each space-time point such that f (x, p)∆ 3 x∆ 3 p gives the average number of particles at a given time t in the volume element ∆ 3 x at point x with momenta in the range ( p, p + ∆ p). However, this definition of the single-particle phase-space distribution function f (x, p) assumes that, while on one hand, the number of particles contained in ∆ 3 x is large, on the other hand, ∆ 3 x is small compared to macroscopic point of view.
The particle density n(x) is introduced to describe, in general, a non-uniform system, such that n(x)∆ 3 x is the average number of particles in volume ∆ 3 x at ( x, t). Similarly, particle flow j(x) is defined as the particle current. With the help of the distribution function, the particle density and particle flow are given by where v = p/p 0 is the particle velocity. These two local quantities, particle density and particle flow constitute a four-vector field N µ = (n, j), called particle four-flow, and can be written in a unified way as Note that since d 3 p/p 0 is a Lorentz invariant quantity, f (x, p) should be a scalar in order that N µ transforms as a four-vector. Since the energy per particle is p 0 , the average energy density and the energy flow can be written in terms of the distribution function as The momentum density is defined as the average value of particle momenta p i , and, the momentum flow or pressure tensor is defined as the flow in direction j of momentum in direction i. For these two quantities, we have Combining all these in a compact covariant form using v i = p i /p 0 , we obtain the energy-momentum tensor of a macroscopic system Observe that the above definition of the energy momentum tensor corresponds to second moment of the distribution function, and hence, it is a symmetric quantity.
The H-function introduced by Boltzmann implies that the nonequilibrium local entropy density of a system can be written as The entropy flow corresponding to the above entropy density is These two local quantities, entropy density and entropy flow constitute a four-vector field S µ = (s, S), called entropy four-flow, and can be written in a unified way as The above definition of entropy four-current is valid for a system comprised of Maxwell-Boltzmann gas. This expression can also be extended to a system consisting of particles obeying Fermi-Dirac statistics (r = 1), or Bose-Einstein statistics (r = −1) as wheref ≡ 1 − rf . The expressions for the entropy fourcurrent given in Eqs. (77) and (78) can be used to formulate the generalized second law of thermodynamics (entropy law), and, define thermodynamic equilibrium.
For small departures from equilibrium, f (x, p) can be written as f = f 0 + δf . The equilibrium distribution function f 0 is defined as where the scalar product is defined as u · p ≡ u µ p µ and r = 0 for Maxwell-Boltzmann statistics. Note that in equilibrium, i.e., for f (x, p) = f 0 (x, p), the particle fourflow and energy momentum tensor given in Eqs. (71) and (74) reduce to that of ideal hydrodynamics N µ (0) and T µν (0) . Therefore using Eq. (51), the dissipative quantities, viz., the bulk viscous pressure Π, the particle diffusion current n µ , and the shear stress tensor π µν can be written as The evolution equations for the dissipative quantities expressed in terms of the non-equilibrium distribution function, Eqs. (80)- (82), can be obtained provided the evolution of distribution function is specified from some microscopic considerations. Boltzmann equation governs the evolution of the phase-space distribution function which provides a reliably accurate description of the microscopic dynamics. Relativistic Boltzmann equation can be written as where dp ≡ d 3 p/p 0 and C[f ] is the collision functional. For microscopic interactions restricted to 2 ↔ 2 elastic collisions, the form of the collision functional is given by where W pp →kk is the collisional transition rate. The first and second terms within the integral of Eq. (84) refer to the processes kk → pp and pp → kk , respectively. In the relaxation-time approximation, where it is assumed that the effect of the collisions is to restore the distribution function to its local equilibrium value exponentially, the collision integral reduces to [56] where τ R is the relaxation time.

F. Dissipative fluid dynamics from kinetic theory
The derivation of a causal theory of relativistic dissipative hydrodynamics by Israel and Stewart [40] proceeds by invoking the second law of thermodynamics, viz., ∂ µ S µ ≥ 0, from the algebraic form of the entropy four-current given in Eq. (61). As noted earlier, the new parameters, β 0 , β 1 , β 2 , α 0 and α 1 , cannot be determined within the framework of thermodynamics alone and microscopic theories, such as kinetic theory, have to be invoked in order to determine these coefficients. On the other hand, one may demand the second law of thermodynamics from the definition of the entropy four-current, given in Eqs. (77) and (78), in order to obtain the dissipative equations [57]. This essentially ensures that the nonequilibrium corrections to the distribution function, δf , does not violate the second law of thermodynamics. In Ref. [57], the generalized method of moments developed by Denicol et al. [58] was used to quantify the dissipative corrections to the distribution function. The form of the resultant dissipative equations, obtained in Ref. [57], are identical to Eqs. (66)- (68), with the welcome exception that all the transport coefficients are now determined in terms of the thermodynamical quantities.
The moment method, originally proposed by Grad [49], has been used quite extensively to quantify the dissipative corrections to the distribution function [57][58][59][60][61][62][63][64][65]. In this method, the distribution function is Taylor expanded in powers of four-momenta around its local equilibrium value. Truncating the Taylor expansion at second-order in momenta results in 14 unknowns that have to be determined to describe the distribution function. This expansion implicitly assumes a converging series in powers of momenta. An alternative derivation of causal dissipative equations, which do not make use of the moment method, was proposed in Ref. [66]. In this method, which is based on a Chapman-Enskog like expansion, the Boltzmann equation in the relaxation time approximation is solved iteratively to obtain δf up to any arbitrary order in derivatives. To first and second-order in gradients, one obtains This method of obtaining the form of the nonequilibrium distribution function is consistent with dissipative hydrodynamics, which is also formulated as a gradient expansion.
Since it is well established that QGP formed in high energy heavy-ion collisions is strongly coupled, it is of interest to compare the transport coefficients obtained from kinetic theory with that of a strongly coupled system [86]. In contrast to kinetic theory, strongly coupled quantum systems, in general, does not allow for a quasiparticle interpretation. This can be attributed to the fact that the quasiparticle notion hinges on the presence of a well-defined peak in the spectral density, which may not exist at strong coupling. Therefore it is interesting to study the hydrodynamic limit of an infinitely strongly coupled system, which are different than sys-τπT /(η/s) λ1T /(η/s) λ2T /(η/s) λ3 ADS/CFT 2(2 − ln 2) 2η 4η ln 2 0 KT 5 (25/7)η 10η 0 tems described by kinetic theory. In the following we discuss the evolution equation for shear stress tensor for a strongly coupled conformal system which is equivalent to a system of massless particles in kinetic theory. For such a system, the evolution of shear stress tensor is governed by the equation For a system of massless particles, the kinetic theory results for second-order transport coefficients agree in the case of both moment method [58] as well as the the Chapman-Enskog like iterative solution of the Boltzmann equation [66]. In Table I, we compare the transport coefficients obtained from kinetic theory and from calculations employing ADS/CFT correspondence of strongly coupled N = 4 SYM and its supergravity dual [87,88]. We see that the second-order transport coefficients obtained from Kinetic theory are, in general, larger than those obtained from the ADS/CFT calculations.

III. INITIAL CONDITIONS
In order to apply hydrodynamics to study the collective phenomena observed at relativistic heavy-ion collisions, one needs to first characterize the system. To this end, we shall discuss here, and in the next few sections about initial conditions, Equation of State (EoS), and freeze-out procedure as used in state of the art relativistic hydrodynamics simulations. However, we note that the following discussions are in no way complete but we will try to provide appropriate references wherever possible. Most of the following discussions can be also found in more details in Ref. [32,[89][90][91][92][93][94].
In high energy heavy-ion collisions, bunch of nucleus of heavy elements are accelerated inside the beam pipes and in the final state (after the collisions) we have hundred or thousands of newly created particles coming out from the collision point in all directions. The underlying processes of the collisions between the constituent partons of the colliding nucleus and the conversion of initial momentum along the beam direction to the (almost) isotropic particle production is still not very well understood. Particularly the state just after the collisions when the longitudinal momentum distribution of the partons started to become isotropic and subsequently achieve the local thermal equilibrium state is poorly understood. But the precise knowledge of this so called pre-equilibrium stage is essential input in the viscous hydrodynamics models. The knowledge of distribution of energy/entropy density and the thermalisation time is one of the uncertainty present in the current hydrodynamics model studies. Below we discuss four most popular initial condition models used in hydrodynamics simulation of heavy-ion collisions.

A. Glauber model
The Glauber model of nuclear collisions is based on the original idea of Roy J. Glauber to describe the quantum mechanical scattering of proton-nucleus and nucleusnucleus collisions at low energies. The original idea of Glauber was further modified by Bialas et. al. [95] to explain inelastic nuclear collisions. For a nice and more complete review of Glauber model see Ref. [96] and references therein for more details. We discuss below the very essential part of this model as used in heavy-ion collisions, particularly in the context of relativistic hydrodynamics.
At present, there are two main variation of Glauber model in use. One of them is based on the optical limit approximation for nuclear scattering, where the nuclear scattering amplitude can be described by an eikonal approach. In this limit each of the colliding nucleons see a smooth density of nucleon distribution in the other nucleus. This variation of Glauber model, also known as optical Glauber model, uses the Wood-Saxon nuclear density distribution for a nucleus with mass number A as where R 0 , a 0 are the nuclear radius and skin thickness parameter of the nucleus, and ρ 0 is an overall constant that is determined by requiring d 3 xρ A (x, y, z) = A. One additionally defines the "thickness function" [96] T which indicates the Lorentz contraction in the laboratory frame. The Wood-Saxon nuclear density distribution is used along with the experimentally measured inelastic nuclear cross section to calculate the number of participating nucleons (N part ) and number of binary collisions (N coll ) for the two colliding nucleus. In order to calculate the N coll (x, y, b) and N part (x, y, b) from Glauber model, one can choose the X axis along the impact parameter vector b. The N part and N coll distribution are functions of impact parameter, the inelastic Nucleon-Nucleon cross-section and the nuclear density distribution function. For a collision of two spherical nuclei with different mass number 'A' and 'B', the transverse density of binary collision and wounded nucleon profile is given by [96] where Here σ in is the inelastic nucleon-nucleon cross section whose value depends on the √ s N N and is obtained from the experimental data. The distribution of N part and N coll in the transverse plane as obtained from Glauber model is used to calculate the initial energy/entropy density for the hydrodynamics simulation. The exact form for calculation of the energy density in the transverse plane using optical Glauber model is given by where 0 is a multiplicative constant used to fix the charged hadron multiplicity, α is the fraction of hard scattering [97]. The energy density corresponds to the MC-Glauber model is obtained with similar contribution from number of binary collisions and number of participants.
In the second variation, the distribution of nucleons inside the colliding nucleus are sampled according to the nuclear density distribution by using statistical Monte-Carlo (MC) method. The collisions between two nucleons occurs when the distance between them becomes equal or smaller than the radius obtained from the inelastic nucleon-nucleon cross section. This is also known as MC-Glauber model.
In MC-Glauber model the positions of binary collisions and participating nucleons are random and they are delta function in configuration space. These delta functions cannot be used in the numerical simulation of hydrodynamics. The usual practice is to use two-dimensional Gaussian profile to make a smooth profile of initial energy density as given by [98] where σ is a free parameter controlling the width of the Gaussian. Typical values for this fluctuation size parameter are of the order of 0.5 fm, WN is the abbreviation for wounded nucleons which is same as number of participant N part and BC represent number of binary collisions. With this short discussion we now move on to the next topic.

B. Color-Glass-Condensate
The Color-Glass-Condensate (CGC) model takes into account the non-linear nature of the QCD interactions. Due to Lorentz contraction at relativistic energies, the nucleus in the laboratory frame is contracted into a sheet and therefore one only needs to consider the transverse plane. The density of partons inside such a highly Lorentz contracted nucleus is dominated by gluons. According to the uncertainty principle, the radius, r gl , of a gluon is related to its momentum, Q, via |r gl |×|Q| ∼ = 1. Therefore the cross-section of gluon-gluon interactions is where α s is the strong coupling constant. The total number of gluons in a nucleus can be considered to be approximately proportional to the number of partons, and therefore also to its mass number A. The density of gluons in the transverse plane is then given by A/(πR 2 0 ), where R 0 is the radius of the nucleus. Gluons starts interacting with each other when the scattering probability becomes of the order unity, This indicates that there exists a typical momentum scale Q 2 s = α s A/R 2 0 separating perturbative (Q 2 Q 2 s ) and non-perturbative (Q 2 Q 2 s ) regimes. Classical chromodynamics is a good approximation at low momenta due to the high occupation number ("saturation"). The CGC model was developed to incorporate the saturation physics at low momenta Q 2 in relativistic heavy-ion collisions [99,100].
The presence of non-abelian plasma instabilities [101][102][103] makes it difficult to determine the energy density distribution in the transverse plane. As a result, one has to resort to phenomenological models for the transverse energy density distribution in the CGC model [104] Here N g is the number of gluons produced in the collision whose momentum distribution is given by where x = p T / √ s. It is important to note that the Glauber and CGC models lead to different values of spatial eccentricity, defined by where represents averaging over the transverse plane with weight (x ⊥ , b). It has been observed that the CGC model typically has a larger eccentricity than the Glauber model which means that the anisotropy in fluid velocities is larger for the CGC model.
In a variant of the CGC model, also known as the Kharzeev-Levin-Nardi (KLN) model [105][106][107], the entropy production is determined by the initial gluon multiplicity. A Monte-Carlo version of KLN model (MC-KLN) has also been proposed to incorporate event-by-event fluctuations in the nucleon positions [108,109]. In these models, the initial gluon production is calculated using the perturbative merging of two gluons from the target and projectile nuclei. The gluon structure functions are parametrized by a position-dependent gluon saturation momentum, Q s , which is computed from the longitudinally projected density of wounded nucleons. The positions of the wounded nucleons are sampled according to Eq. (97) using the MC-Glauber model. However, one should keep in mind that the MC-Glauber and MC-KLN models are unable to account for fluctuations of the gluon fields inside the colliding nucleons.

IV. PRE-EQUILIBRIUM DYNAMICS
The initial condition models described in the previous section are static models because after the collisions the energy/entropy remains constant in space-time until the initial time τ 0 when the hydrodynamics evolution starts. More realistic condition should include dynamical evolution of the constituent partons in the preequilibrium phase. The simplest choice for the dynamical evolution is the free-streaming of the produced partons in the pre-equilibrium phase, but this is in contrary to the assumption of local thermal equilibrium which needs multiple collisions among the constituent to achieve the local thermal equilibrium. In the following, we describe a few state-of-the-art models which takes into account the pre-equilibrium dynamics, until hydrodynamics sets in.

A. IP-Glasma
In the IP-Glasma model, the initial conditions is determined within the CGC framework by combining the impact parameter dependent saturation model (IP-Sat). In addition to fluctuations of nucleon positions within a nucleus, the IP-Glasma description also incorporates quantum fluctuations of color charges on the length-scale determined by the inverse saturation scale, 1/Q s . The initial Glasma fields are then evolved using the classical Yang-Mills (CYM) equation. One of the most important feature of this model is that long-range rapidity correlations from the initial state wavefunctions are efficiently converted into hydrodynamic flow of the final state quark-gluon matter [110,111]. Moreover, initial energy fluctuations produced within this model naturally follows a negative binomial distribution.
The color charges, ρ a (x − , x ⊥ ), in the IP-Sat model behaves as local sources for small-x classical gluon Glasma fields. The classical gluon fields are then determined by solving the classical Yang-Mills equations, The color current in the above equation, generated by a nucleus A (B) moving along the x + (x − ) direction, is given by where the upper indices are for nucleus A.
It is easy to solve Eq. (108) in Lorentz gauge, ∂ µ A µ = 0, where the equation transforms into a two-dimensional Poisson equation The solution of the above equation can be written as Using the path-ordered exponential one can gauge transform the results of Lorentz gauge to light-cone gauge, A + (A − ) = 0. The pure gauge fields are then given by [112][113][114] The discontinuity in the fields on the light-cone corresponds to the localized valence charge source [115]. The initial condition for a heavy-ion collision, at time τ = 0, is determined by the solution of the CYM equations in Fock-Schwinger gauge A τ = (x + A − + x − A + )/τ = 0, where the τ, η coordinates are defined as τ = √ 2x + x − and η = 0.5 ln(x + /x − ). The Fock-Schwinger gauge is a natural gauge choice because it interpolates between the light-cone gauge conditions of the incoming nuclei. In terms of the gauge fields of the colliding nuclei, one obtains [115,116]: In the limit τ → 0, A η = −E η /2, where E η is the longitudinal component of the electric field. At τ = 0, one can non-perturbatively calculate the longitudinal magnetic and electric fields, which are the only non-vanishing components of the field strength tensor. These fields determine the energy density of the Glasma at each transverse position in a single event [117][118][119]. The Glasma fields are then evolved in time numerically according to Eq. (108), up to a proper time τ switch , which is the switching time from classical Yang-Mills dynamics to hydrodynamics [120]. At the switching time, one can construct the fluid's initial energy momentum tensor T µν fluid = ( + P)u µ u ν − Pg µν + Π µν from the energy density in the fluid's rest frame ε and the flow velocity u µ . The local pressure P at each transverse position is obtained using an equation of state. The hydrodynamic quantities ε and u µ are obtained by solving the Landau frame condition, u µ T µν CYM = εu ν .

B. Transport: AMPT and UrQMD
In Refs. [121][122][123] a different approach was taken in order to incorporate the pre-equilibrium dynamics for obtaining the initial condition of hydrodynamics evolution. While the authors of Ref. [121] employ ultrarelativistic quantum molecular dynamics (UrQMD) string dynamics model, A Multi Phase Transport Model (AMPT) was used in Refs. [122,123] to simulate the pre-equilibrium dynamics. In these studies the partons produced in the collisions were evolved until the initial time τ 0 according to a simplified version of Boltzmann transport equation. We shall discuss here the particular procedure used in Ref. [122] for calculating initial conditions for a (3+1)D hydrodynamics evolution with the parton transport in the pre-equilibrium phase. Additional benefit for choosing this type of initial condition is that one naturally incorporate the fluctuating energy density in the longitudinal direction due to the discrete nature of partons, details of which will be discussed in a later section.
In Ref. [122], A Multi Phase Transport Model (AMPT) [124,125] was used to obtain the local initial energymomentum tensor in each computational cell. The AMPT model uses the Heavy-Ion Jet INteraction Generator (HIJING) model [126,127] to generate initial partons from hard and semi-hard scatterings and excited strings from soft interactions. The number of excited strings in each event is equal to that of participant nucleons. The number of mini-jets per binary nucleon-nucleon collision follows a Poisson distribution with the average number given by the mini-jet cross section, which depends on both the colliding energy and the impact parameter via an impact-parameter dependent parton shadowing in a nucleus. The total energy-momentum density of parton depends on the number of participants, number of binary collisions, multiplicity of mini jets in each nucleon-nucleon collisions and the fragmentation of excited strings. HIJING uses MC-Glauber model to cal-culate number of participant and binary collisions with the Wood-Saxon nuclear density distribution function.
After the production of partons from hard collisions and from the melting strings, they are evolved within a parton cascade model, where only two parton collisions are considered. The positions and momentum of each partons are then recorded and used to calculate the initial energy-momentum tensor using a Gaussian smearing at time τ 0 as where p τ i = m iT cosh (Y i − η is ), and p x,y i = p x,yi , p η i = m iT sinh (Y i − η is ) /τ 0 are the four-momenta of the ith parton and Y i , η is , and m iT are the momentum rapidity, the spatial rapidity, and the transverse mass of the ith parton, respectively. Unless otherwise stated, the smearing parameters are taken as: σ r = 0.6 fm and σ ηs = 0.6 from Refs. [122] where the soft hadron spectra, rapidity distribution and elliptic flow can be well described. The sum index i runs over all produced partons (N ) in a given nucleus-nucleus collision. The scale factor K and the initial proper time τ 0 are the two free parameters that we adjust to reproduce the experimental measurements of hadron spectra for central Pb+Pb collisions at mid-rapidity [122]. The initial energy density and the local fluid velocity in each cell is obtained from the calculated T µν via a root finding method which is used as an input to the subsequent hydrodynamics evolution, see Ref. [122] for further details.

C. Numerical relativity: AdS/CFT
Another method to simulate the pre-equilibrium stage is via numerical relativity solutions to AdS/CFT [128]. In this method, one employs the dynamics of the energymomentum tensor of the strongly coupled Conformal Field Theory (CFT) on the boundary using the gravitational field in the bulk of AdS 5 . Therefore a relativistic nucleus may be described using a gravitational shockwave in AdS, whereby the energy-momentum tensor of a nucleus can be exactly matched [129]. For a central collision, the dynamics of the colliding shockwaves has been solved near the boundary of AdS in Ref. [130], resulting in the energy-momentum tensor at early times. The starting point of this simulation is the energy density of a highly boosted and Lorentz contracted nucleus, T tt = δ(t + z)T A (x, y). Here the thickness function, T A (x, y), is the same as defined in Eq. (95) but with an extra normalization, 0 , which is used to match the experimentally observed particle multiplicity, dN/dY .
In terms of the polar Milne coordinates τ, ξ, ρ, θ with t = τ cosh ξ, z = τ sinh ξ, ρ 2 = x 2 + y 2 , tan θ = y/x, the energy density, fluid velocity and pressure anisotropy was found up to leading order in t [130] where in the local rest frame T µ ν = diag(− , P T , P T , P L ) [131][132][133][134]. One finds that the corresponding line-element ds 2 turns out to be ξ-independent (boost-invariant), up to leading order in τ , and can be written as Here all functions depend on τ , ρ and the fifth AdS space dimension r only. In this scenario, the space boundary is located at r → ∞ where the induced metric is given by g µν = diag(g τ τ , g ρρ , g θθ , g ξξ ) = diag(−1, 1, ρ 2 , τ 2 ). The metric is then expanded near the boundary, where B 0 is given by the vacuum value. In order to have a stable time evolution, a function with one bulk parameter σ has been introduced to extend the metric functions to arbitrary r. An analogous expansion is also made for C. Using Eq. (120) to fix the near-boundary coefficients at a time τ init , and choosing a value for σ, the time evolution of the metric can be determined by solving the Einstein equations. This is done numerically by adopting a pseudo-spectral method based on Refs. [135][136][137]. At a proper time τ hydro , which is the switching time from AdS/CFT to hydrodynamics, the evolution using Einstein equations is stopped and hydrodynamic quantities such as , u µ , π µν are extracted from the metric using Eq. (120). These quantities are then used to create the energy-momentum tensor which provides the initial conditions for the subsequent relativistic viscous hydrodynamic evolution. The initial conditions for hydrodynamic evolution is therefore determined using an early-time, far-from-equilibrium dynamics, modeled as a strongly coupled CFT described by gravity in AdS. Recently, a non-conformal extension has also been studied in order to incorporate bulk viscosity [138].

V. EQUATION OF STATE
Equation of State (EoS) is the functional relationship between thermodynamic variables pressure (P) and number density (n) to the energy density ( ). The conservation equations, ∂ µ T νµ = 0, contains one additional variable than the number of equations. EoS closes the system of equations by providing another functional relationship and it is one of the important input to hydrodynamics. For a relativistic simple fluid the acceleration under a given pressure gradient ∇P is governed by the following relationship where D = u µ ∂ µ is the covariant derivative, and P are the energy density and pressure respectively. Clearly the fluid expansion is governed by the gradient of pressure as well as the combined value of pressure and energy density. The pressure for a given energy density is defined via the EoS and hence the EoS governs the rate of change of fluid expansion. At present, the most reliable calculation of EoS for nuclear matter at high temperature (> 100 MeV) is obtained from lattice QCD (lQCD) calculations. However, at present the lQCD calculations are not reliable at lower temperatures (because of the large grid size needed at lower temperatures) and at higher baryon densities (due to the so called sign problem for finite chemical potential). The usual practice in the heavy-ion community is to use lQCD calculation at high temperature and a hadron resonance gas (HRG) model at lower temperature to construct the equation of state for vanishing baryon chemical potential (µ b ). The EoS for finite µ b is usually obtained by employing some approximation such as Taylor series expansion around µ b = 0. For more details about the nuclear EoS relevant to the heavy-ion collisions see Ref. [139] and references therein. Here we briefly outline the procedure used to calculate the lQCD+HRG equation of state for vanishing baryon chemical potential.
Usual lQCD calculations for the thermodynamical variables assume that the system has infinite extent (volume V → ∞) and it is homogeneous [140]. All thermodynamic quantities can be derived from the partition function Z(T, V ). The energy density and pressure are derivatives of the partition function with respect to T and V respectively The pressure for a homogeneous system of infinite extent can be simply expressed in terms of f as Using the above relations one can arrive at the following relationships where the trace anomaly, Θ(T ) = ( − 3P )/T 4 . In the high temperature, Θ(T ) can be reliably calculated from lQCD. On the other hand, at lower temperature the lQCD results are affected by possibly large discretisation effect. Therefore the usual practice to construct realistic EoS is to use the lattice data for the trace anomaly in the high temperature region (T > 250M eV ), and use HRG model in the low temperature region (T < 180M eV ). In the intermediate temperature Θ(T ) is obtained by joining the parametrised high temperature and low temperature values smoothly (continuous first and second derivative). Once Θ(T ) is known, pressure can be calculated by using Eq. (128). The energy density then can be readily obtained from Eq. (127). The top panel of Fig. 3 shows the trace anomaly calculated in lattice QCD with p4 and asqtad actions on N τ = 6 and 8 lattices [141] compared with various parametrizations given by the solid, dotted and dashed lines [139]. The bottom panel shows the trace anomaly calculated in lattice QCD compared with the HRG model given by the solid and dashed lines [139].

VI. FREEZE-OUT: SPECTRA AND FLOW
In the late stage of hydrodynamics evolution of hot and dense nuclear matter created in high energy heavyion collisions, the density and the temperature reaches a critical value when the constituents no longer collides among themselves and thereafter they move in a straight trajectory towards the detectors. This phenomenon is known as freeze-out, more precisely the kinetic freezeout. There is another chemical freeze-out when the particle number changing processes ceases. The chemical freeze-out temperature so far known to be higher than the kinetic freeze-out temperature.
In relativistic hydrodynamics simulations one also needs to stop the hydrodynamics evolution when the system reaches the kinetic freeze-out criterion. For this one needs to impose some physical constraints to calculate the freeze-out hyper-surface. This can be done by more than one way, we discuss here only the most popular choices used to calculate the freeze-out hyper-surface namely (i) the constant temperature freeze-out (ii) the constant energy density freeze-out (iii) the dynamical freeze-out. Among the three choices, the first two are based on the general idea that the pion (or other hadron) cross section is very sensitive to the temperature/energy density of the system, thus within short interval of temperature/energy density the condition for kinetic freezeout is achieved. For the computational purpose this is realised by choosing a constant temperature/energy density surface.
The dynamical freeze-out is based on the idea that the ratio of expansion rate (θ) to the collision rate (Γ) should be much much less than unity ( θ Γ 1) in order to maintain the local thermal equilibrium essential for the applicability of the hydrodynamics evolution. One can then define the freeze-out criterion based on some predefined value of θ Γ smaller than 1. Though the idea of dynamical freeze-out sounds more realistic it is not easy to implement in numerical calculation see Ref. [139] for details.
The thermodynamical quantities of the fluid such as energy density, pressure and the fluid velocity obtained from the hydrodynamics simulation on the freeze-out surface are used to evaluate momentum distributions of the identified hadrons. The conversion of fluid to hadrons is done by using the Cooper-Frye procedure, see Ref. [142] for details. In Cooper-Frye procedure the momentum distribution (or invariant yield) of hadrons are calculated as [142] where E, N, and p µ are energy, number and fourmomentum of hadrons, dΣ µ is the differential freezeout hyper-surface element. The distribution function, f (x, p), consists of an equilibrium part, f 0 (x, p), and dissipative corrections, δf (x, p). While the equilibrium distribution corresponding to the local thermodynamic quantities, as given in Eq. (79), is taken to be either Bose-Einstein or Fermi-Dirac distribution depending on the spin of the hadronic species, the dissipative correction is not unique, and will be explained in the following. In the simple case when the dissipation is only due to the shear viscosity, leading-order moment method, also known as the Grad's 14-moment approximation, leads to the well-known form of the viscous correction [40,49] δf (x, p) = f 0f0 2( + P )T 2 p α p β π αβ , wheref 0 ≡ 1 − rf 0 , with r = 1, −1, 0 for Fermi, Bose, and Boltzmann gases, respectively. Note that the viscous correction in this case increases with quadratic power of momenta. On the other hand, the Chapman-Enskog like iterative solution of the Boltzmann equation, Eq. (87), leads to a viscous correction which is effectively linear in momenta [69], It has been shown that, in contrast to Eq. (130) obtained using moment method, Eq. (131) leads to phenomenologically consistent corrections to the equilibrium distribution function, and is therefore a better alternative for hydrodynamic modeling of relativistic heavy-ion collisions [69]. We note here that the calculation of four-dimensional freeze-out hyper-surface and the numerical evaluation of it is not trivial, for example see Ref. [143] for more details. Once we know the invariant momentum distribution the "n-th" order Fourier coefficient the flow harmonics v n can be readily obtained as These above mentioned quantities are directly compared to the corresponding experimental data in order to obtain information about the transport coefficients such as shear and bulk viscosity of the QGP.

VII. RESONANCE DECAY AND HADRONIC RESCATTERING
In high energy nuclear collisions various hadronic resonances are formed. The life time of most of the resonance particles are of the order of the expansion life time of the nuclear matter. The end product for the most of the decay channels involve pions. The decay of hadron resonances to pion enhances the pion yield specially at low transverse momentum, p T . One can use the formalism given in [144] to calculate the relative contribution of the resonance decay to thermal pion spectra.
According to the formalism given in [144], to calculate the pion contribution from resonances, one need to provide the source temperature. The parametric fit to the ratio of the total pion to the thermal pion for the calculation at two different freeze-out temperature T fo = 130 MeV and 150 MeV are approximately given by [145] π ± total π ± thermal T fo =130MeV = 1.0121 + 1.4028 π ± total π ± thermal T fo =150MeV = 1.0252 + 3.0495 where m π = 139 MeV is the pion mass. Note that about ∼ 50% of the total pion yield come from resonance decay at LHC energy ( √ s N N = 2.76 TeV), whereas for RHIC energy ( √ s N N = 200 GeV) the resonance contribution to total pion yield is ∼ 30% for T fo = 130 MeV. The sudden conversion of fluid to non-interacting hadrons at the freeze-out hyper-surface in the fluid dynamical evolution is hard to happen in practice. In reality the hydrodynamical picture should work fine for the early hot and dense phase of the QGP evolution when the scattering rate is comparatively large compared to the expansion rate. As the system grows in size and cools down with time the scattering rate goes down compared to the expansion rate. At some point of space-time, particularly in the late hadronic phase it is expected that the dynamical evolution most probably be governed by the microscopic Boltzmann equations considering multiple hadronic species and their collisions rather than the simplified macroscopic hydrodynamics evolution. Thus a complete dynamical evolution of high energy heavyion collisions contains simpler hydrodynamics evolution in the early time and a much computational expensive hadronic transport evolution in the late stage with the additional complexity of transforming fluid variables to position and momentum of hadrons.
For the hadronic rescattering phase several microscopic algorithms that solve coupled Boltzmann equations for a hadronic gas were developed in the 1980s and 1990s [21,124,[146][147][148][149][150]. Hybrid codes that coupled an ideal fluid dynamical description of an expanding QGP to hadronic rescattering codes and compared the results with purely fluid dynamical calculations began to appear around 2000 [151][152][153][154]. One of the first numerical code VISHNU that couples (2+1)D viscous hydro with a late hadronic Boltzmann cascade appeared in 2011 [155]. The use of these more sophisticated hybrid models are believed to reduce the uncertainty in the extracted value of η/s of QGP, since the late hadronic stage is known to have larger shear viscosity which in usual viscous hydrodynamics simulations is not taken into account properly. We shall not go into the details of the hadronic transport model nor to the technical details of various techniques and uncertainties arising due to the matching of viscous hydrodynamics to the hadronic transport, details of which can be found in Ref. [155] and references therein. Before finishing this section we should point out one of the major findings of Ref. [155], the η/s of hadronic matter is found to be quite sensitive to the details of preceding hydrodynamics phase and on the switching temperature when the viscous hydrodynamics is switched to the hadronic transport evolution. The effort to better constraint the η/s of QGP by using such sophisticated numerical models is a topic of current ongoing research.

VIII. TRANSPORT COEFFICIENTS
Determination of transport coefficients of the hot and dense QCD matter is one of the primary goal of theoretical simulations of relativistic heavy-ion collisions; see [156] for a recent review. Ideal hydrodynamics has been proved to be quite successfully in the past to describe the spectra of produced particles in relativistic heavy-ion collisions. The presence of dissipation leads to dissipative entropy generation via Eq. (62), which results in the increase of total particle multiplicity for a fixed initial entropy. Shear viscosity, in particular, also leads to stronger radial flow leading to an increase in the mean transverse momentum of particles. However, the most important effect of shear viscosity is to suppress the elliptic flow coefficient, v 2 , defined in Eq. (132) strongly. Therefore, in order to estimate the viscosity of the QCD matter within a hydrodynamic simulation, one has to tune the value of the specific shear viscosity, η/s, in order to fit the experimental data for v 2 . One of the first estimates of η/s was made within a hydrodynamics inspired blast-wave model [157]. Since then there has been a lot of activity in this field, which is briefly reviewed in the following.  TeV [167][168][169][170]. All the model calculations indicates that the values of η/s of the QCD matter formed in heavy-ion collision at LHC lies between 1-4 ×(1/4π). The specific shear viscosity was obtained in reference [167] by using A Multi Phase Transport model (AMPT). Bozek [168] has estimated the specific shear viscosity of the fluid for LHC energy by using a (2+1)D viscous hydrodynamics model. In addition to shear viscosity, bulk viscosity (ζ/s = 0.04) in the hadronic phase was considered. Freeze-out and resonance decay was based on THERMINATOR event generator [171]. Experimental data are best fitted with η/s ∼ 0.08. A (3+1)D viscous hydrodynamics calculation with fluctuating initial conditions was done by Schenke et al. [169]. They explain the v 2 (p T ) and p T integrated v 2 for different centralities. Their calculation shows that the experimental data measured at LHC by the ALICE collaboration are best described for η/s value 0.08 or smaller. Luzum et al. [170] have estimated η/s by using a (2+1)D viscous hydrodynamics simulation with smooth initial conditions for LHC energy to be same as at RHIC, η/s = 0.1 ± 0.1(theory) ± 0.08(experiment). Comparison of experimentally measured integrated and differential v 2 , the charged hadron p T spectra and multiplicity in the mid-rapidity and their global fit by minimising χ 2 was done in Ref. [172] by using a (2+1)D viscous hydrodynamics simulation, the extracted value of η/s is ∼ 0.07 ± 0.01. The effects of bulk viscosity in hydrodynamic simulations of relativistic heavy-ion collisions have not been investigated as thoroughly as that of shear viscosity. In principle, the bulk viscosity of the QCD matter should not be zero for the range of temperatures achieved at the RHIC and the LHC, and it may become large enough to significantly affect the evolution of the medium [173,174]. There has been several simulations of heavy-ion collisions that include the effect of bulk viscosity where it has been demonstrated that bulk viscosity can have a non-negligible effect on heavy-ion observables [175][176][177][178][179][180]. However, there are various uncertainties in the extraction of bulk viscosity from the anisotropic flow data of heavy-ion collisions. For example, the theoretical uncertainties arising due to the ambiguities in the form of the specific bulk viscosity, ζ/s, its relaxation time and the bulk viscous corrections to the freeze-out process, makes it difficult to study the effect of bulk viscosity on the evolution of QCD matter. Unlike shear viscosity, the extraction of bulk viscosity from hydrodynamic simulations is still unresolved and is currently an active research area.

IX. RECENT DEVELOPMENTS
A. Flow in small systems: proton-proton and proton-nucleus collisions As mentioned in the introduction, among the recent developments in the field of high energy heavy-ion collisions the most striking observation is the existence of radial flow like pattern in high multiplicity proton-proton (p-p) and proton-lead (p-Pb)collisions, for example see Ref. [181] for a summary of recent experimental results. At this point it needs some explanation why the observation of flow in small system is remarkable. One of the fundamental assumption when applying hydrodynamics to high energy nuclear collisions is that the system reaches a state of "local thermal equilibrium" very quickly because of the strong interactions among the quarks and gluons. In heavy-ion collisions such as Pb-Pb or Au-Au the number of participating nucleons are large, for example for a head-on Au-Au or Pb-Pb collisions there are 197+197 and 208+208 participating nucleons. Each of these col- 5: (Color online) Elliptic flow, v2, of identified particles from the hydrodynamic model compared to the ALICE data [186]. The figure is taken from Ref. [187].
liding nucleons on average produce more than one particles in each collisions, thus the total number of degrees of freedom in the system created just after the collisions are large. They collide among themselves through strong interaction and subsequently reaches local thermal equilibrium (we note here that the full mechanism by which the system reaches local thermal equilibrium within such a short period of time is not fully understood yet).
The situation is very different in smaller colliding systems such as p-p or p-Pb where the number of participating nucleons and the numbers of produced particles are comparatively smaller in numbers. It is very counterintuitive that with such few number of particles the system reaches local thermal equilibrium within a very short time-period. On the other hand, event generators based on perturbative QCD such as PYTHIA and HI-JING have successfully described various observables associated with particle production in p-p collisions. Thus p-p as well as p-nucleus collisions have been for long time considered qualitatively different from heavy-ion collisions, for which the hydrodynamic description became a mainstream since its successful explanation of RHIC data. It is interesting to note that the p-p collisions were used as a benchmark for studying the existence of QGP in larger systems where a thermalised medium is believed to be created.
This situation changed recently as the CMS and AT-LAS collaboration observed a "ridge" like correlation in the azimuthal distribution of charged hadrons produced in high multiplicity p-p or p-Pb collisions. In those experiments a mass dependence of the slope of identified hadron's m T spectrum in high multiplicity p-Pb and pp collisions were also observed. All these phenomenon are known to be the most significant indication for existence of hydrodynamic flow in larger colliding systems such as Au-Au or Pb-Pb. Like in heavy-ion collisions, the PYTHIA model failed to describe these observed experimental measurement for high multiplicity p-p or p-Pb events unless it employs some special mechanism like Color Reconnection (CR) and Multi Parton Interac- tion (MPI) with an additional free parameter to explain the experimental data [182]. On the other hand, the relativistic hydrodynamic models with large radial velocity have been proved to be quite successful in describing the same experimental data. It is also worthwhile to mention that there are some other theoretical conjectures about these recent observation which does not incorporate this hydrodynamics like flow, but till now those studies lack detailed numerical calculation in order to compare it with the experimental data Ref. [183]. TeV is similar to those in peripheral Pb-Pb collisions Ref. [184]. Considering the fact that final charged multiplicity is proportional to the initial energy/entropy density it is clear that the initial energy density in most violent p-Pb collisions is similar as in heavy-ion collisions. In fact the collision zone in p-Pb collisions is expected to be smaller than the peripheral Pb-Pb collisions, consequently the energy density is higher in high multiplicity p-Pb events than the peripheral Pb-Pb events. The initial high energy density within small volume in p-Pb collisions creates favourable condition for the subsequent hydrodynamics evolution. Another strong evidence of hydrodynamical flow in p-Pb collisions came from the observation of mass dependence of slope of identified hadrons p T spectra, measured in experiment Ref. [185]. The experimentally measured v n and v 2 (p T ) data for identified hadrons in p-Pb collisions by CMS [184] and ALICE [186] collaboration is nicely explained by a (3+1)D viscous hydrodynamics model study by Bozek et al. in Ref. [187] (see Fig. 5). In addition to Open symbols correspond to the CMS data [190] for |η| < 2.4 and √ s = 7 TeV, while the solid ones are obtained from the best one-parameter fit of the Gubser's flow. The figure is taken from Ref. [188].
that the mean transverse momentum of identified hadrons is also explained within the same (3+1)D hydrodynamic model, whereas, the Monte-Carlo event generator model HIJING which is based on the perturbative QCD processes relevant to the collisions fails to explain the same experimental data as can be seen in Fig. 6. This already gives the indication that for the high multiplicity p-Pb collisions QGP is produced and it flows like fluid before freezing out to hadrons.
• p-p collisions: Qualitatively similar signature of collective behaviour is also observed in high multiplicity p-p collisions like high multiplicity p-Pb collisions. However, the initial measurement shows that the p T integrated v 2 and v 3 is 30% and 50% smaller than in p-Pb at similar multiplicity. Like heavy-ion and p-Pb collisions a simple hydrodynamic inspired model with large radial velocity has successfully explained the experimental observation of mass dependence of slope in p-p collisions; see Fig. 7 which is taken from Ref. [188]. In a similar effort a blast wave model fit was shown to be inconsistent with the experimental data, see Ref. [189] for details. There are some studies of viscous hydrodynamics for p-p collisions for example see Ref. [191][192][193], more extensive study is needed for the comparison of all experimentally available data. We note here that it is still an open question whether the small system created in p-p collisions are big enough or live long enough for hydrodynamics to be applicable, detailed discussion of which is out of the scope of the present review. We refer to see the Ref. [194] for a detailed discussion about the applicability of hydrodynamics in small systems.  [195].

B. Flow in ultra central collisions
As mentioned earlier, the hydrodynamic response of the anisotropy in the initial overlap geometry in the configuration space transforms to the final momentum space anisotropy giving rise to non-zero values of flow harmonics v n . The most prominent flow harmonics v 2 originates as a hydrodynamic expansion of the initial elliptic shape of the fireball. The conversion efficiency of the spatial deformation into the momentum space anisotropy is very sensitive to the shear viscosity over entropy density (η/s) and the initial configuration of the system. The extraction of η/s of QGP by comparing hydrodynamic simulation results to the corresponding experimental data is riddled with large uncertainties in our understanding of the initial-state conditions of heavy-ion collisions. For example, viscous hydrodynamic simulation with MC-Glauber initial condition gives very different values of η/s compared to the same simulation with different initial condition such as MC-KLN. This uncertainty due to the poorly known initial condition can be minimised in case of ultra central collisions. In ultra-central collisions v 2 and other higher flow harmonics solely originate from the ini- For ultra central collisions the initial collision geometry is predominantly generated by fluctuations such that various orders of eccentricities predicted by different models tend to converge. Therefore, studies of v n in ultra-central heavy-ion collisions can help to reduce the systematic uncertainties of initial-state modelling in extracting the η/s value of the system. Let us first discuss the recent experimental results for ultra central Pb-Pb collisions at LHC, after that we shall also discuss the corresponding results from viscous hydrodynamics simulations. Top panel of Fig. 8 shows the experimentally measured differential flow coefficients v n as a function of p T for 0-0.2% centrality Pb-Pb collisions at √ s N N = 2.76 TeV. The v n 's were calculated using 2 particle correlation method with large pseudo-rapidity gap |∆η| > 2 between the two hadrons. The bottom panel of Fig. 8 shows the p T integrated v n (n=2-7) in ultra central Pb-Pb collisions for five different collisions centrality. The experimental data and the figure are taken from Ref. [195]. Before we proceed any further we note the following experimental observation from the CMS paper.
• At higher transverse momentum (p T ≥ 2 GeV), v 2 becomes even smaller than the higher-order v 3 , v 4 , and at much higher values of p T it becomes smaller than other higher order v n .
• The p T averaged v 2 and v 3 are found to be equal within 2%, while other higher-order v n decrease as n increases.  [196].
The evolution of the QGP according to relativistic hydrodynamics simulations have been able to consistently explain experimentally measured v n 's for different centralities and for different colliding energies, it is natural to expect that it should also explain the measured v n in the ultra central collisions. Before we discuss the results of hydrodynamic simulations, we note that one need to carefully select events into centrality classes since the integrated v n 's are quite sensitive on the selection of centrality class as can be seen from the bottom panel of Fig. 8. We also note that it is computationally expensive to simulate such ultra-central collisions since the number of events within the given centrality class is significantly small compared to the total number of minimum bias events. Although the essential p T dependence of charged hadrons v n and their observed ordering for ultra-central Pb-Pb collisions was nicely explained by a viscous hydrodynamic simulation using initial conditions from AMPT model [123]; see Fig. 9. However on careful observation we notice that at low p T < 1.5 GeV, the splitting from hydrodynamics simulation is larger than the corresponding experimental measurement. Similar disagreements are also evident for p T > 1.5 GeV in Fig. 10, which is taken from Ref. [196]. This can be seen more clearly from the p T integrated v 2 and v 3 in Fig. 11 which is also taken from Ref. [196]. In Ref. [196], the p T integrated v n was studied using (2+1)D viscous hydrodynamics model with MC-Glauber and MC-KLN initial conditions.
The nucleon-nucleon correlations in the colliding nucleus were also considered as a potential cause behind the experimentally measured v 2 ∼ v 3 . However, none of the initial condition model has so far been able to simultaneously explain the experimentally measured v n 's, as can be seen in Figs. 9, 10 and 11. In this regard we note that Denicol et al.,Ref. [197], have considered bulk viscosity along with the shear viscosity and the nucleon-nucleon correlations in order to explain this apparent discrepancy between the experimental data and corresponding theo- retical results, although there was some improvement but so far the effort remains unfruitful.

C. Longitudinal fluctuations and correlations
In relativistic heavy-ion collision experiments, a fraction of the incoming kinetic energy is converted into new matter deposited in the collision zone. The distribution of this matter in the plane transverse to the colliding beams is inhomogeneous and fluctuates from collision to collision. The lumpy initial energy density distribution and its event-by-event fluctuations lead to anisotropic flows of final hadrons through collective expansion in high energy heavy-ion collisions. The first numerical demonstration of the role of lumpy initial energy density (or event-by-event fluctuation) in the transverse plane (plane defined by the impact parameter vector and one of the perpendicular axis to the beam direction) to the experimentally observed non-zero odd flow harmonics (particularly third harmonics v 3 ) in heavy-ion collision was made by Alver and Roland [198]. From then on experimentally measured flow harmonics for all order (even and odd) has been successfully explained by viscous hydrodynamics model studies with fluctuating initial conditions such as Monte-Carlo (from now on we denote it by MC) Glauber [199,200], MC-CGC [201], URQMD [202], EPOS [203], AMPT [122], and IP-Glasma [204]. Fluctuations in the transverse plane not only give rise to odd flow harmonics but also significant even and odd v n in ultra central collisions [205]. They also result in p T dependent event planes, which break down the flow factorization v n,n (p T 1 , p T 2 ) = v n (p T 1 )v n (p T 2 ) [206]. Like the lumpy initial energy density in the transverse plane, it is also expected (the reason for which will be discussed shortly) that the energy density is lumpy in the longitudinal (space rapidity) direction.
Recent measurement of decorrelation of anisotropic flow along longitudinal direction by CMS collaboration has corroborated the above expectation. Studies of fluctuations along the longitudinal direction and their effects on anisotropic flows of final charged hadrons have only recently been started. At present the current understanding of longitudinal correlation (or decorrelation) of flow harmonics is as follows • The fluctuations of energy density along the longitudinal direction due to the fragmentation and different lengths of the coloured string produced in the scattering of nucleons [207][208][209].
• A gradual twist of the fireball (or more specifically the event plane) along the longitudinal direction Ref. [210,211].
Let us discuss each of them separately. Regarding the contribution of colour string we shall particularly discuss here a recent study Ref. [207] where AMPT transport model is used to evaluate the initial conditions for (3+1)D hydrodynamic model. AMPT uses HIJING to generate initial partons from hard and semi-hard scatterings and excited strings from soft interactions. The number of mini-jet partons per binary nucleon-nucleon collision in hard and semi-hard scatterings follow a Poisson distribution with the mean value given by the jet cross-section. The number of excited strings is equal to the number of participant nucleons in each event. Besides random fluctuations from mini-jet partons, the parton density fluctuates along longitudinal direction according to the length of strings. There are basically three types of strings: 1. Strings associated with each wounded nucleon (between a valence quark and a diquark), 2. Single strings between q-q pairs from quark annihilation and gluon fusion processes, 3. Strings between one hard parton from parton scatterings and valence quark or diquark in wounded nucleons.
These strings finally fragment into the partons along the longitudinal direction and giving rise to fluctuating energy density distribution in η s , see Fig. 12 for the distribution of coloured strings in the longitudinal direction for a typical Pb-Pb collision at √ s N N = 2.76 TeV. For more details about the longitudinal fluctuations and the visualization of parton density distribution in η s see Ref. [207]. The idea of a gradual twist of the fireball (or torqued fireball) along the longitudinal direction is due to Piotr Bozek et al. [210]. According to Bozek et al. the following ingredients are responsible for the appearance of the torque effect: 1. Statistical fluctuations of the transverse density of the sources (wounded nucleons), and 2. The asymmetric shape of the particle emission function, peaked in the forward (backward) rapidity for the forward (backward) moving wounded nucleons.
We note that in this model, the initial energy density profile is parametrized in such a way that after the hydrodynamics evolution and the freeze-out the hadronic spectra produced at different rapidities match with the corresponding experimental data. Whereas in the case of AMPT initial condition we do not need to use such procedure in order to explain the corresponding experimental data, for example in Fig. 13 we show the comparison of experimental data of longitudinal correlation for Pb-Pb collisions from CMS collaboration [212] and a (3+1)D hydrodynamics simulation result with AMPT initial condition. Note that with AMPT initial condition the experimental data is quite well described by the (3+1)D ideal hydrodynamics simulation. In the AMPT initial condition both longitudinal fluctuations and torque effects are present, the interplay of twist and fluctuation and the relative contribution of this two effects in heavy-ion collisions was studied within (3+1)D hydrodynamics model and AMPT in Ref. [208]. Many techniques have been proposed to study the longitudinal structure of final hadron production in heavyion collisions and the underlying mechanisms. For example, three-particle correlations were suggested to measure the twist effect [213] in heavy-ion collisions at RHIC. One can also characterise the longitudinal fluctuations in terms of coefficients in the Chebyshev polynomials [214] and the Legendre polynomial expansion of two-particle correlations in pseudorapidity [215,216]. The most intuitive method is to measure the forward-backward (FB) event plane angles or anisotropic flow differences [208] with varying pseudorapidity gaps. These methods are used within the torqued fireball model [210], (3+1)D hydrodynamics model and the AMPT model [208,217] to study the decorrelation of event plane angles or anisotropic flow along the pseudorapidity direction. Jia et al. [211] also proposed an "event-shape twist" technique to study the event plane decorrelation due to the twist in initial energy density distributions by selecting events with big FB event plane angle differences. Alternatively by selecting events with vanishingly small FB event plane angle differences, one can then eliminate the twist effect and the measured decorrelation of anisotropic flow with finite pseudorapidity gaps should be caused only by random fluctuations of event plane angles as was done in Ref. [208]. Before ending this section we note that the experimentally observed difference in the longitudinal correlation (r 2 and r 3 ) for different reference rapidity bin [212] is not yet understood within theoretical model studies [207]. We need further studies in order to understand those finer details.

D. Flow in intense magnetic field
The most strongest known magnetic field (| B| ∼ 10 18 − 10 19 Gauss) in the universe is produced in laboratory experiments of Au-Au or Pb-Pb collisions such as at RHIC and at LHC. Previous theoretical studies show that the intensity of the produced magnetic field rises approximately linearly with the centre of mass energy ( √ s NN ) of the colliding nucleons [218,219]. The corresponding electric fields in such collisions also becomes very strong which is same order of magnitude as the magnetic field (e B ≈ e E ∼ 10m 2 π for a typical Au-Au collision at top RHIC energy √ s NN = 200 GeV) [220], where m π is the pion mass. Such intense electric and magnetic fields are strong enough to initiate the particle production from vacuum via Schwinger mechanism [221].
The origin of such large electric and magnetic field is the relativistic velocities of the positive charge nucleus.
Within a MC Glauber model the electric ( E) and magnetic ( B) field at position r and at time t for a nucleus of charge Ze moving with velocity v in straight line is given by Here R i = x − x i (t) is the distance from a proton at position x i to x where the field is evaluated. In the above expression the summation index i denotes the contribution of all protons inside the colliding nucleus, for example Fig. 14 shows the positions of nucleons inside the two Au nucleus for a typical peripheral collision calculated in MC-Glauber model. Due to the fluctuating proton position from event to event the electric and magnetic field becomes irregular both in direction and in magnitude in the transverse plane. Moreover, the magnetic field in the central collisions becomes non-zero for such initial random proton positions. This can be seen from Fig. 15 where the event averaged value of B y as a function of impact parameter b is shown. The black dashed-dotted line corresponds to the average of the absolute magnitude of B y which is clearly non-zero even for b = 0 fm collisions. Note that we have used the natural unit where = c = k b = ε 0 = µ 0 = 1, with this choice the electric charge e = » 4π 137 becomes a dimensionless number. In the limit v ∼ c, the denominator in Eq. (135) and (136) becomes very small and we have large E and B.
There are large number of theoretical predictions based on the expectation of the creation of large magnetic field in heavy-ion collisions such as chiral magnetic effect, chiral electric effect, chiral magnetic waves etc., the discussion of which is beyond the scope of this review. For more details, we refer the reader to the following references [222][223][224][225][226]. Here we will concentrate on the possible effect of this large electro-magnetic field on the initial energy density and the subsequent hydrodynamics evolution of QGP produced in RHIC or LHC experiments. To the best of our knowledge, one of the first numerical study of the effect of magnetic fields on the hydrodynamics evolution in heavy-ion collisions are by Gursoy et al. [227] and by Hirono et al. [228]. However, those studies were based on several assumptions and none of them have considered the full magneto-hydrodynamics solution for the QGP evolution. We note here that the electric and magnetic field might affect the initial energy density, subsequent hydrodynamics stage, and the freeze-out distribution functions provided that the field is strong enough and lives until freeze-out.
One of us has recently studied the importance of electromagnetic field energy density compared to the energy density of the QGP fluid for Au-Au collisions at √ s N N = 200 GeV in Ref. [229]. The ratio (σ) of the magnetic field energy density to the fluid energy density was found to be ∼ 1 for peripheral collisions, but in central collisions σ << 1. It was also found that electric field also contributes similar energy density as magnetic field. Recent study by Tuchin [230] shows that the decay of the initial magnetic field can be substantially delayed in the case of finite electrical conductivity of QGP. Thus it becomes increasingly important to consider the electromagnetic field in the hydrodynamic evolution of heavy-ion collisions [231]. In Refs. [232][233][234] analytic solution of relativistic hydrodynamics for simplified cases was obtained. Finding analytic solution for general initial conditions is very difficult and there are very few analytical solution that exists for relativistic magnetohydrodynamics. The only possible way is to use numerical methods to solve magnetohydrodynamic equations relevant for heavy-ion collisions. This is not an easy task to accomplish. Initial effort in this direction can be found in Ref. [235]. However, we note that the authors of Ref. [235] have solved usual hydrodynamics conservation equations (without magnetic field) by considering an external force originating due to the paramagnetic interaction of QGP with the magnetic field.
One of us has also recently solved the hydrodynamics equations where magnetic field is taken into account in the energy-momentum tensor of the fluid. A realistic space and time dependence of magnetic field is considered for an ideal fluid evolution in Au-Au collisions at √ s N N = 200 GeV. It is found that in the presence of a finite electrical conductivity of QGP, the elliptic flow of π − increases noticeably, depending on the details of the magnitude of the magnetic field and the subsequent time evolution of the field. This is still a very new field of study and at present more detailed investigations are underway. In Fig. 16 we show the v 2 of π − obtained for impact parameter b = 10 fm Au-Au collisions at √ s N N = 200 GeV collisions. Initial value of magnetic field is taken to be 10m 2 π and the time variation of the magnetic field is obtained by parametrizing the results from Ref. [230]. A realistic spatial profile for y-component of magnetic field was considered for the simulation. From Fig. 16 one can see that v 2 of π − is noticeably enhanced in the presence of magnetic field (blue dashed line) compared to the case of no magnetic field (red line).

X. OUTLOOK
In this review article, we have discussed various aspects of relativistic dissipative hydrodynamics and its application to high energy heavy-ion collisions. While considerable success has been achieved in explaining many experimental observations, there are several issues that still needs further investigation. For example the experimentally measured longitudinal correlation of flow harmonics shows a splitting in the quantities corresponds to the correlation measure namely r 2 (η a , η b ) and r 3 (η a , η b ) for two different reference rapidity windows, which cannot be explained within a (3+1)D ideal hydrodynamics model with initial condition obtained from HIJING model. The reason behind this experimentally observed difference in r 2 (η a , η b ) and r 3 (η a , η b ) is still poorly understood.
It was also argued in the present review that the effect of magnetic field might not be negligible on the hydrodynamics evolution of QGP produced in heavy-ion collisions. Particularly it may be important for the reason that the elliptic and higher-order flow harmonics might be affected under such strong magnetic field. However, at present there are some open issues in this regard: the electrical conductivity of the QGP might play an essential role in the temporal decay of the magnetic field. A poor knowledge of the temperature dependent electrical conductivity is one of the major source of uncertainties. In addition to that, the magnetic susceptibility of the QGP and the hadron resonance gas should be included for a realistic calculation.
Another unsolved problem is the experimentally measured v 2 , v 3 in ultra-central collisions. Within error bars the magnitude of v 2 and v 3 are observed to be same for ultra-central collisions (0-0.2%). Although viscous hydrodynamics model with MC-KLN initial condition considering nucleon-nucleon correlations produces quite close result to the experimental measurement, it is still not fully explained within the given error. Another puzzling aspect of high energy collisions is the observation of flow like behaviour in small systems. Initial study shows that the experimentally observed flow harmonics and other bulk observables for high multiplicity p-p and p-Pb collisions can be well described within viscous hydrodynamics model simulation. A detailed theoretical explanation of how such small system behave collectively is still not well understood. There are also possibility for non-hydrodynamical origin of this observed flow in such small system. This is a topic of current research and we hope we will have more clear theoretical understanding within next few years when further studies including alternative possibilities will be available.
Finally we note that the field of high energy collisions is a very active area of research, we have not covered all aspects of the recent developments in the field related to hydrodynamics/collectivity of the QGP. For example the event-by-event distribution of flow harmonics [236] and their correlations [237,238] emerges as a promising observable to better constrain the initial conditions and the shear viscosity of the QGP. The event-by-event study of photon and dilepton production within viscous hydrodynamics provide us another window to look at the early stages of the heavy-ion collisions [239,240]. We hope that through all these ongoing experimental and theoretical/phenomenological studies we will have much refined understanding about the collective behaviour and the transport properties of the QGP in the near future.