#### 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.

#### 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].

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 Figure 2 [20].

Assuming that thermalization is achieved in the early stages of heavy-ion collisions and that the interaction between the quarks is strong enough to maintain local thermodynamic equilibrium during the subsequent expansion, the time evolution of the QGP and hadronic matter can be described by the laws of fluid dynamics [21–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 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 . Kovtun, Son, and Starinets (KSS) conjectured this strong coupling result to be the absolute lower bound for all fluids; that is, [29, 30]. This specific combination of hydrodynamic quantities, , 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 remains finite [31]. Similar result for a lower bound also follows from the quantum mechanical uncertainty principle. The kinetic theory prediction for viscosity, , 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, . For a weakly interacting relativistic Bose gas, the entropy per particle is given by . This leads to which is very close to the lower KSS bound.

Indeed, the estimated 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 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 Section 2 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 , with , where is the Planck constant, is the velocity of light in vacuum, and is the Boltzmann constant. Unless stated otherwise, the space-time is always taken to be Minkowskian, where the metric tensor is given by . Apart from Minkowskian coordinates , we will also regularly employ Milne coordinate system or , with proper time , the radial coordinate , the azimuthal angle , and space-time rapidity . Hence, , , , and . For the coordinate system , the metric becomes , whereas for , the metric is . 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 and , respectively. The notation represents scalar product of a covariant and a contravariant four-vector. Tensors without indices shall always correspond to Lorentz scalars. We follow Einstein summation convention, which states that repeated indices in a single term are implicitly summed over all the values of that index.

We denote the fluid four-velocity by and the Lorentz contraction factor by . The projector onto the space orthogonal to is defined as . Hence, satisfies the conditions with trace . The partial derivative can then be decomposed asIn the fluid rest frame, reduces to the time derivative and reduces to the spacial gradient. Hence, the notation is also commonly used. We also frequently use the symmetric, antisymmetric, and angular brackets notations defined as where is the traceless symmetric projection operator orthogonal to satisfying the conditions .

#### 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 , the total energy , entropy , and number of particles , 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 operations

The first law of thermodynamics postulates that the changes in the total energy of a thermodynamic system must result from (1) heat exchange, (2) the mechanical work done by an external force, and (3) particle exchange with an external medium. Hence, the conservation law relating the small changes in state variables, , , and iswhere and are the pressure and chemical potential, respectively, and 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, , and is defined in terms of the heat exchange as , with the temperature being the proportionality constant. Then, when considering variations between equilibrium states that are infinitesimally close to each other, it is possible to write the first law of thermodynamics in terms of differentials of the state variables:Hence, using (5), the intensive quantities, , and , can be obtained in terms of partial derivatives of the entropy as

The entropy is mathematically defined as an extensive and additive function of the state variables, which means thatDifferentiating both sides with respect to , we obtainwhich holds for any arbitrary value of . Setting and using (6), we obtain the so-called Euler’s relation:Using Euler’s relation, (9), along with the first law of thermodynamics, (5), we arrive at the Gibbs-Duhem relation:

In terms of energy, entropy, and number of densities defined as , , and , respectively, Euler’s relation, (9), and Gibbs-Duhem relation, (10), reduce toDifferentiating (11) and using (12), we obtain the relation analogous to first law of thermodynamics:It is important to note that all the densities defined above 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 , there can be assigned a temperature , a chemical potential , and a collective four-velocity field:The proper time increment is given by the line elementwhere . This implies thatwhere . In the nonrelativistic limit, we obtain . It is important to note that the four-vector only contains three independent components, since it obeys the relationThe quantities , , and 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, , and the (net) particle four-current, . 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, , and the energy-momentum tensor, , the (net) particle four-current, , and the entropy four-current, , should have the characteristic form of a system in static equilibrium. In other words, in local rest frame, there is no flow of energy , the force per unit surface element is isotropic , and there is no particle and entropy flow ( and ). Consequently, the energy-momentum tensor, particle four-current, and entropy four-current in this frame take the following simple forms:

For an ideal relativistic fluid, the general form of the energy-momentum tensor, , (net) particle four-current, , and the entropy four-current, , has to be built out of the hydrodynamic tensor degrees of freedom, namely, the vector, , and the metric tensor, . Since should be symmetric and transform as a tensor and and should transform as a vector, under Lorentz transformations, the most general form allowed is thereforeIn the local rest frame, . Hence, in this frame, (19) takes the formBy comparing the above equation with the corresponding general expressions in the local rest frame, (18), one obtains the following expressions for the coefficients:The conserved currents of an ideal fluid can then be expressed aswhere is the projection operator on the three-space orthogonal to and satisfies the following properties of an orthogonal projector:

The dynamical description of an ideal fluid is obtained using the conservation laws of energy, momentum, and (net) particle number. These conservation laws can be mathematically expressed using the four-divergences of energy-momentum tensor and particle four-current which leads to the following equations:where the partial derivative transforms as a covariant vector under Lorentz transformations. Using the four-velocity, , and the projection operator, , the derivative, , can be projected along and orthogonal to :Projection of energy-momentum conservation equation along and orthogonal to together with the conservation law for particle number leads to the equations of motion of ideal fluid dynamics:where . It is important to note that an ideal fluid is described by four fields, , , , and , 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, , 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 Section 2.1, (11), (12), and (13), in a covariant form [39, 45]. For this purpose, it is convenient to introduce the following notations:In these notations, the covariant version of Euler’s relation, (11), and the Gibbs-Duhem relation, (12), can be postulated asrespectively. The above equations can then be used to derive a covariant form of the first law of thermodynamics, (13):

The covariant thermodynamic relations were constructed in such a way that when (30), (31), and (32) are contracted with , we obtain the usual thermodynamic relations, (11), (12), and (13). Here we have used the property of the fluid four-velocity, . The projection of (30), (31), and (32) on the three-space orthogonal to just leads to trivial identities:From the above equations, we conclude that the covariant thermodynamic relations do not contain more information than the usual thermodynamic relations.

The first law of thermodynamics, (32), leads to the following expression for the entropy four-current divergence:After employing the conservation of energy-momentum and net particle number, (24), the above equation leads to the conservation of entropy, . 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 asWe 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 () 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: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 are added to the previously derived ideal currents, and :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 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]:These matching conditions enforce the following constraints on the dissipative currents:Subsequently, using 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 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 ,with the following properties:The parentheses in the above equation denote symmetrization of the Lorentz indices; that is, . The dissipative current then can be tensor decomposed in its irreducible form by using , , and aswhere we have definedThe scalar is the bulk viscous pressure, the vector 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 and are orthogonal to and, additionally, is traceless. Armed with these definitions, all the irreducible hydrodynamic fields are expressed in terms of and aswhere the angular bracket notations are defined as and .

We observe that is a symmetric second-rank tensor with ten independent components and is a four-vector; overall they have fourteen independent components. Next we count the number of independent components in the tensor decompositions of and . Since and are orthogonal to , they can have only three independent components each. The shear stress tensor is symmetric, traceless, and orthogonal to and hence can have only five independent components. Together with , , , and , which have in total six independent components ( 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 was introduced as a general normalized four-vector and was not specified. Hence 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 fix . In the case of ideal fluids, the local rest frame was defined as the frame in which there is, simultaneously, no net energy and particle flow. While the definition of local rest frame was unambiguous for ideal fluids, this definition is no longer possible in the case of dissipative fluids due to the presence of both energy and particle diffusion. From a mathematical perspective, the fluid velocity can be defined in numerous ways. However, from the physical perspective, there are two natural choices: the* Eckart definition* [34], in which the velocity is defined by the flow of particles,and the* Landau definition* [35], in which the velocity is specified by the flow of the total energy,

We note that the above two definitions of 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:

As done for ideal fluids, the energy-momentum conservation equation in (37) is decomposed parallel and orthogonal to . Using (48) together with the conservation law for particle number in (37) leads to the equations of motion for dissipative fluids. For , , and , one obtainsrespectively. Here , and the shear tensor .

We observe that while there are fourteen total independent components of and , (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 , , and .

###### 2.4.4. Relativistic Navier-Stokes Theory

In the presence of dissipative currents, the entropy is no longer a conserved quantity; that is, . Since the form of the entropy four-current for a dissipative fluid is not known* a priori*, it is not trivial to obtain its equation. We proceed by recalling the form of the entropy four-current for ideal fluids, (30), and extending it for dissipative fluids:The above extension remains valid because an artificial equilibrium state was constructed using the matching conditions to satisfy the thermodynamic relations as if in equilibrium. This was the key step proposed by Eckart, Landau, and Lifshitz in order to derive the relativistic Navier-Stokes theory [34, 35]. The next step is to calculate the entropy generation, , in dissipative fluids. To this end, we substitute the form of and for dissipative fluids from (48) in (52). Taking the divergence and using (49)–(51), we obtain

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 must always be positive,The above inequality can be satisfied for all possible fluid configurations if one assumes that the bulk viscous pressure , the particle-diffusion four-current , and the shear stress tensor are linearly proportional to , , and , respectively. This leads towhere 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:In the above equation, since is orthogonal to the timelike four-vector , it is spacelike and hence . Moreover, is symmetric in its Lorentz indices and in the local rest frame . Since the trace of the square of a symmetric matrix is always positive, . Hence, as long as , 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 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 bywhere is a function of deviations from local equilibrium, , and . Using (48) and Taylor-expanding to second order in dissipative fluxes, we obtainwhere denotes third-order terms in the dissipative currents and , and 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, , using the more general entropy four-current given in (58):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, , is guaranteed to be satisfied if we impose linear relationships between thermodynamical fluxes and extended thermodynamic forces, leading to the following evolution equations for bulk pressure, particle-diffusion four-current, and shear stress tensor:where . This implies that the dissipative currents must satisfy the dynamical equations:The above equations for the dissipative quantities are relaxation-type equations with the relaxation times defined asSince the relaxation times must be positive, the Taylor expansion coefficients , , and 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, , , , , and , 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 mass , momentum , and energy . Therefore, from relativity, we have . For a large number of particles, we introduce a single-particle distribution function which gives the distribution of the four-momentum at each space-time point such that gives the average number of particles at a given time in the volume element at point with momenta in the range . However, this definition of the single-particle phase-space distribution function assumes that while, on one hand, the number of particles contained in is large, on the other hand, is small compared to macroscopic point of view.

The particle density is introduced to describe, in general, a nonuniform system, such that is the average number of particles in volume at . Similarly, particle flow is defined as the particle current. With the help of the distribution function, the particle density and particle flow are given bywhere is the particle velocity. These two local quantities, particle density and particle flow, constitute a four-vector field , called particle four-flow, and can be written in a unified way asNote that since is a Lorentz invariant quantity, should be a scalar in order that transforms as a four-vector.

Since the energy per particle is , the average energy density and the energy flow can be written in terms of the distribution function asThe momentum density is defined as the average value of particle momenta , and the momentum flow or pressure tensor is defined as the flow in direction of momentum in direction . For these two quantities, we haveCombining all these in a compact covariant form using , we obtain the energy-momentum tensor of a macroscopic system:Observe that the above definition of the energy-momentum tensor corresponds to second moment of the distribution function, and, hence, it is a symmetric quantity.

-function introduced by Boltzmann implies that the nonequilibrium local entropy density of a system can be written asThe entropy flow corresponding to the above entropy density isThese two local quantities, entropy density and entropy flow, constitute a four-vector field , called entropy four-flow, and can be written in a unified way asThe 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 () or Bose-Einstein statistics () aswhere . 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, can be written as . The equilibrium distribution function is defined aswhere the scalar product is defined as and for Maxwell-Boltzmann statistics. Note that, in equilibrium, that is, for , the particle four-flow and energy-momentum tensor given in (64) and (67) reduce to those of ideal hydrodynamics and . Therefore, using (48), the dissipative quantities, namely, the bulk viscous pressure , the particle diffusion current , and the shear stress tensor , can be written as

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 aswhere and is the collision functional. For microscopic interactions restricted to elastic collisions, the form of the collision functional is given bywhere is the collisional transition rate. The first and second terms within the integral of (75) refer to the processes and , 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]where 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, , from the algebraic form of the entropy four-current given in (58). As noted earlier, the new parameters, , , , , and , 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, , 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,is solved iteratively to obtain up to any arbitrary order in derivatives. To first order and second order in gradients, one obtainsThis 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 from (78) and (79) in (73):where 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]: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 equationFor 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 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.

#### 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 number as where and are the nuclear radius and skin thickness parameter of the nucleus and is an overall constant that is determined by requiring . One additionally defines the “thickness function” [96],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 ( and number of binary collisions () for the two colliding nuclei.

In order to calculate and from Glauber model, one can choose -axis along the impact parameter vector . and 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 “” and “,” the transverse density of binary collision and wounded nucleon profile is given by [96]where Here is the inelastic nucleon-nucleon cross section whose value depends on and is obtained from the experimental data.

The distribution of and 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 bywhere 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]where is a free parameter controlling the width of the Gaussian. Typical values for this fluctuation size parameter are of the order of fm, is the abbreviation for wounded nucleons which is same as number of participants , and 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, , of a gluon is related to its momentum, , via ~. Therefore the cross section of gluon-gluon interactions is where 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 . The density of gluons in the transverse plane is then given by , where is the radius of the nucleus. Gluons start interacting with each other when the scattering probability becomes of the order unity: This indicates that there exists a typical momentum scale separating perturbative () and nonperturbative () 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 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]:Here is the number of gluons produced in the collision whose momentum distribution is given by where . It is important to note that the Glauber and CGC models lead to different values of spatial eccentricity, defined bywhere represents averaging over the transverse plane with weight . 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, , 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 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, . 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, , in the IP-Sat model behave as local sources for small- classical gluon Glasma fields. The classical gluon fields are then determined by solving the classical Yang-Mills equations:The color current in the above equation, generated by a nucleus moving along direction, is given bywhere the upper indices are for nucleus .

It is easy to solve (95) in Lorentz gauge, , where the equation transforms into a two-dimensional Poisson equation: The solution of the above equation can be written as Using the path-ordered exponentialone can gauge transform the results of Lorentz gauge to light-cone gauge, . The pure gauge fields are then given by [112–114]The discontinuity in the fields on the light cone corresponds to the localized valence charge source [115].

The initial condition for a heavy-ion collision, at time , is determined by the solution of the CYM equations in Fock-Schwinger gauge , where coordinates are defined as and . The Fock-Schwinger gauge is a natural gauge choice because it interpolates between the light-cone gauge conditions of the incoming nuclei. In terms of the gauge fields of the colliding nuclei, one obtains [115, 116]In the limit , , where is the longitudinal component of the electric field. At , 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 , 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 from the energy density in the fluid’s rest frame and the flow velocity . The local pressure at each transverse position is obtained using an equation of state. The hydrodynamic quantities and are obtained by solving the Landau frame condition, .

##### 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 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 aswhere and ; are the four-momenta of the th parton and , , and are the momentum rapidity, the spatial rapidity, and the transverse mass of the th parton, respectively. Unless otherwise stated, the smearing parameters are taken as fm and from [122], where the soft hadron spectra, rapidity distribution, and elliptic flow can be well described. The sum index runs over all produced partons () in a given nucleus-nucleus collision. The scale factor and the initial proper time 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 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 AdS_{5}. Therefore, a relativistic nucleus may be described using a gravitational shockwave in AdS, whereby the energy-momentum tensor of a nucleus can be exactly matched [129]. For a central collision, the dynamics of the colliding shockwaves has been solved near the boundary of AdS in [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, . Here the thickness function, , is the same as defined in (84) but with extra normalization, , which is used to match the experimentally observed particle multiplicity, .

In terms of the polar Milne coordinates with , , , and , the energy density, fluid velocity, and pressure anisotropy were found up to leading order in [130]:where in the local rest frame [131–134]. One finds that the corresponding line element turns out to be -independent (boost-invariant), up to leading order in , and can be written as Here all functions depend on , , and the fifth AdS space dimension only. In this scenario, the space boundary is located at , where the induced metric is given by .

The metric is then expanded near the boundary,