Causal Viscous Hydrodynamics for Relativistic Heavy Ion Collisions
Abstract
\startonehalfspaceThe viscosity of the QGP is a presently hotly debated subject. Since its computation from first principles is difficult, it is desirable to try to extract it from experimental data. Viscous hydrodynamics provides a tool that can attack this problem and which may work in regions where ideal hydrodynamics begins to fail.
\startonehalfspaceThis thesis focuses on viscous hydrodynamics for relativistic heavy ion collisions. We first review the 2nd order viscous equations obtained from different approaches, and then report on the work of the Ohio State University group on setting up the equations for causal viscous hydrodynamics in 2+1 dimensions and solving them numerically for central and noncentral Cu+Cu and Au+Au collisions at RHIC energies and above. We discuss shear and bulk viscous effects on the hydrodynamic evolution of entropy density, temperature, collective flow, and flow anisotropies, and on the hadron multiplicity, single particle spectra and elliptic flow. Viscous entropy production and its influence on the centrality dependence of hadron multiplicities and the multiplicity scaling of eccentricityscaled elliptic flow are studied in viscous hydrodynamics and compared with experimental data. The dynamical effects of using different versions of the IsraelStewart second order formalism for causal viscous fluid dynamics are discussed, resolving some of the apparent discrepancies between early results reported by different groups. Finally, we assess the present status of constraining the shear viscosity to entropy ratio of the hot and dense matter created at RHIC.
Graduate Program in Physics \advisornameUlrich W. Heinz \memberRichard J. Furnstahl \memberMichael A. Lisa \memberJunko Shigemitsu
Acknowledgements.
\startonehalfspaceFirst I would like to express my deepest gratitude to my advisor Prof. Ulrich Heinz, for his thoughtful guidance during my graduate study in the past five years. I enjoy to work with him very much. He led me to a fantastic research area, which made me busy and fruitful; he provided me with all kinds of opportunities to broaden my horizon and allow me to grow up as a scientist. Especially, I thank him for his patience, encouragement and continuous support even during the hardest times of my research without which my early research would not have found its proper direction. I also enjoy the research atmosphere he created for us, through which I learned a lot from other group members. Here I especially thank Evan for many technical discussions on hydrodynamics, and GuangYou and Abhijit for discussions beyond hydrodynamics. \startonehalfspaceFor work covered in this thesis, I gratefully acknowledgement the fruitful and enlightening discussions with S. Bass, K. Dusling, R. Fries, P. Huovinen, P. Kolb, M. Lisa, D. Molnar, P. Petreczky, S. Pratt, P. Romatschke and D. Teaney. I especially thank K. Dusling and P. Romatschke for their close collaborations during the code verification. I am also grateful for many thoughtstimulating discussions and debates with N. Demir, S. Jeon, J. Y. Jia, V. Koch, Y. Kovchegov, R. Lacey, R. Neufeld, J. NoronhaHostler, J. Randrup, R. Snellings, P. Steinberg, A. Tang, I. Vitev, X. N. Wang, X. Zhe, and many others (This includes numerous people in the RHIC community, whom I met at conferences, workshops and seminars – I apologize for not listing all of them here individually). \startonehalfspaceI also thank my B. S. and M. S. thesis advisors Prof. GuoMo Zeng and Prof. YuXin Liu, who gave me early scientific training and encouraged me to switch to new research areas to further explore my potentials. \startonehalfspaceLast but not least, I thank my husband Yulu Che for the life I could enjoy with him together and for a warm home where I could rest with peace of mind. I thank my parents for their unconditional love and support.Contents:

\thechapter Introduction
 1 The Quark Gluon Plasma and Relativistic Heavy Ion Collisions
 2 Transport coefficients for nonrelativistic fluids
 3 Transport coefficients for relativistic fluids: weak vs. strong coupling
 4 Shear and bulk viscosity for QCD matter
 5 Shear viscosity and bulk viscosity from SYM and AdS/CFT
 6 Notation and outline of this thesis

\thechapter Relativistic Viscous Hydrodynamics: the Formalism
 7 From ideal to viscous hydrodynamics
 8 NavierStokes formalism
 9 IsraelStewart formalism from the 2nd law of thermodynamics
 10 IsraelStewart formalism from kinetic theory
 11 IsraelStewart formalism from conformal symmetry constraints
 12 ÖttingerGrmela formalism
 13 Final comments on the state of the formalism
 \thechapter Causal Viscous Hydrodynamics in 2+1 Dimensions
 \thechapter Generic Viscous Effects: Shear Viscosity
 \thechapter Multiplicity Scaling of and Shear Viscosity Effects
 \thechapter Generic Viscous Effects: Bulk Viscosity
 \thechapter Recent Developments in Viscous Hydrodynamics
 \thechapter Extracting from Experimental Data: Uncertainty Analysis
 \thechapter from Other Considerations and Extraction Methods
 \thechapter Summary and Concluding Remarks
 \thechapter
 A Coordinates and Transformations
 B Details of the viscous hydro code VISH2+1
 C Tests of the viscous hydro code VISH2+1
 D Hydrodynamics vs. blast wave model
 E Cooling rates for ideal and viscous heavyion fireballs
 \thechapter Glossary
Chapter \thechapter Introduction
1 The Quark Gluon Plasma and Relativistic Heavy Ion Collisions
Shortly after the discovery of the asymptotic freedom of QCD, people realized that common nuclear matter, confining quarks and gluons within individual protons and neutrons, could transform into a new deconfinement phase (at high energy densities) – the quark gluon plasma (QGP) [1, 2], a form of matter that once existed in the very early universe or some variant of which possibly still exists in the inner core of a neutron star. It was theoretically conjectured that such extreme conditions can also be realized on earth through colliding two heavy nuclei with ultrarelativistic energies [3], which transform a fraction of the kinetic energies of the two nuclei into heating the QCD vacuum within an extremely small volume. However, it actually took more than 25 years of efforts to reach the threshold for the QGP phase transition, from the first relativistic heavy ion program at the Bevalac at LBL (with beam energies of up to ) in the early 1980s [4], to the AGS at BNL (with center of mass energies per nucleon pair ) and the SPS at CERN (with ) in the late 1980s [5], to the RHIC at BNL (with ), starting in 2000 [6, 7, 8, 9], and to the inpending heavy ion program at the LHC at CERN (with ), which will probably start towards the end of next year [10].
At AGS and SPS energies and below, there was no unambiguous evidence
for QGP formation, although a number of signals found at the SPS
strongly suggested the formation of a “new state of
matter” [11]
1.1 From weakly coupled QGP to strongly coupled QGP
The quark gluon plasma is defined as a thermalized state of quarks and gluons without color confinement. It was originally conceived as a weakly coupled gas, motivated by the asymptotic freedom of QCD at high energies. Using the statistics of relativistic massless fermions and bosons, one obtains the equation of state (EOS) for a free massless QGP gas [15]:
(1) 
Here, is the total degeneracy of the QGP. For gluons with 8 colors and 2 polarizations, . For quarks with 3 colors, 2 spins and flavors, . The factor comes from FermiDirec statistics. The above equations give the ideal EOS for a noninteraction massless QGP, . If one assumes that the phase transition happens at , roughly seven times that of normal nuclear matter, then one finds a critical temperature of for a massless 2 flavor QGP.
In principle, the free energy and the EOS of the weakly coupled QGP can be perturbatively calculated order by order in thermal QCD. However, the QCD free energy shows bad convergence of the perturbative expansion in powers of the strong coupling constant [16, 17]. To solve this problem, one needs to reorganize the perturbation theory. This led to the recent developments of hard thermal loop (HTL) perturbation theory, the so called derivable approach, dimensionally reduced screened perturbation theory, and others (see the review article [18] for details), all of which yield much improved convergence for the QCD free energy and an improved EOS at high temperature. Still, the calculation needs nonperturbative input, due to the existence of a nonperturbative magnetic mass scale , which enters the EOS at order [19, 20]. On the other hand, even the most sophisticated weakly coupled QCD methods fail near the phase transition where nonperturbative effects become dominant.
Lattice QCD offers a nonperturbative approach to study the QCD
properties at finite temperature. Recent developments in this field
have made available quite precise lattice calculations of the EOS
with almost physical quark
masses [21, 22].
Fig 1. shows the Lattice EOS (at zero chemical
potential) from Bazavov et al., using two different improved
staggered fermion actions (the asqtad and p4 action,
respectively) [21]. They found that both
deconfinement and chiral symmetry restoration happen in a narrow
temperature region:
1.2 Stages of a heavy ion collision and theoretical tools
Stages of a heavy ion collision
In this subsection, we introduce different stages of a heavy ion collision and the main theoretical tools to describe or simulate these stages. After a discussion of these theoretical tools, we will list several experimental probes and signatures to detect the transient QGP phase.
The different stages for an ultrarelativistic heavy ion collision () are schematically illustrated in Fig 2 and 3. The initial stage before the collision is followed after impact by a preequilibrium stage, an expanding QGP and hadron resonance gas (HRG) stage and a final freezeout/decoupling stage. In more detail, these stages are characterized as follows [32]:

Initial stage at : Two Lorentz contracted heavy nuclei approach each other with more than 99.9% of the speed of light. At sufficiently high collision energy (possibly reached at RHIC), this initial stage can be described by dense gluon walls known as the Color Glass Condensate (CGC) (this is further explained in the theoretical tools below).

Preequilibrium stage and thermalization: The energetic collision of the two heavy nuclei excites the QCD vacuum and produces a dense preequilibrium matter consisting of quarks, antiquarks and gluons. It takes around for the preequilibrium bulk matter to achieve local thermalization and form the quarkgluon plasma. In this very early preequilibrium stage, the primary collisions between fast partons inside the colliding nuclei also generate “hard probes” with either large mass or large transverse momentum, such as heavy quark pairs ( and ), preequilibrium real or virtual photons, and very energetic quarks and gluons with large transverse momentum (from which jets are formed after hadronization).

QGP expansion and hadronization: After thermalization, the quarkgluon plasma (QGP), driven by thermal pressure gradients, expands and cools down very quickly. After reaching the critical temperature , it hadronizes and turns into hadronic matter, consisting of a mixture of stable and unstable hadrons and hadron resonances. Hadronization happens continuously at the edge of the QGP fireball during the whole QGP expansion period. In central Au+Au collisions at RHIC, it takes around for the QGP fireball to expand and fully convert to hadronic matter.

Hadronic expansion and decoupling: The hadronic matter continues to expand until the system becomes very dilute. Then individual hadrons decouple from the system (kinetic freezeout) and freestream to the detector. Like the QGP hadronization process, the hadronic decoupling happens continuously at the edge of the fireball, where the density is low. After complete hadronization, it takes another for the hadronic matter to completely freezeout.
Theoretical tools for different stages of a heavy
collisions
There does not exists a unique theoretical tool to describe the
whole heavy ion collision process from the very beginning till the
end. The different energy and time scales during different collision
stages imply dramatic changes of the effective physical degrees of
freedom and their interactions, which thus require different tools
for their description. A complete description requires matching
these tools to each other, generating socalled hybrid approaches.
For example, the results from the Color Glass Condensate (CGC)
theory can be used as initial conditions for dynamical evolution
models such as hydrodynamics or the parton cascade. Hydrodynamics or
parton cascade models are then connected with a hadron cascade model
for a description of the late hadronic stage. The hadron cascade
model is equipped with an afterburner to generate quantum
statistical or finalstate interaction induced 2particle
correlations. Fig 4 shows the different theoretical
approaches and their ranges of applicability. Although plotted by S.
Bass as early as 2001, it still very nicely illustrates our present
understanding, in spite of several new theoretical developments in
the past years
Color Glass Condensate: At very high collision energies,
particle production at mid rapidity probes the nuclear structure
functions in the small regime ( is the fraction between
parton momentum and beam momentum). It is well known that the gluon
distribution function increases dramatically with
decreasing of , for large enough probe resolution . When the
gluon density becomes high enough, two gluons start to recombine to
one. This leads to gluon saturation below some momentum scale
. At this momentum scale, each gluon mode has macroscopic
occupation number , which is why this state
has been called a condensate. These small gluons are generated
by radiation corrections from gluons at large , whose natural
evolution time scale is Lorentz dilated. This time dilation is
transferred to the small degrees of freedoms, making them evolve
very slowly compared to other natural time scales: their behavior is
“glassy”. Considering these factors, together with the color nature
of gluons, this initial stage has been named Color Glass Condensate
(CGC) [40, 41, 49]. The
production of the preequilibrium secondary soft gluons is related
to breaking the phase coherence in this initial CGC wave
function [50].
Parton cascade model: After the initial parton production,
one needs to describe the preequilibrium stage and its
thermalization. For dilute systems with incoherent parton
configurations, the classical motion of onshell partons can be
described by the Parton Cascade Model
(PCM) [44, 51, 52], which
solves the Boltzmann equation with a leading order PQCD collision
term. Initially, the PCM was assumed to work for the preequilibrium
stage, the subsequent thermalization period and the succeeding QGP
expansion stage. However, the typical thermalization time obtained
from PCM with 2body () PQCD scattering cross
sections is of the order of 5 fm/c [52], which is
too long for the fast thermalization ()
required by RHIC data [24, 53]. The recent
development of a PCM with () processes leads to faster thermalization and
indeed reproduces the large observed elliptic
flow [54, 55]. However, it is
applied to a dense system, which creates tension with the
assumptions under which the Boltzmann equation is valid. In a very
dense system (such as the very early preequilibrium stage at RHIC
and LHC energies) the partons scatter so frequently that they are no
longer on shell. To treat partons with offshell energies requires
the use of quantum transport theory. The theoretical framework for
quarkgluon quantum transport was developed by Heinz more than 20
years
ago [56, 57, 58, 59],
but its numerical demands are exorbitant, and it has not yet been
implemented numerically.
Thermalization: The formation of the quark gluon plasma
requires two aspects: local equilibrium (thermalization) and local
momentum isotropy. Generally, it is believed that the system
achieves momentum isotropy before completing thermalization, which
is assumed to happen at a time scale of
at RHIC energies as indicated by the validity of ideal hydrodynamic
simulations [24, 53]. Elucidating the
mechanism for fast isotropization and thermalization has been a
theoretical challenge. Although much progress has been made during
the past years, the thermalization mechanism is still not fully
understood. Recently, people realized that the Weibel instability,
which is known from traditional electromagnetic
plasmas [60], might help to understand the fast
isotropization and thermalization of the
QGP [61, 62, 35, 36].
The presence of a ChromoWeibel instability in the early parton
plasma, which is characterized by strong momentum anisotropies, has
been confirmed by numerical simulations in
1+1dimensional [35, 36, 63] and
3+1dimensional [64, 65] hard thermal loop
effective theory and in a 3+1dimensional expanding
Glasma
Hydrodynamics: If thermalization is achieved and can be
locally maintained during the subsequent expansion, the further
evolution of the QGP and hadronic matter can be described by
hydrodynamics [45]. Hydrodynamics is a
macroscopic approach which describes the system by macroscopic
variables, such as local energy density, pressure, temperature and
flow velocity. It requires knowledge of the equation of state, which
gives a relation between pressure, energy and baryon density, but no
detailed knowledge of the microscopic dynamics. The simplest version
is ideal hydrodynamics [26, 27], which
totally neglects viscous effects and assumes that local equilibrium
is always perfectly maintained during the fireball expansion.
Microscopically, this requires that the microscopic scattering time
is very much shorter than the macroscopic expansion (evolution) time
and that the mean free path is much smaller than the system size. If
this is not satisfied, viscous effects come in, and one can apply
viscous hydrodynamics as long as the deviation from local
equilibrium remains small [68, 69]. If the
system is far away from equilibrium, one has to switch to a kinetic
theory
approach, such as parton [44] or hadron [47] cascade models.
Hadron cascade model and hybrid approaches: The hadron
cascade model [46, 47, 48], which
solves the Boltzmann equation for a variety of hadron species with
flavordependent crosssections, is a successful tool to describe
the hadronic matter created at AGS and SPS energies. At these
collision energies, the hadron cascade model is initialized by a
superposition of hadrons and hadronic strings, produced in the
primary nucleonnucleon collisions. At RHIC energies and above,
hybrid approaches that combine a parton cascade model or
hydrodynamics with a hadron cascade, provide a “unified”
description of the evolution of the QGP and the succeeding hadronic
matter. However, some caution must be taken for the transition
between the models. Parton + hadron cascade hybrids must deal with
the problem of converting partons to hadrons without violating the
second law of thermodynamics (i.e. without losing entropy), and they
have difficulties to incorporate the change in the structure of the
QCD vacuum during the phase transition
1.3 QGP signatures
Possible signals and probes for the quarkgluon plasma have been
investigated for around 30 years since the birth of the field. Such
signatures include: collective
flow [72, 73, 26, 27],
strangeness enhancement [74, 75],
charmonium suppression [76, 77], thermal
photon and dilepton emission [78, 79], jet
quenching [80, 81, 82, 83],
critical fluctuations [84], and others. Some of the
predicted signature were already found in earlier heavy ion
experiments at AGS and SPS energies [11].
However, none of these signatures allow individually to prove QGP
formation, as they are contaminated by the dynamical evolution of
the fireball through various stages, usually from the very early
preequilibrium stage through (perhaps) a QGP phase to the late
hadronic stage. The combination of three observations at RHIC, that
finally convinced the community that the QGP has been successfully
created, were the measurements of strong anisotropic collective
flow, valence quark number scaling of the elliptic flow , and
jet quenching [85, 86]. This led to the
announcement in 2005 that the QGP had been discovered at
RHIC [12, 13, 85, 87].
Collective flow:
The hadron momentum spectra, their angular distribution (flow
patterns) and the particle yield ratios are primary observables for
the bulk medium created in heavy ion collisions. One of the main
discoveries of RHIC is that the medium displays strong collective
dynamics [88, 89, 90, 91],
which, for the first time in the history of particle and nuclear
physics, could be quantitatively well described by ideal
hydrodynamics
The elliptic flow is defined as the 2nd Fourier coefficient of the azimuthal distribution of hadron spectra:
(2) 
where is the angular distribution of the transverse momentum () dependent spectra, and is the rapidity of the particles.
Fig 5 compares the experimental elliptic flow data with ideal hydrodynamic predictions [94, 90, 95, 96]. For , where most (more than 98%) of the particles are produced, ideal hydrodynamics shows excellent agreement with the experimental data and correctly predicts the observed splitting for different hadron species. This strongly indicates that the bulk of the matter is strongly coupled and behaves like an almost perfect fluid [85].
For a successful description of the RHIC data, especially for the
elliptic flow, ideal hydrodynamic requires a fast thermalization of
the system, which must happen on a time scale of about [24, 53]. Around that time,
the early matter has a peak temperature of
– about twice the QGP phase transition temperature. Although this
gives indirect evidence for QGP formation, the success of ideal
hydrodynamics is primarily evidence for the formation of a
thermalized
new form of matter at [24, 12, 13].
Quark degrees of freedom and partonic collectivity:
The success of the statistical model in analyzing of the measured hadron abundances [98, 99] gives additional evidence for the QGP phase transition. Particle ratios including a large number of hadron species are well described by a thermal model with just two parameters, and , yielding a chemical freezeout temperature , which is approximately equal to the QCD phase transition temperature as determined by lattice QCD simulations. This strongly indicates that hadrons are born during the QGP to hadron phase transition, and that flavor changing interactions (from inelastic collisions) quickly cease right after the phase transition [32].
Direct evidence for quark degrees of freedom in the newly formed matter can be extracted from detailed measurements of the differential elliptic flow for a large variety of mesons and baryons. This observable, shown in Fig. 5, shows a characteristic splitting between mesons and baryons at intermediate transverse momentum . Even though this is above the range where hydrodynamics is valid, the underlying collective flow of the fireball still affects hadron production at this [32]. After replacing the transverse momentum by the transverse kinetic energy , the mass splitting at disappears. Fig. 6 left shows, as a result of this procedure, two universal scaling curves for baryons and mesons, respectively [100]. Baryons, which contain three valance quarks, show stronger at intermediate transverse momentum than mesons, containing only two valance quarks. This phenomenon can be well explained by the quark recombination model [101, 102, 103, 104], in which collectively flowing baryons and mesons are generated by the coalescence of quarks that collectively flow with the medium. According to the recombination model, the baryon and meson elliptic flow coefficients are expressed in terms of the quark elliptic flow in the following way [105]:
(3) 
As the right panel in Fig. 6 shows, a corresponding
rescaling of both and by the number of valence quark (2
for mesons and 3 for baryons) leads to a universal scaling
curve for all hadron species. In short, the universal
valence quark number scaling suggests that the collectively
flowing matter directly involves quarks, and that the quark
collective properties are
transferred to those of hadrons by quark recombination.
Jet quenching:
At the very beginning of a heavy ion collision, hard scatterings of incoming quarks and gluons every now and then create a pair of energetic fast partons with large transverse momentum. Each of the two fast partons will finally fragment into a spray of hadrons, forming what is called a jet. The rates for such hard process are small but grow rapidly with increasing collision energy. At RHIC energies, the fast partons for the first time became sufficiently abundant as useful probes for the hot medium. Their production rates are well understood, both experimentally from pp collision at the same energy or theoretically via PQCD. What makes them useful for heavy ion collisions is that, before hadronization to jets, they have to travel through the hot fireball medium formed in the collision, which modifies their initial energy and momentum.
One of the most exciting results from RHIC, right after the discovery of strong elliptic flow, was the observation of a strong suppression of high hadrons in central Au+Au collision, compared with the scaled results from the pp collisions [109, 110, 111]. The experimentally observed suppression by a factor 45 agrees qualitatively with theoretical predictions for jet quenching [80, 112, 82, 83], which argued that the QGP formation could lead to large energy loss of fast partons by collision induced gluon radiation and thus suppresses the production of high hadrons fragmented from such partons. The discovery of jet quenching is illustrated in Fig 7. The strong suppression of pions and mesons is distinctly in contrast to that of direct photons, which show no suppression since they interact only electrodynamically and thus directly escape from the medium without further interactions.
Angular correlations between a high leading (trigger) hadron
with other energetic hadrons provide additional strong support for
the picture of significant parton energy loss in the QGP
medium [81]. Energy momentum conservation requires
that fast partons are always generated in pairs, moving along
opposite directions in the pair center of mass frame. This leads to
backtoback correlations between the resulting jets, as shown by
two peaks in the azimuthal angle correlation separated by .
Such a backtoback correlation is clearly seen in p+p and d+Au
collisions (see left panel of
Fig. 8 [107, 113]). In central
Au+Au collisions, however, one sees only one such peak in the
direction of fast trigger hadron (blue stars in the left panel of
Fig. 8). This can be understood if one assumes that the
energetic parton pair is created near the surface of fireball: the
nearside outgoing parton quickly escapes from the medium and
fragments into the leading hadron and other softer hadrons
correlated in angle with the leading hadron, while its
inwardtraveling partner at loses most of its energy through
interactions with the QGP medium and no longer contributes to
energetic hadrons in the awayside (recoiling)
direction [81]. The energy carried by the awayside
fast parton is deposited in the medium, which leads to an
enhancement of soft hadron production in the awayside hemisphere,
as shown in the right panel of Fig 8 [108].
Soon, people realized that the much broadened awayside correlations
may be related to a possible collective “hydrodynamic” response of
the medium to the energy and momentum deposition from the fast
parton, called Mach
Cone [114, 115, 116].
1.4 “RHIC scientists serve up the perfect liquid ”
The 2007 NSAC Long Range Plane for Nuclear Physics States that “The experiments performed at RHIC since it began operation in 2000 have provided spectacular evidence that the QGP does exist” [85], but its properties are quite different from the earlier expectation based on PQCD, which had suggested that the QGP should behave like a dilute gas. In fact, the strong collective flow together with the very good description from ideal hydrodynamics indicates that the quark gluon plasma created at RHIC is strongly coupled and behaves like a nearly perfect liquid with very low viscosity [85].
The discovery of the strongly coupled QGP is intellectually exciting since it surprisingly connects superstring theory with relativistic heavy ion experiments. Mapping strongly coupled quantum field theories to weakly coupled gravity by the AdS/CFT correspondence [117, 118], string theorists have developed tools for gaining (at least) qualitative insights into the strongly coupled QCDlike systems [119, 120] where traditional PQCD methods fail to apply. Meanwhile, insights from other strongly coupled systems may help us to understand the nature of the strongly coupled QGP. An example comes from strongly interacting fermion systems, created in optical traps at extremely low temperatures, which show similar hydrodynamic behavior [121, 122]. The advantage of such cold atom experiments lies in the tunable coupling strength, via Feshbach resonances, which allows to switch between strongly and weakly coupled fermion systems and thus may help us to understand the transition between strongly and weakly coupled QGP.
With these exciting developments in the past years, the next phase of the RHIC physics program and of the incoming heavy ion program at the LHC will focus on detailed investigations of the QGP, “both to quantify its properties and to understand precisely how they emerge from the fundamental properties of QCD” [85]. Fundamental questions that need to be addressed include (see Ref [123] for details):

What is the mechanism of the unexpectedly fast thermal equilibrium?

What is the initial temperature and thermal evolution of the produced matter?

What is the energy density and equation of state of the medium?

What is the viscosity of the produced matter?

Is there direct evidence for deconfinement,color screening, and a partonic nature of the hot dense medium? What is the screening length?

Is the chiral symmetry restored by QCD?

How does the new form of matter hadronize at the phase transition?
In this thesis, I will concentrate on one of the above questions,
the viscosity of the QGP. To answer it requires continuous
progresses on both theoretical and experimental sides, especially
the development of viscous hydrodynamics and future high statistics
flow measurements for a variety of identified hadrons species. This
thesis will focus on viscous hydrodynamics for relativistic heavy
collisions
2 Transport coefficients for nonrelativistic fluids
There are several transport coefficients that characterize the
internal ”friction” in a fluid. For example, the shear viscosity
(dynamic viscosity) measures the fluid’s resistance to flow,
the bulk viscosity (volume viscosity) measures the fluid’s
resistance to expansion, and the thermal conductivity (heat
conductivity) measures the fluid’s ability to conduct
heat. In this section, we will discuss shear and bulk viscosity in
some detail, since they will be used in later parts
of this thesis.
Shear viscosity
Classically, the shear viscosity is defined as the ratio between the friction force per area and the transverse flow gradient ,
(4) 
Microscopically, shear viscosity is associated with momentum transfer between particles in different fluid cells. For a dilute gas of nonrelativistic particles, the shear viscosity and shear viscosity to entropy density ratio are estimated as [147]:
(5) 
where is the particle density, is the particle mass,
is the mean velocity,
(where is the transport cross section) is the mean free
path, and the entropy density is proportional to the particle
density . The shear viscosity of a nonrelativistic dilute
gas increases with temperature (taking
) and is approximately independent of
density
(6) 
in the high temperature gas phase.
For a liquid, the momentum transfer between different fluid cells involves the motion of voids. The density of voids decreases with temperature, which leads to an increase of shear viscosity due to the growth of mean free path. Detailed calculation [148] shows that
(7) 
where is an activation energy and is the Planck constant, and again . This result shows that for a liquid increases with decreasing temperature.
Taking the liquid phase and gas phase together, Eq.(6) and
Eq.(7) imply that is likely to reach a minimum
during the liquidgas phase transition. This is confirmed by
experimental results shown in Fig 9 which plots
in the right panel as a function of temperature for three
different fluids (helium, nitrogen and water) [146],
and for water at different pressures in the left
panel [145]. In all of these cases,
reaches a minimum near the phase transition.
Bulk viscosity
In nonrelativistic fluid dynamics, the bulk viscosity is
defined as a combination of shear viscosity and volume
viscosity
Experimentally, bulk viscosity is generally determined from measuring the sound absorption coefficient. In contrast to shear viscosity, there are no comprehensive bulk viscosity data sets for many varieties of fluids over a wide range of temperature and density, and the existing bulk viscosity data have large error bars.
The phase transition behavior of the bulk viscosity of spherical
molecules was studied by Meier, Laesecke and Kabelac using
moleculardynamics simulations with LennardJones potentials. They
found that the bulk viscosity shows a peak in the vicinity of the
gasliquid phase transition [152] (see also Fig.6 and
related comments in Ref. [153]).
In short, the shear (bulk) viscosity to entropy ratio reaches a
minimum (maximum) near phase transitions for common nonrelativistic
fluids. As we will see in Sec. 4, the same holds for
relativistic QCD matter. This is not a rigorous statement from first
principles calculations, but appears to be supported by specific
examples of both
experimental data and theoretical results.
3 Transport coefficients for relativistic fluids: weak coupling vs. strong coupling
Relativistic hydrodynamics can be viewed as an effective theory for quantum field systems at large distance and time scales. The dynamics of long wavelength, low frequency fluctuations are characterized by transport coefficients (shear viscosity , bulk viscosity and heat conductivity , electric conductivity etc.). For a relativistic fluid, the transport coefficients define the leading order corrections of the energy momentum tensor and charge current from their local equilibrium forms. For example, in the local rest frame the stress tensor reads at leading order in gradients: . In this section, we will briefly review methods for calculating these transport coefficients (especially shear and bulk viscosity) in both weak and strong coupling regimes.
For a weakly coupled system, the transport coefficients can be calculated from kinetic theory or from linear response theory. The kinetic theory approach [154, 155, 156, 157, 158, 159] starts from local equilibrium distribution functions and expands the full distribution function order by order in terms of gradients of the fourvelocity, temperature and chemical potential. The viscosity coefficients are associated with the first order terms in velocity gradients, and can be determined from the Boltzmann equations if the collision term is known explicitly (they depend on the transport cross section in the collision term).
In linear response theory, the transport coefficients can be rigorously expressed by Kubo formulae [160], that relate the transport coefficients to the slope of associated spectral functions at zero frequency. For example,
(8) 
where the spectral functions are given by the imaginary part of the retarded green functions, , for certain components of the energymomentum tensor:
The Kubo formulae can be applied to both weakly and strongly coupled
systems (see Chap. 1.5) and can be generalized for lattice
simulations. In lattice QCD, the spectral function is
obtained from the Matsubara (imaginary time) Green function ,
instead of the retarded Green function . (Both of them have
spectral representation in terms of
(9) 
In principle, one can obtain the spectral function by
inverting eq.(9), and then use
eq.(8) to calculate the transport coefficients. In
practice, one can only compute a finite number of points for
on the lattice, which makes the unique extraction of an
analytical spectral function impossible. Generally,
one makes an ansatz for with a small number of
parameters, motivated by (usually somewhat modeldependent)
considerations on the expected behavior of at small
and large , and then fits these parameters to the Monte
Carlo
data [161, 162, 163, 164].
4 Shear and bulk viscosity for QCD matter
Shear viscosity
Attempts to calculate the QGP shear viscosity using weakly coupled QCD started more than 20 years ago [165, 166]. Using kinetic theory in the relaxation time approximation and using a simple perturbative estimate for the latter, one found that the shear viscosity of the QCD matter behaves as [165, 166]. A first calculation of the leading logarithmic contribution from the Boltzmann equation was performed in [167], with a complete result finally published in [156]. A full leading order calculation that also computed the coefficient under logarithm was performed by Arnold, Moore and Yaffe using effective kinetic theory in the hard thermal loop approximation [157]. For a weakly coupled QGP with three massless quark flavors, the shear viscosity to entropy ratio is [157]:
(10) 
This result is shown as the solid line in the left panel of Fig. 10, using the twoloop renormalization group expression for the running coupling [145]. (Note that using Eq.10 with a running coupling is phenomenological and beyond the order of accuracy at which Eq.10 was derived.)
In principle, the QGP shear viscosity can also be nonperturbatively calculated by lattice QCD using Kubo formulae. However, such calculations are highly nontrivial in the standard MonteCarlo simulations due to the large noise to signal ratio for the relevant operators and the illposed inversion problem for the spectral function , given the finite number of data points for the Euclidean correlators [163]. A first attempt to calculate the shear viscosity of SU(3) gluonic matter on an lattice comes from the pioneering work of Karsch and Wyld [161], using a threeparameter ansatz for the spectral function . This method was later implemented by Nakamura and Sakai to calculate both shear and bulk viscosity on larger lattices ( and ) with improved Iwasaki action [162]. They found that the shear viscosity to entropy density ratio is less than one, , and that the bulk viscosity is zero within errors for the temperature range . A new evaluation of the shear viscosity for SU(3) gluon dynamics was performed by Meyer [163] using a twolevel algorithm [168], which dramatically improves the statistical accuracy of the relevant Euclidean correlators. He derived a robust upper bound for the shear viscosity entropy ratio and estimated that:
(11) 
Early calculations of the shear viscosity for the hadronic matter started from a theory of massless pions in the lowenergy chiral limit, giving [170]
(12) 
where is the pion decay constant. The ratio diverges at zero temperature; together with the logarithmic perturbative results discussed above, is seen to reach a minimum near the QGP phase transition. A more detailed calculation of for hadronic matter with both pions and kaons can be found in Ref. [170]; it goes beyond the chiral approximation and includes intermediate resonances such as the meson (see also [171]). The left panel of Fig. 10 shows these two results, which qualitatively agree with each other at low temperature but gradually deviate from each other above due to kaon excitations.
Recently, for a hadron resonance gas including a large set
of hadron species up to 2 GeV was extracted from
UrQMD
It is worthwhile to point out that Ref. [172]
argued that including Hagedorn states
Bulk viscosity
At sufficiently high temperature, the equation of state of a massless QGP satisfies which, on the classical level, leads to a vanishing bulk viscosity for the weakly coupled QGP. However, quantum corrections break the conformal symmetry and result in a nonzero bulk viscosity. Detailed calculations from Arnold, Dogan and Moore [175] showed:
(13) 
Here, , and depend on the number of quark flavors in QCD. For QCD with three massless quark flavors, , and . Relating Eq. (13) with the shear viscosity calculated within weakly coupled QCD [156, 157] one finds an approximate, but simple relationship between shear and bulk viscosity:
(14) 
While this relation does not hold exactly [175], it gives the right order of magnitude of .
The bulk viscosity of interacting massless pions was calculated by Chen and Wang in the framework of kinetic theory, combined with chiral perturbation theory. They found that monotonically increases with temperature below : . The bulk viscosity for massive pions was calculated by Prakash et al. in the 1990’s, using kinetic theory with experimental elastic scattering cross sections as inputs [170, 171]. In contrast to the case of massless pions, one finds that is a decreasing function of temperature. The bulk viscosity to entropy density ratio for massless and massive poins, together with the result for a weakly coupled QGP, are illustrated in Fig. 11.
The behavior of bulk viscosity near the QCD phase transition was
first investigated by Paech and Pratt within the linear sigma
model [176]. Their work showed qualitatively that
reaches a maximum near the phase transition. This behavior
was confirmed in later research by Kharzeev et
al. [177, 177] and
Meyer [164], respectively. Using low energy theorems,
Kharzeev et al. connected the bulk viscosity with the lattice
interaction measure and extracted a temperature dependent
from lattice data for pure SU(3) gluon
dynamics [177] and for full
QCD [178]. Meyer directly extracted from
lattice QCD simulations by calculating the corresponding correlation
function for a pure SU(3)gluon plasma. He showed that the peak value
of near can reach as high as
0.73 [164] (a value that is around 10 times larger
than the string theory estimates from holographical
models [179, 180]). Both of these methods
involve an ansatz for the spectral function. However, the
parameterized spectral functions used by these authors were
challenged by Moore and Saremi [181], who examined the
behavior of the spectral function with analytical methods both near
the QCD phase transition
and in the weak coupling regime.
5 Shear viscosity and bulk viscosity from SYM and the AdS/CFT correspondence
The transport coefficients calculated from weakly coupled QCD are valid at sufficiently high temperature where the running coupling constant becomes small. However, the temperature reached in relativistic heavy ion collisions is not very high (top RHIC energies: , and LHC energies: ). It is thus questionable to directly apply the weakly coupled QCD results at LHC/RHIC temperatures. Some insights for this problem can be gained from studies of supersymmetric YangMills theory (SYM), where the shear viscosity can be calculated in both the strong coupling [182, 183] and weak coupling regimes [184]. Without too much effort, one can parameterize the behavior at intermediate coupling by interpolating between the two limits [184].
For SYM theory with infinite and large, but finite t’Hooft coupling (strongly coupling regime), was calculated by Buchel, Liu and Starinets (using the gauge/gravity correspondence, see below). They found [183]
(15) 
The of SYM theory in the weakly coupled regime was investigated by the McGill group using the kinetic theory approach. They found a much smaller than what was previously found in weakly coupled QCD at the same value of the coupling constant. However, good agreement can be achieved after rescaling each result by a combination of Debye screening mass and the number of degrees of freedom in the theory. Extending this mapping between weakly coupled QCD and weakly coupled N=4 SYM to lower temperature or into the strongly coupled regime, one finds that in QCD corresponds to in N=4 SYM, where from N=4 SYM is still relatively large, of the order of 1 [184].
The finite temperature version of the gauge/gravity (AdS/CFT)
correspondence [185, 186], which maps finite
temperature gauge field theory at strong coupling onto weakly
coupled gravity in a curved space with a black whole, provides a new
method to study strongly coupled systems and helps us to gain
insights for the properties of a strongly coupled
QGP
(16) 
is a universal value for a large class of black holes or branes.
Given that the corresponding dual gauge field theories are very
different, they conjectured that is an absolute low bound
for , now known as the KSS bound
The bulk viscosity vanishes in conformal field theories such as SYM. To apply string theory techniques to the calculation of the bulk viscosity in the strongly coupled regime thus requires investigating nonconformal deformations of the original AdS/CFT correspondence. Some recent developments can be found in Ref [179, 191, 192, 193]. Based on various holographical model computations, Buchel proposed a lower bound for the bulk viscosity to entropy ratio [194]:
(17) 
(Note the different dependence on from the
weakly coupled result (14)). Using a scalar potential
tuned to reproduce the lattice equation of state, Gubser et al.
calculated the bulk viscosity via gravity duals of nonconformal
gauge theories [179, 180]. They found that
rises sharply near , with
.
6 Notation and outline of this thesis
Throughout thesis we adopt units in which . The metric tensor is always taken to be . is the projector onto the space transverse to the fluid velocity , . The partial derivative can then be decomposed as:
where, is the transverse component of the partial derivative , and is the corresponding longitudinal component. In the fluid rest frame, reduces to the time derivation and reduces to the spacial gradient. One therefore also uses the denotation .
We also frequently use the symmetric, antisymmetric and the brackets in Chap. 2, following Ref. [195, 196, 197, 198]:
Using the above notations, the commonly used local fluid rest frame variables in dissipative viscous hydrodynamics are expressed in terms of the energy momentum tensor , charge current and entropy current as following (see Chap. 2.1 and Ref. [195, 196, 197] for details):
p: thermal pressure, : bulk pressure;  
In the 1st and 2nd order formalisms for viscous hydrodynamics (Chap. 2) one also frequently encounters the following scalar and tensors constructed from the gradients of the four velocity :
Note that the above notations are for Cartesian coordinates . In Chap. 3, we change to curvilinear coordinates
for the convenience of describing a
longitudinally boost invariant system (i.e. viscous hydrodynamics in
2+1 dimensions). To express the above variables in
coordinates, one can notationally replace by , but
must also replace the partial derivative by the
covariant derivative , see Chap. 3.1 and Appendix A. 1 for detail.
Outlines:
This thesis focuses on viscous hydrodynamics for relativistic heavy ion collisions. Chapter 2 outlines the 2nd order formalism for viscous hydrodynamics, obtained from different approaches. Chapter 3 sets us up for numerical calculations in (2+1)dimension, assuming exact longitudinal boost invariance. The discussion there includes explicit transport equations in (2+1)dimension, the initial conditions, the EOS, final conditions and the transport coefficients used in viscous hydrodynamic simulations. Generic shear and bulk viscosity effects are studied in Chapter 4 and Chapter 6, respectively, by comparing runs with ideal and viscous hydrodynamics using identical initial and final conditions. Chapter 5 focuses on the system size dependence of the shear viscous effects and investigates the multiplicity scaling of the elliptic flow coefficient for both ideal and viscous hydrodynamics. In Chapter 7, we report on some of the most recent developments in viscous hydrodynamics, including the search for the optimized IS equations for numerical implementations and recent efforts on code verifications among different groups. Chapter 8 assesses the current uncertainties for extracting the QGP viscosity from experimental data and briefly comments on some future directions for viscous hydrodynamics. Chapter 9 briefly outlines other methods for estimating the QGP shear viscosity. Our conclusions are summarized in Chapter 10.
Chapter \thechapter Relativistic Viscous Hydrodynamics – the Formalism
In this chapter we introduce the
formalism for viscous (dissipative) hydrodynamics. Initial attempts
to formulate relativistic dissipative fluid dynamics started as a
relativistic generalization of the NavierStokes (NS)
equations [199, 200]. Unfortunately, it turned out that
the relativistic NS equations are unsuited for numerical
implementation since they developed instabilities due to
exponentially growing modes, in particular at high frequencies, and
violate causality (the NS formalism will be briefly described in
Chap. 2.2.). These difficulties are largely avoided in the 2nd order
formalism developed 30 years ago by Israel and
Stewart [201, 198, 202],
which goes beyond the socalled 1st order NavierStokes approach
by expanding the entropy current to 2nd order in
dissipative flows, replacing the instantaneous identification of the
dissipative flows with their driving forces multiplied by some
transport coefficient (as is done in NavierStokes theory) by a
kinetic equation that evolves the dissipative flows rapidly but
smoothly towards their NavierStokes limit. The deduction of the IS
equation from the entropy current expansion and the 2nd law of
thermodynamics will be introduced in Chap. 2.3. The IS formalism
can also be derived from the Boltzmann
equation [131, 203, 159], by expanding
the distribution function around local equilibrium. The most general
form of 2nd order viscous hydrodynamics was derived
in [204, 159] for systems with conformal
symmetry and in [203] for systems without such symmetry
(i.e. which also feature bulk viscosity and heat conductivity). Some
of this will be discussed in Chap. 2.5. All versions of the IS
formalism directly solve evolution equations for the dissipative
flows. In contrast, another type of 2ndorder formalism, developed
by the Ottinger and Grmela (OG formalism) [205, 206, 207], uses
evolution equations for auxiliary fields, rather than the
dissipative flows (see Chap. 2.6). The standard dissipative flows
can be approximately reconstructed from the auxiliary fields [208].
7 From ideal to viscous hydrodynamics
Hydrodynamics is a macroscopic tool to describe the expansion of the QGP and hadronic matter. It starts from the conservation laws for the conserved charge currents and the energy momentum tensor ( denotes the 4dimensional spacetime coordinates ) [209]:
(18a)  
(18b) 
For simplicity, one sets and only considers the conserved net baryon number current. This leaves independent variables: from the symmetric energy momentum tensor and from the charge current . However, this system of equations is unclosed since (18) only offers independent equations. To solve this problem, one needs either to reduce the number of independent variables through physically assumptions, or to provide more equations.
Ideal hydrodynamics [26, 27] solves this problem via the first route. By assuming perfect local thermal equilibrium, one can decompose the charge current and energy momentum tensor as follows:
(19a)  
(19b) 
where , , are the local net baryon density, energy density and pressure, respectively, and is the 4flow velocity which is time like and normalized: . The above ideal fluid decomposition reduces the 14 independent variables into 6 (1 each for , and and 3 independent components in ). With one additional input – the equation of state (EOS) – the system of equations is closed and can be solved to simulate the evolution of the system. Using the fundamental thermodynamic identity , it is easy to show that, in the absence of shocks (i.e. discontinuity in , or ), Eqs. (18) conserve entropy, , with .
In classical kinetic theory, and are associated with the microscopic phasespace distribution function as follows [209, 210]:
(20a)  
(20b) 
Here are the baryon charges per particle of species (particles and antiparticles count as sperate species). Plugging in the equilibrium distribution function (where and are the local chemical potential of particle species and the temperature, respectively) into eq.(20), one directly obtains the ideal fluid decomposition (19). Local thermal equilibrium is thus the basic assumption behind ideal hydrodynamics. In kinetic theory, it requires that the microscopic mean free path is much smaller than the system size and the microscopic collision time scale is much shorter than the macroscopic evolution time scale. The concept of local thermal equilibrium is, however, more general and can be formulated without recourse to kinetic theory. Ideal fluid dynamics always applies if local thermal equilibrium is ensured.
If microscopic processes are not fast enough to satisfy the above conditions, the system is no longer in local equilibrium. For a nearequilibrium system, the distribution function can be decomposed into an equilibrium part plus a small deviation:
(21) 
Putting eq. (21) into eq. (20) one finds
(22a)  
(22b) 
where the dissipative flows and are generated by the nonequilibrium contribution . By imposing the “Landau matching conditions” [200] for an arbitrary frame : and , one associates and with the equilibrium definitions of net baryon density and energy density and pressure (, ). In other words, these matching conditions fix the temperature and chemical potential of the equilibrium distribution in (21) such that and are defined in terms of in the standard way. Then , can be fully decomposed as [209]:
(23a)  
(23b) 
where , , and are called dissipative or viscous flows. More specifically, describes a baryon flow in the local rest frame, and (where is the ”heat flow vector”) describes an energy flow in the local rest frame. , and are all transverse to the frame : , and , so each of them has 3 independent components. Since is defined through and , the 3 vectors leave 6 independent components altogether. is the viscous bulk pressure, which contributes to the trace of energy momentum tensor. (where the expression is the shorthand for traceless and transverse to and as defined by the projector in square brackets) is called the shear stress tensor. is traceless and transverse to , and , and it is symmetric (like ); it thus has 5 independent components.
The Landau matching conditions, however, leave the choice of frame unconstrained. Two commonly used frames are the Eckart frame [199] and the Landau frame [200]. The Eckart frame sets parallel to the charge flow . As a consequence, disappears and the energy flow vector reduces to the heat flow vector . However, for a system with very small or vanishing net baryon number, which is approximately realized in nuclear collisions at top RHIC energy, this frame is ill defined [211]. Therefore, this thesis and other theoretical viscous hydrodynamics frameworks aimed at RHIC physics [195, 196, 197, 124, 212, 213] adopt the second choice – the Landau frame. In the Landau frame, is taken parallel to the 4velocity of energy flow (), which leads to a zero value for . Then the tensor decomposition shown in eq. (23) simplifies as follows:
(24a)  
(24b) 
For the full evolution of all components of and
, we need in addition to the 5 conservation law equations
(18) additional equations for ,
and .
8 NavierStokes formalism
The NavierStokes (NS) formalism comes from the relativistic generalization of the nonrelativistic NS equations, which impose linear relationships between the dissipative flows and the corresponding thermodynamic forces [199, 200]:
(25a)  
(25b)  
(25c) 
The (positive) “transport” coefficients , and are the bulk viscosity, heat conductivity and shear viscosity, respectively. The scalar, vector and tensor thermodynamical forces on the right hand side are explicitly given by , and .
Unfortunately the NS equations violate
causality [214, 215, 216] and
lead to infinite speed of signal propagation: if the thermodynamic
forces turn off suddenly, the dissipative flows will also disappear
instantaneously. This conflicts with the fact that a macroscopic
process caused by microscopic scattering is delayed by a relaxation
time comparable to the kinetic scattering time scale. Connected
with the acausality problem in the NS formalism are numerical
instabilities [214, 215, 216],
which render it useless for numerical simulations. Both of these two
problems are avoided in the 2nd order formalism which
is introduced in the following sections.
9 IsraelStewart formalism from the 2nd law of thermodynamics
The 2nd order formalism for relativistic dissipative fluids was first obtained by Israel and Stewart in the late 1970’s, generalizing earlier nonrelativistic work by I. [217]. They gave two derivations, a macroscopic one based on the 2nd law of thermodynamics, derived by expanding the entropy current up to the 2nd order in deviations from local equilibrium, and a microscopic one based on a nearequilibrium expansion of the Boltzmannn equation. Detailed presentations can be found in the original articles of Israel and Stewart [201, 202] and in recent articles by Muronga [195, 196, 197, 212, 213], who first brought this work to the attention of the RHIC community.
In a relativistic equilibrium system, the entropy current is written as:
(26) 
where and are related to the local temperature and chemical potential by , . For a nearequilibrium system, one can generalize the entropy current by the following offequilibrium expansion:
(27) 
where, in addition to the first order corrections and appearing in the middle terms, the last term includes second and higher order corrections in terms of and . After some algebra, the entropy production rate can be written as:
(28) 
After rewriting and in terms of the corresponding dissipative flows , and , expressing and through the thermodynamics forces and defined in eq. (25), this can be further recast into:
(29) 
The inequality on the right implements the 2nd law of thermodynamics . Note that the first three terms on the r.h.s. are first order, while the last term is higher order in the dissipative flows.
The first order theory [199, 200] neglects the second order and higher order contributions and sets . By postulating linear relationships between the dissipative flows and the thermodynamic forces with nonnegative coefficients described by the NS equations (25), the inequality is automatically satisfied:
(30) 
Note that, since , is spacelike, . The NS equations thus corresponds to a 1st order expansion of the entropy current, which is why the NavierStokes formalism is known as 1st order viscous hydrodynamics.