\thechapter Introduction

Causal Viscous Hydrodynamics for Relativistic Heavy Ion Collisions



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


This 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 eccentricity-scaled elliptic flow are studied in viscous hydrodynamics and compared with experimental data. The dynamical effects of using different versions of the Israel-Stewart 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

\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 Guang-You 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 thought-stimulating discussions and debates with N. Demir, S. Jeon, J. -Y. Jia, V. Koch, Y. Kovchegov, R. Lacey, R. Neufeld, J. Noronha-Hostler, 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. Guo-Mo Zeng and Prof. Yu-Xin 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.

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 de-confinement 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 ultra-relativistic 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. Only after the beginning of the RHIC program during the summer of 2000, more and more evidence showed that the QGP had been discovered [6, 7, 8, 9, 12, 13, 14]. Before discussing these QGP signatures, we will first explain the concept of the QGP and how its features have become richer during the past few years.

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]:


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 Fermi-Direc statistics. The above equations give the ideal EOS for a non-interaction 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 non-perturbative input, due to the existence of a non-perturbative 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 non-perturbative effects become dominant.

Figure 1: Energy density and three times the pressure as a function of temperature, calculated from lattice QCD with asqtad and p4 action [21].

Lattice QCD offers a non-perturbative 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: 2, which is indicated by the narrow band in Fig 1. At higher temperature , one finds that the lattice EOS approaches 85-90% of the Stefan-Boltzmann limit for an ideal non-interacting QGP gas. But this does not necessarily mean that the QGP is already weakly coupled there. From the EOS alone, lattice QCD simulations can not distinguish a weakly coupled QGP from a strongly coupled QGP. Indeed, the concept of a strongly coupled QGP, which behaves like an almost perfect liquid with very low viscosity [24, 12, 13, 25], came from the strong collective flow observed in RHIC experiments and its very successful descriptions by ideal hydrodynamics [26, 27]. However, why the QGP is strongly coupled at RHIC energies is still a theoretical challenge3, and it is also unknown whether the QGP created in the future LHC experiments will continue to be a strongly coupled liquid or behave more like a weakly coupled gas.

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 ultra-relativistic heavy ion collision () are schematically illustrated in Fig 2 and 3. The initial stage before the collision is followed after impact by a pre-equilibrium stage, an expanding QGP and hadron resonance gas (HRG) stage and a final freeze-out/decoupling stage. In more detail, these stages are characterized as follows [32]:

Figure 2: Different stages for relativistic heavy ion collisions [33].
  1. 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).

  2. Pre-equilibrium stage and thermalization: The energetic collision of the two heavy nuclei excites the QCD vacuum and produces a dense pre-equilibrium matter consisting of quarks, anti-quarks and gluons. It takes around for the pre-equilibrium bulk matter to achieve local thermalization and form the quark-gluon plasma. In this very early pre-equilibrium 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 ), pre-equilibrium real or virtual photons, and very energetic quarks and gluons with large transverse momentum (from which jets are formed after hadronization).

    Figure 3: The space-time picture of a heavy ion collision [34].
  3. QGP expansion and hadronization: After thermalization, the quark-gluon 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.

  4. Hadronic expansion and decoupling: The hadronic matter continues to expand until the system becomes very dilute. Then individual hadrons decouple from the system (kinetic freeze-out) and free-stream 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 freeze-out.

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 so-called 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 final-state interaction induced 2-particle 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 years4.

Figure 4: Theory tools for RHIC and their range of applicability [39]. Solid bands denotes the safe range to apply respective theoretical tools, while the dashed and dotted bands refer to the region where the approach is still applied, but may be questionable or unsafe. CYM: Classical Yang-Mills Theory [40, 41], LGT: Lattice Gauge Transport [42, 43], PCM: Parton Cascade Model [44], NFD: Nuclear Fluid Dynamics or Hydrodynamics [45, 26, 27]. String and Hadronic TM: String and Hadronic Transport Model [46, 47, 48].

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 pre-equilibrium 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 pre-equilibrium stage and its thermalization. For dilute systems with incoherent parton configurations, the classical motion of on-shell 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 pre-equilibrium stage, the subsequent thermalization period and the succeeding QGP expansion stage. However, the typical thermalization time obtained from PCM with 2-body () 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 pre-equilibrium stage at RHIC and LHC energies) the partons scatter so frequently that they are no longer on shell. To treat partons with off-shell energies requires the use of quantum transport theory. The theoretical framework for quark-gluon 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 Chromo-Weibel instability in the early parton plasma, which is characterized by strong momentum anisotropies, has been confirmed by numerical simulations in 1+1-dimensional [35, 36, 63] and 3+1-dimensional [64, 65] hard thermal loop effective theory and in a 3+1-dimensional expanding Glasma5  [66, 67]. However, non-Abelian saturation effects temper the exponential growth of the Chromo-Weibel instability, leading to significantly larger time scales even for isotropization than those required for the fast thermalization approximately observed at RHIC [64, 66].

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 flavor-dependent cross-sections, 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 nucleon-nucleon 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 transition6. Hydrodynamics + hadron cascade hybrids [37, 38] can more easily accommodate these, by employing a realistic EOS from lattice QCD. One generally stops hydrodynamics at a switching temperature slightly below , converting the fluid to hadrons using a Cooper-Frye prescription (see Chap. 3.4) for phase-space distributions to generate (via Monte-Carlo) initial momentum and spectra profiles for the hadron cascade simulations. However, this procedure can not deal with a potential feed-back from cascade hadrons to the hydrodynamics fluid along space-like parts of the matching freeze-out hypersurface.

1.3 QGP signatures

Possible signals and probes for the quark-gluon 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 pre-equilibrium 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 7. In the language of hydrodynamics, the collective flow is driven by pressure gradients, thus providing access to the equation of state (EOS) of the medium. Whereas the azimuthally averaged radial flow receives contributions from all expansion stages, the anisotropic “elliptic” flow seen in non-central collisions is generated mostly during the hot early stage, and thus provides the information about the QGP phase, namely its thermalization and its EOS  [26, 27, 93].

Figure 5: Elliptic flow for different hadron species [94, 90, 95], plotted as a function of transverse momentum, compared with ideal hydrodynamics predictions [96, 97].

The elliptic flow is defined as the 2nd Fourier coefficient of the azimuthal distribution of hadron spectra:


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 freeze-out 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].

Figure 6: left: The transverse kinetic energy dependence of for various hadron species. Mesons and baryons fail to different universal curves respectively. right: and are scaled by the number of valence quarks. All of the hadron species fail to the same universal curve: the differential elliptical flow per quark. [100]

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]:


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.

Figure 7: The nuclear modification factor factor as a function of transverse momentum for direct photons , as well as and mesons in central Au + Au Collisions [106]. is defined as the ratio of the cross section per nucleon-nucleon collision measured in a heavy ion collision divided by the cross section measured in collisions. if the heavy collisions can be viewed as a simple superposition of collisions.
Figure 8: Azimuthal angular correlations between high momentum hadrons in , d + Au, and central Au + Au collisions. In both cases, the trigge particle of the pair has high-momentum (). Left: Associated particles recoiling with high momenta () exhibit strong suppression in Au + Au [107]. Right: Associated particles with low recoil momenta () are strongly enhanced in Au + Au collisions [108].

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 4-5 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 back-to-back correlations between the resulting jets, as shown by two peaks in the azimuthal angle correlation separated by . Such a back-to-back 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 near-side 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 inward-traveling partner at loses most of its energy through interactions with the QGP medium and no longer contributes to energetic hadrons in the away-side (recoiling) direction [81]. The energy carried by the away-side fast parton is deposited in the medium, which leads to an enhancement of soft hadron production in the away-side hemisphere, as shown in the right panel of Fig 8 [108]. Soon, people realized that the much broadened away-side 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 super-string 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 QCD-like 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 collisions8. The transport coefficients (such as shear viscosity and bulk viscosity) are free parameters in viscous hydrodynamic calculations. The hope is that by tuning these parameters for best fits of a sufficiently large set of sensitive experimental observables, the QGP viscosity can be extracted phenomenologically from experimental data. This requires not only the development of a practical and accurate viscous hydrodynamics code, but a careful investigation of other ingredients (EOS, initial and final conditions, etc.), which may affect the viscosity-sensitive observables. Some of these issues will be addressed in this thesis. (Of course, the QGP transport coefficients are not just simple numbers, but depend on the thermodynamics properties of the fireball, which evolve with time). We should therefore use, as much theoretical knowledge as available, to constrain the the temperature dependence of these variables, at least at the qualitative level. In the rest of this chapter we will briefly review the present status of theoretical understanding and knowledge of the transport coefficients for QCD matter and other strongly coupled systems.

2 Transport coefficients for non-relativistic 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 ,

Figure 9: left: as a function of temperature for water at different pressures  [145]. right: for helium, nitrogen and water  [146].

Microscopically, shear viscosity is associated with momentum transfer between particles in different fluid cells. For a dilute gas of non-relativistic particles, the shear viscosity and shear viscosity to entropy density ratio are estimated as [147]:


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 non-relativistic dilute gas increases with temperature (taking ) and is approximately independent of density9. The density dependence of translates into a pressure dependence of , as shown in the left panel of Fig 9. Using the ideal gas EOS , one finds that, for fixed pressure, increases with temperature approximately as


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


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 liquid-gas 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 non-relativistic fluid dynamics, the bulk viscosity is defined as a combination of shear viscosity and volume viscosity10 : . Bulk viscosity vanishes for a scale invariant system. For a dilute monatomic gas, experimental data and kinetic theory show that bulk viscosity is vanishingly small or zero  [149]. The bulk viscosity of a dilute diatomic gas, however, has an appreciable value due to the exchange of energy between translational and rotational degrees of freedom during collisions  [150]. In a dense gas, bulk viscosity is associated with the internal friction force arising from the change of volume at constant shape [151]. In a liquid, bulk viscosity is related to the rearrangement of molecules during acoustic compression and rarefaction [151].

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 molecular-dynamics simulations with Lennard-Jones potentials. They found that the bulk viscosity shows a peak in the vicinity of the gas-liquid 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 non-relativistic 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 four-velocity, 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,


where the spectral functions are given by the imaginary part of the retarded green functions, , for certain components of the energy-momentum 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  11. ) The calculation of shear and bulk viscosity on the lattice [161, 162, 163, 164] starts with the calculation of the corresponding imaginary time correlators of the energy momentum tensor in Euclidean time , which is related to the spectral function by


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 model-dependent) 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]:


This result is shown as the solid line in the left panel of Fig. 10, using the two-loop 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 non-perturbatively calculated by lattice QCD using Kubo formulae. However, such calculations are highly non-trivial in the standard Monte-Carlo simulations due to the large noise to signal ratio for the relevant operators and the ill-posed 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 three-parameter 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 two-level 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:

Figure 10: left: for low temperature hadronic phase and high temperature QGP phase [145]. right: for hadron resonance gas in chemical and kinetic equilibrium [169].

Early calculations of the shear viscosity for the hadronic matter started from a theory of massless pions in the low-energy chiral limit, giving [170]


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 UrQMD12 simulations, using the Kubo formula  [169]. Fig. 10 shows the extracted as a function of temperature for a chemically equilibrated hadron resonance gas with zero chemical potential. Below , quickly rises with decreasing temperature, in qualitative agreement with the chiral pion result [170]. Between 100 MeV and 160 MeV, saturates at a value . For a hadron resonance gas out of chemical equilibrium, a similar tendency was demonstrated in Fig. 4 of Ref. [169], where it is shown that at fixed temperature , a non-zero baryon chemical potential reduces . For (an unrealistically large chemical potential for RHIC energies), can decrease to 0.3-0.4 near the phase transition, and the plateau structure between and 160 MeV is replaced by a tendency that continues to decrease with increasing .

It is worthwhile to point out that Ref. [172] argued that including Hagedorn states13 can significantly reduce in the hadronic phase near since the highly degenerate Hagedorn states could dramatically increase the entropy density. The authors of Ref. [172] investigated Hagedorn state effects on for a gas of pions and nucleons [173] and for a hadron-resonance gas with excluded volume corrections [174] and found that near , can be significantly reduced to values close to the KSS bound (see Chap. 1.5).

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 non-zero bulk viscosity. Detailed calculations from Arnold, Dogan and Moore [175] showed:


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:


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.

Figure 11: left: for high temperature QGP and for massless (left) and massive (right) pions [153].

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 Yang-Mills 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]


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 re-scaling 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 QGP14. One example is the shear viscosity in the strong coupling limit. In the language of black holes gravity, the shear viscosity, defined by the Kubo formula, is connected with the absorption cross section of the black hole: , where the last equality comes from the general theorem on black holes [187] which states that the cross section is equal to the horizon area for a broad class of black holes. On the other hand, the horizon area also represents the entropy of the black hole: . This shows that the shear viscosity to entropy ratio is a constant [188]. The constant was found by Kovtun, Son, and Starinets [146], who showed that


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

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 non-conformal 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]:


(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 non-conformal 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, anti-symmetric 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.


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 I-S 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 Navier-Stokes (N-S) equations [199, 200]. Unfortunately, it turned out that the relativistic N-S equations are unsuited for numerical implementation since they developed instabilities due to exponentially growing modes, in particular at high frequencies, and violate causality (the N-S 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 so-called 1st order Navier-Stokes 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 Navier-Stokes theory) by a kinetic equation that evolves the dissipative flows rapidly but smoothly towards their Navier-Stokes limit. The deduction of the I-S equation from the entropy current expansion and the 2nd law of thermodynamics will be introduced in Chap. 2.3. The I-S 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 I-S formalism directly solve evolution equations for the dissipative flows. In contrast, another type of 2nd-order formalism, developed by the Ottinger and Grmela (O-G 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 4-dimensional space-time coordinates [209]:


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:


where , , are the local net baryon density, energy density and pressure, respectively, and is the 4-flow 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 phase-space distribution function as follows [209, 210]:


Here are the baryon charges per particle of species (particles and anti-particles 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 near-equilibrium system, the distribution function can be decomposed into an equilibrium part plus a small deviation:


Putting eq. (21) into eq. (20) one finds


where the dissipative flows and are generated by the non-equilibrium 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]:


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 4-velocity of energy flow (), which leads to a zero value for . Then the tensor decomposition shown in eq. (23) simplifies as follows:


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 Navier-Stokes formalism

The Navier-Stokes (N-S) formalism comes from the relativistic generalization of the non-relativistic N-S equations, which impose linear relationships between the dissipative flows and the corresponding thermodynamic forces [199, 200]:


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 N-S 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 N-S 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 Israel-Stewart 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 non-relativistic 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 near-equilibrium 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:


where and are related to the local temperature and chemical potential by , . For a near-equilibrium system, one can generalize the entropy current by the following off-equilibrium expansion:


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:


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:


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 non-negative coefficients described by the N-S equations (25), the inequality is automatically satisfied:


Note that, since , is space-like, . The N-S equations thus corresponds to a 1st order expansion of the entropy current, which is why the Navier-Stokes formalism is known as 1st order viscous hydrodynamics.

The 2nd order theory of Israel and Stewart [201, 202, 195, 196] keeps , more precisely, it keeps all terms that are second order in the dissipative flows: