# Bulk viscosity, particle spectra and flow

in heavy-ion collisions

###### Abstract

We study the effects of bulk viscosity on spectra and elliptic flow in heavy ion collisions. For this purpose we compute the dissipative correction to the single particle distribution functions in leading-log QCD, and in several simplified models. We consider, in particular, the relaxation time approximation and a kinetic model for the hadron resonance gas. We implement these distribution functions in a hydrodynamic simulation of collisions at RHIC. We find significant corrections due to bulk viscosity in hadron spectra and the differential elliptic flow parameter . We observe that bulk viscosity scales as the second power of conformality breaking, , whereas scales as the first power. Corrections to the spectra are therefore dominated by viscous corrections to the distribution function, and reliable bounds on the bulk viscosity require accurate calculations of in the hadronic resonance phase. Based on viscous hydrodynamic simulations and a simple kinetic model of the resonance phase which correctly extrapolates to the kinetic description of a dilute pion gas we conclude that it is difficult to describe the spectra at RHIC unless near freeze–out. We also find that effects of the bulk viscosity on the integrated are small.

## 1 Introduction

One of the fascinating discoveries of the Relativistic Heavy Ion Collider (RHIC) program is the near ideal nature of the fluid produced in the collision of two heavy nuclei [1, 2, 3, 4, 5]. There is a general consensus in the community that the ratio of the shear viscosity to the entropy density of the system is no more than a few times the bound conjectured by Kovtun, Son, and Starinets [6]. However, it is difficult to determine the level of accuracy that can be obtained when extracting the transport properties. To date, the best estimate of the shear viscosity comes from a detailed comparison of particle spectra and elliptic flow with viscous hydrodynamic simulations [7]. But within these state of the art calculations there are many systematic uncertainties which are not fully under control. Some of these include the precise form of the initial condition, the details of the equation of state, the handling of the freeze–out dynamics, and the role of bulk viscosity. Irrespective of its role in constraining shear viscosity, the bulk viscosity of the matter produced at RHIC and the LHC is clearly an interesting quantity in itself. In this work we will study the effects of bulk viscosity on the spectra and the elliptic flow parameter. Our goal is to assess the uncertainty in the extraction of due to the bulk viscosity, and to identify observables that constrain the bulk viscosity.

The earliest viscous hydrodynamic simulations only included corrections due to shear viscosity. One could argue that this may be a safe assumption as there are a number of physical systems, possibly relevant to heavy–ion collisions, where the bulk viscosity is zero or negligible. For example, it is well known that bulk viscosity vanishes in both the non–relativistic and ultra–relativistic limits of a gas when the number of particles are conserved [10]. In a weakly coupled quark–gluon plasma, it was found that the bulk viscosity is on the order of times smaller than the shear viscosity [11]. Finally, in the simplest kinetic model, the relaxation time approximation, one finds that the bulk viscosity goes as the square of the deviation from conformality,

(1) |

The above relation was first found by Weinberg for a photon gas coupled to matter [12]. It also happens to give parametrically correct results for weakly coupled QCD but not for a scalar field theory. In the context of AdS/CFT an analogous relationship [13] has been found,

(2) |

In this case the bulk viscosity is proportional to the first power of conformal breaking. Based on these above examples, it is clear that for a system which is nearly conformally invariant (such as weakly coupled QCD) the bulk viscosity will be small. However, lattice QCD computations [14] have shown that the equation of state differs strongly from the conformal limit at temperatures relevant to heavy–ion collisions (see fig. 1). For example, if the speed of sound approaches near the phase transition we find using either of the expressions (1) or (2) given above. Even larger values have been obtained in direct lattice studies of the bulk viscosity in the regime [15]. It is therefore important to study how bulk viscosity modifies hadronic observables, such as spectra and elliptic flow. Previous studies of this type can be found in [16, 17, 18, 19, 20, 21, 22, 23].

We begin by reminding the reader how shear viscosity manifests itself in the spectra of produced particles. The equation of hydrodynamics express the conservation of the energy momentum tensor,

(3) |

which is given as a sum of ideal and dissipative parts,

(4) |

In the above expression for the stress–energy tensor we have used the definition of the three–frame projector . In the first–order (or Navier–Stokes) approximation the dissipative parts of the stress–energy tensor can be written in the local rest frame as

(5) | |||||

(6) |

where () is the shear (bulk) viscosity and indicates that the bracketed tensor should be symmetrized and made traceless. In principle, it would be satisfactory to solve the relativistic Navier–Stokes equations in order to compute the first–order viscous correction to particle spectra. However, the first order theory is plagued with difficulties such as instabilities and violations of causality. In order to circumvent these difficulties it is necessary to use a second order theory, like the one proposed by Israel and Stewart [24, 25] or Öttinger and Grmela [26, 27]. The two theories are qualitatively the same in that they both approach the first order theory for small relaxation times. In this work we will not be interested in the higher–order corrections arising from the second order theory. Instead we use second order hydrodynamics as a practical way to obtain the lowest order correction in going from ideal to Navier–Stokes hydrodynamics.

The solution to the Navier–Stokes equations will lead to viscous corrections to the resulting temperature and flow profiles. Particle spectra are then computed using the Cooper–Frye [28] formula

(7) |

where is the freeze–out hypersurface taken as a surface of constant energy density in this work. For a system out of equilibrium is not the equilibrium distribution function but also contains viscous corrections

(8) |

where is the usual equilibrium Bose/Fermi distribution function. The only constraint on is that the stress–energy tensor remains continuous across the freeze–out hypersurface;

(9) |

As shown in [29] this constraint still leaves a lot of freedom in the form of for shear viscosity. It was argued that the functional form of could fall anywhere between a linearly increasing function of momentum to a quadratically increasing function of momentum. These two forms of the distribution function lead to qualitatively different behavior for as demonstrated by the right plot of fig. 2. By definition is given by

(10) |

where is short for and is the first viscous correction to this. If, as a pedagogical exercise we neglect the viscous correction to the distribution function all together (which violates energy–momentum conservation across the freeze–out surface), would follow the curve labeled ‘’ as shown in the left plot of fig. 2. Clearly, the form of the viscous correction to the distribution function will play an important role in extracting the shear viscosity.

There is an analogous viscous correction to the distribution function coming from bulk viscosity as well. The main goal of this work is to characterize the functional form of due to bulk viscosity for various theories and models. We will also show how bulk viscous corrections exhibit themselves in spectra as well as some phenomenological consequences.

## 2 The Boltzmann transport equation

Let us first start by setting up the notation that will be used throughout this work. The equilibrium distribution functions for bosons and fermions are

(11) |

where the upper (minus) sign is for bosons and the lower (plus) sign is for fermions. We will use capital letters to label 4–vectors and bold–type for their corresponding 3–vector components having energy . The magnitude of the three–momentum will be written as . The sign convention for the metric tensor is and therefore the hydrodynamic fluid four–velocity obeys the normalization condition . We also use the notation for the quasi–particle’s energy in the laboratory frame having four momentum in the local rest frame.

The starting point for our analysis will always be the Boltzmann transport equation

(12) |

where is the particle’s velocity and is the external force on the particle,

(13) |

In this work we will consider only small deviations from local thermal equilibrium and therefore expand the Boltzmann equation around the local thermal equilibrium solution

(14) |

This procedure is known as the Chapman-Enskog expansion. In the
Chapman-Enskog procedure we expand the left hand side of the
Boltzmann equation in gradients of the thermodynamic variables
and linearize the collision operator in .
Using the following relations^{1}^{1}1Two useful identities are
and .

(15) | |||||

(16) |

the left--hand side of the Boltzmann equation can be written as^{2}^{2}2Even
though we are working in the local rest frame, gradients that are acting on
the flow velocity are still non–vanishing. For example, but since .

(17) |

Let us now assume that the quasi–particles in our system have a dispersion relation of the form

(18) |

where we have implicitly included a mass that may be a function of temperature. With this dispersion relation the following identities hold

(19) |

Making use of the above relations the left--hand side of the Boltzmann
equation can be rewritten as^{3}^{3}3In deriving this expression
we have used the two equilibrium identities and .

(20) |

where we have defined

(21) |

In order to match the kinetic description to hydrodynamics we need to define a covariantly conserved energy–momentum tensor in the kinetic theory. There is a subtlety that comes about due to the space–time dependence of the mass in the dispersion relation. In order to see this, let us first start with the canonical form of the stress–energy tensor which is typically used in kinetic theory

(22) |

For situations where the dispersion relation is independent of the medium
this form is satisfactory as one can show that energy and momentum is
covariantly conserved^{4}^{4}4This can be seen by using the definition
of the stress–energy tensor given in eq. (22) and
differentiating both sides. For the specific case where the dispersion
relation is independent of space–time we find

(25) |

In the case where we have a non–trivial dispersion relation the partial integration can not pass through the integration measure. Instead we find that

(26) |

where

(27) |

We would like to modify the stress–energy tensor such that the above source term vanishes. This can be achieved by using the definition

(28) |

Throughout this work we will always use this modified form of the stress–energy tensor when matching from the kinetic theory to the macroscopic hydrodynamic fields. We stress that if the quasi–particle’s mass is space–time independent the above two definitions of the stress–energy tensor coincide. We also note that these observations are not new. The modified form of the stress–energy tensor was used in studies of the bulk viscosity of a hadronic gas [30, 31, 32, 33] and of scalar field theory [34, 35].

## 3 Relaxation time approximation

In this section we consider the simplest form of the collision kernel, which is known as the relaxation time approximation (RTA) or Bhatnagar-Gross-Krook (BGK) approximation. In this model the collision term has the simple form

(29) |

If we define the deviation from equilibrium as and use the linearized form of the streaming operator given in eq. (20) we find that

(30) |

We would now like to identify the relaxation time encoded in with the transport coefficients and . First we start with the shear viscosity. Looking at any of the off-diagonal components of the stress–energy tensor given in eqs. (4) and (28) we find in the local rest frame

(31) |

and the shear viscosity can be identified as

(32) |

If we take a relaxation time of the form^{5}^{5}5We follow the
notation of [29] whereby taking corresponds
to the usual quadratic ansatz. In this case the relaxation time grows
linearly with momentum, , and . The other
extreme case follows from where now the relaxation is independent
of momentum, , and . For leading
order QCD one numerically finds and .

(33) |

we find the following relation between the shear viscosity and relaxation time

(34) |

where the dimensionless phase space integral is worked out in appendix B.1.

We now come to bulk viscosity, which characterizes the deviation of the
pressure from its equilibrium value as the fluid expands or contracts more
quickly than the time it takes the pressure to relax back to its equilibrium
value. The bulk viscous pressure, , is therefore related to the extra
pressure from the departure from equilibrium . However,
the departure from equilibrium can not only shift the pressure but
also the energy density by an amount . This shift in
energy density will also lead to a shift in pressure, which should
not be included in the bulk viscous pressure. This is because the
bulk viscous pressure should only include the difference between the
actual pressure and the pressure determined by thermodynamics
[11] which in our case will be . This additional pressure shift must therefore
be subtracted when defining the bulk viscous pressure^{6}^{6}6We have
used

(37) |

Making use of the form of the dispersion relation in eq. (18) it will be convenient to define the quantities and via

(38) |

The following relation between the relaxation time and bulk viscosity coefficient then holds,

(39) |

where the dimensionless phase space integral depends on both the thermal mass and the shifted mass . This phase space integral is discussed at length in appendix B.1. In the high temperature limit, , one finds

(40) |

where the function , defined in appendix B.1, depends on the statistics of the particles. For classical statistics is the usual Gamma function. From the above formulas we can recover the well–known relationship [36] between shear and bulk viscosity,

(41) |

We note that this relation is independent of the momentum dependence of the relaxation time.

### 3.1 Landau matching in the relaxation time approximation

Landau matching is a way to uniquely specify the energy density and fluid four velocity in terms of four components of . If we use the Landau–Lifshitz convention

(42) | |||||

(43) |

then the other six independent components of are given by a non-equilibrium stress tensor satisfying . In order that the stress–energy tensor remains continuous across the freeze–out surface the functional form of must be such that the Landau matching condition is satisfied; . From eq. (28) the matching condition is

(44) |

It is sufficient for the above matching condition to be satisfied in the local rest frame. This corresponds to the condition that the shift in energy density stemming from vanishes,

(45) |

Let us now look at the energy density shift coming from the off–equilibrium distribution given in eq. (30)

(46) |

The above expression simplifies considerably when there are no mean fields, ,

(47) |

The above integral vanishes only for , which is the case where
the relaxation time is momentum--independent^{7}^{7}7This
is easily seen by using the definition of the sound speed,

In the presence of mean–fields (i.e. the quasi–particle’s mass is temperature dependent) we can write eq. (46) as

(49) | |||||

In this case taking makes the first term vanish, but the second term remains finite (even though it may be parametrically small since it is proportional to the coupling). It is possible, however, to use the relaxation time approximation consistent with Landau matching by a fine–tuning of the parameter .

## 4 Scalar field theory

The case of a weakly coupled scalar field theory was studied by Jeon
[34] where the Boltzmann equation and collision kernel
were derived from first principles. While the full computation of the
transport coefficients are numerically intensive a lot can be said
about the form of the off–equilibrium distribution function from
certain general considerations. As shown in [35] one
can compute the transport coefficients in
theory at weak coupling by solving Boltzmann equation^{8}^{8}8For
our discussion it will be sufficient to look at a pure
theory.,

(50) |

where the collision operator has been split into a term containing processes and a second term involving number changing processes. While the number changing processes are higher order in the coupling constant (), they are required in order for a system undergoing a uniform expansion or contraction to equilibrate. If number changing processes were not included the above Boltzmann equation would have no solution. Formally, this is due to the presence of a (spurious) zero mode associated with particle number conservation in the processes. This zero–mode is not orthogonal to the source term and subsequently renders the linearized Boltzmann equation non–invertible. We should also point out that there is a zero mode corresponding to energy conservation. This zero–mode is not problematic since it is orthogonal to the source.

It is precisely the above behavior of a scalar field theory that allows one to obtain the approximate form of the off–equilibrium distribution function. In order to see how this works out let us start by linearizing the above Boltzmann equation around its equilibrium solution

(51) |

This equation for follows from the Chapman-Enskog expansion eq. 17. The equations in the shear and bulk channels can be separated. In the spin 0 (bulk) channel we find

(52) |

where we have written to make it explicit that the collision term should be linearized around the equilibrium solution. The resulting operators (including the final state symmetry factors) are

(53) | |||||

(54) | |||||

where we have used the shorthand . The transition rates are given as

(55) |

(56) |

Formally, we can solve eq. 52 by inverting the collision operator. Lu and Moore observed that the largest contribution will come from the near–zero mode [37] which has the form

(57) |

where are constants to be determined. Substituting the above form of into the spin 0 channel of the linearized Boltzmann equation, eq. (52), and integrating both sides over all phase space we obtain

(58) |

where

(59) |

and we have defined the function

(60) |

which characterizes the deviation of the theory from conformality. The total inelastic cross–section given in eq. (59) can be computed by doing the phase space integrals numerically. From a phenomenological perspective this is not necessary. Instead, the total inelastic cross–section can be related to the bulk viscosity coefficient by using eq. (37). This identification leads to

(61) |

The constant is undetermined by the Boltzmann equation. Instead it is constrained by requiring that the deviation from equilibrium does not bring about a shift in the energy density,

(62) |

We therefore find the following form for the off–equilibrium distribution function

(63) |

where has been defined in eq. (60) and

(64) |

For completeness, it is worth discussing the parametric behavior of the bulk viscosity at high temperature. The bulk viscosity coefficient is given by

(65) |

In the high temperature limit we can evaluate semi–analytically (see appendix B.2). In this regime we can ignore the bare and thermal mass of the scalar quasi–particles (up to logarithms). The deviation from conformality contained in is controlled by the running of the coupling. For a scalar field theory we have

(66) |

and using we find that

(67) |

Naively the total inelastic rate would go as . However, there is a soft enhancement which leads to [35]. We therefore find that

(68) |

## 5 Leading log treatment in QCD

In this section we will use the Boltzmann equation in the leading approximation. In this approximation the dynamics can be summarized by a Fokker–Plank equation which describes the momentum diffusion of the quasi–particles. The functional form of can be found by solving a simple ordinary differential equation. We start by discussing the pure glue theory and then consider a multi–component QGP.

### 5.1 Pure Glue

In a leading log approximation, is considered to be parametrically large. The resulting dynamics describes Coulomb scattering with a small momentum transfer of order but with a rapid collision rate of (up to logarithms). At leading log order the linearized Boltzmann equation can be recast as a Fokker-Planck equation [38, 39]. This equation allows us to determine in a suitable limit (absence of “gain” terms) by solving a differential equation rather than an integral equation. The Fokker-Planck equation is

(69) | |||||

where is the drag coefficient in the leading log approximation

(70) |

The Debye mass is given by . Eq. (69) without the gain terms is a Fokker–Planck equation for a hard particle undergoing drag and diffusion in a thermal bath. In order to conserve energy and momentum the gain terms must be included. The gain terms can be written as [38] with

(71) |

where and are the energy and momentum transfer to the hard particle from the thermal bath per unit time;

(72) |

where

(73) |

We express the off–equilibrium distribution function in terms of and as in equ. (51). Substituting this expression into the Fokker–Planck equation we find that the shear and bulk contributions decouple. In the shear shear channel the gain terms vanish and we are left with the following ordinary differential equation for

(74) |

At high momentum and we find [29]

(75) |

The above differential equation can also be solved numerically. For this purpose two boundary conditions must be specified. The first boundary condition is that , which implies that in QCD soft gluons equilibrate rapidly. The second boundary condition follows from the structure of the solution at large momentum. In general the differential equation has two independent solutions; one being a polynomial in and the other growing exponentially in . We choose the second boundary condition so that the exponentially growing solution is suppressed. In practice, this can be done using a shooting method on such that , which removes the exponential solution. The result of this procedure is shown in fig. 3. The shear viscosity can be found using the relation

(76) |

and we find in agreement with [39].

In the case of bulk channel, while is zero the gain term is non–vanishing. In order to understand the role of this term we first analyze the Fokker–Planck without the gain term

(77) |

This differential equation has one exact zero mode, , related to particle number conservation in scattering. This zero mode is removed if splitting and joining processes are included. We can take this into account by imposing the boundary condition . The second boundary condition is chosen in order to suppress the exponentially growing solution as discussed in the shear case.

The solution obtained in this way is not physically acceptable because it does not respect energy conservation. The fact that the collision term conserves energy implies that the most general solution of the linearized Boltzmann equation must be of the form , where is a constant and we have used the fact that the leading log collision integral is computed using . It is easy to see that this is a property of the Fokker–Planck equation in the bulk channel with the gain term included, but not without it. We find that restoring the zero mode is the dominant effect of the gain term, and that is very well approximated by the solution of the ordinary differential equation (77).

The freedom in adding the zero mode has no effect on the calculation of the bulk viscous pressure via eq. (37), because any shift in the pressure due to a shift in the energy density is projected out. However, in this work we are also interested in the correction to the single particle spectra, and in that context the linear term in matters. We therefore fix by the requirement that does not contribute to the energy density as required by the Landau matching conditions

(78) |

In this case there is no need to remove the shift in pressure due to the shift in energy density when computing the bulk viscosity

(79) |

The numerical solution of eq. (77) is shown in fig. (3). We observe that changes sign at , and that for large values of the momentum, , the non-equilibrium distribution function in the bulk channels scales as the distribution function in the shear channel multiplied by one power of the conformal symmetry breaking parameter

(80) |

Integrating the solution gives , in agreement with the result in [11]. The bulk viscosity scales as the second power of the conformal symmetry parameter,

(81) |

This result has the same structure as the relation obtained in the relaxation time approximation, eq. (41), but with a larger numerical coefficient.