# Relativistic Hydrodynamics in Heavy-Ion Collisions: General Aspects and Recent Developments

Authors: Amaresh Jaiswal; Victor Roy
Journal: Advances in High Energy Physics
 (2016)

---

## Abstract

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 ultracentral collisions, longitudinal fluctuations, and correlations and flow in intense magnetic field are also discussed.


---

## Body

## 1. 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–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 realised on earth, in a controlled experimental environment, by colliding two heavy nuclei with ultrarelativistic 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 experiments 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–14] and recently at LHC energies [15–18]. The regime with relatively large baryon chemical potential will be probed by the upcoming experimental facilities like Facility for Antiproton 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 Figure 1 [19].Figure 1
Schematic phase diagram of the QCD matter. The net baryon density onx-axis is normalized to that of the normal nuclear matter [19].It is possible to create hot and dense nuclear matter with very high energy densities in relatively large volumes by colliding ultrarelativistic 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 noninteracting. 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 modeled: starting from the stage where two highly Lorentz contracted heavy nuclei collide with each other, the formation and thermalization of the QGP or deconfined phase in the initial stages of the collision, its subsequent space-time evolution, the phase transition to the hadronic or 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 ultrarelativistic heavy-ion collisions are schematically illustrated in Figure2 [20].Figure 2
Various stages of ultrarelativistic heavy-ion collisions [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–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, that is, 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, and so forth. 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; that is, η/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)nlmfpp¯, 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; that is, lmfpp¯≳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 organized as follows. In Section2 we discuss the general formalism of a causal theory of relativistic dissipative fluid dynamics. Section 3 deals with the initial conditions necessary for the modeling of relativistic heavy-ion, p-p, and p-Pb collisions. In Section 4, we discuss some models of preequilibrium dynamics employed before hydrodynamic evolution. Section 5 briefly covers various equations of state, necessary to close the hydrodynamic equations. In Section 6, particle production mechanism via Cooper-Frye freeze-out and anisotropic flow generation is discussed. In Section 7, we discuss models of hadronic rescattering after freeze-out and contribution to particle spectra and flow from resonance decays. Section 8 deals with the extraction of transport coefficients from hydrodynamic analysis of flow data. Finally, in Section 9, 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=kB=1, with ħ=h/2π, where h is the Planck constant, c is the velocity of light in vacuum, and kB is the Boltzmann constant. Unless stated otherwise, the space-time 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 τ=t2-z2, the radial coordinate r=x2+y2, the azimuthal angle φ=tan-1(y/x), and space-time rapidity ηs=tanh-1(z/t). Hence, t=τcosh⁡ηs, x=rcos⁡φ, y=rsin⁡φ, 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,-r2,-τ2). Roman letters are used to indicate indices that vary from 1 to 3 and Greek letters are used for indices that vary from 0 to 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 byuμ 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(1)∂μ=∇μ+uμD,where∇μ≡Δμν∂ν,D≡uμ∂μ.In the fluid rest frame, D reduces to the time derivative and ∇μ reduces to the spacial gradient. Hence, the notation f˙≡Df is also commonly used. We also frequently use the symmetric, antisymmetric, and angular brackets notations defined as (2)AμBν≡12AμBν+AνBμ,AμBν≡12AμBν-AνBμ,AμBν≡ΔμναβAαBβ,where (3)Δμναβ≡12ΔμαΔνβ+ΔναΔμβ-23ΔαβΔμν is the traceless symmetric projection operator orthogonal to uμ satisfying the conditions ΔμναβΔαβ=ΔμναβΔμν=0.


## 2. Relativistic Fluid Dynamics

The physical characterization of a system consisting of many degrees of freedom is in general a nontrivial 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 which 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; that is, any spatial gradients can be ignored, and it 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 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. 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, 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 ultrarelativistic heavy-ion 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 [33], 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; that is, 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 [34] and Landau and Lifshitz [35]. 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 [36–38].The second-order Israel-Stewart (IS) theory [39] restores causality but may not guarantee stability [40]. 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 which 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 ultrarelativistic 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 modeling 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.


### 2.1. 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 [41–43]:(i)
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(ii)
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(iii)
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(iv)
Third Law: it is also known as Nernst’s heat theorem and 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 operationsThe 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) particle exchange with an external medium. Hence, the conservation law relating the small changes in state variables,E, V, and N is(4)δE=δQ-PδV+μδN,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 which 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:(5)dE=TdS-PdV+μdN.Hence, using (5), the intensive quantities, T, μ and P, can be obtained in terms of partial derivatives of the entropy as(6)∂S∂EN,V=1T,∂S∂VN,E=PT,∂S∂NE,V=-μT.The entropy is mathematically defined as an extensive and additive function of the state variables, which means that(7)SλE,λV,λN=λSE,V,N.Differentiating both sides with respect to λ, we obtain(8)S=E∂S∂λEλN,λV+V∂S∂λVλN,λE+N∂S∂λNλE,λV,which holds for any arbitrary value of λ. Setting λ=1 and using (6), we obtain the so-called Euler’s relation:(9)E=-PV+TS+μN.Using Euler’s relation, (9), along with the first law of thermodynamics, (5), we arrive at the Gibbs-Duhem relation:(10)VdP=SdT+Ndμ.In terms of energy, entropy, and number of densities defined asϵ≡E/V, s≡S/V, and n≡N/V, respectively, Euler’s relation, (9), and Gibbs-Duhem relation, (10), reduce to(11)ϵ=-P+Ts+μn,(12)dP=sdT+ndμ.Differentiating (11) and using (12), we obtain the relation analogous to first law of thermodynamics:(13)dϵ=Tds+μdn⟹ds=1Tdϵ-μTdn.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 [43]. In the next section, we introduce and derive the equations of relativistic ideal fluid dynamics.


### 2.2. Relativistic Ideal Fluid Dynamics

An ideal fluid is defined by the assumption of local thermal equilibrium; that is, all fluid elements must be exactly in thermodynamic equilibrium [35, 44]. 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:(14)uμx≡dxμdτ.The proper time increment dτ is given by the line element(15)dτ2=gμνdxμdxν=dt2-dx→2=dt21-v→2,where v→≡dx→/dt. This implies that(16)uμx=dtdτdxμdt=γv→1v→,where γ(v→)=1/1-v→2. In the nonrelativistic 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(17)u2≡uμxgμνuνx=γ2v→1-v→2=1.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, that is, 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, TLRFμν, the (net) particle four-current, NLRFμ, and the entropy four-current, SLRFμ, should have the characteristic form of a system in static equilibrium. In other words, in local rest frame, there is no flow of energy (TLRFi0=0), the force per unit surface element is isotropic (TLRFij=δijP), and there is no particle and entropy flow (N→=0 and S→=0). Consequently, the energy-momentum tensor, particle four-current, and entropy four-current in this frame take the following simple forms:(18)TLRFμν=ϵ0000P0000P0000P,NLRFμ=n000,SLRFμ=s000.For an ideal relativistic fluid, the general form of the energy-momentum tensor,T(0)μν, (net) particle four-current, 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(19)T0μν=c1uμuν+c2gμν,N0μ=c3uμ,S0μ=c4uμ.In the local rest frame, v→=0⇒uμ=(1,0→). Hence, in this frame, (19) takes the form(20)T0LRFμν=c1+c20000-c20000-c20000-c2,N0LRFμ=c3000,S0LRFμ=c4000.By comparing the above equation with the corresponding general expressions in the local rest frame, (18), one obtains the following expressions for the coefficients:(21)c1=ϵ+P,c2=-P,c3=n,c4=s.The conserved currents of an ideal fluid can then be expressed as(22)T0μν=ϵuμuν-PΔμν,N0μ=nuμ,S0μ=suμ,where Δμν=gμν-uμuν is the projection operator on the three-space orthogonal to uμ and satisfies the following properties of an orthogonal projector:(23)uμΔμν=Δμνuν=0,ΔρμΔρν=Δμν,Δμμ=3.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:(24)∂μT0μν=0,∂μN0μ=0,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μ:(25)D≡uμ∂μ,∇μ≡Δμρ∂ρ,⟹∂μ=uμD+∇μ.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:(26)uμ∂νT0μν=0⟹Dϵ+ϵ+Pθ=0,(27)Δμα∂νT0μν=0⟹ϵ+PDuα-∇αP=0,(28)∂μN0μ=0⟹Dn+nθ=0,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.


### 2.3. Covariant Thermodynamics

In the following, we rewrite the equilibrium thermodynamic relations derived in Section2.1, (11), (12), and (13), in a covariant form [39, 45]. For this purpose, it is convenient to introduce the following notations:(29)β≡1T,α≡μT,βμ≡uμT.In these notations, the covariant version of Euler’s relation, (11), and the Gibbs-Duhem relation, (12), can be postulated as(30)S0μ=Pβμ+βνT0μν-αN0μ,(31)dPβμ=N0μdα-T0μνdβν,respectively. The above equations can then be used to derive a covariant form of the first law of thermodynamics, (13):(32)dS0μ=βνdT0μν-αdN0μ.The covariant thermodynamic relations were constructed in such a way that when (30), (31), and (32) are contracted with uμ,(33)uμS0μ-Pβμ-βνT0μν+αN0μ=0⟹s+αn-βϵ+P=0,uμdPβμ-N0μdα+T0μνdβν=0⟹dβP-ndα+ϵdβ=0,uμdS0μ-βνdT0μν+αdN0μ=0⟹ds-βdϵ+αdn=0, we obtain the usual thermodynamic relations, (11), (12), and (13). Here we have used the property of the fluid four-velocity, uμuμ=1⇒uμduμ=0. The projection of (30), (31), and (32) on the three-space orthogonal to uμ just leads to trivial identities:(34)ΔμαS0μ-Pβμ-βνT0μν+αN0μ=0⟹0=0,ΔμαdPβμ-N0μdα+T0μνdβν=0⟹0=0,ΔμαdS0μ-βνdT0μν+αdN0μ=0⟹0=0.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, (32), leads to the following expression for the entropy four-current divergence:(35)∂μS0μ=βμ∂νT0μν-α∂μN0μ.After employing the conservation of energy-momentum and net particle number, (24), 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(36)∂μS0μ=0⟹Ds+sθ=0.We observe that the rate equation of the entropy density in the above equation is identical to that of the net particle number, (28). Therefore, we conclude that, for ideal hydrodynamics, the ratio of entropy density to number density (s/n) is a constant of motion.


### 2.4. Relativistic Dissipative Fluid Dynamics

The derivation of relativistic ideal fluid dynamics proceeds by employing the properties of the Lorentz transformation and 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 [33], the assumption of local thermodynamic equilibrium is never strictly realised in practice. In the following, we consider a more general theory of fluid dynamics which 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 was due to Eckart [34], in 1940, and, later, by Landau and Lifshitz [35], in 1987. 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 nonrelativistic fluids from weakly coupled gases such as air to strongly coupled fluids such as water. Hence, a relativistic generalization 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 generalization of Navier-Stokes theory is intrinsically unstable [36–38]. The source of such instability is attributed to the inherent acausal behaviour of this theory [46, 47]. A straightforward relativistic generalization of Navier-Stokes theory allows signals to propagate with infinite speed in a medium. While in nonrelativistic 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:(37)∂μTμν=0,∂μNμ=0.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)μ:(38)Tμν=T0μν+τμν=ϵuμuν-PΔμν+τμν,Nμ=N0μ+nμ=nuμ+nμ,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.


#### 2.4.1. 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 was 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 [39]:(39)ϵ≡uμuνTμν,n≡uμNμ.These matching conditions enforce the following constraints on the dissipative currents:(40)uμuντμν=0,uμnμ=0.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.


#### 2.4.2. Tensor Decompositions of Dissipative Quantities

To proceed further, it is convenient to decomposeτμν in terms of its irreducible components, that is, 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, (40), satisfied by τμν. To this end, we introduce another projection operator, the double symmetric, traceless projector orthogonal to uμ,(41)Δαβμν≡12ΔαμΔβν+ΔβμΔαν-23ΔμνΔαβ,with the following properties:(42)Δμναβ=Δαβμν,ΔρσμνΔαβρσ=Δαβμν,uμΔαβμν=gμνΔαβμν=0,Δμνμν=5.The parentheses in the above equation denote symmetrization of the Lorentz indices; that is, A(μν)≡(Aμν+Aνμ)/2. The dissipative current τμν then can be tensor decomposed in its irreducible form by using uμ, Δμν, and Δαβμν as(43)τμν≡-ΠΔμν+2uμhν+πμν,where we have defined(44)Π≡-13Δαβταβ,hμ≡Δαμuβταβ,πμν≡Δαβμνταβ.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(45)ϵ=uαuβTαβ,n=uαNα,Π=-P-13ΔαβTαβ,hμ=uαTμα,nμ=Nμ,πμν=Tμν,where the angular bracket notations are defined as A〈μ〉≡ΔαμAα and B〈μν〉≡ΔαβμνBαβ.We observe thatTμν is a symmetric second-rank tensor with ten independent components and Nμ is a four-vector; 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 is 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.


#### 2.4.3. Definition of the Velocity Field

In the process of formulating the theory of dissipative fluid dynamics, the next important step is to fixuμ. 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: theEckart definition [34], in which the velocity is defined by the flow of particles,(46)Nμ=nuμ⟹nμ=0,and theLandau definition [35], in which the velocity is specified by the flow of the total energy,(47)uνTμν=ϵuμ⟹hμ=0.We note that the above two definitions ofuμ 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, (47). The conserved currents in this frame take the following form:(48)Tμν=ϵuμuν-P+ΠΔμν+πμν,Nμ=nuμ+nμ.As done for ideal fluids, the energy-momentum conservation equation in (37) is decomposed parallel and orthogonal to uμ. Using (48) together with the conservation law for particle number in (37) leads to the equations of motion for dissipative fluids. For uμ∂νTμν=0, Δμα∂νTμν=0, and ∂μNμ=0, one obtains(49)ϵ˙+ϵ+P+Πθ-πμνσμν=0,(50)ϵ+P+Πu˙α-∇αP+Π+Δμα∂νπμν=0,(51)n˙+nθ+∂μnμ=0,respectively. Here A˙≡DA=uμ∂μA, and the shear tensor σμν≡∇μuν=Δαβμν∇αuβ.We observe that while there are fourteen total independent components ofTμν and Nμ, (49)–(51) 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 (49)–(51). Eventually, this corresponds to finding the closed dynamical or constitutive relations satisfied by the dissipative tensors Π, nμ, and πμν.


#### 2.4.4. Relativistic Navier-Stokes Theory

In the presence of dissipative currents, the entropy is no longer a conserved quantity; that is,∂μSμ≠0. Since the form of the entropy four-current for a dissipative fluid is not knowna priori, it is not trivial to obtain its equation. We proceed by recalling the form of the entropy four-current for ideal fluids, (30), and extending it for dissipative fluids:(52)Sμ=Pβμ+βνTμν-αNμ.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 [34, 35]. 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 (48) in (52). Taking the divergence and using (49)–(51), we obtain(53)∂μSμ=-βΠθ-nμ∇μα+βπμνσμν.The relativistic Navier-Stokes theory can then be obtained by applying the second law of thermodynamics to each fluid element; that is, by requiring that the entropy production∂μSμ must always be positive,(54)-βΠθ-nμ∇μα+βπμνσμν≥0.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(55)Π=-ζθ,nμ=κ∇μα,πμν=2ησμν,where the proportionality coefficients ζ, κ, and η refer to the bulk viscosity, the particle diffusion, and the shear viscosity, respectively. Substituting the above equation in (53), we observe that the source term for entropy production becomes a quadratic function of the dissipative currents:(56)∂μSμ=βζΠ2-1κnμnμ+β2ηπμνπμν.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, πμνπμν>0. Hence, as long as ζ,κ,η≥0, the entropy production is always positive. Constitutive relations for the dissipative quantities, (55), along with (49)–(51) are known as the relativistic Navier-Stokes equations.The relativistic Navier-Stokes theory in this form was obtained originally by Landau and Lifshitz [35]. A similar theory was derived independently by Eckart, using a different definition of the fluid four-velocity [34]. 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, (55). 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 [39, 45, 48], the divergence-type theory [49, 50], extended irreversible thermodynamics [51], Carter’s theory [52], Grmela and Ottinger theory [53], 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.


#### 2.4.5. 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 nonequilibrium entropy four-current [39, 45, 48]. In equilibrium, the entropy four-current was expressed exactly in terms of the primary fluid dynamical variables, (30). Strictly speaking, the nonequilibrium entropy four-current should depend on a larger number of independent dynamical variables, and a direct extension of (30) to (52) 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(57)Sμ=Pβμ+βνTμν-αNμ-QμδNμ,δTμν,where Qμ is a function of deviations from local equilibrium, δNμ≡Nμ-N(0)μ, and δTμν≡Tμν-T(0)μν. Using (48) and Taylor-expanding Qμ to second order in dissipative fluxes, we obtain(58)Sμ=suμ-αnμ-β0Π2-β1nνnν+β2πρσπρσuμ2T-α0ΠΔμν+α1πμνnνT+Oδ3,where O(δ3) denotes third-order terms in the dissipative currents and β0,β1,β2,α0, and α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 (58) 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 recalculate the entropy production, ∂μSμ, using the more general entropy four-current given in (58):(59)∂μSμ=-βΠθ+β0Π˙+βΠΠΠθ+ψαnΠnμu˙μ+α0∇μnμ+ψαΠnnμ∇μα-βnμT∇μα-β1n˙μ-βnnnμθ+α0∇μΠ+α1∇νπμν+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπμν∇να+χ~αnππμνu˙ν+βπμνσμν-β2π˙μν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν.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:(60)Π=-ζθ+β0Π˙+βΠΠΠθ+α0∇μnμ+ψαnΠnμu˙μ+ψαΠnnμ∇μα,nμ=λT∇μα-β1n˙μ-βnnnμθ+α0∇μΠ+α1Δρμ∇νπρν+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπνμ∇να+χ~αnππνμu˙ν,πμν=2ησμν-β2π˙μν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν,where λ≡κ/T. This implies that the dissipative currents must satisfy the dynamical equations:(61)Π˙+ΠτΠ=-1β0θ+βΠΠΠθ+ψαnΠnμu˙μ+α0∇μnμ+ψαΠnnμ∇μα,n˙μ+nμτn=1β1T∇μα-βnnnμθ+α1Δρμ∇νπρν+α0∇μΠ+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπνμ∇να+χ~αnππνμu˙ν,π˙μν+πμντπ=1β2σμν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν.The above equations for the dissipative quantities are relaxation-type equations with the relaxation times defined as(62)τΠ≡ζβ0,τn≡λβ1=κβ1T,τπ≡2ηβ2.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, that is, 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.


### 2.5. 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 single-particle 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 [54].Let us consider a system of relativistic particles, each having rest massm, momentum p→, and energy p0. Therefore, from relativity, we have p0=(p→)2+m2. 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μ=(p0,p→) at each space-time point such that f(x,p)Δ3xΔ3p gives the average number of particles at a given time t in the volume element Δ3x 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 Δ3x is large, on the other hand, Δ3x is small compared to macroscopic point of view.The particle densityn(x) is introduced to describe, in general, a nonuniform system, such that n(x)Δ3x is the average number of particles in volume Δ3x 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(63)nx=∫d3pfx,p,j→x=∫d3pv→fx,p,where v→=p→/p0 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(64)Nμx=∫d3pp0pμfx,p.Note that since d3p/p0 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 isp0, the average energy density and the energy flow can be written in terms of the distribution function as(65)T00x=∫d3pp0fx,p,T0ix=∫d3pp0vifx,p.The momentum density is defined as the average value of particle momenta pi, 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(66)Ti0x=∫d3ppifx,p,Tijx=∫d3ppivjfx,p.Combining all these in a compact covariant form using vi=pi/p0, we obtain the energy-momentum tensor of a macroscopic system:(67)Tμνx=∫d3pp0pμpνfx,p.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.H-function introduced by Boltzmann implies that the nonequilibrium local entropy density of a system can be written as(68)sx=-∫d3pfx,plnfx,p-1.The entropy flow corresponding to the above entropy density is(69)S→x=-∫d3pv→fx,plnfx,p-1.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(70)Sμx=-∫d3pp0pμfx,plnfx,p-1.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(71)Sμx=-∫d3pp0pμfx,plnfx,p+rf~x,plnf~x,p,where f~≡1-rf. The expressions for the entropy four-current given in (70) and (71) 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=f0+δf. The equilibrium distribution function f0 is defined as(72)f0x,p=1expβu·p-α+r,where the scalar product is defined as u·p≡uμpμ and r=0 for Maxwell-Boltzmann statistics. Note that, in equilibrium, that is, for f(x,p)=f0(x,p), the particle four-flow and energy-momentum tensor given in (64) and (67) reduce to those of ideal hydrodynamics N(0)μ and T(0)μν. Therefore, using (48), the dissipative quantities, namely, the bulk viscous pressure Π, the particle diffusion current nμ, and the shear stress tensor πμν, can be written as(73)Π=-13Δαβ∫d3pp0pαpβδf,nμ=Δμν∫d3pp0pνδf,πμν=Δαβμν∫d3pp0pαpβδf.The evolution equations for the dissipative quantities expressed in terms of the nonequilibrium distribution function, (73), 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(74)pμ∂μf=Cf,where dp≡d3p/p0 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(75)Cf=12∫dp′dkdk′Wpp′→kk′fkfk′f~pf~p′-fpfp′f~kf~k′,where Wpp′→kk′ is the collisional transition rate. The first and second terms within the integral of (75) 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 [55](76)Cf=-u·pδfτR,where τR is the relaxation time.


### 2.6. Dissipative Fluid Dynamics from Kinetic Theory

The derivation of a causal theory of relativistic dissipative hydrodynamics by Israel and Stewart [39] proceeds by invoking the second law of thermodynamics, namely, ∂μSμ≥0, from the algebraic form of the entropy four-current given in (58). 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 (70) and (71), in order to obtain the dissipative equations [56]. This essentially ensures that the nonequilibrium corrections to the distribution function, δf, do not violate the second law of thermodynamics. In [56], the generalized method of moments developed by Denicol et al. [57] was used to quantify the dissipative corrections to the distribution function. The form of the resultant dissipative equations, obtained in [56], is identical to (61), 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 [48], has been used quite extensively to quantify the dissipative corrections to the distribution function [56–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 [66]. In this method, which is based on a Chapman-Enskog like expansion, the Boltzmann equation in the relaxation time approximation,(77)pμ∂μf=-u·pτRf-f0,is solved iteratively to obtain δf up to any arbitrary order in derivatives. To first order and second order in gradients, one obtains(78)δf1=-τRu·ppμ∂μf0,(79)δf2=τRu·ppμpν∂μτRu·p∂νf0.This method of obtaining the form of the nonequilibrium distribution function is consistent with dissipative hydrodynamics, which is also formulated as a gradient expansion.The second-order evolution equations for the dissipative quantities are then obtained by substitutingδf=δf(1)+δf(2) from (78) and (79) in (73):(80)ΠτΠ=-Π˙-βΠθ-δΠΠΠθ+λΠππμνσμν-τΠnn·u˙-λΠnn·∇α-lΠn∂·n,nμτn=-n˙μ+βn∇μα-nνωνμ-λnnnνσνμ-δnnnμθ+λnΠΠ∇μα-λnππμν∇να-τnππνμu˙ν+τnΠΠu˙μ+lnπΔμν∂γπνγ-lnΠ∇μΠ,πμντπ=-π˙μν+2βπσμν+2πγμωνγ-τπππγμσνγ-δπππμνθ+λπΠΠσμν-τπnnμu˙ν+λπnnμ∇να+lπn∇μnν,where ωμν=(∇μuν-∇νuμ)/2 is the vorticity tensor. It is interesting to note that although the form of the evolution equations for dissipative quantities in (80) is identical to those obtained in [58] using the moment method, the transport coefficients are, in general, different [67, 68]. Moreover, it was shown that the above described method, based on iterative solution of Boltzmann equation, leads to phenomenologically consistent corrections to the distribution function [69] and the transport coefficients exhibit intriguing similarities with strongly coupled conformal field theory [70, 71].Proceeding in a similar way, a third-order dissipative evolution equation can also be obtained [72–74]:(81)π˙μν=-πμντπ+2βπσμν+2πγμωνγ-107πγμσνγ-43πμνθ+257βππρμωνγπργ-13βππγμπνγθ-38245βππμνπργσργ-2249βππρμπνγσργ-2435∇μπνγu˙γτπ+435∇μτπ∇γπνγ-27∇γτπ∇μπνγ+127∇γτπu˙μπνγ-17∇γτπ∇γπμν+67∇γτπu˙γπμν-27τπωρμωνγπργ-27τππρμωνγωργ-1063τππμνθ2+2621τππγμωνγθ.It is reassuring that the results obtained using third-order evolution equation indicate convergence of the gradient expansion and show improvement over second order, when compared to the direct solutions of the Boltzmann equation [72–75].Apart from these standard formulations, there are several other formulations of relativistic dissipative hydrodynamics from kinetic theory. Among them, the ones which have gained widespread interest are anisotropic hydrodynamics and derivations based on renormalization group method. Anisotropic hydrodynamics is a nonperturbative reorganization of the standard relativistic hydrodynamics which takes into account the large momentum space anisotropies generated in ultrarelativistic heavy-ion collisions [76–81]. On the other hand, the derivation based on renormalization group method attempts to solve the Boltzmann equation, as faithfully as possible, in an organized perturbation scheme and resum the possible secular terms by a suitable setting of the initial value of the distribution function [82–85]. Since it is widely accepted that the QGP is momentum space anisotropic, application of anisotropic hydrodynamics to high energy heavy-ion collisions has phenomenological implications. Nevertheless, the dissipative hydrodynamic formulation based on renormalization group method is important in order to accurately determine the higher-order transport coefficients.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, do 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 is different from systems 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(82)πμν=2ησμν-τππ˙μν+43πμνθ-λ1η2πγμπνγ+λ2ηπγμωνγ+λ3ωγμωνγ.For a system of massless particles, the kinetic theory results for second-order transport coefficients agree in the case of both moment method [57] and the Chapman-Enskog like iterative solution of the Boltzmann equation [66]. In Table 1, 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.Table 1
Comparison of transport coefficients.
τ π T / ( η / s ) λ 1 T / ( η / s ) λ 2 T / ( η / s ) λ 3 ADS/CFT 2 ( 2 - ln ⁡ 2 ) 2 η 4 η ln ⁡ 2 0 KT 5 ( 25 / 7 ) η 10 η 0


## 2.1. 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 [41–43]:(i)
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(ii)
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(iii)
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(iv)
Third Law: it is also known as Nernst’s heat theorem and 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 operationsThe 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) particle exchange with an external medium. Hence, the conservation law relating the small changes in state variables,E, V, and N is(4)δE=δQ-PδV+μδN,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 which 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:(5)dE=TdS-PdV+μdN.Hence, using (5), the intensive quantities, T, μ and P, can be obtained in terms of partial derivatives of the entropy as(6)∂S∂EN,V=1T,∂S∂VN,E=PT,∂S∂NE,V=-μT.The entropy is mathematically defined as an extensive and additive function of the state variables, which means that(7)SλE,λV,λN=λSE,V,N.Differentiating both sides with respect to λ, we obtain(8)S=E∂S∂λEλN,λV+V∂S∂λVλN,λE+N∂S∂λNλE,λV,which holds for any arbitrary value of λ. Setting λ=1 and using (6), we obtain the so-called Euler’s relation:(9)E=-PV+TS+μN.Using Euler’s relation, (9), along with the first law of thermodynamics, (5), we arrive at the Gibbs-Duhem relation:(10)VdP=SdT+Ndμ.In terms of energy, entropy, and number of densities defined asϵ≡E/V, s≡S/V, and n≡N/V, respectively, Euler’s relation, (9), and Gibbs-Duhem relation, (10), reduce to(11)ϵ=-P+Ts+μn,(12)dP=sdT+ndμ.Differentiating (11) and using (12), we obtain the relation analogous to first law of thermodynamics:(13)dϵ=Tds+μdn⟹ds=1Tdϵ-μTdn.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 [43]. In the next section, we introduce and derive the equations of relativistic ideal fluid dynamics.


## 2.2. Relativistic Ideal Fluid Dynamics

An ideal fluid is defined by the assumption of local thermal equilibrium; that is, all fluid elements must be exactly in thermodynamic equilibrium [35, 44]. 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:(14)uμx≡dxμdτ.The proper time increment dτ is given by the line element(15)dτ2=gμνdxμdxν=dt2-dx→2=dt21-v→2,where v→≡dx→/dt. This implies that(16)uμx=dtdτdxμdt=γv→1v→,where γ(v→)=1/1-v→2. In the nonrelativistic 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(17)u2≡uμxgμνuνx=γ2v→1-v→2=1.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, that is, 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, TLRFμν, the (net) particle four-current, NLRFμ, and the entropy four-current, SLRFμ, should have the characteristic form of a system in static equilibrium. In other words, in local rest frame, there is no flow of energy (TLRFi0=0), the force per unit surface element is isotropic (TLRFij=δijP), and there is no particle and entropy flow (N→=0 and S→=0). Consequently, the energy-momentum tensor, particle four-current, and entropy four-current in this frame take the following simple forms:(18)TLRFμν=ϵ0000P0000P0000P,NLRFμ=n000,SLRFμ=s000.For an ideal relativistic fluid, the general form of the energy-momentum tensor,T(0)μν, (net) particle four-current, 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(19)T0μν=c1uμuν+c2gμν,N0μ=c3uμ,S0μ=c4uμ.In the local rest frame, v→=0⇒uμ=(1,0→). Hence, in this frame, (19) takes the form(20)T0LRFμν=c1+c20000-c20000-c20000-c2,N0LRFμ=c3000,S0LRFμ=c4000.By comparing the above equation with the corresponding general expressions in the local rest frame, (18), one obtains the following expressions for the coefficients:(21)c1=ϵ+P,c2=-P,c3=n,c4=s.The conserved currents of an ideal fluid can then be expressed as(22)T0μν=ϵuμuν-PΔμν,N0μ=nuμ,S0μ=suμ,where Δμν=gμν-uμuν is the projection operator on the three-space orthogonal to uμ and satisfies the following properties of an orthogonal projector:(23)uμΔμν=Δμνuν=0,ΔρμΔρν=Δμν,Δμμ=3.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:(24)∂μT0μν=0,∂μN0μ=0,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μ:(25)D≡uμ∂μ,∇μ≡Δμρ∂ρ,⟹∂μ=uμD+∇μ.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:(26)uμ∂νT0μν=0⟹Dϵ+ϵ+Pθ=0,(27)Δμα∂νT0μν=0⟹ϵ+PDuα-∇αP=0,(28)∂μN0μ=0⟹Dn+nθ=0,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.


## 2.3. Covariant Thermodynamics

In the following, we rewrite the equilibrium thermodynamic relations derived in Section2.1, (11), (12), and (13), in a covariant form [39, 45]. For this purpose, it is convenient to introduce the following notations:(29)β≡1T,α≡μT,βμ≡uμT.In these notations, the covariant version of Euler’s relation, (11), and the Gibbs-Duhem relation, (12), can be postulated as(30)S0μ=Pβμ+βνT0μν-αN0μ,(31)dPβμ=N0μdα-T0μνdβν,respectively. The above equations can then be used to derive a covariant form of the first law of thermodynamics, (13):(32)dS0μ=βνdT0μν-αdN0μ.The covariant thermodynamic relations were constructed in such a way that when (30), (31), and (32) are contracted with uμ,(33)uμS0μ-Pβμ-βνT0μν+αN0μ=0⟹s+αn-βϵ+P=0,uμdPβμ-N0μdα+T0μνdβν=0⟹dβP-ndα+ϵdβ=0,uμdS0μ-βνdT0μν+αdN0μ=0⟹ds-βdϵ+αdn=0, we obtain the usual thermodynamic relations, (11), (12), and (13). Here we have used the property of the fluid four-velocity, uμuμ=1⇒uμduμ=0. The projection of (30), (31), and (32) on the three-space orthogonal to uμ just leads to trivial identities:(34)ΔμαS0μ-Pβμ-βνT0μν+αN0μ=0⟹0=0,ΔμαdPβμ-N0μdα+T0μνdβν=0⟹0=0,ΔμαdS0μ-βνdT0μν+αdN0μ=0⟹0=0.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, (32), leads to the following expression for the entropy four-current divergence:(35)∂μS0μ=βμ∂νT0μν-α∂μN0μ.After employing the conservation of energy-momentum and net particle number, (24), 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(36)∂μS0μ=0⟹Ds+sθ=0.We observe that the rate equation of the entropy density in the above equation is identical to that of the net particle number, (28). Therefore, we conclude that, for ideal hydrodynamics, the ratio of entropy density to number density (s/n) is a constant of motion.


## 2.4. Relativistic Dissipative Fluid Dynamics

The derivation of relativistic ideal fluid dynamics proceeds by employing the properties of the Lorentz transformation and 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 [33], the assumption of local thermodynamic equilibrium is never strictly realised in practice. In the following, we consider a more general theory of fluid dynamics which 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 was due to Eckart [34], in 1940, and, later, by Landau and Lifshitz [35], in 1987. 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 nonrelativistic fluids from weakly coupled gases such as air to strongly coupled fluids such as water. Hence, a relativistic generalization 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 generalization of Navier-Stokes theory is intrinsically unstable [36–38]. The source of such instability is attributed to the inherent acausal behaviour of this theory [46, 47]. A straightforward relativistic generalization of Navier-Stokes theory allows signals to propagate with infinite speed in a medium. While in nonrelativistic 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:(37)∂μTμν=0,∂μNμ=0.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)μ:(38)Tμν=T0μν+τμν=ϵuμuν-PΔμν+τμν,Nμ=N0μ+nμ=nuμ+nμ,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.


### 2.4.1. 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 was 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 [39]:(39)ϵ≡uμuνTμν,n≡uμNμ.These matching conditions enforce the following constraints on the dissipative currents:(40)uμuντμν=0,uμnμ=0.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.


### 2.4.2. Tensor Decompositions of Dissipative Quantities

To proceed further, it is convenient to decomposeτμν in terms of its irreducible components, that is, 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, (40), satisfied by τμν. To this end, we introduce another projection operator, the double symmetric, traceless projector orthogonal to uμ,(41)Δαβμν≡12ΔαμΔβν+ΔβμΔαν-23ΔμνΔαβ,with the following properties:(42)Δμναβ=Δαβμν,ΔρσμνΔαβρσ=Δαβμν,uμΔαβμν=gμνΔαβμν=0,Δμνμν=5.The parentheses in the above equation denote symmetrization of the Lorentz indices; that is, A(μν)≡(Aμν+Aνμ)/2. The dissipative current τμν then can be tensor decomposed in its irreducible form by using uμ, Δμν, and Δαβμν as(43)τμν≡-ΠΔμν+2uμhν+πμν,where we have defined(44)Π≡-13Δαβταβ,hμ≡Δαμuβταβ,πμν≡Δαβμνταβ.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(45)ϵ=uαuβTαβ,n=uαNα,Π=-P-13ΔαβTαβ,hμ=uαTμα,nμ=Nμ,πμν=Tμν,where the angular bracket notations are defined as A〈μ〉≡ΔαμAα and B〈μν〉≡ΔαβμνBαβ.We observe thatTμν is a symmetric second-rank tensor with ten independent components and Nμ is a four-vector; 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 is 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.


### 2.4.3. Definition of the Velocity Field

In the process of formulating the theory of dissipative fluid dynamics, the next important step is to fixuμ. 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: theEckart definition [34], in which the velocity is defined by the flow of particles,(46)Nμ=nuμ⟹nμ=0,and theLandau definition [35], in which the velocity is specified by the flow of the total energy,(47)uνTμν=ϵuμ⟹hμ=0.We note that the above two definitions ofuμ 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, (47). The conserved currents in this frame take the following form:(48)Tμν=ϵuμuν-P+ΠΔμν+πμν,Nμ=nuμ+nμ.As done for ideal fluids, the energy-momentum conservation equation in (37) is decomposed parallel and orthogonal to uμ. Using (48) together with the conservation law for particle number in (37) leads to the equations of motion for dissipative fluids. For uμ∂νTμν=0, Δμα∂νTμν=0, and ∂μNμ=0, one obtains(49)ϵ˙+ϵ+P+Πθ-πμνσμν=0,(50)ϵ+P+Πu˙α-∇αP+Π+Δμα∂νπμν=0,(51)n˙+nθ+∂μnμ=0,respectively. Here A˙≡DA=uμ∂μA, and the shear tensor σμν≡∇μuν=Δαβμν∇αuβ.We observe that while there are fourteen total independent components ofTμν and Nμ, (49)–(51) 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 (49)–(51). Eventually, this corresponds to finding the closed dynamical or constitutive relations satisfied by the dissipative tensors Π, nμ, and πμν.


### 2.4.4. Relativistic Navier-Stokes Theory

In the presence of dissipative currents, the entropy is no longer a conserved quantity; that is,∂μSμ≠0. Since the form of the entropy four-current for a dissipative fluid is not knowna priori, it is not trivial to obtain its equation. We proceed by recalling the form of the entropy four-current for ideal fluids, (30), and extending it for dissipative fluids:(52)Sμ=Pβμ+βνTμν-αNμ.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 [34, 35]. 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 (48) in (52). Taking the divergence and using (49)–(51), we obtain(53)∂μSμ=-βΠθ-nμ∇μα+βπμνσμν.The relativistic Navier-Stokes theory can then be obtained by applying the second law of thermodynamics to each fluid element; that is, by requiring that the entropy production∂μSμ must always be positive,(54)-βΠθ-nμ∇μα+βπμνσμν≥0.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(55)Π=-ζθ,nμ=κ∇μα,πμν=2ησμν,where the proportionality coefficients ζ, κ, and η refer to the bulk viscosity, the particle diffusion, and the shear viscosity, respectively. Substituting the above equation in (53), we observe that the source term for entropy production becomes a quadratic function of the dissipative currents:(56)∂μSμ=βζΠ2-1κnμnμ+β2ηπμνπμν.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, πμνπμν>0. Hence, as long as ζ,κ,η≥0, the entropy production is always positive. Constitutive relations for the dissipative quantities, (55), along with (49)–(51) are known as the relativistic Navier-Stokes equations.The relativistic Navier-Stokes theory in this form was obtained originally by Landau and Lifshitz [35]. A similar theory was derived independently by Eckart, using a different definition of the fluid four-velocity [34]. 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, (55). 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 [39, 45, 48], the divergence-type theory [49, 50], extended irreversible thermodynamics [51], Carter’s theory [52], Grmela and Ottinger theory [53], 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.


### 2.4.5. 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 nonequilibrium entropy four-current [39, 45, 48]. In equilibrium, the entropy four-current was expressed exactly in terms of the primary fluid dynamical variables, (30). Strictly speaking, the nonequilibrium entropy four-current should depend on a larger number of independent dynamical variables, and a direct extension of (30) to (52) 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(57)Sμ=Pβμ+βνTμν-αNμ-QμδNμ,δTμν,where Qμ is a function of deviations from local equilibrium, δNμ≡Nμ-N(0)μ, and δTμν≡Tμν-T(0)μν. Using (48) and Taylor-expanding Qμ to second order in dissipative fluxes, we obtain(58)Sμ=suμ-αnμ-β0Π2-β1nνnν+β2πρσπρσuμ2T-α0ΠΔμν+α1πμνnνT+Oδ3,where O(δ3) denotes third-order terms in the dissipative currents and β0,β1,β2,α0, and α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 (58) 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 recalculate the entropy production, ∂μSμ, using the more general entropy four-current given in (58):(59)∂μSμ=-βΠθ+β0Π˙+βΠΠΠθ+ψαnΠnμu˙μ+α0∇μnμ+ψαΠnnμ∇μα-βnμT∇μα-β1n˙μ-βnnnμθ+α0∇μΠ+α1∇νπμν+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπμν∇να+χ~αnππμνu˙ν+βπμνσμν-β2π˙μν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν.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:(60)Π=-ζθ+β0Π˙+βΠΠΠθ+α0∇μnμ+ψαnΠnμu˙μ+ψαΠnnμ∇μα,nμ=λT∇μα-β1n˙μ-βnnnμθ+α0∇μΠ+α1Δρμ∇νπρν+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπνμ∇να+χ~αnππνμu˙ν,πμν=2ησμν-β2π˙μν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν,where λ≡κ/T. This implies that the dissipative currents must satisfy the dynamical equations:(61)Π˙+ΠτΠ=-1β0θ+βΠΠΠθ+ψαnΠnμu˙μ+α0∇μnμ+ψαΠnnμ∇μα,n˙μ+nμτn=1β1T∇μα-βnnnμθ+α1Δρμ∇νπρν+α0∇μΠ+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπνμ∇να+χ~αnππνμu˙ν,π˙μν+πμντπ=1β2σμν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν.The above equations for the dissipative quantities are relaxation-type equations with the relaxation times defined as(62)τΠ≡ζβ0,τn≡λβ1=κβ1T,τπ≡2ηβ2.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, that is, 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.


## 2.4.1. 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 was 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 [39]:(39)ϵ≡uμuνTμν,n≡uμNμ.These matching conditions enforce the following constraints on the dissipative currents:(40)uμuντμν=0,uμnμ=0.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.


## 2.4.2. Tensor Decompositions of Dissipative Quantities

To proceed further, it is convenient to decomposeτμν in terms of its irreducible components, that is, 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, (40), satisfied by τμν. To this end, we introduce another projection operator, the double symmetric, traceless projector orthogonal to uμ,(41)Δαβμν≡12ΔαμΔβν+ΔβμΔαν-23ΔμνΔαβ,with the following properties:(42)Δμναβ=Δαβμν,ΔρσμνΔαβρσ=Δαβμν,uμΔαβμν=gμνΔαβμν=0,Δμνμν=5.The parentheses in the above equation denote symmetrization of the Lorentz indices; that is, A(μν)≡(Aμν+Aνμ)/2. The dissipative current τμν then can be tensor decomposed in its irreducible form by using uμ, Δμν, and Δαβμν as(43)τμν≡-ΠΔμν+2uμhν+πμν,where we have defined(44)Π≡-13Δαβταβ,hμ≡Δαμuβταβ,πμν≡Δαβμνταβ.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(45)ϵ=uαuβTαβ,n=uαNα,Π=-P-13ΔαβTαβ,hμ=uαTμα,nμ=Nμ,πμν=Tμν,where the angular bracket notations are defined as A〈μ〉≡ΔαμAα and B〈μν〉≡ΔαβμνBαβ.We observe thatTμν is a symmetric second-rank tensor with ten independent components and Nμ is a four-vector; 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 is 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.


## 2.4.3. Definition of the Velocity Field

In the process of formulating the theory of dissipative fluid dynamics, the next important step is to fixuμ. 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: theEckart definition [34], in which the velocity is defined by the flow of particles,(46)Nμ=nuμ⟹nμ=0,and theLandau definition [35], in which the velocity is specified by the flow of the total energy,(47)uνTμν=ϵuμ⟹hμ=0.We note that the above two definitions ofuμ 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, (47). The conserved currents in this frame take the following form:(48)Tμν=ϵuμuν-P+ΠΔμν+πμν,Nμ=nuμ+nμ.As done for ideal fluids, the energy-momentum conservation equation in (37) is decomposed parallel and orthogonal to uμ. Using (48) together with the conservation law for particle number in (37) leads to the equations of motion for dissipative fluids. For uμ∂νTμν=0, Δμα∂νTμν=0, and ∂μNμ=0, one obtains(49)ϵ˙+ϵ+P+Πθ-πμνσμν=0,(50)ϵ+P+Πu˙α-∇αP+Π+Δμα∂νπμν=0,(51)n˙+nθ+∂μnμ=0,respectively. Here A˙≡DA=uμ∂μA, and the shear tensor σμν≡∇μuν=Δαβμν∇αuβ.We observe that while there are fourteen total independent components ofTμν and Nμ, (49)–(51) 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 (49)–(51). Eventually, this corresponds to finding the closed dynamical or constitutive relations satisfied by the dissipative tensors Π, nμ, and πμν.


## 2.4.4. Relativistic Navier-Stokes Theory

In the presence of dissipative currents, the entropy is no longer a conserved quantity; that is,∂μSμ≠0. Since the form of the entropy four-current for a dissipative fluid is not knowna priori, it is not trivial to obtain its equation. We proceed by recalling the form of the entropy four-current for ideal fluids, (30), and extending it for dissipative fluids:(52)Sμ=Pβμ+βνTμν-αNμ.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 [34, 35]. 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 (48) in (52). Taking the divergence and using (49)–(51), we obtain(53)∂μSμ=-βΠθ-nμ∇μα+βπμνσμν.The relativistic Navier-Stokes theory can then be obtained by applying the second law of thermodynamics to each fluid element; that is, by requiring that the entropy production∂μSμ must always be positive,(54)-βΠθ-nμ∇μα+βπμνσμν≥0.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(55)Π=-ζθ,nμ=κ∇μα,πμν=2ησμν,where the proportionality coefficients ζ, κ, and η refer to the bulk viscosity, the particle diffusion, and the shear viscosity, respectively. Substituting the above equation in (53), we observe that the source term for entropy production becomes a quadratic function of the dissipative currents:(56)∂μSμ=βζΠ2-1κnμnμ+β2ηπμνπμν.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, πμνπμν>0. Hence, as long as ζ,κ,η≥0, the entropy production is always positive. Constitutive relations for the dissipative quantities, (55), along with (49)–(51) are known as the relativistic Navier-Stokes equations.The relativistic Navier-Stokes theory in this form was obtained originally by Landau and Lifshitz [35]. A similar theory was derived independently by Eckart, using a different definition of the fluid four-velocity [34]. 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, (55). 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 [39, 45, 48], the divergence-type theory [49, 50], extended irreversible thermodynamics [51], Carter’s theory [52], Grmela and Ottinger theory [53], 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.


## 2.4.5. 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 nonequilibrium entropy four-current [39, 45, 48]. In equilibrium, the entropy four-current was expressed exactly in terms of the primary fluid dynamical variables, (30). Strictly speaking, the nonequilibrium entropy four-current should depend on a larger number of independent dynamical variables, and a direct extension of (30) to (52) 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(57)Sμ=Pβμ+βνTμν-αNμ-QμδNμ,δTμν,where Qμ is a function of deviations from local equilibrium, δNμ≡Nμ-N(0)μ, and δTμν≡Tμν-T(0)μν. Using (48) and Taylor-expanding Qμ to second order in dissipative fluxes, we obtain(58)Sμ=suμ-αnμ-β0Π2-β1nνnν+β2πρσπρσuμ2T-α0ΠΔμν+α1πμνnνT+Oδ3,where O(δ3) denotes third-order terms in the dissipative currents and β0,β1,β2,α0, and α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 (58) 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 recalculate the entropy production, ∂μSμ, using the more general entropy four-current given in (58):(59)∂μSμ=-βΠθ+β0Π˙+βΠΠΠθ+ψαnΠnμu˙μ+α0∇μnμ+ψαΠnnμ∇μα-βnμT∇μα-β1n˙μ-βnnnμθ+α0∇μΠ+α1∇νπμν+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπμν∇να+χ~αnππμνu˙ν+βπμνσμν-β2π˙μν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν.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:(60)Π=-ζθ+β0Π˙+βΠΠΠθ+α0∇μnμ+ψαnΠnμu˙μ+ψαΠnnμ∇μα,nμ=λT∇μα-β1n˙μ-βnnnμθ+α0∇μΠ+α1Δρμ∇νπρν+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπνμ∇να+χ~αnππνμu˙ν,πμν=2ησμν-β2π˙μν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν,where λ≡κ/T. This implies that the dissipative currents must satisfy the dynamical equations:(61)Π˙+ΠτΠ=-1β0θ+βΠΠΠθ+ψαnΠnμu˙μ+α0∇μnμ+ψαΠnnμ∇μα,n˙μ+nμτn=1β1T∇μα-βnnnμθ+α1Δρμ∇νπρν+α0∇μΠ+ψ~αnΠΠu˙μ+ψ~αΠnΠ∇μα+χ~απnπνμ∇να+χ~αnππνμu˙ν,π˙μν+πμντπ=1β2σμν-βππθπμν-α1∇μnν-χαπnnμ∇να-χαnπnμu˙ν.The above equations for the dissipative quantities are relaxation-type equations with the relaxation times defined as(62)τΠ≡ζβ0,τn≡λβ1=κβ1T,τπ≡2ηβ2.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, that is, 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.


## 2.5. 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 single-particle 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 [54].Let us consider a system of relativistic particles, each having rest massm, momentum p→, and energy p0. Therefore, from relativity, we have p0=(p→)2+m2. 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μ=(p0,p→) at each space-time point such that f(x,p)Δ3xΔ3p gives the average number of particles at a given time t in the volume element Δ3x 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 Δ3x is large, on the other hand, Δ3x is small compared to macroscopic point of view.The particle densityn(x) is introduced to describe, in general, a nonuniform system, such that n(x)Δ3x is the average number of particles in volume Δ3x 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(63)nx=∫d3pfx,p,j→x=∫d3pv→fx,p,where v→=p→/p0 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(64)Nμx=∫d3pp0pμfx,p.Note that since d3p/p0 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 isp0, the average energy density and the energy flow can be written in terms of the distribution function as(65)T00x=∫d3pp0fx,p,T0ix=∫d3pp0vifx,p.The momentum density is defined as the average value of particle momenta pi, 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(66)Ti0x=∫d3ppifx,p,Tijx=∫d3ppivjfx,p.Combining all these in a compact covariant form using vi=pi/p0, we obtain the energy-momentum tensor of a macroscopic system:(67)Tμνx=∫d3pp0pμpνfx,p.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.H-function introduced by Boltzmann implies that the nonequilibrium local entropy density of a system can be written as(68)sx=-∫d3pfx,plnfx,p-1.The entropy flow corresponding to the above entropy density is(69)S→x=-∫d3pv→fx,plnfx,p-1.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(70)Sμx=-∫d3pp0pμfx,plnfx,p-1.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(71)Sμx=-∫d3pp0pμfx,plnfx,p+rf~x,plnf~x,p,where f~≡1-rf. The expressions for the entropy four-current given in (70) and (71) 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=f0+δf. The equilibrium distribution function f0 is defined as(72)f0x,p=1expβu·p-α+r,where the scalar product is defined as u·p≡uμpμ and r=0 for Maxwell-Boltzmann statistics. Note that, in equilibrium, that is, for f(x,p)=f0(x,p), the particle four-flow and energy-momentum tensor given in (64) and (67) reduce to those of ideal hydrodynamics N(0)μ and T(0)μν. Therefore, using (48), the dissipative quantities, namely, the bulk viscous pressure Π, the particle diffusion current nμ, and the shear stress tensor πμν, can be written as(73)Π=-13Δαβ∫d3pp0pαpβδf,nμ=Δμν∫d3pp0pνδf,πμν=Δαβμν∫d3pp0pαpβδf.The evolution equations for the dissipative quantities expressed in terms of the nonequilibrium distribution function, (73), 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(74)pμ∂μf=Cf,where dp≡d3p/p0 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(75)Cf=12∫dp′dkdk′Wpp′→kk′fkfk′f~pf~p′-fpfp′f~kf~k′,where Wpp′→kk′ is the collisional transition rate. The first and second terms within the integral of (75) 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 [55](76)Cf=-u·pδfτR,where τR is the relaxation time.


## 2.6. Dissipative Fluid Dynamics from Kinetic Theory

The derivation of a causal theory of relativistic dissipative hydrodynamics by Israel and Stewart [39] proceeds by invoking the second law of thermodynamics, namely, ∂μSμ≥0, from the algebraic form of the entropy four-current given in (58). 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 (70) and (71), in order to obtain the dissipative equations [56]. This essentially ensures that the nonequilibrium corrections to the distribution function, δf, do not violate the second law of thermodynamics. In [56], the generalized method of moments developed by Denicol et al. [57] was used to quantify the dissipative corrections to the distribution function. The form of the resultant dissipative equations, obtained in [56], is identical to (61), 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 [48], has been used quite extensively to quantify the dissipative corrections to the distribution function [56–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 [66]. In this method, which is based on a Chapman-Enskog like expansion, the Boltzmann equation in the relaxation time approximation,(77)pμ∂μf=-u·pτRf-f0,is solved iteratively to obtain δf up to any arbitrary order in derivatives. To first order and second order in gradients, one obtains(78)δf1=-τRu·ppμ∂μf0,(79)δf2=τRu·ppμpν∂μτRu·p∂νf0.This method of obtaining the form of the nonequilibrium distribution function is consistent with dissipative hydrodynamics, which is also formulated as a gradient expansion.The second-order evolution equations for the dissipative quantities are then obtained by substitutingδf=δf(1)+δf(2) from (78) and (79) in (73):(80)ΠτΠ=-Π˙-βΠθ-δΠΠΠθ+λΠππμνσμν-τΠnn·u˙-λΠnn·∇α-lΠn∂·n,nμτn=-n˙μ+βn∇μα-nνωνμ-λnnnνσνμ-δnnnμθ+λnΠΠ∇μα-λnππμν∇να-τnππνμu˙ν+τnΠΠu˙μ+lnπΔμν∂γπνγ-lnΠ∇μΠ,πμντπ=-π˙μν+2βπσμν+2πγμωνγ-τπππγμσνγ-δπππμνθ+λπΠΠσμν-τπnnμu˙ν+λπnnμ∇να+lπn∇μnν,where ωμν=(∇μuν-∇νuμ)/2 is the vorticity tensor. It is interesting to note that although the form of the evolution equations for dissipative quantities in (80) is identical to those obtained in [58] using the moment method, the transport coefficients are, in general, different [67, 68]. Moreover, it was shown that the above described method, based on iterative solution of Boltzmann equation, leads to phenomenologically consistent corrections to the distribution function [69] and the transport coefficients exhibit intriguing similarities with strongly coupled conformal field theory [70, 71].Proceeding in a similar way, a third-order dissipative evolution equation can also be obtained [72–74]:(81)π˙μν=-πμντπ+2βπσμν+2πγμωνγ-107πγμσνγ-43πμνθ+257βππρμωνγπργ-13βππγμπνγθ-38245βππμνπργσργ-2249βππρμπνγσργ-2435∇μπνγu˙γτπ+435∇μτπ∇γπνγ-27∇γτπ∇μπνγ+127∇γτπu˙μπνγ-17∇γτπ∇γπμν+67∇γτπu˙γπμν-27τπωρμωνγπργ-27τππρμωνγωργ-1063τππμνθ2+2621τππγμωνγθ.It is reassuring that the results obtained using third-order evolution equation indicate convergence of the gradient expansion and show improvement over second order, when compared to the direct solutions of the Boltzmann equation [72–75].Apart from these standard formulations, there are several other formulations of relativistic dissipative hydrodynamics from kinetic theory. Among them, the ones which have gained widespread interest are anisotropic hydrodynamics and derivations based on renormalization group method. Anisotropic hydrodynamics is a nonperturbative reorganization of the standard relativistic hydrodynamics which takes into account the large momentum space anisotropies generated in ultrarelativistic heavy-ion collisions [76–81]. On the other hand, the derivation based on renormalization group method attempts to solve the Boltzmann equation, as faithfully as possible, in an organized perturbation scheme and resum the possible secular terms by a suitable setting of the initial value of the distribution function [82–85]. Since it is widely accepted that the QGP is momentum space anisotropic, application of anisotropic hydrodynamics to high energy heavy-ion collisions has phenomenological implications. Nevertheless, the dissipative hydrodynamic formulation based on renormalization group method is important in order to accurately determine the higher-order transport coefficients.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, do 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 is different from systems 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(82)πμν=2ησμν-τππ˙μν+43πμνθ-λ1η2πγμπνγ+λ2ηπγμωνγ+λ3ωγμωνγ.For a system of massless particles, the kinetic theory results for second-order transport coefficients agree in the case of both moment method [57] and the Chapman-Enskog like iterative solution of the Boltzmann equation [66]. In Table 1, 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.Table 1
Comparison of transport coefficients.
τ π T / ( η / s ) λ 1 T / ( η / s ) λ 2 T / ( η / s ) λ 3 ADS/CFT 2 ( 2 - ln ⁡ 2 ) 2 η 4 η ln ⁡ 2 0 KT 5 ( 25 / 7 ) η 10 η 0


## 3. 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 [32, 89–94].In high energy heavy-ion collisions, bunch of nuclei of heavy elements are accelerated inside the beam pipes and in the final state (after the collisions) we have hundreds 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 are 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 preequilibrium stage is essential input in the viscous hydrodynamics models. The knowledge of distribution of energy/entropy density and the thermalization time is one of the uncertainties present in the current hydrodynamics model studies. Below we discuss four most popular initial condition models used in hydrodynamics simulation of heavy-ion collisions.


### 3.1. 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 nucleus-nucleus collisions at low energies. The original idea of Glauber was further modified by Białłas et al. [95] to explain inelastic nuclear collisions. For a nice and more complete review of Glauber model see [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 variations 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 nuclei sees 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 numberA as(83)ρAx,y,z=ρ01+expr-R0/a0, where R0 and a0 are the nuclear radius and skin thickness parameter of the nucleus and ρ0 is an overall constant that is determined by requiring ∫d3xρA(x,y,z)=A. One additionally defines the “thickness function” [96],(84)TAx,y=∫-∞∞dzρAx,y,z,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 (Npart) and number of binary collisions (Ncoll) for the two colliding nuclei.In order to calculateNcoll(x,y,b→) and Npart(x,y,b→) from Glauber model, one can choose x-axis along the impact parameter vector b→. Npart and Ncoll distributions 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](85)Ncollx,y;b=σinTAx+b2,yTBx-b2,y,(86)Npartx,y;b=TAx+b2,yFTBx-b2,y,B+TBx-b2,yFTAx+b2,y,A,where(87)FTAx,y,A≡1-1-σinTAx,yAA. Here σin is the inelastic nucleon-nucleon cross section whose value depends on sNN and is obtained from the experimental data.The distribution ofNpart and Ncoll 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(88)ϵx,y=ϵ0αNcollx,y+1-αNpartx,y2,where ϵ0 is a multiplicative constant used to fix the charged hadron multiplicity; α is the fraction of hard scattering [97]. The energy density which 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 is sampled according to the nuclear density distribution by using statistical Monte Carlo (MC) method. The collisions between two nucleons occur when the distance between them becomes equal to 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](89)ϵx,y=K∑WN,BC12πσ2exp-x-xi2+y-yi22σ2,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 participants Npart, and BC represents number of binary collisions. With this short discussion we now move on to the next topic.


### 3.2. Color-Glass-Condensate

The Color-Glass-Condensate (CGC) model takes into account the nonlinear 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,rgl, of a gluon is related to its momentum, Q, via |rgl|×|Q|~ħ=1. Therefore the cross section of gluon-gluon interactions is(90)σ~αsQ2πrgl2~αsQ2πQ2, 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/(πR02), where R0 is the radius of the nucleus. Gluons start interacting with each other when the scattering probability becomes of the order unity:(91)AπR02σ=αsQ2AR02Q2~1. This indicates that there exists a typical momentum scale Qs2=αsA/R02 separating perturbative (Q2≫Qs2) and nonperturbative (Q2≪Qs2) 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 Q2 in relativistic heavy-ion collisions [99, 100].The presence of nonabelian plasma instabilities [101–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]:(92)ϵx⊥,b=const×dNgd2xTdYxT,b4/3.Here Ng is the number of gluons produced in the collision whose momentum distribution is given by(93)dNgd2xTdY~∫d2pTpT2∫pTd2kTαskTϕ+pT+kT24;xTϕ-pT-kT24;xT,ϕ±kT2;xT=Qs21-x4αsQs2maxQs2,kT2npartAx⊥,±bTAx±b/2,y,Qs2x,x⊥=2TA2x±b/2,yGeV2npartAx⊥,±bfm21.530.01x0.288, where x=pT/s. It is important to note that the Glauber and CGC models lead to different values of spatial eccentricity, defined by(94)exb=y2-x2ϵy2+x2ϵ,where 〈〉ϵ represents averaging over the transverse plane with weight ϵ(x⊥,b). It has been observed that the CGC model typically has 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–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, Qs, which is computed from the longitudinally projected density of wounded nucleons. The positions of the wounded nucleons are sampled according to (86) 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.


## 3.1. 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 nucleus-nucleus collisions at low energies. The original idea of Glauber was further modified by Białłas et al. [95] to explain inelastic nuclear collisions. For a nice and more complete review of Glauber model see [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 variations 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 nuclei sees 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 numberA as(83)ρAx,y,z=ρ01+expr-R0/a0, where R0 and a0 are the nuclear radius and skin thickness parameter of the nucleus and ρ0 is an overall constant that is determined by requiring ∫d3xρA(x,y,z)=A. One additionally defines the “thickness function” [96],(84)TAx,y=∫-∞∞dzρAx,y,z,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 (Npart) and number of binary collisions (Ncoll) for the two colliding nuclei.In order to calculateNcoll(x,y,b→) and Npart(x,y,b→) from Glauber model, one can choose x-axis along the impact parameter vector b→. Npart and Ncoll distributions 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](85)Ncollx,y;b=σinTAx+b2,yTBx-b2,y,(86)Npartx,y;b=TAx+b2,yFTBx-b2,y,B+TBx-b2,yFTAx+b2,y,A,where(87)FTAx,y,A≡1-1-σinTAx,yAA. Here σin is the inelastic nucleon-nucleon cross section whose value depends on sNN and is obtained from the experimental data.The distribution ofNpart and Ncoll 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(88)ϵx,y=ϵ0αNcollx,y+1-αNpartx,y2,where ϵ0 is a multiplicative constant used to fix the charged hadron multiplicity; α is the fraction of hard scattering [97]. The energy density which 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 is sampled according to the nuclear density distribution by using statistical Monte Carlo (MC) method. The collisions between two nucleons occur when the distance between them becomes equal to 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](89)ϵx,y=K∑WN,BC12πσ2exp-x-xi2+y-yi22σ2,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 participants Npart, and BC represents number of binary collisions. With this short discussion we now move on to the next topic.


## 3.2. Color-Glass-Condensate

The Color-Glass-Condensate (CGC) model takes into account the nonlinear 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,rgl, of a gluon is related to its momentum, Q, via |rgl|×|Q|~ħ=1. Therefore the cross section of gluon-gluon interactions is(90)σ~αsQ2πrgl2~αsQ2πQ2, 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/(πR02), where R0 is the radius of the nucleus. Gluons start interacting with each other when the scattering probability becomes of the order unity:(91)AπR02σ=αsQ2AR02Q2~1. This indicates that there exists a typical momentum scale Qs2=αsA/R02 separating perturbative (Q2≫Qs2) and nonperturbative (Q2≪Qs2) 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 Q2 in relativistic heavy-ion collisions [99, 100].The presence of nonabelian plasma instabilities [101–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]:(92)ϵx⊥,b=const×dNgd2xTdYxT,b4/3.Here Ng is the number of gluons produced in the collision whose momentum distribution is given by(93)dNgd2xTdY~∫d2pTpT2∫pTd2kTαskTϕ+pT+kT24;xTϕ-pT-kT24;xT,ϕ±kT2;xT=Qs21-x4αsQs2maxQs2,kT2npartAx⊥,±bTAx±b/2,y,Qs2x,x⊥=2TA2x±b/2,yGeV2npartAx⊥,±bfm21.530.01x0.288, where x=pT/s. It is important to note that the Glauber and CGC models lead to different values of spatial eccentricity, defined by(94)exb=y2-x2ϵy2+x2ϵ,where 〈〉ϵ represents averaging over the transverse plane with weight ϵ(x⊥,b). It has been observed that the CGC model typically has 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–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, Qs, which is computed from the longitudinally projected density of wounded nucleons. The positions of the wounded nucleons are sampled according to (86) 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.


## 4. Preequilibrium 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 preequilibrium phase, but this is in contrast 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 take into account the preequilibrium dynamics until hydrodynamics sets in.


### 4.1. IP-Glasma

In the IP-Glasma model, the initial conditions are 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/Qs. The initial Glasma fields are then evolved using the classical Yang-Mills (CYM) equation. One of the most important features 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 follow a negative binomial distribution.The color charges,ρa(x-,x⊥), in the IP-Sat model behave as local sources for small-x classical gluon Glasma fields. The classical gluon fields are then determined by solving the classical Yang-Mills equations:(95)Dμ,Fμν=Jν.The color current in the above equation, generated by a nucleus A(B) moving along x+(x-) direction, is given by(96)Jν=δν±ρABx∓,x⊥,where the upper indices are for nucleus A.It is easy to solve (95) in Lorentz gauge, ∂μAμ=0, where the equation transforms into a two-dimensional Poisson equation:(97)-∇⊥2AAB±=ρABx∓,x⊥. The solution of the above equation can be written as(98)AAB±=-ρABx∓,x⊥/∇⊥2. Using the path-ordered exponential(99)VABx⊥=Pexp-ig∫dx-ρABx-,x⊥∇T2+m2,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–114](100)AABix⊥=θx-x+igVABx⊥∂iVAB†x⊥,A-A+=0.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 τ,η coordinates are defined as τ=2x+x- and η=0.5ln(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](101)Ai=AAi+ABi,Aη=ig2AAi,ABi,∂τAi=0,∂τAη=0.In the limit τ→0, Aη=-Eη/2, where Eη is the longitudinal component of the electric field. At τ=0, one can nonperturbatively calculate the longitudinal magnetic and electric fields, which are the only nonvanishing components of the field strength tensor. These fields determine the energy density of the Glasma at each transverse position in a single event [117–119].The Glasma fields are then evolved in time numerically according to (95), 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 Tfluidμν=(ϵ+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μTCYMμν=εuν.


### 4.2. Transport: AMPT and UrQMD

In [121–123], a different approach was taken in order to incorporate the preequilibrium dynamics for obtaining the initial condition of hydrodynamics evolution. While the authors of [121] employ ultrarelativistic quantum molecular dynamics (UrQMD) string dynamics model, a multiphase transport model (AMPT) was used in [122, 123] to simulate the preequilibrium 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 [122] for calculating initial conditions for a (3+1)D hydrodynamics evolution with the parton transport in the preequilibrium phase. Additional benefit for choosing this type of initial condition is that one naturally incorporates 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 [122], a multiphase transport model (AMPT) [124, 125] was used to obtain the local initial energy-momentum 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 semihard 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 calculate number of participants 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 evolve 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(102)Tμντ0,x,y,ηs=K∑i=1Npiμpiνpiτ12πτ0σr22πσηs2exp-x-xi2+y-yi22σr2-ηs-ηis22σηs2,where piτ=miTcosh⁡Yi-ηis and pix,y=px,yi; piη=miTsinh⁡Yi-ηis/τ0 are the four-momenta of the ith parton and Yi, ηis, and miT 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 [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 are obtained from the calculated Tμν via a root finding method which is used as an input to the subsequent hydrodynamics evolution; see [122] for further details.


### 4.3. Numerical Relativity: AdS/CFT

Another method to simulate the preequilibrium stage is via numerical relativity solutions to AdS/CFT [128]. In this method, one employs the dynamics of the energy-momentum tensor of the strongly coupled conformal field theory (CFT) on the boundary using the gravitational field in the bulk of AdS5. 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 [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, Ttt=δ(t+z)TA(x,y). Here the thickness function, TA(x,y), is the same as defined in (84) but with 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=x2+y2, and tan⁡θ=y/x, the energy density, fluid velocity, and pressure anisotropy were found up to leading order in t [130]:(103)ϵ=2TA2ρτ2,uρ=-TA′ρ3TAρτ,PLPT=-32,where in the local rest frame Tνμ=diag(-ϵ,PT,PT,PL) [131–134]. One finds that the corresponding line element ds2 turns out to be ξ-independent (boost-invariant), up to leading order in τ, and can be written as(104)ds2=-Adτ2+Σ2e-B-Cdξ2+eBdρ2+eCdθ2+2drdτ+2Fdρdτ. 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,(105)Br,τ,ρ⟶B0r,τ,ρ+∑i=06biτ,ρr-i1+σ7r-7,where B0 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 (103) 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 pseudospectral method based on [135–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 (103). 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 are therefore determined using an early-time, far-from-equilibrium dynamics, modeled as a strongly coupled CFT described by gravity in AdS. Recently, a nonconformal extension has also been studied in order to incorporate bulk viscosity [138].


## 4.1. IP-Glasma

In the IP-Glasma model, the initial conditions are 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/Qs. The initial Glasma fields are then evolved using the classical Yang-Mills (CYM) equation. One of the most important features 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 follow a negative binomial distribution.The color charges,ρa(x-,x⊥), in the IP-Sat model behave as local sources for small-x classical gluon Glasma fields. The classical gluon fields are then determined by solving the classical Yang-Mills equations:(95)Dμ,Fμν=Jν.The color current in the above equation, generated by a nucleus A(B) moving along x+(x-) direction, is given by(96)Jν=δν±ρABx∓,x⊥,where the upper indices are for nucleus A.It is easy to solve (95) in Lorentz gauge, ∂μAμ=0, where the equation transforms into a two-dimensional Poisson equation:(97)-∇⊥2AAB±=ρABx∓,x⊥. The solution of the above equation can be written as(98)AAB±=-ρABx∓,x⊥/∇⊥2. Using the path-ordered exponential(99)VABx⊥=Pexp-ig∫dx-ρABx-,x⊥∇T2+m2,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–114](100)AABix⊥=θx-x+igVABx⊥∂iVAB†x⊥,A-A+=0.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 τ,η coordinates are defined as τ=2x+x- and η=0.5ln(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](101)Ai=AAi+ABi,Aη=ig2AAi,ABi,∂τAi=0,∂τAη=0.In the limit τ→0, Aη=-Eη/2, where Eη is the longitudinal component of the electric field. At τ=0, one can nonperturbatively calculate the longitudinal magnetic and electric fields, which are the only nonvanishing components of the field strength tensor. These fields determine the energy density of the Glasma at each transverse position in a single event [117–119].The Glasma fields are then evolved in time numerically according to (95), 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 Tfluidμν=(ϵ+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μTCYMμν=εuν.


## 4.2. Transport: AMPT and UrQMD

In [121–123], a different approach was taken in order to incorporate the preequilibrium dynamics for obtaining the initial condition of hydrodynamics evolution. While the authors of [121] employ ultrarelativistic quantum molecular dynamics (UrQMD) string dynamics model, a multiphase transport model (AMPT) was used in [122, 123] to simulate the preequilibrium 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 [122] for calculating initial conditions for a (3+1)D hydrodynamics evolution with the parton transport in the preequilibrium phase. Additional benefit for choosing this type of initial condition is that one naturally incorporates 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 [122], a multiphase transport model (AMPT) [124, 125] was used to obtain the local initial energy-momentum 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 semihard 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 calculate number of participants 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 evolve 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(102)Tμντ0,x,y,ηs=K∑i=1Npiμpiνpiτ12πτ0σr22πσηs2exp-x-xi2+y-yi22σr2-ηs-ηis22σηs2,where piτ=miTcosh⁡Yi-ηis and pix,y=px,yi; piη=miTsinh⁡Yi-ηis/τ0 are the four-momenta of the ith parton and Yi, ηis, and miT 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 [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 are obtained from the calculated Tμν via a root finding method which is used as an input to the subsequent hydrodynamics evolution; see [122] for further details.


## 4.3. Numerical Relativity: AdS/CFT

Another method to simulate the preequilibrium stage is via numerical relativity solutions to AdS/CFT [128]. In this method, one employs the dynamics of the energy-momentum tensor of the strongly coupled conformal field theory (CFT) on the boundary using the gravitational field in the bulk of AdS5. 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 [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, Ttt=δ(t+z)TA(x,y). Here the thickness function, TA(x,y), is the same as defined in (84) but with 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=x2+y2, and tan⁡θ=y/x, the energy density, fluid velocity, and pressure anisotropy were found up to leading order in t [130]:(103)ϵ=2TA2ρτ2,uρ=-TA′ρ3TAρτ,PLPT=-32,where in the local rest frame Tνμ=diag(-ϵ,PT,PT,PL) [131–134]. One finds that the corresponding line element ds2 turns out to be ξ-independent (boost-invariant), up to leading order in τ, and can be written as(104)ds2=-Adτ2+Σ2e-B-Cdξ2+eBdρ2+eCdθ2+2drdτ+2Fdρdτ. 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,(105)Br,τ,ρ⟶B0r,τ,ρ+∑i=06biτ,ρr-i1+σ7r-7,where B0 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 (103) 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 pseudospectral method based on [135–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 (103). 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 are therefore determined using an early-time, far-from-equilibrium dynamics, modeled as a strongly coupled CFT described by gravity in AdS. Recently, a nonconformal extension has also been studied in order to incorporate bulk viscosity [138].


## 5. 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, contain one additional variable compared to the number of equations. EoS closes the system of equations by providing another functional relationship and it is one of the important inputs to hydrodynamics. For a relativistic simple fluid the acceleration under a given pressure gradient ∇P is governed by the following relationship:(106)Duμ=-1ϵ+P∇P, 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 [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 (volumeV→∞) 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:(107)ϵ=T2V∂∂Tln⁡ZT,V,P=T∂∂Vln⁡ZT,V. The pressure for a homogeneous system of infinite extent can be simply expressed in terms of f as(108)P=TVln⁡ZT,V. Using the above relations one can arrive at the following relationships:(109)ϵ=T∂P∂T-P,(110)ΘT=T∂∂TPT4,where the trace anomaly Θ(T)=(ϵ-3P)/T4.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>250 MeV) and use HRG model in the low temperature region (T<180 MeV). In the intermediate temperature Θ(T) is obtained by joining the parametrized high temperature and low temperature values smoothly (continuous first and second derivatives). Once Θ(T) is known, pressure can be calculated by using (110). The energy density then can be readily obtained from (109). Figure 3(a) 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]. Figure 3(b) shows the trace anomaly calculated in lattice QCD compared with the HRG model given by the solid and dashed lines [139].Figure 3
The trace anomaly calculated in lattice QCD with p4 and asqtad actions onNτ=6 and 8 lattices compared with various parametrizations given by the solid, dotted, and dashed lines (a) and the trace anomaly calculated in lattice QCD compared with the HRG model given by the solid and dashed lines (b). The figures are taken from [139].
(a) (b)


## 6. Freeze-Out: Spectra and Flow

In the late stage of hydrodynamics evolution of hot and dense nuclear matter created in high energy heavy-ion collisions, the density and the temperature reach a critical value when the constituents no longer collide among themselves and thereafter they move in a straight trajectory towards the detectors. This phenomenon is known as freeze-out, more precisely the kinetic freeze-out. When the particle number changing processes cease, the system is said to have attained chemical freeze-out. The chemical freeze-out temperature is 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 hypersurface. This can be done by more than one way; we discuss here only the most popular choices used to calculate the freeze-out hypersurface, namely, (i) the constant temperature freeze-out, (ii) the constant energy density freeze-out, and (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 freeze-out 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 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 [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 [142] for details. In Cooper-Frye procedure the momentum distribution (or invariant yield) of hadrons is calculated as [142](111)Ed3Nd3p=d3Nd2pTdy=∫Σfx,ppμdΣμ,where E, N, and pμ are energy, number, and four-momentum of hadrons; dΣμ is the differential freeze-out hypersurface element. The distribution function, f(x,p), consists of an equilibrium part, f0(x,p), and dissipative corrections, δf(x,p). While the equilibrium distribution corresponding to the local thermodynamic quantities, as given in (72), 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 Grad’s 14-moment approximation, leads to the well-known form of the viscous correction [39, 48]:(112)δfx,p=f0f~02ϵ+PT2pαpβπαβ,where f~0≡1-rf0, 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, (78), leads to a viscous correction which is effectively linear in momenta [69]:(113)δfx,p=5f0f~02ϵ+PT1u·ppαpβπαβ.It has been shown that, in contrast to (112) obtained using moment method, (113) 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 hypersurface and the numerical evaluation of it is not trivial; for example, see [143] for more details. Once we know the invariant momentum distribution the “nth” order Fourier coefficient the flow harmonics vn can be readily obtained as(114)vn=∫ydy∫pTd2pTd3N/d2pTdycos⁡nϕn-ψn∫ydy∫pTd2pTd3N/d2pTdy.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.


## 7. Resonance Decay and Hadronic Rescattering

In high energy nuclear collisions various hadronic resonances are formed. The lifetime of most of the resonance particles is of the order of the expansion lifetime of the nuclear matter. The end product for the most of the decay channels involves pions. The decay of hadron resonances to pion enhances the pion yield especially at low transverse momentum,pT. One can use the formalism given in [144] to calculate the relative contribution of the resonance decay to thermal pion spectra. The relative contribution of the resonance decay to pion spectra is a function of both the freeze-out temperature, Tfo, and pT. Thus the final pT spectra of π are obtained by adding the contribution from resonance decay to the thermal pT spectra calculated from Cooper-Frye formula. The most dominant hadronic decay channels contributing to pion yield are ρ±→π±π0, ρ0→π-π+, K∗±→π±K0, K∗0→π-K+, Δ→π±,0N, ω→π+π-π0, and η→π+π-π0, which should be considered with their corresponding branching ratios [144].According to the formalism given in [144], to calculate the pion contribution from resonances, one needs 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 temperatures Tfo=130 MeV and 150 MeV is approximately given by [145] (115)πtotal±πthermal±Tfo=130MeV=1.0121+1.40281+pT/mπ-0.0964/3.6662,πtotal±πthermal±Tfo=150MeV=1.0252+3.04951+pT/mπ-0.2302/2.7922,where mπ=139 MeV is the pion mass. Note that about ~50% of the total pion yield comes from resonance decay at LHC energy (sNN=2.76 TeV), whereas for RHIC energy (sNN=200 GeV) the resonance contribution to total pion yield is ~30% for Tfo=130 MeV.The sudden conversion of fluid to noninteracting hadrons at the freeze-out hypersurface 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 is 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 heavy-ion 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–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–154]. One of the first numerical codes, VISHNU, which couples (2+1)D viscous hydro with a late hadronic Boltzmann cascade appeared in 2011 [155]. The use of these more sophisticated hybrid models is 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 [155] and references therein. Before finishing this section we should point out one of the major findings of [155]; η/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 constrain η/s of QGP by using such sophisticated numerical models is a topic of current ongoing research.


## 8. Transport Coefficients

Determination of transport coefficients of the hot and dense QCD matter is one of the primary goals of theoretical simulations of relativistic heavy-ion collisions; see [156] for a recent review. Ideal hydrodynamics has been proven to be quite successful 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 (59), 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, v2, defined in (114) 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 v2. 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.Figure4(a) shows the extracted values of η/s in different model calculations for Au-Au collisions at sNN=200 GeV [25, 26, 158–166]. Most of the estimates are obtained by comparing experimental data for elliptic flow with model calculations. Some of the estimates used pT correlations and heavy meson v2 data. The theoretical calculations include simulations with transport based approach as well as (2+1)D and (3+1)D viscous hydrodynamics with various initial conditions. Also shown are the results from lattice QCD calculation. All these results indicate that η/s value of the QGP fluid produced at top RHIC energies lies within 1–5 ×  1/4π and is below η/s value of helium (blue dashed line) at Tc. The spread in the estimated values of η/s reflects the current uncertainties associated with the theoretical calculations.Figure 4
Extracted values ofη/s for Au-Au collision at sNN=200 GeV (a) and for Pb-Pb collision at sNN=2.76 TeV (b) by different model calculations using different experimental observables. The solid vertical line at the left shows the lower limit of η/s in unit of 1/(4π) [30]. For comparison we have also shown η/s of Helium at Tc (dashed vertical line on the right).
(a) (b)Figure4(b) shows η/s estimated in various model calculations for Pb-Pb collision at sNN=2.76 TeV [167–170]. All the model calculations indicate that the value of η/s of the QCD matter formed in heavy-ion collision at LHC lies within 1–4 ×  (1/4π). The specific shear viscosity was obtained in reference [167] by using a multiphase 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 v2(pT) and pT integrated v2 for different centralities. Their calculation shows that the experimental data measured at LHC by the ALICE collaboration are best described for η/s value of 0.08 or smaller. Luzum and Romatschke [170] have estimated η/s by using a (2+1)D viscous hydrodynamics simulation with smooth initial conditions for LHC energy to be the same as at RHIC, η/s=0.1±0.1(theory)±0.08(experiment). Comparison of experimentally measured integrated and differential v2, the charged hadron pT spectra, and multiplicity in the mid-rapidity and their global fit by minimising χ2 was done in [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 which include the effect of bulk viscosity, where it has been demonstrated that bulk viscosity can have a nonnegligible effect on heavy-ion observables [175–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 make 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.


## 9. Recent Developments

### 9.1. Flow in Small Systems: Proton-Proton and Proton-Nucleus Collisions

As mentioned in 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 [181] for a summary of recent experimental results. At this point, why the observation of flow in small system is remarkable needs some explanation. One of the fundamental assumptions 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 is large; for example, for head-on Au-Au or Pb-Pb collisions, there are 197 + 197 and 208 + 208 participating nucleons. Each of these colliding nucleons on average produces more than one particle in each collision; thus the total number of degrees of freedom in the system created just after the collisions is large. They collide among themselves through strong interaction and subsequently reach 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. 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 HIJING 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 because of 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 thermalized medium is believed to be created.This situation changed recently as the CMS and ATLAS 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, mass dependence of the slope of identified hadron’smT spectrum in high multiplicity p-Pb and p-p collisions was also observed. All these phenomena 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 this observed experimental measurement for high multiplicity p-p or p-Pb events unless it employs some special mechanism like Color Reconnection (CR) and Multiparton Interaction (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 proven 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 observations which do not incorporate this hydrodynamics like flow, but till now those studies lack detailed numerical calculation in order to compare it with the experimental data [183].As for (i) p-Pb collisions, recent experimental measurement shows that the number of charged particles produced in p-Pb collisions atsNN=5.02 TeV is similar to that in peripheral Pb-Pb collisions [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 to that 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 in 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 pT spectra, measured in experiment [185]. The experimentally measured vn and v2(pT) data for identified hadrons in p-Pb collisions by CMS [184] and ALICE [186] collaboration are nicely explained by a (3+1)D viscous hydrodynamics model study by Bożek and Broniowski in [187] (see Figure 5). In addition to 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 Figure 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.Figure 5
Elliptic flow,v2, of identified particles from the hydrodynamic model compared to the ALICE data [186]. The figure is taken from [187].Figure 6
Average transverse momentum of identified particles in p-Pb collisions atsNN=5.02 TeV; the experimental data (various symbols) from ALICE collaboration [186], compared to the results of the HIJING model and of the viscous hydrodynamics. The figure is taken from [187].
(a)
3+1D hydro (b)
HIJINGRegarding (ii) 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 thatpT integrated v2 and v3 are 30% and 50% smaller compared to 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 Figure 7 which is taken from [188]. In a similar effort a blast wave model fit was shown to be inconsistent with the experimental data; see [189] for details. There are some studies of viscous hydrodynamics for p-p collisions; for example, see [190–192]; 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 systems 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 [193] for a detailed discussion about the applicability of hydrodynamics in small systems.Figure 7
Normalized spectra of pions (squares), kaons (triangles), and protons (discs) for p-p collisions. Open symbols correspond to the CMS data [194] for |η|<2.4 and s=7 TeV, while the solid ones are obtained from the best one-parameter fit of Gubser’s flow. The figure is taken from [188].


### 9.2. 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 nonzero values of flow harmonicsvn. The most prominent flow harmonics v2 originate 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 ultracentral collisions. In ultracentral collisions v2 and other higher flow harmonics solely originate from the initial fluctuating energy density, since the overlap zone in ultracentral collisions is almost circular.For ultracentral 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 ofvn in ultracentral heavy-ion collisions can help reduce the systematic uncertainties of initial-state modeling in extracting η/s value of the system. Let us first discuss the recent experimental results for ultracentral Pb-Pb collisions at LHC; after that we shall also discuss the corresponding results from viscous hydrodynamics simulations. Figure 8(a) shows the experimentally measured differential flow coefficients vn as a function of pT for 0–0.2% centrality Pb-Pb collisions at sNN=2.76 TeV. vn’s were calculated using 2 particle correlation methods with large pseudorapidity gap |Δη|>2 between the two hadrons. Figure 8(b) shows pT integrated vn (n= 2–7) in ultracentral Pb-Pb collisions for five different collisions’ centrality. The experimental data and the figure are taken from [195]. Before we proceed any further we note the following experimental observation from the CMS paper:(i)
At higher transverse momentum (pT≥2 GeV), v2 becomes even smaller than the higher-order v3,v4, and at much higher values of pT it becomes smaller than other higher-order vn.(ii)
pT averaged v2 and v3 are found to be equal within 2%, while other higher-order vn decrease as n increases.Figure 8
(a)vn (n=2–6) values as a function of pT in 0.0–0.2% central Pb+Pb collisions at sNN=2.76 TeV. (b) Experimental measurement of pT averaged (0.3–3.0 GeV) vn as a function of “n" in five centrality classes (2.5–5.0%, 0–2.5%, 0-1%, 0–0.2%, and 0–0.02%) for Pb+Pb sNN=2.76 TeV collisions. Error bars denote the statistical uncertainties, while the shaded color boxes correspond to the systematic uncertainties. The data was measured by the CMS collaboration and the figures are taken from [195].
(a) (b)The evolution of the QGP according to relativistic hydrodynamics simulations has been able to consistently explain experimentally measuredvn’s for different centralities and for different colliding energies, it is natural to expect that it should also explain measured vn in the ultracentral collisions. Before we discuss the results of hydrodynamic simulations, we note that one needs to carefully select events into centrality classes, since integrated vn’s are quite sensitive on the selection of centrality class as can be seen from Figure 8(b). We also note that it is computationally expensive to simulate such ultracentral 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 pT dependence of charged hadrons vn and their observed ordering for ultracentral Pb-Pb collisions was nicely explained by a viscous hydrodynamic simulation using initial conditions from AMPT model [123] (see Figure 9), on careful observation we notice that, at low pT<1.5 GeV, the splitting from hydrodynamics simulation is larger than the corresponding experimental measurement. Similar disagreements are also evident for pT>1.5 GeV in Figure 10, which is taken from [196]. This can be seen more clearly from pT integrated v2 and v3 in Figure 11 which is also taken from [196]. In [196], pT integrated vn was studied using (2+1)D viscous hydrodynamics model with MC-Glauber and MC-KLN initial conditions.Figure 9
Comparison of a (2+1)D viscous hydrodynamics simulation with the initial condition from AMPT model to corresponding experimental data for root mean squared values ofvn for n = 2–5; figure is from [123].Figure 10
Comparison ofpT dependent vn2 of charged hadrons in 2.76 A TeV Pb+Pb collisions at 0–0.2% centrality for viscous hydrodynamics simulations (various lines) with the corresponding experimental results (solid squares); the figure is from [196].
(a) (b) (c) (d) (e) (f)Figure 11
Comparison ofpT integrated vn2 of charged hadrons in 2.76 A TeV Pb+Pb collisions at 0–0.2% centrality for viscous hydrodynamics simulations (various lines) with the corresponding experimental results (solid squares). Plots (a) and (c) are for MC-Glauber initial conditions, and (b) and (d) correspond to MC-KLN initial conditions. Results in (a) and (b) were obtained by considering nucleons with a repulsive hard core, whereas the results in (c) and (d) are obtained for the initial conditions with finite nucleon-nucleon correlations. The figure is taken from [196].
(a) (b) (c) (d)The nucleon-nucleon correlations in the colliding nucleus were also considered as a potential cause behind the experimentally measuredv2~v3. However, none of the initial condition models has so far been able to simultaneously explain the experimentally measured vn’s, as can be seen in Figures 9, 10, and 11. In this regard, we note that Denicol et al. [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 theoretical results, although there was some improvement, but so far the effort remains unfruitful.


### 9.3. 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 axes to the beam direction) to the experimentally observed nonzero odd flow harmonics (particularly third harmonicsv3) in heavy-ion collision was made by Alver and Roland [198]. From then on experimentally measured flow harmonics for all order (even and odd) have 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 give rise to not only odd flow harmonics but also significant even and odd vn in ultracentral collisions [205]. They also result in pT dependent event planes, which break down the flow factorization vn,n(pT1,pT2)=vn(pT1)vn(pT2) [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:(i)
The fluctuations of energy density along the longitudinal direction due to the fragmentation and different lengths of the colored string produced in the scattering of nucleons [207–209](ii)
A gradual twist of the fireball (or more specifically the event plane) along the longitudinal direction [210, 211] Let us discuss each of them separately. Regarding the contribution of color string, we shall particularly discuss here a recent study [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 semihard scatterings and excited strings from soft interactions. The number of mini-jet partons per binary nucleon-nucleon collision in hard and semihard scatterings follows 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 giving rise to fluctuating energy density distribution in ηs; see Figure 12 for the distribution of colored strings in the longitudinal direction for a typical Pb-Pb collision at sNN=2.76 TeV. For more details about the longitudinal fluctuations and the visualization of parton density distribution in ηs, see [207].Figure 12
Distribution of strings created between the partons of two colliding Pb nuclei as a function of space-time rapidity (ηs) at sNN=2.76 TeV for 0-1% (a) and 40–50% (b) collision centrality. The figure is taken from [207].
(a) (b)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)(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 rapidity 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 Figure 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 are 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 these two effects in heavy-ion collisions was studied within (3+1)D hydrodynamics model and AMPT in [208].Figure 13
(a) The factorization ratior2 as a function of space-time rapidity ηa for two different reference rapidity bins 3.0<ηb<4.0 and 4.4<ηb<5.0 in Pb-Pb collisions at sNN=2.76 TeV (open and solid diamonds) and for 2.5<ηb<3.0 and 3.0<ηb<4.0 in Au-Au collisions at sNN=200 GeV (open and solid circles) from event-by-event (3+1)D ideal hydrodynamics simulations compared with experimental data from CMS collaboration [212] for Pb-Pb collisions at sNN=2.76 TeV (empty and solid squares). (b) shows the factorization ratio r3. The figures are from [207].
(a) (b)Many techniques have been proposed to study the longitudinal structure of final hadron production in heavy-ion 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 characterize 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 and Huo [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 [208]. Before ending this section we note that the experimentally observed difference in the longitudinal correlation (r2 and r3) 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.


### 9.4. Flow in Intense Magnetic Field

The strongest known magnetic field (|B→|~1018–1019 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 (sNN) of the colliding nucleons [218, 219]. The corresponding electric fields in such collisions also become very strong which is same order of magnitude as the magnetic field (eB→≈eE→~10mπ2 for a typical Au-Au collision at top RHIC energy sNN=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(116)E→r→,t=e4π∑i=1NprotonZiR→i-Riv→iRi-R→i·v→i31-vi2,B→r→,t=e4π∑i=1NprotonZiv→i×R→iRi-R→i·v→i31-vi2.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, Figure 14 shows the positions of nucleons inside the two Au nuclei 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 nonzero for such initial random proton positions. This can be seen from Figure 15, where the event averaged value of By as a function of impact parameter b is shown. The black dashed-dotted line corresponds to the average of the absolute magnitude of By which is clearly nonzero even for b=0 fm collisions. Note that we have used the natural unit, where ħ=c=kb=ε0=μ0=1; with this choice the electric charge e=4π/137 becomes a dimensionless number. In the limit v~c, the denominator in (116) becomes very small and we have large E→ and B→.Figure 14
Distribution of nucleons inside target and projectile nuclei in a typical Au-Au collision atsNN=200 GeV for b=12 fm. Figure is from [229].Figure 15
Impact parameter dependence of event averaged magnetic and electric fields at the centre of the fireball for Au-Au collisions atsNN=200 GeV. Figure is from [229].There are a 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, and chiral magnetic waves, the discussion of which is beyond the scope of this review. For more details, we refer the reader to the following references: [222–226]. Here we will concentrate on the possible effect of this large electromagnetic 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, two of the first numerical studies of the effect of magnetic fields on the hydrodynamics evolution in heavy-ion collisions are by Gürsoy 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 magnetohydrodynamics 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 atsNN=200 GeV in [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 [232–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 solutions that exist for relativistic magnetohydrodynamics. The only possible way is to use numerical methods to solve magnetohydrodynamic equations relevant to heavy-ion collisions. This is not an easy task to accomplish. Initial effort in this direction can be found in [235]. However, we note that the authors of [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 atsNN=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 Figure 16, we show v2 of π- obtained for impact parameter b=10 fm Au-Au collisions at sNN=200 GeV. 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 [230]. A realistic spatial profile for y-component of magnetic field was considered for the simulation. From Figure 16, one can see that v2 of π- is noticeably enhanced in the presence of magnetic field (blue dashed line) compared to the case of no magnetic field (red line).Figure 16
Elliptic flow ofπ- as a function of pT for b=10 fm Au-Au collisions at sNN=200 GeV. Red line corresponds to (2+1)D ideal hydrodynamics result without magnetic field; the blue dashed line corresponds to v2 with external magnetic field.


## 9.1. Flow in Small Systems: Proton-Proton and Proton-Nucleus Collisions

As mentioned in 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 [181] for a summary of recent experimental results. At this point, why the observation of flow in small system is remarkable needs some explanation. One of the fundamental assumptions 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 is large; for example, for head-on Au-Au or Pb-Pb collisions, there are 197 + 197 and 208 + 208 participating nucleons. Each of these colliding nucleons on average produces more than one particle in each collision; thus the total number of degrees of freedom in the system created just after the collisions is large. They collide among themselves through strong interaction and subsequently reach 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. 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 HIJING 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 because of 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 thermalized medium is believed to be created.This situation changed recently as the CMS and ATLAS 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, mass dependence of the slope of identified hadron’smT spectrum in high multiplicity p-Pb and p-p collisions was also observed. All these phenomena 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 this observed experimental measurement for high multiplicity p-p or p-Pb events unless it employs some special mechanism like Color Reconnection (CR) and Multiparton Interaction (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 proven 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 observations which do not incorporate this hydrodynamics like flow, but till now those studies lack detailed numerical calculation in order to compare it with the experimental data [183].As for (i) p-Pb collisions, recent experimental measurement shows that the number of charged particles produced in p-Pb collisions atsNN=5.02 TeV is similar to that in peripheral Pb-Pb collisions [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 to that 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 in 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 pT spectra, measured in experiment [185]. The experimentally measured vn and v2(pT) data for identified hadrons in p-Pb collisions by CMS [184] and ALICE [186] collaboration are nicely explained by a (3+1)D viscous hydrodynamics model study by Bożek and Broniowski in [187] (see Figure 5). In addition to 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 Figure 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.Figure 5
Elliptic flow,v2, of identified particles from the hydrodynamic model compared to the ALICE data [186]. The figure is taken from [187].Figure 6
Average transverse momentum of identified particles in p-Pb collisions atsNN=5.02 TeV; the experimental data (various symbols) from ALICE collaboration [186], compared to the results of the HIJING model and of the viscous hydrodynamics. The figure is taken from [187].
(a)
3+1D hydro (b)
HIJINGRegarding (ii) 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 thatpT integrated v2 and v3 are 30% and 50% smaller compared to 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 Figure 7 which is taken from [188]. In a similar effort a blast wave model fit was shown to be inconsistent with the experimental data; see [189] for details. There are some studies of viscous hydrodynamics for p-p collisions; for example, see [190–192]; 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 systems 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 [193] for a detailed discussion about the applicability of hydrodynamics in small systems.Figure 7
Normalized spectra of pions (squares), kaons (triangles), and protons (discs) for p-p collisions. Open symbols correspond to the CMS data [194] for |η|<2.4 and s=7 TeV, while the solid ones are obtained from the best one-parameter fit of Gubser’s flow. The figure is taken from [188].


## 9.2. 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 nonzero values of flow harmonicsvn. The most prominent flow harmonics v2 originate 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 ultracentral collisions. In ultracentral collisions v2 and other higher flow harmonics solely originate from the initial fluctuating energy density, since the overlap zone in ultracentral collisions is almost circular.For ultracentral 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 ofvn in ultracentral heavy-ion collisions can help reduce the systematic uncertainties of initial-state modeling in extracting η/s value of the system. Let us first discuss the recent experimental results for ultracentral Pb-Pb collisions at LHC; after that we shall also discuss the corresponding results from viscous hydrodynamics simulations. Figure 8(a) shows the experimentally measured differential flow coefficients vn as a function of pT for 0–0.2% centrality Pb-Pb collisions at sNN=2.76 TeV. vn’s were calculated using 2 particle correlation methods with large pseudorapidity gap |Δη|>2 between the two hadrons. Figure 8(b) shows pT integrated vn (n= 2–7) in ultracentral Pb-Pb collisions for five different collisions’ centrality. The experimental data and the figure are taken from [195]. Before we proceed any further we note the following experimental observation from the CMS paper:(i)
At higher transverse momentum (pT≥2 GeV), v2 becomes even smaller than the higher-order v3,v4, and at much higher values of pT it becomes smaller than other higher-order vn.(ii)
pT averaged v2 and v3 are found to be equal within 2%, while other higher-order vn decrease as n increases.Figure 8
(a)vn (n=2–6) values as a function of pT in 0.0–0.2% central Pb+Pb collisions at sNN=2.76 TeV. (b) Experimental measurement of pT averaged (0.3–3.0 GeV) vn as a function of “n" in five centrality classes (2.5–5.0%, 0–2.5%, 0-1%, 0–0.2%, and 0–0.02%) for Pb+Pb sNN=2.76 TeV collisions. Error bars denote the statistical uncertainties, while the shaded color boxes correspond to the systematic uncertainties. The data was measured by the CMS collaboration and the figures are taken from [195].
(a) (b)The evolution of the QGP according to relativistic hydrodynamics simulations has been able to consistently explain experimentally measuredvn’s for different centralities and for different colliding energies, it is natural to expect that it should also explain measured vn in the ultracentral collisions. Before we discuss the results of hydrodynamic simulations, we note that one needs to carefully select events into centrality classes, since integrated vn’s are quite sensitive on the selection of centrality class as can be seen from Figure 8(b). We also note that it is computationally expensive to simulate such ultracentral 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 pT dependence of charged hadrons vn and their observed ordering for ultracentral Pb-Pb collisions was nicely explained by a viscous hydrodynamic simulation using initial conditions from AMPT model [123] (see Figure 9), on careful observation we notice that, at low pT<1.5 GeV, the splitting from hydrodynamics simulation is larger than the corresponding experimental measurement. Similar disagreements are also evident for pT>1.5 GeV in Figure 10, which is taken from [196]. This can be seen more clearly from pT integrated v2 and v3 in Figure 11 which is also taken from [196]. In [196], pT integrated vn was studied using (2+1)D viscous hydrodynamics model with MC-Glauber and MC-KLN initial conditions.Figure 9
Comparison of a (2+1)D viscous hydrodynamics simulation with the initial condition from AMPT model to corresponding experimental data for root mean squared values ofvn for n = 2–5; figure is from [123].Figure 10
Comparison ofpT dependent vn2 of charged hadrons in 2.76 A TeV Pb+Pb collisions at 0–0.2% centrality for viscous hydrodynamics simulations (various lines) with the corresponding experimental results (solid squares); the figure is from [196].
(a) (b) (c) (d) (e) (f)Figure 11
Comparison ofpT integrated vn2 of charged hadrons in 2.76 A TeV Pb+Pb collisions at 0–0.2% centrality for viscous hydrodynamics simulations (various lines) with the corresponding experimental results (solid squares). Plots (a) and (c) are for MC-Glauber initial conditions, and (b) and (d) correspond to MC-KLN initial conditions. Results in (a) and (b) were obtained by considering nucleons with a repulsive hard core, whereas the results in (c) and (d) are obtained for the initial conditions with finite nucleon-nucleon correlations. The figure is taken from [196].
(a) (b) (c) (d)The nucleon-nucleon correlations in the colliding nucleus were also considered as a potential cause behind the experimentally measuredv2~v3. However, none of the initial condition models has so far been able to simultaneously explain the experimentally measured vn’s, as can be seen in Figures 9, 10, and 11. In this regard, we note that Denicol et al. [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 theoretical results, although there was some improvement, but so far the effort remains unfruitful.


## 9.3. 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 axes to the beam direction) to the experimentally observed nonzero odd flow harmonics (particularly third harmonicsv3) in heavy-ion collision was made by Alver and Roland [198]. From then on experimentally measured flow harmonics for all order (even and odd) have 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 give rise to not only odd flow harmonics but also significant even and odd vn in ultracentral collisions [205]. They also result in pT dependent event planes, which break down the flow factorization vn,n(pT1,pT2)=vn(pT1)vn(pT2) [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:(i)
The fluctuations of energy density along the longitudinal direction due to the fragmentation and different lengths of the colored string produced in the scattering of nucleons [207–209](ii)
A gradual twist of the fireball (or more specifically the event plane) along the longitudinal direction [210, 211] Let us discuss each of them separately. Regarding the contribution of color string, we shall particularly discuss here a recent study [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 semihard scatterings and excited strings from soft interactions. The number of mini-jet partons per binary nucleon-nucleon collision in hard and semihard scatterings follows 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 giving rise to fluctuating energy density distribution in ηs; see Figure 12 for the distribution of colored strings in the longitudinal direction for a typical Pb-Pb collision at sNN=2.76 TeV. For more details about the longitudinal fluctuations and the visualization of parton density distribution in ηs, see [207].Figure 12
Distribution of strings created between the partons of two colliding Pb nuclei as a function of space-time rapidity (ηs) at sNN=2.76 TeV for 0-1% (a) and 40–50% (b) collision centrality. The figure is taken from [207].
(a) (b)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)(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 rapidity 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 Figure 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 are 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 these two effects in heavy-ion collisions was studied within (3+1)D hydrodynamics model and AMPT in [208].Figure 13
(a) The factorization ratior2 as a function of space-time rapidity ηa for two different reference rapidity bins 3.0<ηb<4.0 and 4.4<ηb<5.0 in Pb-Pb collisions at sNN=2.76 TeV (open and solid diamonds) and for 2.5<ηb<3.0 and 3.0<ηb<4.0 in Au-Au collisions at sNN=200 GeV (open and solid circles) from event-by-event (3+1)D ideal hydrodynamics simulations compared with experimental data from CMS collaboration [212] for Pb-Pb collisions at sNN=2.76 TeV (empty and solid squares). (b) shows the factorization ratio r3. The figures are from [207].
(a) (b)Many techniques have been proposed to study the longitudinal structure of final hadron production in heavy-ion 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 characterize 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 and Huo [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 [208]. Before ending this section we note that the experimentally observed difference in the longitudinal correlation (r2 and r3) 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.


## 9.4. Flow in Intense Magnetic Field

The strongest known magnetic field (|B→|~1018–1019 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 (sNN) of the colliding nucleons [218, 219]. The corresponding electric fields in such collisions also become very strong which is same order of magnitude as the magnetic field (eB→≈eE→~10mπ2 for a typical Au-Au collision at top RHIC energy sNN=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(116)E→r→,t=e4π∑i=1NprotonZiR→i-Riv→iRi-R→i·v→i31-vi2,B→r→,t=e4π∑i=1NprotonZiv→i×R→iRi-R→i·v→i31-vi2.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, Figure 14 shows the positions of nucleons inside the two Au nuclei 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 nonzero for such initial random proton positions. This can be seen from Figure 15, where the event averaged value of By as a function of impact parameter b is shown. The black dashed-dotted line corresponds to the average of the absolute magnitude of By which is clearly nonzero even for b=0 fm collisions. Note that we have used the natural unit, where ħ=c=kb=ε0=μ0=1; with this choice the electric charge e=4π/137 becomes a dimensionless number. In the limit v~c, the denominator in (116) becomes very small and we have large E→ and B→.Figure 14
Distribution of nucleons inside target and projectile nuclei in a typical Au-Au collision atsNN=200 GeV for b=12 fm. Figure is from [229].Figure 15
Impact parameter dependence of event averaged magnetic and electric fields at the centre of the fireball for Au-Au collisions atsNN=200 GeV. Figure is from [229].There are a 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, and chiral magnetic waves, the discussion of which is beyond the scope of this review. For more details, we refer the reader to the following references: [222–226]. Here we will concentrate on the possible effect of this large electromagnetic 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, two of the first numerical studies of the effect of magnetic fields on the hydrodynamics evolution in heavy-ion collisions are by Gürsoy 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 magnetohydrodynamics 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 atsNN=200 GeV in [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 [232–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 solutions that exist for relativistic magnetohydrodynamics. The only possible way is to use numerical methods to solve magnetohydrodynamic equations relevant to heavy-ion collisions. This is not an easy task to accomplish. Initial effort in this direction can be found in [235]. However, we note that the authors of [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 atsNN=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 Figure 16, we show v2 of π- obtained for impact parameter b=10 fm Au-Au collisions at sNN=200 GeV. 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 [230]. A realistic spatial profile for y-component of magnetic field was considered for the simulation. From Figure 16, one can see that v2 of π- is noticeably enhanced in the presence of magnetic field (blue dashed line) compared to the case of no magnetic field (red line).Figure 16
Elliptic flow ofπ- as a function of pT for b=10 fm Au-Au collisions at sNN=200 GeV. Red line corresponds to (2+1)D ideal hydrodynamics result without magnetic field; the blue dashed line corresponds to v2 with external magnetic field.


## 10. 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 need further investigation. For example, the experimentally measured longitudinal correlation of flow harmonics shows that splitting in the quantities corresponds to the correlation measure, namely,r2ηa,ηb and r3η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 r2ηa,ηb and r3η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 sources 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 measuredv2 and v3 in ultracentral collisions. Within error bars the magnitudes of v2 and v3 are observed to be same for ultracentral 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 behaves collectively is still not well understood. There is also possibility for nonhydrodynamical 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 provides us with 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.



---

*Source: 9623034-2016-10-27.xml*