Causal dissipative hydrodynamics for QGP fluid in 2+1 dimensions
A. K. Chaudhuri∗
Variable Energy Cyclotron Centre, 1/AF, Bidhan Nagar, Kolkata 700 064, India
(Dated: February 21, 2013)
In 2nd order causal dissipative theory, space-time evolution of QGP fluid is studied in 2+1 di-
mensions. Relaxation equations for shear stress tensors are solved simultaneously with the energy-
momentum conservation equations. Comparison of evolution of ideal and viscous QGP fluid, ini-
tialized under the same conditions, e.g. same equilibration time, energy density and velocity profile,
indicate that in a viscous dynamics, energy density or temperature of the fluid evolve slowly, than
in an ideal fluid. Cooling gets slower as viscosity increases. Transverse expansion also increases in
a viscous dynamics. For the first time we have also studied elliptic flow of ’quarks’ in causal viscous
dynamics. It is shown that elliptic flow of quarks saturates due to non-equilibrium correction to
equilibrium distribution function, and can not be mimicked by an ideal hydrodynamics.
PACS numbers: 47.75.+f, 25.75.-q, 25.75.Ld
I. INTRODUCTION
One of the most important discoveries in Relativistic
Heavy ion collider (RHIC) at Brokhaven National Lab-
oratory is the large elliptic flow in non-central Au+Au
collisions [1, 2, 3, 4] . Elliptic flow measures the momen-
tum anisotropy of produced particles and is quantified by
the 2nd harmonic of the azimuthal distribution,
v2(pT ) =< cos(2φ) >=
dyd2pT
cos(2φ)dφ
dyd2pT
(1.1)
Elliptic flow is naturally explained in hydrodynamics.
Hydrodynamic pressure is built up from rescattering of
secondaries, and pressure gradients drive the subsequent
collective motion. In non-central Au+Au collisions, ini-
tially, the reaction zone is asymmetric (almond shaped).
The pressure gradient is large in one direction and small
in the other. The asymmetric pressure gradients gener-
ates the elliptic flow. Naturally, in a central collision,
reaction zone is symmetric and elliptic flow vanishes.
Observed elliptic flow then give the strongest indication
that in non-central Au+Au collisions, a collective QCD
matter is produced. Whether the formed matter can be
identified as the much sought after Quark-Gluon Plasma
(QGP) as predicted in Lattice QCD simulations [5] is
presently debatable.
Ideal hydrodynamics has been partly successful in ex-
plaining the observed elliptic flow, quantitatively [6]. El-
liptic flow of identified particles, up to pT ∼1.5 GeV are
well reproduced in ideal hydrodynamics. Ideal hydrody-
namics also explains the transverse momentum spectra
of identified particles (up to pT ∼ 1.5 GeV). Success of
ideal hydrodynamics in explaining bulk of the data [6],
together with the string theory motivated lower limit of
∗E-mail:akc@veccal.ernet.in
shear viscosity η/s ≥ 1/4π [7, 8] has led to a paradigm
that in Au+Au collisions, a nearly perfect fluid is created.
However, the paradigm of ”perfect fluid” produced in
Au+Au collisions at RHIC need to be clarified. As in-
dicated above, the ideal hydrodynamics is only partially
successful and in a limited pT range (pT ≤1.5 GeV) [9].
The transverse momentum spectra of identified particles
also starts to deviate form ideal fluid dynamics prediction
beyond pT ≈ 1.5 GeV. Experimentally determined HBT
radii are not reproduced in the ideal fluid dynamic mod-
els, the famous ”HBT puzzle” [10]. It also do not repro-
duce the experimental trend that elliptic flow saturates
at large transverse momentum. These shortcomings of
ideal fluid dynamics indicate greater importance of dis-
sipative effects in the pT ranges greater than 1.5 GeV
or in more peripheral collisions. Indeed, ideal fluid is a
concept, never realized in nature. As suggested in string
theory motivated models [7, 8], QGP viscosity could be
small, η/s ≥ 1/4π, nevertheless it is non-zero. It is im-
portant to study the effect of viscosity, even if small, on
space-time evolution of QGP fluid and quantify its effect.
This requires a numerical implementation of relativistic
dissipative fluid dynamics. Furthermore, if QGP fluid
is formed in heavy ion collisions, it has to be charac-
terized by measuring its transport coefficients, e.g. heat
conductivity, bulk and shear viscosity. Theoretically, it
is possible to obtain those transport coefficients in a ki-
netic theory model. However, in the present status of
theory, the goal can not be achieved immediately, even
more so for a strongly interacting QGP (sQGP). Alter-
natively, one can use the experimental data to obtain
a ”phenomenological” limit of transport coefficients of
sQGP. It will also require a numerical implementation of
relativistic dissipative fluid dynamics. There is another
incentive to study dissipative hydrodynamics. Ideal hy-
drodynamics depends on the assumption of local equilib-
rium. Before local equilibrium is attained, the system has
to pass through a non-equilibrium stage, where (if non-
equilibrium effects are small) dissipative hydrodynamics
may be applicable. Indeed, we can explore early times of
fluid evolution better in a dissipative hydrodynamics.
http://arxiv.org/abs/0704.0134v2
mailto:akc@veccal.ernet.in
Theory of dissipative relativistic fluid has been formu-
lated quite early. The original dissipative relativistic fluid
equations were given by Eckart [11] and Landau and Lif-
shitz [12]. They are called 1st order theories. Formally,
relativistic dissipative hydrodynamics are obtained from
an expansion of entropy 4-current, in terms of dissipative
fluxes. In 1st order theories, entropy 4-current contains
terms linear in dissipative quantities. 1st order theory
of dissipative hydrodynamics suffer from the problem of
causality violation. Signal can travel faster than light.
Causality violation is unwarranted in any theory, even
more in a relativistic theory. The problem of causal-
ity violation is removed in the Israel-Stewart’s 2nd order
theory of dissipative fluid [13]. In 2nd order theory, ex-
pansion of entropy 4-current contains terms 2nd order in
dissipative fluxes. However, these leads to complications
that dissipative fluxes are no longer function of the state
variables only. They become dynamic. The space of ther-
modynamic variables has to be extended to include the
dissipative fluxes (e.g. heat conductivity, bulk and shear
viscosity).
Even though 2nd order theory was formulated some
30 years back, significant progress towards its numer-
ical implementation has only been made very recently
[14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. At the Cyclotron
Centre, Kolkata, we have developed a code ”AZHYDRO-
KOLKATA” to simulate the hydrodynamic evolution of
dissipative QGP fluid. Presently only dissipative ef-
fect included is the shear viscosity. Some results of
AZHYDRO-KOLKATA, for first order dissipative hydro-
dynamics have been published earlier [19, 20, 21]. In the
present paper, for the first time, we will present some
results for 2nd order dissipative hydrodynamics in 2+1
dimensions. In the present paper, we will consider effect
of dissipation in the QGP phase only. Effect of phase
transition will be studied in a later publication.
The paper is organized as follows: In section II we
briefly review relativistic dissipative fluid dynamics. In
section III we derive the relevant equations in 2+1 di-
mension (assuming boost-invariance). Required inputs
e.g. the equation of state, viscosity coefficient and initial
conditions are discussed in section IV. Simulation re-
sults from AZHYDRO-KOLKATA are shown in section
V. In section VII we compare the transverse momentum
spectra and elliptic flow of quarks in ideal and viscous
dynamics. The concluding section IX summarizes our
results.
II. DISSIPATIVE FLUID DYNAMICS
In this section, I briefly discuss the phenomenological
theory of dissipative hydrodynamics. More detailed ex-
position can be found in [13].
A simple fluid, in an arbitrary state, is fully speci-
fied by primary variables: particle current (Nµ), energy-
momentum tensor (T µν) and entropy current (Sµ) and
a number of additional (unknown) variables. Primary
variables satisfies the conservation laws;
µ = 0, (2.1)
µν = 0, (2.2)
and the 2nd law of thermodynamics,
µ ≥ 0. (2.3)
In relativistic fluid dynamics, one defines a time-like
hydrodynamic 4-velocity, uµ (normalized as u2 = 1).
One also define a projector, ∆µν = gµν − uµuν , orthog-
onal to the 4-velocity (∆µνuν = 0). In equilibrium, an
unique 4-velocity (uµ) exists such that the particle den-
sity (n), energy density (ε) and the entropy density (s)
can be obtained from,
Nµeq = nuµ (2.4)
T µνeq = εu
µuν − p∆µν (2.5)
Sµeq = suµ (2.6)
An equilibrium state is assumed to be fully specified
by 5-parameters, (n, ε, uµ) or equivalently by the thermal
potential, α = µ/T (µ being the chemical potential) and
inverse 4-temperature, βµ = uµ/T . Given a equation of
state, s = s(ε, n), pressure p can be obtained from the
generalized thermodynamic relation,
Sµeq = pβ
µ − αNµeq + βλT λµeq (2.7)
Using the Gibbs-Duhem relation, d(pβµ) = Nµeqdα −
T λµeq dβλ, following relations can be established on the
equilibrium hyper-surface Σeq(α, β
dSµeq = −αdNµeq + βλdT λµeq (2.8)
In a non-equilibrium system, no 4-velocity can be
found such that Eqs.2.4,2.5,2.6 remain valid. Tensor de-
composition leads to additional terms,
Nµ = Nµeq + δN
µ = nuµ + V µ (2.9)
T µν = T µνeq + δT
= [εuµuν − p∆µν ] + Π∆µν + πµν
+(Wµuν +W νuµ) (2.10)
Sµ = Sµeq + δS
µ = suµ + Φµ (2.11)
The new terms describe a net flow of charge V µ =
∆µνNν , heat flow, W
µ = (ε + p)/nV µ + qµ (where qµ
is the heat flow vector), and entropy flow Φµ. Π =
µν − p is the bulk viscous pressure and πµν =
(∆µσ∆ντ +∆νσ∆µτ − 1
∆µν∆στ ]Tστ is the shear stress
tensor. Hydrodynamic 4-velocity can be chosen to elimi-
nate either V µ (the Eckart frame, uµ is parallel to particle
flow) or the heat flow qµ (the Landau frame, uµ is par-
allel to energy flow). In relativistic heavy ion collisions,
central rapidity region is nearly baryon free and Lan-
dau’s frame is more appropriate than the Eckart’s frame.
Dissipative flows are transverse to uµ and additionally,
shear stress tensor is traceless. Thus a non-equilibrium
state require 1+3+5=9 additional quantities, the dissi-
pative flows Π, qµ (or V µ) and πµν . In kinetic theory,
Nµ and T µν are the 1st and 2nd moment of the distri-
bution function. Unless the function is known a-priori,
two moments do not furnish enough information to enu-
merate the microscopic states required to determine Sµ,
and in an arbitrary non-equilibrium state, no relation
exists between, Nν , T µν and Sµ. Only in a state, close
to a equilibrium one, such a relation can be established.
Assuming that the equilibrium relation Eq.2.8 remains
valid in a ”near equilibrium state” also, the entropy cur-
rent can be generalized as,
Sµ = Sµeq + dS
µ = pβµ − αNµ + βλT λµ +Qµ (2.12)
where Qµ is an undetermined quantity in 2nd order in
deviations, δNµ = Nµ − Nµeq and δT µν = T µν − T µνeq .
Detail form of Qµ is constrained by the 2nd law ∂µS
0. With the help of conservation laws and Gibbs-Duhem
relation, entropy production rate can be written as,
µ = −δNµ∂µα+ δT µν∂µβν + ∂µQµ (2.13)
Choice of Qµ leads to 1st order or 2nd order theories
of dissipative hydrodynamics. In 1st order theories the
simplest choice is made, Qµ = 0, entropy current con-
tains terms up to 1st order in deviations, δNµ and δT µν.
Entropy production rate can be written as,
µ = ΠX − qµXµ + πµνXµν (2.14)
where, X = −∇.u; Xµ = ∇
− uν∂νuµ and Xµν =
∇<µuν>.
The 2nd law, ∂µS
µ ≥ 0 can be satisfied by postulat-
ing a linear relation between the dissipative flows and
thermodynamic forces,
Π = −ζθ, (2.15)
qµ = −λ nT
∇µ(µ/T ), (2.16)
πµν = 2η∇<µuν> (2.17)
where ζ, λ and η are the positive transport coefficients,
bulk viscosity, heat conductivity and shear viscosity re-
spectively.
In 1st order theories, causality is violated. If, in a given
fluid cell, at a certain time, thermodynamic forces vanish,
corresponding dissipative fluxes also vanish instantly. Vi-
olation of causality is unwanted in any theory, even more
so in relativistic theory. Causality violation of dissipative
hydrodynamics is corrected in 2nd order theories [13]. In
2nd order theories, entropy current contain terms up to
2nd order in the deviations, Qµ 6= 0. The most general
Qµ containing terms up to 2nd order in deviations can
be written as,
Qµ = −(β0Π2−β1qνqν+β2πνλπνλ)
(2.18)
As before, one can cast the entropy production rate
(T∂µS
µ) in the form of Eq.2.14. Neglecting the terms
involving dissipative flows with gradients of equilib-
rium thermodynamic quantities (both are assumed to be
small) and demanding that a linear relation exists be-
tween the dissipative flows and thermodynamic forces,
following relaxation equations for the dissipative flows
can be obtained,
Π = −ζ(θ + β0DΠ) (2.19)
qµ = −λ
) − β1Dqµ
(2.20)
πµν = 2η
∇<µuν> − β2Dπµν
, (2.21)
where D = uµ∂µis the convective time derivative. Unlike
in the 1st order theories, in 2nd order theories, dynamical
equations control the dissipative flows. Even if thermo-
dynamic forces vanish, dissipative flows do not vanish
instantly.
Before we proceed further, it may be mentioned that
the parameters, α and βλ are not connected to the actual
state (Nµ, T µν). The pressure p in Eq.2.12 is also not
the ”actual” thermodynamics pressure, i.e. not the work
done in an isentropic expansion. Chemical potential α
and 4-inverse temperature βλ has meaning only for the
equilibrium state. Their meaning need not be extended
to non-equilibrium states also. However, it is possible to
fit a fictitious ”local equilibrium” state, point by point,
such that pressure p in Eq.2.12 can be identified with
the thermodynamic pressure, at least up to 1st order.
The conditions of fit fixes the underlying non-equilibrium
phase-space distribution.
III. (2+1)-DIMENSIONAL VISCOUS
HYDRODYNAMICS WITH LONGITUDINAL
BOOST INVARIANCE
Complete dissipative hydrodynamics is a numerically
challenging problem. It requires simultaneous solution of
14 partial differential equations (5 conservation equations
and 9 relaxation equations for dissipative flows). We re-
duce the problem to solution of 6 partial differential equa-
tions (3 conservation equations and 3 relaxation equa-
tions). In the following, we will study boost-invariant
evolution of baryon free QGP fluid, including the dissi-
pative effect due to shear viscosity only. Shear viscosity
is the most important dissipative effect. For example, in
a baryon free QGP, heat conduction is zero and we can
disregard Eq.2.20. Bulk viscosity is also zero for the QGP
fluid (point particles) and Eq.2.19 can also be neglected.
Shear pressure tensor has 5 independent components but
the assumption of boost invariance reduces the number
of independent components to three. For a baryon free
fluid, we can also disregard the conservation equation
Eq.2.1. With the assumption of boost-invariance, energy-
momentum conservation equation ∂µT
µη = 0 become re-
dundant and only three energy-momentum conservation
equations are required to be solved.
-10 -5 0 5 10
X (fm)
-10 -5 0 5 10
_______ η/s=0.08
_______ η/s=0.135
_______ η/s=0
-10 -5 0 5 10
energy density in x-y plane, τ =2.6 fm.
FIG. 1: (color online). Constant energy density contours in
x-y plane at τ=2.6 fm. The black lines are for ideal fluid
(η/s=0). The red and blue lines are for viscous fluid with
ADS/CFT and perturbative estimate of viscosity, η/s=0.08
and 0.135.
Heavy ion collisions are best described in (τ, x, y, η)
coordinates, where τ =
t2 − z2 is the longitudinal
proper time and η = 1
ln t+z
is the space-time rapid-
ity. r⊥ = (x, y) are the usual cartisan coordinate in
the plane, transverse to the beam direction. Relevant
equations concerning this coordinate transformations are
given in the appendix A.
Explicit equations for energy-momentum conservation
in (τ ,x,y,η) coordinates are given in the appendix B. We
note that unlike in ideal fluid, in viscous fluid dynam-
ics, conservation equations (see Eqs.B1-B3) contain addi-
tional pressure gradients due to shear viscosity. Both T τx
and T τy components of energy-momentum tensor now
evolve under additional pressure gradients. The right-
most term of Eq.B3 also indicate that in viscous dynam-
ics, longitudinal pressure is effectively reduced (note that
the πηη component is negative). Since pressure can not
be negative, shear viscosity is limited by the condition,
p+ τ2πηη ≥ 0.
As evident from the Eqs.B1-B3, in boost-invariant dis-
sipative hydrodynamics, with shear viscosity taken into
account, fluid evolution depends only on seven compo-
nents of the shear stress tensors. They are πττ , πτx,
πτy, πxx, πyy, πxy and πηη. However, all the seven com-
ponents are not independent. Tracelessness, transver-
sality to uµ and the assumption of boost-invariance re-
duces the independent components to three. Presently,
we choose πxx πyy and πxy as the independent compo-
nents. Relaxation equations for the independent compo-
nents are given in the appendix C (see Eqs.C4-C6). They
are solved simultaneously with the three energy momen-
tum conservation equations Eqs.B1-B3, with inputs as
discussed below.
-10 -5 0 5 10
X (fm)
-10 -5 0 5 10
_______ η/s=0.08
_______ η/s=0.135
_______ η/s=0
-10 -5 0 5 10
energy density in x-y plane, τ =8.6 fm
FIG. 2: (color online). same as Fig.1 but at time τ=8.6 fm.
IV. EQUATION OF STATE, VISCOSITY
COEFFICIENT AND INITIAL CONDITIONS
A. Equation of state
One of the most important inputs of a hydrodynamic
model is the equation of state. Through this input, the
macroscopic hydrodynamic models make contact with
the microscopic world. In the present demonstrative cal-
culation we will show results for the QGP phase only.
In the QGP phase, we use the simple equation of state,
p = 1
ε, with energy density given as,
gqgpT
4 (4.1)
where gqgp = ggluon+
gquark is the degeneracy factor for
QGP. ggluon = 2(helicity) × 8(color) is the degeneracy
factor for gluons and gquark = 2(spin)× 3(color)× 2(q+
q̄) ×Nf is the degeneracy factor for Nf flavored quarks.
For Nf ≈ 2.5, the degeneracy factor is gqgp = 42.25
B. Shear viscosity coefficient
Shear viscosity coefficient (η) of QGP or sQGP is quite
uncertain. In a strongly coupled QGP, shear viscosity
can not be computed. Recently, using the ADS/CFT
correspondence [7, 8] shear viscosity of a strongly coupled
gauze theory, N=4 SUSY YM, has been evaluated, η =
N2c T
3 and the entropy is given by s = π
N2c T
3. Thus
in the strongly coupled field theory,
ADS/CFT
≈ 0.08, (4.2)
Shear viscosity is quite uncertain in perturbative QCD
also. At high temperature, shear viscosity, in leading log,
can be written as [24, 25],
η = κ
g4 ln g−1
, (4.3)
where g is the strong coupling constant. The leading
log shear viscosity coefficient κ depend on the number of
fermion flavors (Nf ). For example, for two flavored QGP,
κ = 86.47 and κ = 106.7 for a three flavored QGP. With
entropy density of QGP, s = π
gqgpT
3. For two flavored
QGP and αs ≈0.5, the ratio of viscosity over the entropy,
in the perturbative regime is estimated as,
≈ 0.135, (4.4)
For lower αs , perturbative estimation of η/s could be
even higher.
Shear viscosity can also be expressed in terms of sound
attenuation length, Γs, defined as,
(4.5)
Γs is equivalent to mean free path and for a valid hy-
drodynamic description Γs/τ << 1, i.e. mean free path
is much less than the system size. Initial conditions of
the fluid must be chosen carefully such that the validity
condition Γs/τ << 1 remains valid initially as well as at
later time also. In the present work, we have treated vis-
cosity as a parameter. To explore the effect of viscosity,
we have used both the ADS/CFT estimate η/s=0.08 and
perturbative estimate η/s=0.135. We have also run the
code with a higher value of viscosity η/s=0.2.
C. Initial conditions
Solution of Eqs.B1-B3 require initial conditions, the
initial time τi, the transverse distribution of energy den-
sity ε(x, y) and the velocities vx(x, y) and vy(x, y). Fol-
lowing [6], initial transverse energy density is parame-
terized geometrically. At an impact parameter ~b, trans-
verse distribution of wounded nucleons NWN (x, y,~b) and
-10 -5 0 5 10
X axis Title
-0.78
-0.56
-0.33
-0.11
-10 -5 0 5 10
contour plot of Vx at τ=8.6 fm
X axis Title
-10 -5 0 5 10
X axis Title
FIG. 3: (color online). contours of constant vx in x-y plane
at τ=8.6 fm. The black lines are for ideal fluid (η/s=0). The
red and blue lines are for viscous fluid with ADS/CFT and
perturbative estimate of viscosity, η/s=0.08 and 0.135.
-10 -5 0 5 10
X axis Title
-0.78
-0.56
-0.33
-0.11
-10 -5 0 5 10
contour plot of Vy at τ=8.6 fm
X axis Title
-10 -5 0 5 10
X axis Title
FIG. 4: (color online). contours of constant vy in x-y plane
at τ=8.6 fm. The black lines are for ideal fluid (η/s=0). The
red and blue lines are for viscous fluid with ADS/CFT and
perturbative estimate of viscosity, η/s=0.08 and 0.135.
of binary NN collisions NBC(x, y,~b) to are calculated in
a Glauber model. A collision at impact parameter ~b is
assumed to contain 25% hard scattering (proportional
to number of binary collisions) and 75% soft scattering
(proportional to number of wounded nucleons). Trans-
verse energy density profile at impact parameter~b is then
obtained as,
ε(x, y,~b) = ε0(0.75×NWN (x, y,~b) + 0.25×NBC(x, y,~b))
(4.6)
with central energy density ε0=30GeV/fm
−3. The equi-
libration time is chosen as τi=0.6 fm [6]. The initial ve-
locities vx and vy are assumed to be zero initially.
______ η/s=0.08
______ η/s=0
______ η/s=0
Temperature in τ-x plane. y=0.
X (fm)0 2 4 6 8 10
FIG. 5: (color online). Constant temperature contour in x −
τ plane, for fixed y=0. The black, red and blue lines are
for ideal, viscous fluid with η/s=0.08 and viscous fluid with
η/s=0.135.
In dissipative hydrodynamics, one requires initial con-
ditions for the viscous pressures also. Due to longitudinal
boost invariance of the problem, we assume that viscous
pressures have attained their boost-invariant values at
the time of equilibration. Boost invariant values of the
three independent shear stress-tensors can be easily ob-
tained from Eqs.C4-C6, σxx = σyy = θ = 1
and σxy = 0
( at the initial time τi, u
µ = (1, 0, 0, 0), Duµ = 0. The
initial distribution of shear pressure tensors are then ob-
tained as,
πxx(x, y,~b) = 2ησxx = 2η/τi (4.7)
πyy(x, y,~b) = 2ησyy = 2η/τi (4.8)
πxy(x, y,~b) = 2ησxy = 0 (4.9)
As stated earlier, the viscous coefficient η is obtained
using the relation, η/s = const, const=0.08, 0.135 or 0.2.
For these values of shear viscosity, the validity condition
Γs/τ << 1 is satisfied initially. The validity condition is
better satisfied at later time.
V. RESULTS
VI. STABILITY OF NUMERICAL SOLUTIONS
A. Evolution of the viscous QGP fluid
The energy-momentum conservation equations B1-B3,
and the relaxation equations C2-C4 are solved simultane-
ously using the code, AZHYDRO-KOLKATA, developed
at the Cyclotron Centre, Kolkata. As mentioned earlier,
we have solved the equations in the QGP phase only and
did not consider any phase transition. In the following
we will show the results for central Au+Au collisions (im-
pact parameter b = 0 fm). To understand the effect of
shear viscosity, with the same initial conditions, we have
solved the energy-momentum conservation equations for
ideal fluid and viscous fluid. As mentioned earlier, we
have considered two values of viscosity, the ADS/CFT
motivated value η/s=0.08 and the perturbative estimate,
η/s=0.135.
0 2 4 6 8 10
X (fm)
_____ η/s=0.08
_____ η/s=0.135
_____ η/s=0
Temperature   in τ-x plane. y=5 fm.
FIG. 6: (color online). same as fig.5 but at y=5.0
In Fig.1, we have shown the contours of constant en-
ergy density in x-y plane, after an evolution of 2.6 fm.
The black lines are for ideal fluid evolution. The red and
blue lines are for viscous fluid with ADS/CFT (η/s=0.08)
and perturbative (η/s=0.135) estimate of viscosity. Con-
stant energy density contours, as depicted in Fig.1, indi-
cate that with viscosity fluid cools slowly. Cooling gets
slower as viscosity increases. Thus at any point in the x-y
plane, energy density of viscous fluid is higher than that
of an ideal fluid. At later time also, compared to an ideal
fluid, viscous fluid evolve slowly. In Fig.2, contours of
constant energy density at time τ=8.6 fm is shown. Here
also we find than at any point energy density of viscous
fluid is higher than its ideal counter part. The result is
in accordance with our expectation. For dissipative fluid,
equation of motion can be written as,
Dε = −(ε+ p)∇µuµ + πµν∇<µuν> (6.1)
Due to viscosity, evolution of energy density (or tem-
perature) is slowed down.
In Fig.3 and 4, we have shown the contour plot of the
fluid velocity, vx and vy, after evolution of 8.6 fm. As
before the black lines are for the ideal fluid evolution.
The red and blue lines are for viscous fluid with η/s=0.08
and 0.135 respectively. Fluid velocities in viscous and
ideal fluid differ very little. Even at late time, as shown
in Fig.3 and 3, we find that for η/s=0.08-0.135, x and y
component of the fluid velocity show marginal difference.
However, there is an indication that in a viscous fluid,
velocity grow faster than in ideal fluid. But as mentioned
earlier, the difference is marginal.
η/s=0.08
η/s=0.135
τ (fm)
0 2 4 6 8 10
FIG. 7: (color online). In the upper panel, temporal evolu-
tion of the shear pressure tensor πxx at the fluid cell x=y=0
is shown. In the lower panel, evolution of πxy at the fluid
cell x=y=5 fm is shown. The black and red lines are for
ADS/CFT motivated viscosity η/s=0.08 and perturbative es-
timate η/s=0.135 respectively.
As seen in Fig.1-2, in viscous dynamics, QGP fluid
evolves slowly. Thus life-time of the QGP phase is en-
hanced in viscous dynamics. To obtain an idea about the
enhanced life-time, in Fig.5, we have shown the constant
temperature contours in τ − x plane , at a fixed value
of y=0 fm. As seen in Fig.5, temperature evolves slowly
in a viscous fluid and life-time of the QGP phase is ex-
tended. For small viscosity η/s=0.08-0.135, the increase
is not large. At the center of the fluid, for η/s=0.135,
QGP life-time is increased approximately by 5% only. It
is even less for the ADS/CFT estimate of viscosity. How-
ever, enhancement of QGP life-time depends on the fluid
cell position. It could be more. In Fig.6, constant tem-
perature contours at y=5 fm is shown. For η/s=0.135, at
x=0,y=5 fm, the QGP life-time is enhanced by ∼ 10%.
We conclude that in a viscous dynamics, with moder-
ate viscosity η/s=0.08-0.135, QGP life-time could be en-
hanced by 5-10%. Enhanced lifetime of QGP in a viscous
fluid can have significant effect on observables produced
early in the collisions e.g. direct photon production or in
J/ψ suppression.
1.01.21.4
0 2 4 6 8 10
0.030
0.180.21
0 2 4 6 8 10
(d) π
 at τ=2.6 fm(c) π
 at τ=0.6 fm
(b) π
 at τ=2.6 fm(a) π
 at τ=0.6 fm
1.01.21.4
0 2 4 6 8 10
0.030
0.150.18
0.210.24
0 2 4 6 8 10
X (fm)
FIG. 8: (color online). In panel (a) and (b), contours of
constant pressure tensor πxx at initial time τi=0.6 fm and at
time τ=2.6 fm is shown. In panel (c) and (d) same results for
shear pressure tensor πyy is shown.
B. Evolution of shear pressure tensors
We have assumed that initially the shear pressure ten-
sors πxx, πyy and πxy attained their longitudinal boost-
invariant values. As the fluid evolve, pressure tensors also
evolve. Here we investigate the evolution of shear pres-
sure tensors with time. In the top panel of Fig.7 evolu-
tion of shear pressure tensor πxx at the fluid cell position
x=y=0 is shown. The black line is for the ADS/CFT
motivated viscosity, η/s=0.08 and the red line is for the
perturbative estimate of viscosity η/s=0.135. Just after
the start of the evolution the shear pressure tensor πxx
increases, but for a short duration and then steadily de-
creases with time. By 4 fm of evolution, πxx at the center
of the fluid reduces to negligibly small values. Identical
behavior is seen for the shear pressure tensor πyy. In
the bottom panel of Fig.7 we have shown the evolution
of the third independent shear pressure tensor πxy. Ini-
tially πxy is zero. As the fluid evolve, it grow in the
negative direction. We find that at the centre of the fluid
(x=y=0), it never grows. In Fig.7, temporal evolution
of πxy at the fluid cell position x = y = 5fm is shown.
From the initial zero value, πxy rapidly increases in the
negative direction. It reaches its maximum around τ ≈1
fm and then decreases again. We also note that πxy never
grows to large values. Compared to πxx or πyy stress ten-
sor πxy is negligible. The results indicate that in a QGP
fluid, viscous effect persist for a short duration (3-4 fm)
only. At late time the fluid evolve essentially as an ideal
fluid. The result is understandable. Shear viscosity de-
pend strongly on temperature (η ∝ T 3). As the fluid
cools, effect of viscosity decreases rapidly.
τ (fm)
0 2 4 6 8 10
η/s=0.08
η/s=0.135
FIG. 9: Evolution of average entropy with time, for two val-
ues of viscosity, the ADS/CFT motivated viscosity η/s=0.08
and perturbative estimate η/s=0.135 are shown.
To show the spatial distribution of the stress tensors,
in Fig.8, πxx and πyy at initial time τi=0.6 fm and after
an evolution of τ = 2.6 fm are shown. As shown earlier,
πxx and also πyy rapidly decreases with time. By 2 fm
of evolution they are reduced by approximately by a fac-
tor 6. It is also interesting to note that the initial x-y
symmetric distribution of πxx and πyy quickly evolves to
asymmetric distribution. With time πxx evolves faster in
the x-direction than in y-direction. Similarly, πyy evolve
faster in the y-direction than in the x-direction. For cen-
tral collisions the asymmetric evolution of πxx and πyy
counter balance each other. As shown in Fig.1 and 2,
the contour plots of energy density do not show any in-
dication of asymmetry even at late time. However, the
asymmetric pressure tensors can have important effects
on elliptic flow of observables produced early in the col-
lisions, say in elliptic flow of direct photons.
C. Entropy generation
In a viscous fluid dynamics, entropy is generated. We
can easily calculate the entropy generated during the evo-
lution,
πµνπµν
[(πττ )2 + (πxx)2 + (πyy)2 + (τ2πηη)2
−2(πτx)2 − 2(πτy)2 + 2(πxy)2] (6.2)
Evolution of spatially averaged entropy is shown in
Fig.9, for the two values of viscosity coefficients η/s=0.08
and 0.135. As expected, entropy generation is more if vis-
cosity is more. For both the values of viscosity, we find
that entropy generation saturates after ≈ 3 fm of evolu-
tion. It is expected also. As shown previously, viscous
fluxes reduces to very small values after τ=3 fm. Natu-
rally, entropy generation is negligible thereafter.
0.190.22
 0 2 4 6 8 10
10 Temperature in τ-x plane, y= 5m.
_____ viscous fluid (2nd order)
_____ ideal fluid
_____ viscous fluid (1st order)
X (fm)
FIG. 10: (color online) constant temperature contours in
x − τ plane at y=5 fm. The black lines are for ideal fluid.
The red and blue lines are for viscous fluid in 1st order and
in 2nd order theory respectively. η/s=0.135.
D. 1st order theory vs. 2nd order theory
As mentioned earlier, 1st order theory of dissipative
hydrodynamics is acausal, signal can travel faster than
light. This is corrected in 2nd order theory, but we
have to pay the price, relaxation equations for dissipa-
tive fluxes are required to be solved. It is interesting to
τ (fm)
0 2 4 6 8 10 12
1st order theory
2nd order theory
η/s=0.135
FIG. 11: Evolution of average entropy production in a 1st
order (solid line) and 2nd order (dashed line) theory. 2nd
order theory generate more entropy.
compare the difference we can expect in a first order the-
ory and in a 2nd order theory of dissipation. In Fig. 10,
we have shown the contours of constant temperature in
x−τ , for a fixed y = 5fm. The black lines are for an ideal
fluid. The red lines are for a viscous fluid treated in the
1st order theory. The blue lines are for viscous fluid in
2nd order theory. In 2nd order theory fluid evolve more
slowly than in a first order theory. Entropy generation is
also more in a 2nd order theory. In Fig.11, average en-
tropy evolution with proper time is shown, both for the
1st order theory (the solid line) and the 2nd order theory.
In 2nd order theory, approximately 80% more entropy is
generated.
VII. TRANSVERSE MOMENTUM AND
ELLIPTIC FLOW OF QUARKS
Presently we can not compare predictions from viscous
hydrodynamics with experimental data. Hadrons are not
included in the model. The initial QGP fluid evolve and
cools but remain in the QGP phase, it did not undergo
a phase transition to hadronic gas. However, from the
momentum distribution of quarks we can get some idea
about the viscous effect on particle production. Viscos-
ity generates entropy, which will be reflected in enhanced
multiplicity. We use the standard Cooper-Frey prescrip-
tion to obtained the transverse momentum distribution
of quarks. In Cooper-Frey prescription, particle distribu-
tion is obtained by convoluting the one body distribution
function over the freeze-out surface,
µf(x, p) (7.1)
where dΣµ is the freeze-out hyper-surface and f(x, p) is
the one-body distribution function. Now in a viscous dy-
namics, the fluid is not in equilibrium and f(x, p) can not
be approximated by the equilibrium distribution func-
tion,
f (0)(x, p) =
exp[β(uµpµ − µ)] ± 1
, (7.2)
with inverse temperature β = 1/T and chemical poten-
tial µ. In a highly non-equilibrium system, distribution
function f(x, p) is unknown. If the system is slightly
off-equilibrium, then it is possible to calculate correction
to equilibrium distribution function due to (small) non-
equilibrium effects. Slightly off-equilibrium distribution
function can be approximated as,
f(x, p) = f (0)(x, p)[1 + φ(x, p)], (7.3)
φ(x, p) is the deviation from equilibrium distribution
function f (0). With shear viscosity as the only dissi-
pative forces, φ(x, p) can be locally approximated by a
quadratic function of 4-momentum,
φ(x, p) = εµνp
µpν . (7.4)
Without any loss of generality εµν can be written as
εµν =
2(ε+ p)T 2
πµν , (7.5)
completely specifying the non-equilibrium distribution
function. As expected, correction factor increases with
increasing viscosity. We also note that non-equilibrium
correction is more on large momentum particles. The
effect of viscosity is more on large momentum parti-
cles. The correction factor reduces if freeze-out occurs
at higher temperature.
With the corrected distribution function, we can cal-
culate the quark momentum spectra at freeze-out sur-
face Σµ. In appendix D, relevant equations are given.
The quark momentum distribution has two parts, (i)
dyd2pT
, obtained by convoluting the equilibrium distribu-
tion function over the freeze-out surface and (ii) dN
dyd2pT
obtained by convoluting the correction to the equilibrium
distribution function over the freeze-out surface. Since
the correction factor is obtained under the assumption
that non-equilibrium effects are small, φ(x, p) << 1, it
necessarily imply that, dN
dyd2pT
<< dN
dyd2pT
. The ratio,
dNneq
dNneq
dyd2pT
dyd2pT
, (7.6)
could at best be unity or less. If the ratio exceeds
unity, it will imply that non-equilibrium effects are large
 (GeV)
0 1 2 3 4 5
η/s=0.2 
η/s=0.135 
η/s=0.08 
Au+Au@b=6.5 fm
τi=0.6fm,Sini=110fm-3,TF=160MeV
FIG. 12: Ratio of quark spectra with non-equilibrium distri-
bution function to that with equilibrium distribution function.
and the distribution function f(x, p) can not be approx-
imated as in Eq.7.3. Using AZHYDRO-KOLKATA, we
have simulated a b=6.5 fm Au+Au collision. dN
dyd2pT
dNneq
dyd2pT
at freeze-out temperature TF =160 MeV are cal-
culated. The ratio dN
for η/s=0.08,0.135 and 0.2, are
shown in Fig.12. With ADS/CFT estimate of viscosity,
η/s=0.08, non-equilibrium correction to particle produc-
tion become comparable to equilibrium contribution only
beyond pT =5 GeV. However, with perturbative estimate,
η/s=0.135, non-equilibrium correction become compara-
ble to or exceeds the equilibrium contribution at pT ∼
4 GeV. pT range is further reduced for higher viscosity
η/s=0.2. Thus with perturbative estimate of viscosity
(η/s = 0.135 − 0.2), hydrodynamic description remain
valid upto transverse momentum pT ∼ 3.5-4 GeV.
In Fig.13, we have compared the transverse momen-
tum spectra of quarks in ideal hydrodynamics with that
in a viscous dynamics. In Fig.13, the dotted line is the
spectra obtained in ideal dynamics (η/s = 0). The pT
spectra in viscous dynamics are shown by black lines.
We have shown the spectra for three values of viscos-
ity η/s=0.08,0.135 and 0.2. Compared to ideal dynam-
ics, quarks yield in viscous dynamics increases. The in-
crease is more at large pT . For low values of viscosity
the increase is modest, a factor of 2 at pT =3 GeV. But
yield increase by a factor or 4(10) if viscosity increases
to η/s=0.135 (2). Please note that even though we have
shown pT spectra upto 5 GeV, for η/s=0.2 and 0.135,
hydrodynamic description fails beyond pT ∼ 3.5 and 4
We have also studied the effect of viscosity on quark
elliptic flow. Effect of viscosity is very prominent on
elliptic flow. In Fig.14, pT dependence of elliptic flow
of quarks, in a b=6.5 fm collision is shown. The black
line is v2 in ideal dynamics. In ideal dynamics, ellip-
pT (GeV)
0 1 2 3 4 5
η/s=0.0(ID.Fluid)
η/s=0.2 
η/s=0.135 
η/s=0.08 
Au+Au@b=6.5 fm
τi=0.6fm,Sini=110fm-3,TF=160MeV
FIG. 13: Qurak transverse momentum spectra at freeze-out
temperature of 160 MeV. The dotted line is the quarks spectra
in ideal hydrodynamics. The solid lines (top to bottom) are
in viscous dynamics with η/s=0.08,0.135 and 0.2.
tic flow continually increases with pT . It well known, in
contrast to experiments, where elliptic flow saturates at
large pT , in ideal hydrodynamics, elliptic flow continue
to increase with pT . Indeed, this is a major problem in
ideal hydrodynamics. The renewed the interest in dis-
sipative hydrodynamics is partly due to the inability of
ideal hydrodynamics to predict the trend of elliptic flow
in Au+Au collisions. In Fig.14, the blue lines are v2
in viscous dynamics with η/s=0.08,0.135 and 0.2 respec-
tively. In a viscous dynamics, pT dependence of v2 is
drastically changed. In contrast to ideal dynamics where
v2 continue to increase with pT , in viscous dynamics, v2
continue to increase only upto pT ∼ 1.5− 2GeV . There-
after v2 decreases. For perturbative estimate of viscos-
ity η/s=0.135 and beyond, v2 even become negative at
large pT . Veering about of v2 after pT ∼1.5-2 GeV is
due to viscous effect only or more explicitly due to the
non-equilibrium correction to the equilibrium distribu-
tion function. This is clearly manifested from the red
lines in Fig.14. The red lines are calculated ignoring the
non-equilibrium corrections to the equilibrium distribu-
tion function. If non-equilibrium correction is ignored, in
viscous dynamics also, v2 continue to increase with pT ,
albeit its magnitude is reduced compared to ideal dy-
namics. The result is very important. It imply that the
experimental trend of elliptic flow (saturation at large
pT ) could only be explained if the QGP fluid is viscous.
An ideal QGP, will not be able to explain the saturation
trend of the experimental data.
As stated earlier, non-equilibrium correction to the
equilibrium distribution function depends on the freeze-
out condition. To show the effect of freeze-out condition,
on v2, in Fig.15 we have shown v2 for a values of freeze-
out temperature TF =160,150,140,130 and 120 MeV. As
pT (GeV)
0 1 2 3 4 5
vis.fluid with Fneq
vis.fluid with Feq
id. fluid
Au+Au@b=6.5 fm
τi=0.6fm,Sini=110fm-3,TF=160MeV
FIG. 14: (color online) Elliptic flow as a function of transverse
momentum. The black line is v2 in ideal hydrodynamics. The
blue lines are v2 in viscous dynamics with viscosity to entropy
ratio η/s=0.08,0.135 and 0.2 (top to bottom) respectively,
including the correction to equilibrium distribution function.
The red lines are same as the blues lines but ignoring the
non-equilibrium correction to the distribution function.
freeze-out occur at higher and higher temperature, the
veering of v2 takes place at larger and larger pT and for
TF =120 MeV, the elliptic flow saturates. The result is
understood easily. With decreasing freeze-out tempera-
ture, the fluid evolves for longer time, the shear stress-
tensor’s at the freeze-out surface is reduced and the non-
equilibrium correction, proportional to shear stress ten-
sors, decreases.
VIII. STABILITY OF NUMERICAL
SOLUTIONS IN AZHYDRO-KOLKATA
Before we summarise our results, we would like to
comment on the stability of numerical solutions in
AZHYDRO-KOLKATA. As indicated above, with shear
viscosity as the only dissipative force, boost-invariant
causal hydrodynamics require simultaneous solution of
six partial differential equations. Numerical solution of
six partial differential equations is non-trivial and it is
important to check for the numerical stability of the
solutions. Analytical solutions of viscous hydrodynam-
ics, even in restrictive conditions are not available, and
we can not check the solutions against analytical re-
sults. However, we can check for the stability of the
numerical solutions. The standard procedure of check-
ing the numerical stability is to change the integration
step lengths and look for the difference in the solution.
In Fig.16, for viscosity η/s=0.135, we have shown the
constant temperature contours in x − τ plane at a fixed
pT (GeV)
0 1 2 3 4 5
Au+Au@b=6.5 fm
τi=0.6fm,η/s=.2,Sini=110fm-3
TF=160,150,140,130 and 120 MeV
(bottom to top)
FIG. 15: Dependence of elliptic flow on the freeze-out tem-
perature. The solid lines (from bottom to top) are elliptic
flow (v2) in viscous dynamics with TF =160,150,140,130 and
120 MeV respectively.
0.250.28
0 2 4 6 8 10
____dx=dy=0.1,dτ=0.01
____dx=dy=0.2,dτ=0.02Au+Au@b=0
X (fm)
FIG. 16: (color online) constant temperature contours in
x − τ plane at a fixed y=0 fm. The black lines are obtained
with integration step lengths dx=dy=0.2 fm and dτ=0.02 fm.
The blue lines are obtained with integration step lengths,
dx=dy=0.1 fm and dτ=0.01 fm. Halving the step lengths do
not change the evolution. The numerical solutions are stable.
y=0 fm. The black and blue lines are obtained when
integration step lengths are dx=dy=0.2fm,dτ=0.02 fm,
and dx=dy=0.1 fm,dτ=0.01fm respectively. Evolution
of QGP fluid donot alter by changing the step lengths,
the solutions are stable against mesh size.
IX. SUMMARY AND CONCLUSIONS
In Israel-Stewart’s 2nd order theory of dissipative rel-
ativistic hydrodynamics, we have studied evolution QGP
fluid. In 2nd order theory, in addition to usual thermo-
dynamic quantities e.g. energy density, pressure, hydro-
dynamic velocities, dissipative flows are treated as ex-
tended thermodynamic variables. Relaxation equations
for dissipative flows are solved, simultaneously with the
energy-momentum conservation equations. This greatly
enhances the complexity of the problem. Altogether 14
partial differential equations are required to be solved.
We simplify the problem to solution of six partial differ-
ential equations by considering the evolution of baryon
free QGP fluid with longitudinal boost-invariance. We
also consider dissipation due to shear viscosity only, dis-
regarding the bulk viscosity and the heat conduction
(for a baryon free QGP fluid they do not contribute).
The six partial differential equations are solved using
the code AZHYDRO-KOLKATA, developed at the Cy-
clotron Centre, Kolkata.
To bring out the effect of viscosity, we have considered
the evolution of ideal as well as viscous QGP fluid. Both
ideal and viscous fluid are initialized similarly, at initial
time τi=0.6 fm, the central entropy density is 110 fm
Viscous dynamics require initial conditions for the shear-
stress tensor components. It is assumed that at the equi-
libration time, the shear stress tensors components have
attained their boost-invariant values.
Explicit simulation of ideal and viscous fluids confirms
that energy density of a viscous fluid, evolve slowly than
its ideal counterpart. Thus in a viscous fluid, lifetime of
the QGP phase will be enhanced. Transverse expansion
is also more in viscous dynamics. For a similar freeze-out
condition freeze-out surface is extended in viscous fluid.
As the fluid evolve, shear pressure tensors also evolve.
Explicit simulations indicate that shear pressure tensors
πxx and πyy which are initially non zero, rapidly de-
creases as the fluid evolve. By 3-4 fm of evolution they
reduced to very small values. The other independent
shear tensor πxx is zero initially. At later time it grow
in the negative direction but never grow to large value
and is always order of magnitude smaller than the stress
tensors (πxx and πyy). Spatial distribution of shear pres-
sure tensors πxx and πyy reveal an interesting feature
of viscous dynamics. Initially πxx and πyy have sym-
metric distribution. As the fluid evolve, pressure tensors
quickly become asymmetric, e.g. πxx evolve faster in the
x-direction than in the y-direction, πyy evolve faster in
y direction than in x-direction. However, in a central
collision, we did not see any effect of asymmetry in the
energy density distribution. In a central b=0 collision,
the two opposite asymmetry cancels each other.
We could not study effect of shear viscosity on par-
ticle production. However, we have explored the effect
of viscosity on parton momentum distribution and ellip-
tic flow. We have simulated b=6.5 fm Au+Au collision.
Using the Cooper-Frey prescription, transverse momen-
tum spectra as well as elliptic flow of quarks at freeze-out
temperature of TF =160 MeV are obtained. Viscous dy-
namics flattens the quark yield at large pT . At pT =3
GeV, even a small viscosity, η/s=0.8, increase the yield
by a factor of 2. The increase is even more if viscosity is
large. Viscous effect is most prominent on elliptic flow.
In ideal hydrodynamics, elliptic flow continue to increase
with pT . But in viscous dynamics v2 veer about around
pT =1.5-2 and even become negative at large pT . With
appropriate choice of viscosity, freeze-out condition, el-
liptic flow show saturation. The saturation effect is es-
sentially due to non-equilibrium correction to the equi-
librium distribution function and can not be mimicked
in an ideal hydrodynamics. Only in viscous dynamics,
saturation of elliptic flow can be explained.
APPENDIX A: COORDINATE
TRANSFORMATIONS
Instead of Cartesian coordinates xµ = (t, x, y, z) we
use curvilinear coordinates in longitudinal proper time
and rapidity, x̄m = (τ, x, y, η):
t = τ cosh η; τ =
t2 − z2 (A1)
z = τ sinh η; η =
. (A2)
The differentials
dt = dτ cosh η + dη τ sinh η, (A3)
dz = dτ sinh η + dη τ cosh η, (A4)
and the metric tensor is easily read off from
ds2 = gµνdx
µdxν = dt2 − dx2 − dy2 − dz2
= ḡmndx̄
mdx̄n = dτ2 − dx2 − y2 − τ2dη2,(A5)
namely
ḡmn =
1 0 0 0
0 −1 0 0
0 0 −1 0
0 0 0 −τ2
, ḡmn =
1 0 0 0
0 −1 0 0
0 0 −1 0
0 0 0 −1/τ2
In curvilinear coordinates we must replace the partial
derivatives with respect to xµ by covariant derivatives
(denoted by a semicolon) with respect to x̄m:
T̄ ik;p =
∂T̄ ik
+ ΓipmT̄
mk + T̄ imΓkmp.
The only non-vanishing Christoffel symbols are
Γτηη = τ ; Γ
τη = Γ
ητ = 1/τ. (A7)
The hydrodynamic 4-velocity uµ = γ(1, vx, vy, vz) is
transformed to ūm = γ(1, vx, vy, 0), with γ⊥ =
1−v2r . From here on, we drop the bars over tensor
components in x̄-coordinates for simplicity.
The projector can be easily calculated,
∆µν = gµν − uµuν
1 − γ2⊥ −γ2⊥vx γ2⊥vy 0
−γ2⊥vx −1 − γ2⊥v2x −γ2⊥vxvy 0
−γ2⊥vy −γ2⊥vxvy −1 − γ2⊥v2y 0
0 0 0 1
.(A8)
In (τ, x, y, η) coordinate system, the convective time
derivative can be obtained as,
D = u · ∂ = γ(∂τ + vx∂x + vy∂y). (A9)
For future reference,we also write down the the scalar
expansion rate
θ = ∂·u = ∂τuτ + ∂xux + ∂yuy +
(A10)
APPENDIX B: ENERGY-MOMENTUM
CONSERVATION
With longitudinal boost-invariance the energy-
momentum conservation equations Tmn;n = 0 yield
∂τ T̃
ττ + ∂x(T̃
ττvx) + ∂y(T̃
ττvy) = − (p+ τ2πηη) (B1)
∂τ T̃
τx + ∂x(T̃
τxvx) + ∂y(T̃
τxvy) = −∂x(p̃+ π̃xx − π̃τxvx) − ∂y(π̃xy − π̃τxvy) (B2)
∂τ T̃
τy + ∂x(T̃
τyvx) + ∂y(T̃
τyvy) = −∂x(π̃xy − π̃τyvx) − ∂y(p̃+ π̃yy − π̃τyvy) (B3)
where Ãmn ≡ τAmn, p̃ ≡ τp, and vx ≡ T τx/T ττ , vy ≡
T τy/T ττ .
The components of the energy momentum tensors, in-
cluding the shear pressure tensor are,
T ττ = (ε+ p)γ2⊥ − p+ πττ (B4)
T τx = (ε+ p)γ2⊥vx + π
τx (B5)
T τy = (ε+ p)γ2⊥vy + π
τy (B6)
In causal dissipative hydrodynamics, energy momen-
tum conservation equations are solved simultaneously
with the relaxation equations. Given an equation of
state, if energy density (ε) and fluid velocity (vx and
vy) distributions, at any time τi are known, Eqs.B1,B2
and B3 can be integrated to obtain ε, vx and vy at the
next time step τi+1. While for ideal hydrodynamics, this
procedure works perfectly, viscous hydrodynamics poses
a problem that shear stress-tensor components contains
time derivatives, ∂τγ⊥, ∂τu
x, ∂τu
x etc. Thus at time step
τi one needs the still unknown time derivatives. Numer-
ically, time derivatives at step τi could be obtained if ve-
locities at time step τi and τi+1 are known. One possible
way to circumvent the problem, is to use time derivatives
of the previous step, i.e. use velocities at time step τi−1
and τi to calculate the derivatives at time step τi [18].
The underlying assumption that fluid velocity changes
slowly with time. In 1st order theories, this problem is
circumvented by calculating the time derivatives from the
ideal equation of motion ,
Duµ =
, (B7)
Dε = −(ε+ p)∇µuµ. (B8)
With the help of these two equations all the time
derivatives can be expressed entirely in terms of spa-
tial gradients [15, 26]. 1st order theories are restricted
to contain terms at most linear in dissipative quantities.
Neglect of viscous terms can contribute only in 2nd or-
der corrections, which are neglected in 1st order theories.
While the procedure is not correct in 2nd order theory,
we still use it in the present calculations. The alternative
procedure of using the derivative of earlier time step is
not correct either.
APPENDIX C: RELAXATION EQUATIONS FOR
THE VISCOUS PRESSURE TENSOR
Being symmetric and traceless, the viscous pressure
tensor πµν has 9 independent components. The as-
sumption of boost invariance reduces this number by
3 (∇〈mu η〉 =0, m 6= η). The transversality condition
mn = 0 eliminates another three components ( uη
vanish and thus yield no constraint). Thus, with boost-
invariance the viscous pressure tensor has only three in-
dependent components. As seen in Eqs.B1,B2 and B3 in
a boost-invariant evolution only seven pressure tensors
πττ , πxx, πyy, πηη, πτx, πτy and πxy are of importance.
Only three of these seven are independent. In an ear-
lier publication [17], we have debated about the choice
of the independent components and suggested use of ei-
ther (πττ , πηη, ∆ = πxx − πyy) or (πττ ,πηη,πτx, πτy)
(which will require solution of an additional relaxation
equation) as choice of independent components. How-
ever, while computing we find that the three pressure
tensors πxx and πyy and πxy as independent components
are computationally more convenient. The choice has
the advantage that the dependent shear stress tensors
can be obtained from the 3 independent stress tensors by
multiplying them by fluid velocity, vx and vy (see Eqs.
C7-C10). In any other choice of independent components
(e.g. πττ ,πηη,∆ = πxx − πyy), the evaluation of depen-
dent stress tensors requires division by fluid velocities.
Since initially, fluid velocities are assumed to be zero and
they grow slowly, these choices will involve division by
very small numbers. Unless proper care is not taken, di-
vision by small numbers can lead to unrealistically large
values for the dependent stress tensors and ruin the com-
putation.
The relaxation equations for the independent shear
stress tensors πxx, πyy and πxy, in (τ ,x,y,η) co-ordinate
can be written as,
xx + vx∂xπ
xx + vy∂yπ
xx = − 1
(πxx − 2ησxx) (C1)
yy + vx∂xπ
yy + vy∂yπ
yy = − 1
(πyy − 2ησyy) (C2)
xy + vx∂xπ
xy + vy∂yπ
xy = − 1
(πxy − 2ησxy) (C3)
where τπ is the relaxation time, τπ = 2ηβ2 (see
Eq.2.21). In ultra-relativistic limit, for a Boltzman gas,
β2 can be evaluated, β2 ≈ 34p where p is the pressure
[13]. In the present paper, we use this limit to obtain the
relaxation time τπ.
The viscous pressure tensor relaxes on a time scale τπ
to 2η times the shear tensor σµν = ∇〈µ u ν〉. The xx,
yy and xy components of the shear tensor σµν can be
written as
σxx = −∂xux − uxDux −
∆xxθ (C4)
σyy = −∂yuy − uyDuy −
∆yyθ (C5)
σxy = −1
y − ∂yux − uxDuy − uyDux]
∆xyθ (C6)
The dependent shear stress tensors can easily be ob-
tained from the independent ones as,
πτx = vxπ
xx + vyπ
xy (C7)
πτy = vxπ
xy + vyπ
yy (C8)
πττ = v2xπ
xx + v2yπ
yy + 2vxvyπ
xy (C9)
τ2πηη = −(1 − v2x)πxx − (1 − v2y)πyy
+2vxvyπ
xy (C10)
The expressions for the convective time derivative D
and expansion scalar θ = ∂u̇, in (τ ,x,y,η) are given in
Eqs. A9 and A10.
APPENDIX D: PARTICLE SPECTRA
With the non-equilibrium distribution function thus
specified, it can be used to calculate the particle spectra
from the freeze-out surface. In the standard Cooper-Frye
prescription, particle distribution is obtained as,
dyd2pT
µf(x, p) (D1)
In (τ, x, y, ηs) coordinate, the freeze-out surface is pa-
rameterised as,
Σµ = (τf (x, y) cosh ηs, x, y, τf (x, y) sinh ηs), (D2)
and the normal vector on the hyper surface is,
dΣµ = (cosh ηs,−
,−∂τf
,− sinh ηs)τfdxdydηs
At the fluid position (τ, x, y, ηs) the particle 4-
momenta are parameterised as,
pµ = (mT cosh(ηs − Y ), px, py,mT sinh(ηs − Y )) (D4)
The volume element pµdΣµ become,
pµdΣµ = (mT cosh(η − Y ) − ~pT .~∇T τf )τfdxdydη (D5)
Equilibrium distribution function involve the term
which can be evaluated as,
γ(mT cosh(η − Y ) − ~vT .~pT − µ/γ)
The non-equilibrium distribution function require the
sum pµpνπµν ,
pµpνπ
µν = a1cosh
2(η − Y ) + a2cosh(η − Y ) + a3 (D7)
a1 = m
ττ + τ2πηη) (D8)
a2 = −2mT (pxπτx + pyπτy) (D9)
a3 = p
xx + p2yπ
yy + 2pxpyπ
xy −m2T τ2πηη(D10)
Inserting all the relevant formulas in Eq.D1 and inte-
grating over spatial rapidity one obtains,
dyd2pT
dyd2pT
dNneq
dyd2pT
(D11)
with,
dyd2pT
(2π)3
dxdyτf [mTK1(nβ) − pT ~∇T τfK0(nβ)] (D12)
dNneq
dyd2pT
(2π)3
dxdyτf [mT {
K3(nβ) +
K2(nβ) + (
+ a3)K1(nβ) +
K0(nβ)}
−~pT .~∇T τf{
K2(nβ) + a2K1(nβ) + (
+ a3)K0(nβ)}] (D13)
where K0, K1, K2 and K3 are the modified Bessel func-
tions.
We will also show results for elliptic flow v2. It is de-
fined as,
dyd2pT
cos(2φ)dφ
dyd2pT
(D14)
Expanding to the 1st order, elliptic flow as a function
of transverse momentum can be obtained as,
v2(pT ) = v
2 (pT )
2Nneq
pT dpT dφ
pT dpT dφ
dφcos(2φ) d
2Nneq
pT dpT dφ
pT dpT dφ
(D15)
where v
2 is the elliptic flow calculated with the equi-
librium distribution feq.
[1] BRAHMS Collaboration, I. Arsene et al., Nucl. Phys. A
757, 1 (2005).
[2] PHOBOS Collaboration, B. B. Back et al., Nucl. Phys.
A 757, 28 (2005).
[3] PHENIX Collaboration, K. Adcox et al., Nucl. Phys. A
757 (2005), in press [arXiv:nucl-ex/0410003].
[4] STAR Collaboration, J. Adams et al., Nucl. Phys. A 757
(2005), in press [arXiv:nucl-ex/0501009].
[5] Karsch F, Laermann E, Petreczky P, Stickan S and Wet-
zorke I, 2001 Proccedings of NIC Symposium (Ed. H.
Rollnik and D. Wolf, John von Neumann Institute for
Computing, Jülich, NIC Series, vol.9, ISBN 3-00-009055-
X, pp.173-82,2002.)
[6] P. F. Kolb and U. Heinz, in Quark-Gluon Plasma 3,
edited by R. C. Hwa and X.-N. Wang (World Scientific,
Singapore, 2004), p. 634.
[7] G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev.
Lett. 87, 081601 (2001) [arXiv:hep-th/0104066].
[8] G. Policastro, D. T. Son and A. O. Starinets, JHEP
0209, 043 (2002) [arXiv:hep-th/0205052].
[9] U. Heinz, J. Phys. G 31, S717 (2005).
[10] U. W. Heinz and P. F. Kolb, arXiv:hep-ph/0204061.
[11] C. Eckart, Phys. Rev. 58, 919 (1940).
[12] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Sect.
127, Pergamon, Oxford, 1963.
[13] W. Israel, Ann. Phys. (N.Y.) 100, 310 (1976); W. Israel
and J. M. Stewart, Ann. Phys. (N.Y.) 118, 349 (1979).
[14] A. Muronga, Phys. Rev. Lett. 88, 062302 (2002) [Er-
ratum ibid. 89, 159901 (2002)]; and Phys. Rev. C 69,
034903 (2004).
[15] D. A. Teaney, J. Phys. G 30, S1247 (2004).
[16] A. Muronga and D. H. Rischke, nucl-th/0407114 (v2).
[17] U. W. Heinz, H. Song and A. K. Chaudhuri, Phys. Rev.
C 73, 034904 (2006) [arXiv:nucl-th/0510014].
[18] A. K. Chaudhuri and U. W. Heinz, J. Phys. Conf. Ser.
50, 251 (2006) [arXiv:nucl-th/0504022].
[19] A. K. Chaudhuri, Phys. Rev. C 74, 044904 (2006)
[arXiv:nucl-th/0604014].
[20] A. K. Chaudhuri, arXiv:nucl-th/0703029.
[21] A. K. Chaudhuri, arXiv:nucl-th/0703027.
[22] T. Koide, G. S. Denicol, Ph. Mota and T. Kodama, Phys.
Rev. C 75, 034909 (2007).
[23] R. Baier and P. Romatschke, arXiv:nucl-th/0610108.
[24] P. Arnold, G. D. Moore and L. G. Yaffe, JHEP 0011,
001 (2000) [arXiv:hep-ph/0010177].
[25] G. Baym, H. Monien, C. J. Pethick and D. G. Ravenhall,
Phys. Rev. Lett. 64, 1867 (1990).
[26] S. R. de Groot, W. A. van Leeuwen and Ch. G. van
Weert, Relativistic Kinetic Theory ( North-Holland, Am-
sterdam, 1980) p.36
http://arxiv.org/abs/nucl-ex/0410003
http://arxiv.org/abs/nucl-ex/0501009
http://arxiv.org/abs/hep-th/0104066
http://arxiv.org/abs/hep-th/0205052
http://arxiv.org/abs/hep-ph/0204061
http://arxiv.org/abs/nucl-th/0407114
http://arxiv.org/abs/nucl-th/0510014
http://arxiv.org/abs/nucl-th/0504022
http://arxiv.org/abs/nucl-th/0604014
http://arxiv.org/abs/nucl-th/0703029
http://arxiv.org/abs/nucl-th/0703027
http://arxiv.org/abs/nucl-th/0610108
http://arxiv.org/abs/hep-ph/0010177

  In 2nd order causal dissipative theory, space-time evolution of QGP fluid is
studied in 2+1 dimensions. Relaxation equations for shear stress tensors are
solved simultaneously with the energy-momentum conservation equations.
Comparison of evolution of ideal and viscous QGP fluid, initialized under the
same conditions, e.g. same equilibration time, energy density and velocity
profile, indicate that in a viscous dynamics, energy density or temperature of
the fluid evolve slowly, than in an ideal fluid. Cooling gets slower as
viscosity increases. Transverse expansion also increases in a viscous dynamics.
For the first time we have also studied elliptic flow of 'quarks' in causal
viscous dynamics. It is shown that elliptic flow of quarks saturates due to
non-equilibrium correction to equilibrium distribution function, and can not be
mimicked by an ideal hydrodynamics.

Causal dissipative hydrodynamics for QGP fluid in 2+1 dimensions
A. K. Chaudhuri∗
Variable Energy Cyclotron Centre, 1/AF, Bidhan Nagar, Kolkata 700 064, India
(Dated: February 21, 2013)
In 2nd order causal dissipative theory, space-time evolution of QGP fluid is studied in 2+1 di-
mensions. Relaxation equations for shear stress tensors are solved simultaneously with the energy-
momentum conservation equations. Comparison of evolution of ideal and viscous QGP fluid, ini-
tialized under the same conditions, e.g. same equilibration time, energy density and velocity profile,
indicate that in a viscous dynamics, energy density or temperature of the fluid evolve slowly, than
in an ideal fluid. Cooling gets slower as viscosity increases. Transverse expansion also increases in
a viscous dynamics. For the first time we have also studied elliptic flow of ’quarks’ in causal viscous
dynamics. It is shown that elliptic flow of quarks saturates due to non-equilibrium correction to
equilibrium distribution function, and can not be mimicked by an ideal hydrodynamics.
PACS numbers: 47.75.+f, 25.75.-q, 25.75.Ld
I. INTRODUCTION
One of the most important discoveries in Relativistic
Heavy ion collider (RHIC) at Brokhaven National Lab-
oratory is the large elliptic flow in non-central Au+Au
collisions [1, 2, 3, 4] . Elliptic flow measures the momen-
tum anisotropy of produced particles and is quantified by
the 2nd harmonic of the azimuthal distribution,
v2(pT ) =< cos(2φ) >=
dyd2pT
cos(2φ)dφ
dyd2pT
(1.1)
Elliptic flow is naturally explained in hydrodynamics.
Hydrodynamic pressure is built up from rescattering of
secondaries, and pressure gradients drive the subsequent
collective motion. In non-central Au+Au collisions, ini-
tially, the reaction zone is asymmetric (almond shaped).
The pressure gradient is large in one direction and small
in the other. The asymmetric pressure gradients gener-
ates the elliptic flow. Naturally, in a central collision,
reaction zone is symmetric and elliptic flow vanishes.
Observed elliptic flow then give the strongest indication
that in non-central Au+Au collisions, a collective QCD
matter is produced. Whether the formed matter can be
identified as the much sought after Quark-Gluon Plasma
(QGP) as predicted in Lattice QCD simulations [5] is
presently debatable.
Ideal hydrodynamics has been partly successful in ex-
plaining the observed elliptic flow, quantitatively [6]. El-
liptic flow of identified particles, up to pT ∼1.5 GeV are
well reproduced in ideal hydrodynamics. Ideal hydrody-
namics also explains the transverse momentum spectra
of identified particles (up to pT ∼ 1.5 GeV). Success of
ideal hydrodynamics in explaining bulk of the data [6],
together with the string theory motivated lower limit of
∗E-mail:akc@veccal.ernet.in
shear viscosity η/s ≥ 1/4π [7, 8] has led to a paradigm
that in Au+Au collisions, a nearly perfect fluid is created.
However, the paradigm of ”perfect fluid” produced in
Au+Au collisions at RHIC need to be clarified. As in-
dicated above, the ideal hydrodynamics is only partially
successful and in a limited pT range (pT ≤1.5 GeV) [9].
The transverse momentum spectra of identified particles
also starts to deviate form ideal fluid dynamics prediction
beyond pT ≈ 1.5 GeV. Experimentally determined HBT
radii are not reproduced in the ideal fluid dynamic mod-
els, the famous ”HBT puzzle” [10]. It also do not repro-
duce the experimental trend that elliptic flow saturates
at large transverse momentum. These shortcomings of
ideal fluid dynamics indicate greater importance of dis-
sipative effects in the pT ranges greater than 1.5 GeV
or in more peripheral collisions. Indeed, ideal fluid is a
concept, never realized in nature. As suggested in string
theory motivated models [7, 8], QGP viscosity could be
small, η/s ≥ 1/4π, nevertheless it is non-zero. It is im-
portant to study the effect of viscosity, even if small, on
space-time evolution of QGP fluid and quantify its effect.
This requires a numerical implementation of relativistic
dissipative fluid dynamics. Furthermore, if QGP fluid
is formed in heavy ion collisions, it has to be charac-
terized by measuring its transport coefficients, e.g. heat
conductivity, bulk and shear viscosity. Theoretically, it
is possible to obtain those transport coefficients in a ki-
netic theory model. However, in the present status of
theory, the goal can not be achieved immediately, even
more so for a strongly interacting QGP (sQGP). Alter-
natively, one can use the experimental data to obtain
a ”phenomenological” limit of transport coefficients of
sQGP. It will also require a numerical implementation of
relativistic dissipative fluid dynamics. There is another
incentive to study dissipative hydrodynamics. Ideal hy-
drodynamics depends on the assumption of local equilib-
rium. Before local equilibrium is attained, the system has
to pass through a non-equilibrium stage, where (if non-
equilibrium effects are small) dissipative hydrodynamics
may be applicable. Indeed, we can explore early times of
fluid evolution better in a dissipative hydrodynamics.
http://arxiv.org/abs/0704.0134v2
mailto:akc@veccal.ernet.in
Theory of dissipative relativistic fluid has been formu-
lated quite early. The original dissipative relativistic fluid
equations were given by Eckart [11] and Landau and Lif-
shitz [12]. They are called 1st order theories. Formally,
relativistic dissipative hydrodynamics are obtained from
an expansion of entropy 4-current, in terms of dissipative
fluxes. In 1st order theories, entropy 4-current contains
terms linear in dissipative quantities. 1st order theory
of dissipative hydrodynamics suffer from the problem of
causality violation. Signal can travel faster than light.
Causality violation is unwarranted in any theory, even
more in a relativistic theory. The problem of causal-
ity violation is removed in the Israel-Stewart’s 2nd order
theory of dissipative fluid [13]. In 2nd order theory, ex-
pansion of entropy 4-current contains terms 2nd order in
dissipative fluxes. However, these leads to complications
that dissipative fluxes are no longer function of the state
variables only. They become dynamic. The space of ther-
modynamic variables has to be extended to include the
dissipative fluxes (e.g. heat conductivity, bulk and shear
viscosity).
Even though 2nd order theory was formulated some
30 years back, significant progress towards its numer-
ical implementation has only been made very recently
[14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. At the Cyclotron
Centre, Kolkata, we have developed a code ”AZHYDRO-
KOLKATA” to simulate the hydrodynamic evolution of
dissipative QGP fluid. Presently only dissipative ef-
fect included is the shear viscosity. Some results of
AZHYDRO-KOLKATA, for first order dissipative hydro-
dynamics have been published earlier [19, 20, 21]. In the
present paper, for the first time, we will present some
results for 2nd order dissipative hydrodynamics in 2+1
dimensions. In the present paper, we will consider effect
of dissipation in the QGP phase only. Effect of phase
transition will be studied in a later publication.
The paper is organized as follows: In section II we
briefly review relativistic dissipative fluid dynamics. In
section III we derive the relevant equations in 2+1 di-
mension (assuming boost-invariance). Required inputs
e.g. the equation of state, viscosity coefficient and initial
conditions are discussed in section IV. Simulation re-
sults from AZHYDRO-KOLKATA are shown in section
V. In section VII we compare the transverse momentum
spectra and elliptic flow of quarks in ideal and viscous
dynamics. The concluding section IX summarizes our
results.
II. DISSIPATIVE FLUID DYNAMICS
In this section, I briefly discuss the phenomenological
theory of dissipative hydrodynamics. More detailed ex-
position can be found in [13].
A simple fluid, in an arbitrary state, is fully speci-
fied by primary variables: particle current (Nµ), energy-
momentum tensor (T µν) and entropy current (Sµ) and
a number of additional (unknown) variables. Primary
variables satisfies the conservation laws;
µ = 0, (2.1)
µν = 0, (2.2)
and the 2nd law of thermodynamics,
µ ≥ 0. (2.3)
In relativistic fluid dynamics, one defines a time-like
hydrodynamic 4-velocity, uµ (normalized as u2 = 1).
One also define a projector, ∆µν = gµν − uµuν , orthog-
onal to the 4-velocity (∆µνuν = 0). In equilibrium, an
unique 4-velocity (uµ) exists such that the particle den-
sity (n), energy density (ε) and the entropy density (s)
can be obtained from,
Nµeq = nuµ (2.4)
T µνeq = εu
µuν − p∆µν (2.5)
Sµeq = suµ (2.6)
An equilibrium state is assumed to be fully specified
by 5-parameters, (n, ε, uµ) or equivalently by the thermal
potential, α = µ/T (µ being the chemical potential) and
inverse 4-temperature, βµ = uµ/T . Given a equation of
state, s = s(ε, n), pressure p can be obtained from the
generalized thermodynamic relation,
Sµeq = pβ
µ − αNµeq + βλT λµeq (2.7)
Using the Gibbs-Duhem relation, d(pβµ) = Nµeqdα −
T λµeq dβλ, following relations can be established on the
equilibrium hyper-surface Σeq(α, β
dSµeq = −αdNµeq + βλdT λµeq (2.8)
In a non-equilibrium system, no 4-velocity can be
found such that Eqs.2.4,2.5,2.6 remain valid. Tensor de-
composition leads to additional terms,
Nµ = Nµeq + δN
µ = nuµ + V µ (2.9)
T µν = T µνeq + δT
= [εuµuν − p∆µν ] + Π∆µν + πµν
+(Wµuν +W νuµ) (2.10)
Sµ = Sµeq + δS
µ = suµ + Φµ (2.11)
The new terms describe a net flow of charge V µ =
∆µνNν , heat flow, W
µ = (ε + p)/nV µ + qµ (where qµ
is the heat flow vector), and entropy flow Φµ. Π =
µν − p is the bulk viscous pressure and πµν =
(∆µσ∆ντ +∆νσ∆µτ − 1
∆µν∆στ ]Tστ is the shear stress
tensor. Hydrodynamic 4-velocity can be chosen to elimi-
nate either V µ (the Eckart frame, uµ is parallel to particle
flow) or the heat flow qµ (the Landau frame, uµ is par-
allel to energy flow). In relativistic heavy ion collisions,
central rapidity region is nearly baryon free and Lan-
dau’s frame is more appropriate than the Eckart’s frame.
Dissipative flows are transverse to uµ and additionally,
shear stress tensor is traceless. Thus a non-equilibrium
state require 1+3+5=9 additional quantities, the dissi-
pative flows Π, qµ (or V µ) and πµν . In kinetic theory,
Nµ and T µν are the 1st and 2nd moment of the distri-
bution function. Unless the function is known a-priori,
two moments do not furnish enough information to enu-
merate the microscopic states required to determine Sµ,
and in an arbitrary non-equilibrium state, no relation
exists between, Nν , T µν and Sµ. Only in a state, close
to a equilibrium one, such a relation can be established.
Assuming that the equilibrium relation Eq.2.8 remains
valid in a ”near equilibrium state” also, the entropy cur-
rent can be generalized as,
Sµ = Sµeq + dS
µ = pβµ − αNµ + βλT λµ +Qµ (2.12)
where Qµ is an undetermined quantity in 2nd order in
deviations, δNµ = Nµ − Nµeq and δT µν = T µν − T µνeq .
Detail form of Qµ is constrained by the 2nd law ∂µS
0. With the help of conservation laws and Gibbs-Duhem
relation, entropy production rate can be written as,
µ = −δNµ∂µα+ δT µν∂µβν + ∂µQµ (2.13)
Choice of Qµ leads to 1st order or 2nd order theories
of dissipative hydrodynamics. In 1st order theories the
simplest choice is made, Qµ = 0, entropy current con-
tains terms up to 1st order in deviations, δNµ and δT µν.
Entropy production rate can be written as,
µ = ΠX − qµXµ + πµνXµν (2.14)
where, X = −∇.u; Xµ = ∇
− uν∂νuµ and Xµν =
∇<µuν>.
The 2nd law, ∂µS
µ ≥ 0 can be satisfied by postulat-
ing a linear relation between the dissipative flows and
thermodynamic forces,
Π = −ζθ, (2.15)
qµ = −λ nT
∇µ(µ/T ), (2.16)
πµν = 2η∇<µuν> (2.17)
where ζ, λ and η are the positive transport coefficients,
bulk viscosity, heat conductivity and shear viscosity re-
spectively.
In 1st order theories, causality is violated. If, in a given
fluid cell, at a certain time, thermodynamic forces vanish,
corresponding dissipative fluxes also vanish instantly. Vi-
olation of causality is unwanted in any theory, even more
so in relativistic theory. Causality violation of dissipative
hydrodynamics is corrected in 2nd order theories [13]. In
2nd order theories, entropy current contain terms up to
2nd order in the deviations, Qµ 6= 0. The most general
Qµ containing terms up to 2nd order in deviations can
be written as,
Qµ = −(β0Π2−β1qνqν+β2πνλπνλ)
(2.18)
As before, one can cast the entropy production rate
(T∂µS
µ) in the form of Eq.2.14. Neglecting the terms
involving dissipative flows with gradients of equilib-
rium thermodynamic quantities (both are assumed to be
small) and demanding that a linear relation exists be-
tween the dissipative flows and thermodynamic forces,
following relaxation equations for the dissipative flows
can be obtained,
Π = −ζ(θ + β0DΠ) (2.19)
qµ = −λ
) − β1Dqµ
(2.20)
πµν = 2η
∇<µuν> − β2Dπµν
, (2.21)
where D = uµ∂µis the convective time derivative. Unlike
in the 1st order theories, in 2nd order theories, dynamical
equations control the dissipative flows. Even if thermo-
dynamic forces vanish, dissipative flows do not vanish
instantly.
Before we proceed further, it may be mentioned that
the parameters, α and βλ are not connected to the actual
state (Nµ, T µν). The pressure p in Eq.2.12 is also not
the ”actual” thermodynamics pressure, i.e. not the work
done in an isentropic expansion. Chemical potential α
and 4-inverse temperature βλ has meaning only for the
equilibrium state. Their meaning need not be extended
to non-equilibrium states also. However, it is possible to
fit a fictitious ”local equilibrium” state, point by point,
such that pressure p in Eq.2.12 can be identified with
the thermodynamic pressure, at least up to 1st order.
The conditions of fit fixes the underlying non-equilibrium
phase-space distribution.
III. (2+1)-DIMENSIONAL VISCOUS
HYDRODYNAMICS WITH LONGITUDINAL
BOOST INVARIANCE
Complete dissipative hydrodynamics is a numerically
challenging problem. It requires simultaneous solution of
14 partial differential equations (5 conservation equations
and 9 relaxation equations for dissipative flows). We re-
duce the problem to solution of 6 partial differential equa-
tions (3 conservation equations and 3 relaxation equa-
tions). In the following, we will study boost-invariant
evolution of baryon free QGP fluid, including the dissi-
pative effect due to shear viscosity only. Shear viscosity
is the most important dissipative effect. For example, in
a baryon free QGP, heat conduction is zero and we can
disregard Eq.2.20. Bulk viscosity is also zero for the QGP
fluid (point particles) and Eq.2.19 can also be neglected.
Shear pressure tensor has 5 independent components but
the assumption of boost invariance reduces the number
of independent components to three. For a baryon free
fluid, we can also disregard the conservation equation
Eq.2.1. With the assumption of boost-invariance, energy-
momentum conservation equation ∂µT
µη = 0 become re-
dundant and only three energy-momentum conservation
equations are required to be solved.
-10 -5 0 5 10
X (fm)
-10 -5 0 5 10
_______ η/s=0.08
_______ η/s=0.135
_______ η/s=0
-10 -5 0 5 10
energy density in x-y plane, τ =2.6 fm.
FIG. 1: (color online). Constant energy density contours in
x-y plane at τ=2.6 fm. The black lines are for ideal fluid
(η/s=0). The red and blue lines are for viscous fluid with
ADS/CFT and perturbative estimate of viscosity, η/s=0.08
and 0.135.
Heavy ion collisions are best described in (τ, x, y, η)
coordinates, where τ =
t2 − z2 is the longitudinal
proper time and η = 1
ln t+z
is the space-time rapid-
ity. r⊥ = (x, y) are the usual cartisan coordinate in
the plane, transverse to the beam direction. Relevant
equations concerning this coordinate transformations are
given in the appendix A.
Explicit equations for energy-momentum conservation
in (τ ,x,y,η) coordinates are given in the appendix B. We
note that unlike in ideal fluid, in viscous fluid dynam-
ics, conservation equations (see Eqs.B1-B3) contain addi-
tional pressure gradients due to shear viscosity. Both T τx
and T τy components of energy-momentum tensor now
evolve under additional pressure gradients. The right-
most term of Eq.B3 also indicate that in viscous dynam-
ics, longitudinal pressure is effectively reduced (note that
the πηη component is negative). Since pressure can not
be negative, shear viscosity is limited by the condition,
p+ τ2πηη ≥ 0.
As evident from the Eqs.B1-B3, in boost-invariant dis-
sipative hydrodynamics, with shear viscosity taken into
account, fluid evolution depends only on seven compo-
nents of the shear stress tensors. They are πττ , πτx,
πτy, πxx, πyy, πxy and πηη. However, all the seven com-
ponents are not independent. Tracelessness, transver-
sality to uµ and the assumption of boost-invariance re-
duces the independent components to three. Presently,
we choose πxx πyy and πxy as the independent compo-
nents. Relaxation equations for the independent compo-
nents are given in the appendix C (see Eqs.C4-C6). They
are solved simultaneously with the three energy momen-
tum conservation equations Eqs.B1-B3, with inputs as
discussed below.
-10 -5 0 5 10
X (fm)
-10 -5 0 5 10
_______ η/s=0.08
_______ η/s=0.135
_______ η/s=0
-10 -5 0 5 10
energy density in x-y plane, τ =8.6 fm
FIG. 2: (color online). same as Fig.1 but at time τ=8.6 fm.
IV. EQUATION OF STATE, VISCOSITY
COEFFICIENT AND INITIAL CONDITIONS
A. Equation of state
One of the most important inputs of a hydrodynamic
model is the equation of state. Through this input, the
macroscopic hydrodynamic models make contact with
the microscopic world. In the present demonstrative cal-
culation we will show results for the QGP phase only.
In the QGP phase, we use the simple equation of state,
p = 1
ε, with energy density given as,
gqgpT
4 (4.1)
where gqgp = ggluon+
gquark is the degeneracy factor for
QGP. ggluon = 2(helicity) × 8(color) is the degeneracy
factor for gluons and gquark = 2(spin)× 3(color)× 2(q+
q̄) ×Nf is the degeneracy factor for Nf flavored quarks.
For Nf ≈ 2.5, the degeneracy factor is gqgp = 42.25
B. Shear viscosity coefficient
Shear viscosity coefficient (η) of QGP or sQGP is quite
uncertain. In a strongly coupled QGP, shear viscosity
can not be computed. Recently, using the ADS/CFT
correspondence [7, 8] shear viscosity of a strongly coupled
gauze theory, N=4 SUSY YM, has been evaluated, η =
N2c T
3 and the entropy is given by s = π
N2c T
3. Thus
in the strongly coupled field theory,
ADS/CFT
≈ 0.08, (4.2)
Shear viscosity is quite uncertain in perturbative QCD
also. At high temperature, shear viscosity, in leading log,
can be written as [24, 25],
η = κ
g4 ln g−1
, (4.3)
where g is the strong coupling constant. The leading
log shear viscosity coefficient κ depend on the number of
fermion flavors (Nf ). For example, for two flavored QGP,
κ = 86.47 and κ = 106.7 for a three flavored QGP. With
entropy density of QGP, s = π
gqgpT
3. For two flavored
QGP and αs ≈0.5, the ratio of viscosity over the entropy,
in the perturbative regime is estimated as,
≈ 0.135, (4.4)
For lower αs , perturbative estimation of η/s could be
even higher.
Shear viscosity can also be expressed in terms of sound
attenuation length, Γs, defined as,
(4.5)
Γs is equivalent to mean free path and for a valid hy-
drodynamic description Γs/τ << 1, i.e. mean free path
is much less than the system size. Initial conditions of
the fluid must be chosen carefully such that the validity
condition Γs/τ << 1 remains valid initially as well as at
later time also. In the present work, we have treated vis-
cosity as a parameter. To explore the effect of viscosity,
we have used both the ADS/CFT estimate η/s=0.08 and
perturbative estimate η/s=0.135. We have also run the
code with a higher value of viscosity η/s=0.2.
C. Initial conditions
Solution of Eqs.B1-B3 require initial conditions, the
initial time τi, the transverse distribution of energy den-
sity ε(x, y) and the velocities vx(x, y) and vy(x, y). Fol-
lowing [6], initial transverse energy density is parame-
terized geometrically. At an impact parameter ~b, trans-
verse distribution of wounded nucleons NWN (x, y,~b) and
-10 -5 0 5 10
X axis Title
-0.78
-0.56
-0.33
-0.11
-10 -5 0 5 10
contour plot of Vx at τ=8.6 fm
X axis Title
-10 -5 0 5 10
X axis Title
FIG. 3: (color online). contours of constant vx in x-y plane
at τ=8.6 fm. The black lines are for ideal fluid (η/s=0). The
red and blue lines are for viscous fluid with ADS/CFT and
perturbative estimate of viscosity, η/s=0.08 and 0.135.
-10 -5 0 5 10
X axis Title
-0.78
-0.56
-0.33
-0.11
-10 -5 0 5 10
contour plot of Vy at τ=8.6 fm
X axis Title
-10 -5 0 5 10
X axis Title
FIG. 4: (color online). contours of constant vy in x-y plane
at τ=8.6 fm. The black lines are for ideal fluid (η/s=0). The
red and blue lines are for viscous fluid with ADS/CFT and
perturbative estimate of viscosity, η/s=0.08 and 0.135.
of binary NN collisions NBC(x, y,~b) to are calculated in
a Glauber model. A collision at impact parameter ~b is
assumed to contain 25% hard scattering (proportional
to number of binary collisions) and 75% soft scattering
(proportional to number of wounded nucleons). Trans-
verse energy density profile at impact parameter~b is then
obtained as,
ε(x, y,~b) = ε0(0.75×NWN (x, y,~b) + 0.25×NBC(x, y,~b))
(4.6)
with central energy density ε0=30GeV/fm
−3. The equi-
libration time is chosen as τi=0.6 fm [6]. The initial ve-
locities vx and vy are assumed to be zero initially.
______ η/s=0.08
______ η/s=0
______ η/s=0
Temperature in τ-x plane. y=0.
X (fm)0 2 4 6 8 10
FIG. 5: (color online). Constant temperature contour in x −
τ plane, for fixed y=0. The black, red and blue lines are
for ideal, viscous fluid with η/s=0.08 and viscous fluid with
η/s=0.135.
In dissipative hydrodynamics, one requires initial con-
ditions for the viscous pressures also. Due to longitudinal
boost invariance of the problem, we assume that viscous
pressures have attained their boost-invariant values at
the time of equilibration. Boost invariant values of the
three independent shear stress-tensors can be easily ob-
tained from Eqs.C4-C6, σxx = σyy = θ = 1
and σxy = 0
( at the initial time τi, u
µ = (1, 0, 0, 0), Duµ = 0. The
initial distribution of shear pressure tensors are then ob-
tained as,
πxx(x, y,~b) = 2ησxx = 2η/τi (4.7)
πyy(x, y,~b) = 2ησyy = 2η/τi (4.8)
πxy(x, y,~b) = 2ησxy = 0 (4.9)
As stated earlier, the viscous coefficient η is obtained
using the relation, η/s = const, const=0.08, 0.135 or 0.2.
For these values of shear viscosity, the validity condition
Γs/τ << 1 is satisfied initially. The validity condition is
better satisfied at later time.
V. RESULTS
VI. STABILITY OF NUMERICAL SOLUTIONS
A. Evolution of the viscous QGP fluid
The energy-momentum conservation equations B1-B3,
and the relaxation equations C2-C4 are solved simultane-
ously using the code, AZHYDRO-KOLKATA, developed
at the Cyclotron Centre, Kolkata. As mentioned earlier,
we have solved the equations in the QGP phase only and
did not consider any phase transition. In the following
we will show the results for central Au+Au collisions (im-
pact parameter b = 0 fm). To understand the effect of
shear viscosity, with the same initial conditions, we have
solved the energy-momentum conservation equations for
ideal fluid and viscous fluid. As mentioned earlier, we
have considered two values of viscosity, the ADS/CFT
motivated value η/s=0.08 and the perturbative estimate,
η/s=0.135.
0 2 4 6 8 10
X (fm)
_____ η/s=0.08
_____ η/s=0.135
_____ η/s=0
Temperature   in τ-x plane. y=5 fm.
FIG. 6: (color online). same as fig.5 but at y=5.0
In Fig.1, we have shown the contours of constant en-
ergy density in x-y plane, after an evolution of 2.6 fm.
The black lines are for ideal fluid evolution. The red and
blue lines are for viscous fluid with ADS/CFT (η/s=0.08)
and perturbative (η/s=0.135) estimate of viscosity. Con-
stant energy density contours, as depicted in Fig.1, indi-
cate that with viscosity fluid cools slowly. Cooling gets
slower as viscosity increases. Thus at any point in the x-y
plane, energy density of viscous fluid is higher than that
of an ideal fluid. At later time also, compared to an ideal
fluid, viscous fluid evolve slowly. In Fig.2, contours of
constant energy density at time τ=8.6 fm is shown. Here
also we find than at any point energy density of viscous
fluid is higher than its ideal counter part. The result is
in accordance with our expectation. For dissipative fluid,
equation of motion can be written as,
Dε = −(ε+ p)∇µuµ + πµν∇<µuν> (6.1)
Due to viscosity, evolution of energy density (or tem-
perature) is slowed down.
In Fig.3 and 4, we have shown the contour plot of the
fluid velocity, vx and vy, after evolution of 8.6 fm. As
before the black lines are for the ideal fluid evolution.
The red and blue lines are for viscous fluid with η/s=0.08
and 0.135 respectively. Fluid velocities in viscous and
ideal fluid differ very little. Even at late time, as shown
in Fig.3 and 3, we find that for η/s=0.08-0.135, x and y
component of the fluid velocity show marginal difference.
However, there is an indication that in a viscous fluid,
velocity grow faster than in ideal fluid. But as mentioned
earlier, the difference is marginal.
η/s=0.08
η/s=0.135
τ (fm)
0 2 4 6 8 10
FIG. 7: (color online). In the upper panel, temporal evolu-
tion of the shear pressure tensor πxx at the fluid cell x=y=0
is shown. In the lower panel, evolution of πxy at the fluid
cell x=y=5 fm is shown. The black and red lines are for
ADS/CFT motivated viscosity η/s=0.08 and perturbative es-
timate η/s=0.135 respectively.
As seen in Fig.1-2, in viscous dynamics, QGP fluid
evolves slowly. Thus life-time of the QGP phase is en-
hanced in viscous dynamics. To obtain an idea about the
enhanced life-time, in Fig.5, we have shown the constant
temperature contours in τ − x plane , at a fixed value
of y=0 fm. As seen in Fig.5, temperature evolves slowly
in a viscous fluid and life-time of the QGP phase is ex-
tended. For small viscosity η/s=0.08-0.135, the increase
is not large. At the center of the fluid, for η/s=0.135,
QGP life-time is increased approximately by 5% only. It
is even less for the ADS/CFT estimate of viscosity. How-
ever, enhancement of QGP life-time depends on the fluid
cell position. It could be more. In Fig.6, constant tem-
perature contours at y=5 fm is shown. For η/s=0.135, at
x=0,y=5 fm, the QGP life-time is enhanced by ∼ 10%.
We conclude that in a viscous dynamics, with moder-
ate viscosity η/s=0.08-0.135, QGP life-time could be en-
hanced by 5-10%. Enhanced lifetime of QGP in a viscous
fluid can have significant effect on observables produced
early in the collisions e.g. direct photon production or in
J/ψ suppression.
1.01.21.4
0 2 4 6 8 10
0.030
0.180.21
0 2 4 6 8 10
(d) π
 at τ=2.6 fm(c) π
 at τ=0.6 fm
(b) π
 at τ=2.6 fm(a) π
 at τ=0.6 fm
1.01.21.4
0 2 4 6 8 10
0.030
0.150.18
0.210.24
0 2 4 6 8 10
X (fm)
FIG. 8: (color online). In panel (a) and (b), contours of
constant pressure tensor πxx at initial time τi=0.6 fm and at
time τ=2.6 fm is shown. In panel (c) and (d) same results for
shear pressure tensor πyy is shown.
B. Evolution of shear pressure tensors
We have assumed that initially the shear pressure ten-
sors πxx, πyy and πxy attained their longitudinal boost-
invariant values. As the fluid evolve, pressure tensors also
evolve. Here we investigate the evolution of shear pres-
sure tensors with time. In the top panel of Fig.7 evolu-
tion of shear pressure tensor πxx at the fluid cell position
x=y=0 is shown. The black line is for the ADS/CFT
motivated viscosity, η/s=0.08 and the red line is for the
perturbative estimate of viscosity η/s=0.135. Just after
the start of the evolution the shear pressure tensor πxx
increases, but for a short duration and then steadily de-
creases with time. By 4 fm of evolution, πxx at the center
of the fluid reduces to negligibly small values. Identical
behavior is seen for the shear pressure tensor πyy. In
the bottom panel of Fig.7 we have shown the evolution
of the third independent shear pressure tensor πxy. Ini-
tially πxy is zero. As the fluid evolve, it grow in the
negative direction. We find that at the centre of the fluid
(x=y=0), it never grows. In Fig.7, temporal evolution
of πxy at the fluid cell position x = y = 5fm is shown.
From the initial zero value, πxy rapidly increases in the
negative direction. It reaches its maximum around τ ≈1
fm and then decreases again. We also note that πxy never
grows to large values. Compared to πxx or πyy stress ten-
sor πxy is negligible. The results indicate that in a QGP
fluid, viscous effect persist for a short duration (3-4 fm)
only. At late time the fluid evolve essentially as an ideal
fluid. The result is understandable. Shear viscosity de-
pend strongly on temperature (η ∝ T 3). As the fluid
cools, effect of viscosity decreases rapidly.
τ (fm)
0 2 4 6 8 10
η/s=0.08
η/s=0.135
FIG. 9: Evolution of average entropy with time, for two val-
ues of viscosity, the ADS/CFT motivated viscosity η/s=0.08
and perturbative estimate η/s=0.135 are shown.
To show the spatial distribution of the stress tensors,
in Fig.8, πxx and πyy at initial time τi=0.6 fm and after
an evolution of τ = 2.6 fm are shown. As shown earlier,
πxx and also πyy rapidly decreases with time. By 2 fm
of evolution they are reduced by approximately by a fac-
tor 6. It is also interesting to note that the initial x-y
symmetric distribution of πxx and πyy quickly evolves to
asymmetric distribution. With time πxx evolves faster in
the x-direction than in y-direction. Similarly, πyy evolve
faster in the y-direction than in the x-direction. For cen-
tral collisions the asymmetric evolution of πxx and πyy
counter balance each other. As shown in Fig.1 and 2,
the contour plots of energy density do not show any in-
dication of asymmetry even at late time. However, the
asymmetric pressure tensors can have important effects
on elliptic flow of observables produced early in the col-
lisions, say in elliptic flow of direct photons.
C. Entropy generation
In a viscous fluid dynamics, entropy is generated. We
can easily calculate the entropy generated during the evo-
lution,
πµνπµν
[(πττ )2 + (πxx)2 + (πyy)2 + (τ2πηη)2
−2(πτx)2 − 2(πτy)2 + 2(πxy)2] (6.2)
Evolution of spatially averaged entropy is shown in
Fig.9, for the two values of viscosity coefficients η/s=0.08
and 0.135. As expected, entropy generation is more if vis-
cosity is more. For both the values of viscosity, we find
that entropy generation saturates after ≈ 3 fm of evolu-
tion. It is expected also. As shown previously, viscous
fluxes reduces to very small values after τ=3 fm. Natu-
rally, entropy generation is negligible thereafter.
0.190.22
 0 2 4 6 8 10
10 Temperature in τ-x plane, y= 5m.
_____ viscous fluid (2nd order)
_____ ideal fluid
_____ viscous fluid (1st order)
X (fm)
FIG. 10: (color online) constant temperature contours in
x − τ plane at y=5 fm. The black lines are for ideal fluid.
The red and blue lines are for viscous fluid in 1st order and
in 2nd order theory respectively. η/s=0.135.
D. 1st order theory vs. 2nd order theory
As mentioned earlier, 1st order theory of dissipative
hydrodynamics is acausal, signal can travel faster than
light. This is corrected in 2nd order theory, but we
have to pay the price, relaxation equations for dissipa-
tive fluxes are required to be solved. It is interesting to
τ (fm)
0 2 4 6 8 10 12
1st order theory
2nd order theory
η/s=0.135
FIG. 11: Evolution of average entropy production in a 1st
order (solid line) and 2nd order (dashed line) theory. 2nd
order theory generate more entropy.
compare the difference we can expect in a first order the-
ory and in a 2nd order theory of dissipation. In Fig. 10,
we have shown the contours of constant temperature in
x−τ , for a fixed y = 5fm. The black lines are for an ideal
fluid. The red lines are for a viscous fluid treated in the
1st order theory. The blue lines are for viscous fluid in
2nd order theory. In 2nd order theory fluid evolve more
slowly than in a first order theory. Entropy generation is
also more in a 2nd order theory. In Fig.11, average en-
tropy evolution with proper time is shown, both for the
1st order theory (the solid line) and the 2nd order theory.
In 2nd order theory, approximately 80% more entropy is
generated.
VII. TRANSVERSE MOMENTUM AND
ELLIPTIC FLOW OF QUARKS
Presently we can not compare predictions from viscous
hydrodynamics with experimental data. Hadrons are not
included in the model. The initial QGP fluid evolve and
cools but remain in the QGP phase, it did not undergo
a phase transition to hadronic gas. However, from the
momentum distribution of quarks we can get some idea
about the viscous effect on particle production. Viscos-
ity generates entropy, which will be reflected in enhanced
multiplicity. We use the standard Cooper-Frey prescrip-
tion to obtained the transverse momentum distribution
of quarks. In Cooper-Frey prescription, particle distribu-
tion is obtained by convoluting the one body distribution
function over the freeze-out surface,
µf(x, p) (7.1)
where dΣµ is the freeze-out hyper-surface and f(x, p) is
the one-body distribution function. Now in a viscous dy-
namics, the fluid is not in equilibrium and f(x, p) can not
be approximated by the equilibrium distribution func-
tion,
f (0)(x, p) =
exp[β(uµpµ − µ)] ± 1
, (7.2)
with inverse temperature β = 1/T and chemical poten-
tial µ. In a highly non-equilibrium system, distribution
function f(x, p) is unknown. If the system is slightly
off-equilibrium, then it is possible to calculate correction
to equilibrium distribution function due to (small) non-
equilibrium effects. Slightly off-equilibrium distribution
function can be approximated as,
f(x, p) = f (0)(x, p)[1 + φ(x, p)], (7.3)
φ(x, p) is the deviation from equilibrium distribution
function f (0). With shear viscosity as the only dissi-
pative forces, φ(x, p) can be locally approximated by a
quadratic function of 4-momentum,
φ(x, p) = εµνp
µpν . (7.4)
Without any loss of generality εµν can be written as
εµν =
2(ε+ p)T 2
πµν , (7.5)
completely specifying the non-equilibrium distribution
function. As expected, correction factor increases with
increasing viscosity. We also note that non-equilibrium
correction is more on large momentum particles. The
effect of viscosity is more on large momentum parti-
cles. The correction factor reduces if freeze-out occurs
at higher temperature.
With the corrected distribution function, we can cal-
culate the quark momentum spectra at freeze-out sur-
face Σµ. In appendix D, relevant equations are given.
The quark momentum distribution has two parts, (i)
dyd2pT
, obtained by convoluting the equilibrium distribu-
tion function over the freeze-out surface and (ii) dN
dyd2pT
obtained by convoluting the correction to the equilibrium
distribution function over the freeze-out surface. Since
the correction factor is obtained under the assumption
that non-equilibrium effects are small, φ(x, p) << 1, it
necessarily imply that, dN
dyd2pT
<< dN
dyd2pT
. The ratio,
dNneq
dNneq
dyd2pT
dyd2pT
, (7.6)
could at best be unity or less. If the ratio exceeds
unity, it will imply that non-equilibrium effects are large
 (GeV)
0 1 2 3 4 5
η/s=0.2 
η/s=0.135 
η/s=0.08 
Au+Au@b=6.5 fm
τi=0.6fm,Sini=110fm-3,TF=160MeV
FIG. 12: Ratio of quark spectra with non-equilibrium distri-
bution function to that with equilibrium distribution function.
and the distribution function f(x, p) can not be approx-
imated as in Eq.7.3. Using AZHYDRO-KOLKATA, we
have simulated a b=6.5 fm Au+Au collision. dN
dyd2pT
dNneq
dyd2pT
at freeze-out temperature TF =160 MeV are cal-
culated. The ratio dN
for η/s=0.08,0.135 and 0.2, are
shown in Fig.12. With ADS/CFT estimate of viscosity,
η/s=0.08, non-equilibrium correction to particle produc-
tion become comparable to equilibrium contribution only
beyond pT =5 GeV. However, with perturbative estimate,
η/s=0.135, non-equilibrium correction become compara-
ble to or exceeds the equilibrium contribution at pT ∼
4 GeV. pT range is further reduced for higher viscosity
η/s=0.2. Thus with perturbative estimate of viscosity
(η/s = 0.135 − 0.2), hydrodynamic description remain
valid upto transverse momentum pT ∼ 3.5-4 GeV.
In Fig.13, we have compared the transverse momen-
tum spectra of quarks in ideal hydrodynamics with that
in a viscous dynamics. In Fig.13, the dotted line is the
spectra obtained in ideal dynamics (η/s = 0). The pT
spectra in viscous dynamics are shown by black lines.
We have shown the spectra for three values of viscos-
ity η/s=0.08,0.135 and 0.2. Compared to ideal dynam-
ics, quarks yield in viscous dynamics increases. The in-
crease is more at large pT . For low values of viscosity
the increase is modest, a factor of 2 at pT =3 GeV. But
yield increase by a factor or 4(10) if viscosity increases
to η/s=0.135 (2). Please note that even though we have
shown pT spectra upto 5 GeV, for η/s=0.2 and 0.135,
hydrodynamic description fails beyond pT ∼ 3.5 and 4
We have also studied the effect of viscosity on quark
elliptic flow. Effect of viscosity is very prominent on
elliptic flow. In Fig.14, pT dependence of elliptic flow
of quarks, in a b=6.5 fm collision is shown. The black
line is v2 in ideal dynamics. In ideal dynamics, ellip-
pT (GeV)
0 1 2 3 4 5
η/s=0.0(ID.Fluid)
η/s=0.2 
η/s=0.135 
η/s=0.08 
Au+Au@b=6.5 fm
τi=0.6fm,Sini=110fm-3,TF=160MeV
FIG. 13: Qurak transverse momentum spectra at freeze-out
temperature of 160 MeV. The dotted line is the quarks spectra
in ideal hydrodynamics. The solid lines (top to bottom) are
in viscous dynamics with η/s=0.08,0.135 and 0.2.
tic flow continually increases with pT . It well known, in
contrast to experiments, where elliptic flow saturates at
large pT , in ideal hydrodynamics, elliptic flow continue
to increase with pT . Indeed, this is a major problem in
ideal hydrodynamics. The renewed the interest in dis-
sipative hydrodynamics is partly due to the inability of
ideal hydrodynamics to predict the trend of elliptic flow
in Au+Au collisions. In Fig.14, the blue lines are v2
in viscous dynamics with η/s=0.08,0.135 and 0.2 respec-
tively. In a viscous dynamics, pT dependence of v2 is
drastically changed. In contrast to ideal dynamics where
v2 continue to increase with pT , in viscous dynamics, v2
continue to increase only upto pT ∼ 1.5− 2GeV . There-
after v2 decreases. For perturbative estimate of viscos-
ity η/s=0.135 and beyond, v2 even become negative at
large pT . Veering about of v2 after pT ∼1.5-2 GeV is
due to viscous effect only or more explicitly due to the
non-equilibrium correction to the equilibrium distribu-
tion function. This is clearly manifested from the red
lines in Fig.14. The red lines are calculated ignoring the
non-equilibrium corrections to the equilibrium distribu-
tion function. If non-equilibrium correction is ignored, in
viscous dynamics also, v2 continue to increase with pT ,
albeit its magnitude is reduced compared to ideal dy-
namics. The result is very important. It imply that the
experimental trend of elliptic flow (saturation at large
pT ) could only be explained if the QGP fluid is viscous.
An ideal QGP, will not be able to explain the saturation
trend of the experimental data.
As stated earlier, non-equilibrium correction to the
equilibrium distribution function depends on the freeze-
out condition. To show the effect of freeze-out condition,
on v2, in Fig.15 we have shown v2 for a values of freeze-
out temperature TF =160,150,140,130 and 120 MeV. As
pT (GeV)
0 1 2 3 4 5
vis.fluid with Fneq
vis.fluid with Feq
id. fluid
Au+Au@b=6.5 fm
τi=0.6fm,Sini=110fm-3,TF=160MeV
FIG. 14: (color online) Elliptic flow as a function of transverse
momentum. The black line is v2 in ideal hydrodynamics. The
blue lines are v2 in viscous dynamics with viscosity to entropy
ratio η/s=0.08,0.135 and 0.2 (top to bottom) respectively,
including the correction to equilibrium distribution function.
The red lines are same as the blues lines but ignoring the
non-equilibrium correction to the distribution function.
freeze-out occur at higher and higher temperature, the
veering of v2 takes place at larger and larger pT and for
TF =120 MeV, the elliptic flow saturates. The result is
understood easily. With decreasing freeze-out tempera-
ture, the fluid evolves for longer time, the shear stress-
tensor’s at the freeze-out surface is reduced and the non-
equilibrium correction, proportional to shear stress ten-
sors, decreases.
VIII. STABILITY OF NUMERICAL
SOLUTIONS IN AZHYDRO-KOLKATA
Before we summarise our results, we would like to
comment on the stability of numerical solutions in
AZHYDRO-KOLKATA. As indicated above, with shear
viscosity as the only dissipative force, boost-invariant
causal hydrodynamics require simultaneous solution of
six partial differential equations. Numerical solution of
six partial differential equations is non-trivial and it is
important to check for the numerical stability of the
solutions. Analytical solutions of viscous hydrodynam-
ics, even in restrictive conditions are not available, and
we can not check the solutions against analytical re-
sults. However, we can check for the stability of the
numerical solutions. The standard procedure of check-
ing the numerical stability is to change the integration
step lengths and look for the difference in the solution.
In Fig.16, for viscosity η/s=0.135, we have shown the
constant temperature contours in x − τ plane at a fixed
pT (GeV)
0 1 2 3 4 5
Au+Au@b=6.5 fm
τi=0.6fm,η/s=.2,Sini=110fm-3
TF=160,150,140,130 and 120 MeV
(bottom to top)
FIG. 15: Dependence of elliptic flow on the freeze-out tem-
perature. The solid lines (from bottom to top) are elliptic
flow (v2) in viscous dynamics with TF =160,150,140,130 and
120 MeV respectively.
0.250.28
0 2 4 6 8 10
____dx=dy=0.1,dτ=0.01
____dx=dy=0.2,dτ=0.02Au+Au@b=0
X (fm)
FIG. 16: (color online) constant temperature contours in
x − τ plane at a fixed y=0 fm. The black lines are obtained
with integration step lengths dx=dy=0.2 fm and dτ=0.02 fm.
The blue lines are obtained with integration step lengths,
dx=dy=0.1 fm and dτ=0.01 fm. Halving the step lengths do
not change the evolution. The numerical solutions are stable.
y=0 fm. The black and blue lines are obtained when
integration step lengths are dx=dy=0.2fm,dτ=0.02 fm,
and dx=dy=0.1 fm,dτ=0.01fm respectively. Evolution
of QGP fluid donot alter by changing the step lengths,
the solutions are stable against mesh size.
IX. SUMMARY AND CONCLUSIONS
In Israel-Stewart’s 2nd order theory of dissipative rel-
ativistic hydrodynamics, we have studied evolution QGP
fluid. In 2nd order theory, in addition to usual thermo-
dynamic quantities e.g. energy density, pressure, hydro-
dynamic velocities, dissipative flows are treated as ex-
tended thermodynamic variables. Relaxation equations
for dissipative flows are solved, simultaneously with the
energy-momentum conservation equations. This greatly
enhances the complexity of the problem. Altogether 14
partial differential equations are required to be solved.
We simplify the problem to solution of six partial differ-
ential equations by considering the evolution of baryon
free QGP fluid with longitudinal boost-invariance. We
also consider dissipation due to shear viscosity only, dis-
regarding the bulk viscosity and the heat conduction
(for a baryon free QGP fluid they do not contribute).
The six partial differential equations are solved using
the code AZHYDRO-KOLKATA, developed at the Cy-
clotron Centre, Kolkata.
To bring out the effect of viscosity, we have considered
the evolution of ideal as well as viscous QGP fluid. Both
ideal and viscous fluid are initialized similarly, at initial
time τi=0.6 fm, the central entropy density is 110 fm
Viscous dynamics require initial conditions for the shear-
stress tensor components. It is assumed that at the equi-
libration time, the shear stress tensors components have
attained their boost-invariant values.
Explicit simulation of ideal and viscous fluids confirms
that energy density of a viscous fluid, evolve slowly than
its ideal counterpart. Thus in a viscous fluid, lifetime of
the QGP phase will be enhanced. Transverse expansion
is also more in viscous dynamics. For a similar freeze-out
condition freeze-out surface is extended in viscous fluid.
As the fluid evolve, shear pressure tensors also evolve.
Explicit simulations indicate that shear pressure tensors
πxx and πyy which are initially non zero, rapidly de-
creases as the fluid evolve. By 3-4 fm of evolution they
reduced to very small values. The other independent
shear tensor πxx is zero initially. At later time it grow
in the negative direction but never grow to large value
and is always order of magnitude smaller than the stress
tensors (πxx and πyy). Spatial distribution of shear pres-
sure tensors πxx and πyy reveal an interesting feature
of viscous dynamics. Initially πxx and πyy have sym-
metric distribution. As the fluid evolve, pressure tensors
quickly become asymmetric, e.g. πxx evolve faster in the
x-direction than in the y-direction, πyy evolve faster in
y direction than in x-direction. However, in a central
collision, we did not see any effect of asymmetry in the
energy density distribution. In a central b=0 collision,
the two opposite asymmetry cancels each other.
We could not study effect of shear viscosity on par-
ticle production. However, we have explored the effect
of viscosity on parton momentum distribution and ellip-
tic flow. We have simulated b=6.5 fm Au+Au collision.
Using the Cooper-Frey prescription, transverse momen-
tum spectra as well as elliptic flow of quarks at freeze-out
temperature of TF =160 MeV are obtained. Viscous dy-
namics flattens the quark yield at large pT . At pT =3
GeV, even a small viscosity, η/s=0.8, increase the yield
by a factor of 2. The increase is even more if viscosity is
large. Viscous effect is most prominent on elliptic flow.
In ideal hydrodynamics, elliptic flow continue to increase
with pT . But in viscous dynamics v2 veer about around
pT =1.5-2 and even become negative at large pT . With
appropriate choice of viscosity, freeze-out condition, el-
liptic flow show saturation. The saturation effect is es-
sentially due to non-equilibrium correction to the equi-
librium distribution function and can not be mimicked
in an ideal hydrodynamics. Only in viscous dynamics,
saturation of elliptic flow can be explained.
APPENDIX A: COORDINATE
TRANSFORMATIONS
Instead of Cartesian coordinates xµ = (t, x, y, z) we
use curvilinear coordinates in longitudinal proper time
and rapidity, x̄m = (τ, x, y, η):
t = τ cosh η; τ =
t2 − z2 (A1)
z = τ sinh η; η =
. (A2)
The differentials
dt = dτ cosh η + dη τ sinh η, (A3)
dz = dτ sinh η + dη τ cosh η, (A4)
and the metric tensor is easily read off from
ds2 = gµνdx
µdxν = dt2 − dx2 − dy2 − dz2
= ḡmndx̄
mdx̄n = dτ2 − dx2 − y2 − τ2dη2,(A5)
namely
ḡmn =
1 0 0 0
0 −1 0 0
0 0 −1 0
0 0 0 −τ2
, ḡmn =
1 0 0 0
0 −1 0 0
0 0 −1 0
0 0 0 −1/τ2
In curvilinear coordinates we must replace the partial
derivatives with respect to xµ by covariant derivatives
(denoted by a semicolon) with respect to x̄m:
T̄ ik;p =
∂T̄ ik
+ ΓipmT̄
mk + T̄ imΓkmp.
The only non-vanishing Christoffel symbols are
Γτηη = τ ; Γ
τη = Γ
ητ = 1/τ. (A7)
The hydrodynamic 4-velocity uµ = γ(1, vx, vy, vz) is
transformed to ūm = γ(1, vx, vy, 0), with γ⊥ =
1−v2r . From here on, we drop the bars over tensor
components in x̄-coordinates for simplicity.
The projector can be easily calculated,
∆µν = gµν − uµuν
1 − γ2⊥ −γ2⊥vx γ2⊥vy 0
−γ2⊥vx −1 − γ2⊥v2x −γ2⊥vxvy 0
−γ2⊥vy −γ2⊥vxvy −1 − γ2⊥v2y 0
0 0 0 1
.(A8)
In (τ, x, y, η) coordinate system, the convective time
derivative can be obtained as,
D = u · ∂ = γ(∂τ + vx∂x + vy∂y). (A9)
For future reference,we also write down the the scalar
expansion rate
θ = ∂·u = ∂τuτ + ∂xux + ∂yuy +
(A10)
APPENDIX B: ENERGY-MOMENTUM
CONSERVATION
With longitudinal boost-invariance the energy-
momentum conservation equations Tmn;n = 0 yield
∂τ T̃
ττ + ∂x(T̃
ττvx) + ∂y(T̃
ττvy) = − (p+ τ2πηη) (B1)
∂τ T̃
τx + ∂x(T̃
τxvx) + ∂y(T̃
τxvy) = −∂x(p̃+ π̃xx − π̃τxvx) − ∂y(π̃xy − π̃τxvy) (B2)
∂τ T̃
τy + ∂x(T̃
τyvx) + ∂y(T̃
τyvy) = −∂x(π̃xy − π̃τyvx) − ∂y(p̃+ π̃yy − π̃τyvy) (B3)
where Ãmn ≡ τAmn, p̃ ≡ τp, and vx ≡ T τx/T ττ , vy ≡
T τy/T ττ .
The components of the energy momentum tensors, in-
cluding the shear pressure tensor are,
T ττ = (ε+ p)γ2⊥ − p+ πττ (B4)
T τx = (ε+ p)γ2⊥vx + π
τx (B5)
T τy = (ε+ p)γ2⊥vy + π
τy (B6)
In causal dissipative hydrodynamics, energy momen-
tum conservation equations are solved simultaneously
with the relaxation equations. Given an equation of
state, if energy density (ε) and fluid velocity (vx and
vy) distributions, at any time τi are known, Eqs.B1,B2
and B3 can be integrated to obtain ε, vx and vy at the
next time step τi+1. While for ideal hydrodynamics, this
procedure works perfectly, viscous hydrodynamics poses
a problem that shear stress-tensor components contains
time derivatives, ∂τγ⊥, ∂τu
x, ∂τu
x etc. Thus at time step
τi one needs the still unknown time derivatives. Numer-
ically, time derivatives at step τi could be obtained if ve-
locities at time step τi and τi+1 are known. One possible
way to circumvent the problem, is to use time derivatives
of the previous step, i.e. use velocities at time step τi−1
and τi to calculate the derivatives at time step τi [18].
The underlying assumption that fluid velocity changes
slowly with time. In 1st order theories, this problem is
circumvented by calculating the time derivatives from the
ideal equation of motion ,
Duµ =
, (B7)
Dε = −(ε+ p)∇µuµ. (B8)
With the help of these two equations all the time
derivatives can be expressed entirely in terms of spa-
tial gradients [15, 26]. 1st order theories are restricted
to contain terms at most linear in dissipative quantities.
Neglect of viscous terms can contribute only in 2nd or-
der corrections, which are neglected in 1st order theories.
While the procedure is not correct in 2nd order theory,
we still use it in the present calculations. The alternative
procedure of using the derivative of earlier time step is
not correct either.
APPENDIX C: RELAXATION EQUATIONS FOR
THE VISCOUS PRESSURE TENSOR
Being symmetric and traceless, the viscous pressure
tensor πµν has 9 independent components. The as-
sumption of boost invariance reduces this number by
3 (∇〈mu η〉 =0, m 6= η). The transversality condition
mn = 0 eliminates another three components ( uη
vanish and thus yield no constraint). Thus, with boost-
invariance the viscous pressure tensor has only three in-
dependent components. As seen in Eqs.B1,B2 and B3 in
a boost-invariant evolution only seven pressure tensors
πττ , πxx, πyy, πηη, πτx, πτy and πxy are of importance.
Only three of these seven are independent. In an ear-
lier publication [17], we have debated about the choice
of the independent components and suggested use of ei-
ther (πττ , πηη, ∆ = πxx − πyy) or (πττ ,πηη,πτx, πτy)
(which will require solution of an additional relaxation
equation) as choice of independent components. How-
ever, while computing we find that the three pressure
tensors πxx and πyy and πxy as independent components
are computationally more convenient. The choice has
the advantage that the dependent shear stress tensors
can be obtained from the 3 independent stress tensors by
multiplying them by fluid velocity, vx and vy (see Eqs.
C7-C10). In any other choice of independent components
(e.g. πττ ,πηη,∆ = πxx − πyy), the evaluation of depen-
dent stress tensors requires division by fluid velocities.
Since initially, fluid velocities are assumed to be zero and
they grow slowly, these choices will involve division by
very small numbers. Unless proper care is not taken, di-
vision by small numbers can lead to unrealistically large
values for the dependent stress tensors and ruin the com-
putation.
The relaxation equations for the independent shear
stress tensors πxx, πyy and πxy, in (τ ,x,y,η) co-ordinate
can be written as,
xx + vx∂xπ
xx + vy∂yπ
xx = − 1
(πxx − 2ησxx) (C1)
yy + vx∂xπ
yy + vy∂yπ
yy = − 1
(πyy − 2ησyy) (C2)
xy + vx∂xπ
xy + vy∂yπ
xy = − 1
(πxy − 2ησxy) (C3)
where τπ is the relaxation time, τπ = 2ηβ2 (see
Eq.2.21). In ultra-relativistic limit, for a Boltzman gas,
β2 can be evaluated, β2 ≈ 34p where p is the pressure
[13]. In the present paper, we use this limit to obtain the
relaxation time τπ.
The viscous pressure tensor relaxes on a time scale τπ
to 2η times the shear tensor σµν = ∇〈µ u ν〉. The xx,
yy and xy components of the shear tensor σµν can be
written as
σxx = −∂xux − uxDux −
∆xxθ (C4)
σyy = −∂yuy − uyDuy −
∆yyθ (C5)
σxy = −1
y − ∂yux − uxDuy − uyDux]
∆xyθ (C6)
The dependent shear stress tensors can easily be ob-
tained from the independent ones as,
πτx = vxπ
xx + vyπ
xy (C7)
πτy = vxπ
xy + vyπ
yy (C8)
πττ = v2xπ
xx + v2yπ
yy + 2vxvyπ
xy (C9)
τ2πηη = −(1 − v2x)πxx − (1 − v2y)πyy
+2vxvyπ
xy (C10)
The expressions for the convective time derivative D
and expansion scalar θ = ∂u̇, in (τ ,x,y,η) are given in
Eqs. A9 and A10.
APPENDIX D: PARTICLE SPECTRA
With the non-equilibrium distribution function thus
specified, it can be used to calculate the particle spectra
from the freeze-out surface. In the standard Cooper-Frye
prescription, particle distribution is obtained as,
dyd2pT
µf(x, p) (D1)
In (τ, x, y, ηs) coordinate, the freeze-out surface is pa-
rameterised as,
Σµ = (τf (x, y) cosh ηs, x, y, τf (x, y) sinh ηs), (D2)
and the normal vector on the hyper surface is,
dΣµ = (cosh ηs,−
,−∂τf
,− sinh ηs)τfdxdydηs
At the fluid position (τ, x, y, ηs) the particle 4-
momenta are parameterised as,
pµ = (mT cosh(ηs − Y ), px, py,mT sinh(ηs − Y )) (D4)
The volume element pµdΣµ become,
pµdΣµ = (mT cosh(η − Y ) − ~pT .~∇T τf )τfdxdydη (D5)
Equilibrium distribution function involve the term
which can be evaluated as,
γ(mT cosh(η − Y ) − ~vT .~pT − µ/γ)
The non-equilibrium distribution function require the
sum pµpνπµν ,
pµpνπ
µν = a1cosh
2(η − Y ) + a2cosh(η − Y ) + a3 (D7)
a1 = m
ττ + τ2πηη) (D8)
a2 = −2mT (pxπτx + pyπτy) (D9)
a3 = p
xx + p2yπ
yy + 2pxpyπ
xy −m2T τ2πηη(D10)
Inserting all the relevant formulas in Eq.D1 and inte-
grating over spatial rapidity one obtains,
dyd2pT
dyd2pT
dNneq
dyd2pT
(D11)
with,
dyd2pT
(2π)3
dxdyτf [mTK1(nβ) − pT ~∇T τfK0(nβ)] (D12)
dNneq
dyd2pT
(2π)3
dxdyτf [mT {
K3(nβ) +
K2(nβ) + (
+ a3)K1(nβ) +
K0(nβ)}
−~pT .~∇T τf{
K2(nβ) + a2K1(nβ) + (
+ a3)K0(nβ)}] (D13)
where K0, K1, K2 and K3 are the modified Bessel func-
tions.
We will also show results for elliptic flow v2. It is de-
fined as,
dyd2pT
cos(2φ)dφ
dyd2pT
(D14)
Expanding to the 1st order, elliptic flow as a function
of transverse momentum can be obtained as,
v2(pT ) = v
2 (pT )
2Nneq
pT dpT dφ
pT dpT dφ
dφcos(2φ) d
2Nneq
pT dpT dφ
pT dpT dφ
(D15)
where v
2 is the elliptic flow calculated with the equi-
librium distribution feq.
[1] BRAHMS Collaboration, I. Arsene et al., Nucl. Phys. A
757, 1 (2005).
[2] PHOBOS Collaboration, B. B. Back et al., Nucl. Phys.
A 757, 28 (2005).
[3] PHENIX Collaboration, K. Adcox et al., Nucl. Phys. A
757 (2005), in press [arXiv:nucl-ex/0410003].
[4] STAR Collaboration, J. Adams et al., Nucl. Phys. A 757
(2005), in press [arXiv:nucl-ex/0501009].
[5] Karsch F, Laermann E, Petreczky P, Stickan S and Wet-
zorke I, 2001 Proccedings of NIC Symposium (Ed. H.
Rollnik and D. Wolf, John von Neumann Institute for
Computing, Jülich, NIC Series, vol.9, ISBN 3-00-009055-
X, pp.173-82,2002.)
[6] P. F. Kolb and U. Heinz, in Quark-Gluon Plasma 3,
edited by R. C. Hwa and X.-N. Wang (World Scientific,
Singapore, 2004), p. 634.
[7] G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev.
Lett. 87, 081601 (2001) [arXiv:hep-th/0104066].
[8] G. Policastro, D. T. Son and A. O. Starinets, JHEP
0209, 043 (2002) [arXiv:hep-th/0205052].
[9] U. Heinz, J. Phys. G 31, S717 (2005).
[10] U. W. Heinz and P. F. Kolb, arXiv:hep-ph/0204061.
[11] C. Eckart, Phys. Rev. 58, 919 (1940).
[12] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, Sect.
127, Pergamon, Oxford, 1963.
[13] W. Israel, Ann. Phys. (N.Y.) 100, 310 (1976); W. Israel
and J. M. Stewart, Ann. Phys. (N.Y.) 118, 349 (1979).
[14] A. Muronga, Phys. Rev. Lett. 88, 062302 (2002) [Er-
ratum ibid. 89, 159901 (2002)]; and Phys. Rev. C 69,
034903 (2004).
[15] D. A. Teaney, J. Phys. G 30, S1247 (2004).
[16] A. Muronga and D. H. Rischke, nucl-th/0407114 (v2).
[17] U. W. Heinz, H. Song and A. K. Chaudhuri, Phys. Rev.
C 73, 034904 (2006) [arXiv:nucl-th/0510014].
[18] A. K. Chaudhuri and U. W. Heinz, J. Phys. Conf. Ser.
50, 251 (2006) [arXiv:nucl-th/0504022].
[19] A. K. Chaudhuri, Phys. Rev. C 74, 044904 (2006)
[arXiv:nucl-th/0604014].
[20] A. K. Chaudhuri, arXiv:nucl-th/0703029.
[21] A. K. Chaudhuri, arXiv:nucl-th/0703027.
[22] T. Koide, G. S. Denicol, Ph. Mota and T. Kodama, Phys.
Rev. C 75, 034909 (2007).
[23] R. Baier and P. Romatschke, arXiv:nucl-th/0610108.
[24] P. Arnold, G. D. Moore and L. G. Yaffe, JHEP 0011,
001 (2000) [arXiv:hep-ph/0010177].
[25] G. Baym, H. Monien, C. J. Pethick and D. G. Ravenhall,
Phys. Rev. Lett. 64, 1867 (1990).
[26] S. R. de Groot, W. A. van Leeuwen and Ch. G. van
Weert, Relativistic Kinetic Theory ( North-Holland, Am-
sterdam, 1980) p.36
http://arxiv.org/abs/nucl-ex/0410003
http://arxiv.org/abs/nucl-ex/0501009
http://arxiv.org/abs/hep-th/0104066
http://arxiv.org/abs/hep-th/0205052
http://arxiv.org/abs/hep-ph/0204061
http://arxiv.org/abs/nucl-th/0407114
http://arxiv.org/abs/nucl-th/0510014
http://arxiv.org/abs/nucl-th/0504022
http://arxiv.org/abs/nucl-th/0604014
http://arxiv.org/abs/nucl-th/0703029
http://arxiv.org/abs/nucl-th/0703027
http://arxiv.org/abs/nucl-th/0610108
http://arxiv.org/abs/hep-ph/0010177
