Quasi-Particle Theory of Shear and Bulk Viscosities of Hadronic Matter

# Quasi-Particle Theory of Shear and Bulk Viscosities of Hadronic Matter

P. Chakraborty and J. I. Kapusta
School of Physics and Astronomy, University of Minnesota
Minneapolis, Minnesota 55455, USA
###### Abstract

A theoretical framework for the calculation of shear and bulk viscosities of hadronic matter at finite temperature is presented. The framework is based on the quasi-particle picture. It allows for an arbitrary number of hadron species with point-like interactions, and allows for both elastic and inelastic collisions. Detailed balance is ensured. The particles have temperature dependent masses arising from mean field or potential effects, which maintains self-consistency between the equation of state and the transport coefficients. As an example, we calculate the shear and bulk viscosity in the linear model. The ratio of shear viscosity to entropy density shows a minimum in the vicinity of a rapid crossover transition, while the ratio of bulk viscosity to entropy density shows a maximum.

PACS numbers: 11.10.Wx, 25.75.Nq, 51.20.+d, 25.75.-q

Keywords: Shear viscosity, bulk viscosity, hadronic matter at finite temperature, quark-gluon plasma, linear sigma model.

## 1 Introduction

One of the amazing experimental discoveries of measurements of heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) is the surprising amount of collective flow exhibited by the outgoing hadrons. Collective flow is observed in both the single-particle transverse momentum distribution [1] (radial flow) and in the asymmetric azimuthal distribution around the beam axis [2] (elliptic flow). It is now generally accepted that collective flow is mostly generated early in the nucleus-nucleus collision and is present before partons fragment or coalesce into hadrons [3]. The quark-gluon matter created in these collisions must be strongly interacting, unlike the type of weakly interacting quark-gluon plasma expected to occur at very high temperatures on the basis of asymptotic freedom [4]. Perfect fluid dynamics with zero viscosity reproduces the measurements of radial and elliptic flow quite well up to transverse momenta on the order 1.5 GeV/c [5]. These results have been interpreted as strong indicators of early thermalization and collective flow on a time scale of several fm/c.

An amazing theoretical discovery was made by Kovtun, Son and Starinets [6]. They showed that certain special field theories (AdS/CFT or Anti-deSitter/Conformal Field Theory) that are dual to black branes in higher space-time dimensions [7]-[9] have the ratio of shear viscosity to entropy density (in natural units with ). The connection between transport coefficients and gravity arises because both involve commutators of the stress-energy-momentum tensor. They conjectured that all substances have this value as a lower limit, and gave as examples various atomic and molecular systems. In fact, it had been argued much earlier that any substance should a lower bound on because of the uncertainty principle [10]. Is the RHIC data telling us that the created matter has a very small viscosity, the minimal value of , that it is a perfect fluid?

The relatively good agreement between perfect fluid calculations and experimental data for hadrons of low to medium transverse momentum at RHIC suggests that the viscosity is small. However, it cannot be zero. Indeed, calculations within AdS/CFT suggest that [11]-[17]. (Whether or not this is a rigorous lower bound is still an open question [18]-[20].) There are strong theoretical arguments, and evidence from atomic and molecular systems, that should be a minimum in the vicinity of the phase transition or rapid crossover between hadronic matter and quark-gluon plasma [21, 22], and that the ratio of bulk viscosity to entropy density should be a maximum there [23]. See also [24]-[27].

It ought to be possible to extract numerical values of the viscosities in heavy ion collisions via scaling violations to perfect fluid flow predictions [28]-[31]. The program is to solve relativistic viscous fluid equations, with appropriate initial conditions and with a hadron cascade afterburner [32], over a range of beam energies and nuclei and extract and from comparison with data. Thus, sufficiently precise calculations and measurements should allow for a determination of the ratio as well as the ratio of bulk viscosity to entropy density as functions of temperature, and that these ratios can pinpoint the location of the phase transition or rapid crossover from hadronic to quark and gluon matter. This is a different method than trying to infer the equation of state of QCD in the form of pressure as a function of temperature or energy density . Because of advances in both theory and computation, vigorous activities are currently underway to determine the dissipative effects in heavy ion collisions [33]-[36].

From the theoretical perspective, it should be possible to compute the shear and bulk viscosities directly from QCD at finite temperature. In practice, this is extremely difficult because QCD is generally a strongly interacting theory. Calculations can and have been done at extremely high temperatures where perturbation theory, applied to quarks and gluons, can be used on account of asymptotic freedom; see [37] for shear viscosity and [38] for bulk viscosity. At extremely low temperatures, perturbation theory can again be used because the matter consists only of a very dilute gas of pions, and low energy pion dynamics is well understood. See [39] for massive pions; for massless pions, see [39] for shear viscosity and [40] for bulk viscosity. There have been a variety of other kinetic theory calculations of the shear and/or bulk viscosities at low to moderate temperatures in the literature in recent years [41]-[45]; these usually include only elastic scattering of one or a few species of hadrons.

In the intermediate region, which may be loosely defined as MeV, neither the low nor high temperature approach is accurate. A few lattice QCD simulations have used the Kubo formulae [46] to compute the shear [47, 48] and bulk [49] viscosities just above the critical temperature of pure gluon/glueball matter. However, accurate lattice QCD simulations of the properties of hadronic matter are extremely time consuming and the final results are still likely to be far in the future. The reason is that the lattice spacing must be small enough to describe the properties of an individual hadron ( fm) while the box size must be large enough to contain many hadrons forming the dilute gas ( fm). Hence the number of spatial lattice sites should be at least 200 in each direction. An interesting alternative approach to the intermediate region is a model of classical, non-relativistic quasi-particles with color charges [50].

Our goal in this paper is to provide a theoretical framework in which to calculate the shear and bulk viscosities of hadronic matter. This framework has the following features.

1. It is relativistic.

2. It allows for an arbitrary number of hadron species.

3. It allows for both elastic and inelastic collisions.

4. It respects detailed balance.

5. It allows for mean fields and temperature-dependent masses.

6. The viscosities and the equation of state are mutually consistent in the sense that the same interactions are used to compute them all.

Obviously some assumptions or approximations must be made for the theory to be applied in practice. The essential assumptions are that quasi-particles are well-defined and that the elementary interactions are local. Thus our proposed theoretical framework goes well beyond the classic works of [51] and [39] which, although relativistic, considered only elastic collisions in dilute gases. The inclusion of not only resonances, but especially inelastic collisions, mean fields, and temperature-dependent masses are essential for an accurate determination of the bulk viscosity [52].

The outline of this paper is as follows. In section 2 we recall the basics of the Boltzmann transport equation. In section 3 we derive the integral equations for the viscosities by using the Boltzmann equation. In section 4 we show how the Landau-Lifshitz condition plays a crucial role for the bulk viscosity. In section 5 we work out formulas for the viscosities in the relaxation time approximation. In section 6 we generalize the previous results to include mean field or potential effects and their significance for the bulk viscosity. In section 7 we apply the framework to the linear model with massive pions. As expected, the ratio has a minimum and the ratio has a maximum near the rapid crossover transition, which is more pronounced for larger vacuum masses. We conclude in section 8. The reader not interested in mathematical details is referred to sections 7 and 8 and to the appendix where the main formulas are summarized.

Since the asymmetry between matter and anti-matter in high energy nuclear collisions at RHIC is very small, so are the baryon and electric charge chemical potentials. Therefore, thermal and electrical conductivity are neglected in this paper. There inclusion is straightforward but tedious, and work to include them is in progress.

For ease and clarity of presentation we will, for the most part, display formulas that include only reactions, both elastic and inelastic, and formation and decay of resonances. At certain key points in the paper we simply write down formulas for the more general cases. In addition, in practice we generally use classical statistics. The lightest hadron for which this would be the most significant approximation is the pion, but even then the difference between Bose-Einstein and classical statistics have been shown to be inconsequential for the transport coefficients for zero chemical potential [53]. While some of the material here is not original, it is basic to developing the theoretical framework. We present it to make the paper self-contained and to define our notation.

## 2 Boltzmann Equation

The rate for processes, that is, the number of reactions per unit time per unit volume of the type is

 rate=11+δab∫d3pa2Ea(2π)3d3pb2Eb(2π)3d3pc2Ec(2π)3d3pd2Ed(2π)3|M(a,b|c,d)|2
 ×(2π)4δ4(pa+pb−pc−pd)fafb(1+(−1)2scfc)(1+(−1)2sdfd) (1)
 =11+δab∫d3pa(2π)3d3pb(2π)3d3pc(2π)3d3pd(2π)3W(a,b|c,d)
 ×fafb(1+(−1)2scfc)(1+(−1)2sdfd) (2)

whence

 W(a,b|c,d)=(2π)4δ4(pa+pb−pc−pd)2Ea2Eb2Ec2Ed|M(a,b|c,d)|2. (3)

Here is the spin of particle , etc. and the factor takes into account the possibility that the incoming particles are identical. The amplitude is dimensionless. There are either Bose-enhancment or Pauli-suppression factors in the final state. In the rest frame of the system the single-particle distributions are normalized such that

 ∫d3p(2π)3fa(x,p,t)=na(x,t) (4)

is the spatial density of particles of type . In thermal equilibrium

 fa(x,p,t)=1e(Ea−μa)/T−(−1)2sa (5)

where is the temperature and is the chemical potential of the particle. It is obvious that the rate is a Lorentz scalar.

The rate for decay processes, that is, the number of decays per unit time per unit volume of the type is

 rate=∫d3pa2Ea(2π)3d3pc2Ec(2π)3d3pd2Ed(2π)3|M(a|c,d)|2
 ×(2π)4δ4(pa−pc−pd)fa(1+(−1)2scfc)(1+(−1)2sdfd)
 =∫d3pa(2π)3d3pc(2π)3d3pd(2π)3W(a|c,d)fa(1+(−1)2scfc)(1+(−1)2sdfd) (6)

whence

 W(a|c,d)=(2π)4δ4(pa−pc−pd)2Ea2Ec2Ed|M(a|c,d)|2. (7)

This has dimension of energy.

How do the relate to cross-sections and decay rates? The relationships for the cross-sections are as follows.

 dσdΩ∗ = 164π2sp∗finalp∗initial|M|2 dσdt = 164πs1(p∗initial)2|M|2 (8)

Here and are the Mandelstam variables (one can tell from the context whether represents a Mandelstam variable or time). The differential cross-section in the center-of-momentum frame is while is usually written as a function of the invariants . Thus

 W(a,b|c,d)=sEaEbEcEdp∗initialp∗finaldσdΩ∗(2π)6δ4(pa+pb−pc−pd). (9)

This agrees with the literature on relativistic Boltzmann equations when only elastic collisions are considered since then . An alternative form, which may be more useful for cross-sections that are not isotropic in the center-of-momentum frame, is

 W(a,b|c,d)=2s(p∗initial)2EaEbEcEddσdt(2π)5δ4(pa+pb−pc−pd) (10)

where

 4s(p∗initial)2 = (s−m2a−m2b)2−4m2am2b, 4s(p∗final)2 = (s−m2c−m2d)2−4m2cm2d. (11)

For the decay consider the particle at rest. It will decay according to the usual exponential law.

 dna(t)dt=−Γa→c+dna(t). (12)

One computes that

 Γa→c+d=p∗final8πm2a|M(a|c,d)|2 (13)

where

 4m2a(p∗final)2=(m2a−m2c−m2d)2−4m2cm2d (14)

so that

 W(a|c,d)=πm2aEaEcEdp∗finalΓa→c+d(2π)4δ4(pa−pc−pd). (15)

Now we consider the Boltzmann equation. Taking into account both gain and loss rates we write it as follows.

 ∂fa∂t+va⋅∇fa=∑bcd∫d3pb(2π)3d3pc(2π)3d3pd(2π)3
 ×{11+δcdW(c,d|a,b)fcfd(1+(−1)2safa)(1+(−1)2sbfb)
 −11+δabW(a,b|c,d)fafb(1+(−1)2scfc)(1+(−1)2sdfd)}
 −W(a|c,d)fa(1+(−1)2scfc)(1+(−1)2sdfd)}
 −11+δabW(a,b|c)fafb(1+(−1)2scfc)} (16)

Due to detailed balance on the microscopic level, energy conservation, and chemical equilibrium as represented by for 2- body reactions and by for 2-body decays, we find relations between the forward and backward going rates.

 (1+δab)W(c,d|a,b) = (1+δcd)W(a,b|c,d) (17) W(c,d|a) = (1+δcd)W(a|c,d) (18)

The Boltzmann equation then becomes

 ∂fa∂t+va⋅∇fa=∑bcd11+δab∫d3pb(2π)3d3pc(2π)3d3pd(2π)3W(a,b|c,d)
 ×{fcfd(1+(−1)2safa)(1+(−1)2sbfb)
 −fafb(1+(−1)2scfc)(1+(−1)2sdfd)}
 +∑cd∫d3pc(2π)3d3pd(2π)3W(a|c,d){fcfd(1+(−1)2safa)
 −fa(1+(−1)2scfc)(1+(−1)2sdfd)}
 +∑bc∫d3pb(2π)3d3pc(2π)3W(c|a,b){fc(1+(−1)2safa)(1+(−1)2sbfb)
 −fafb(1+(−1)2scfc)}. (19)

For classical statistics, where the Bose and Pauli factors are dropped, the equilibrium phase space distribution is

 feqa(x,p,t)=e−(Ea−μa)/T. (20)

The Boltzmann equation then shortens somewhat.

 ∂fa∂t+va⋅∇fa=∑bcd11+δab∫d3pb(2π)3d3pc(2π)3d3pd(2π)3W(a,b|c,d){fcfd−fafb}
 +∑cd∫d3pc(2π)3d3pd(2π)3W(a|c,d){fcfd−fa}
 +∑bc∫d3pb(2π)3d3pc(2π)3W(c|a,b){fc−fafb} (21)

The generalization to arbitrary reactions with particles in the initial state and particles in the final state is now clear. In obvious notation the Boltzmann equation is

 ∂fa∂t+va⋅∇fa=∑{i}{j}1S∫′dPidPjW({i}|{j})F[f], (22)

where the prime indicates that there is no integration over the momentum of . There is a statistical factor for identical particles in the initial state

 S=∏ini! (23)

and products of Bose-Einstein and Fermi-Dirac distributions as appropriate

 F[f]=∏i∏j{fj(1+(−1)sifi)−fi(1+(−1)sjfj)}. (24)

## 3 Viscosities

Now we use the Boltzmann equation to calculate the viscosities. We restrict ourselves to zero chemical potentials. We assume that the system is in approximately local equilibrium, with local temperature and flow velocity . In the Landau-Lifshitz approach, is the velocity of energy transport while in the Eckart approach, would be the velocity of baryon number flow [54, 55]. However, the net baryon number, electric charge, and all other conserved quantum numbers are taken to be zero. Therefore one cannot use the Eckart approach. Another consequence of all conserved quantum numbers being zero is that thermal conductivity has no meaning.

The symmetric energy-momentum tensor is written as

 Tμν=−Pgμν+wUμUν+ΔTμν (25)

where is pressure, is entropy density, is energy density, and is enthalpy density. These are all measured in a frame in which the fluid is instantaneously at rest. The is the dissipative part. It satisfies the condition

 UμΔTμν=0 (26)

on account of the Landau-Lifshitz definition of flow. The entropy current is

 sμ=sUμ (27)

and is conserved if dissipative terms are neglected. The most general form of is given by

 ΔTμν=η(DμUν+DνUμ+23Δμν∂ρUρ)−ζΔμν∂ρUρ. (28)

Here

 Δμν=UμUν−gμν (29)

is a projection tensor normal to , and

 Dμ=∂μ−UμUβ∂β (30)

is a derivative normal to . The is the shear viscosity and the is the bulk viscosity. In the local rest frame of the fluid

 Δ0ν = 0 Δij = δij (31)

and

 D0 = 0 Di = ∂i. (32)

In this frame

 ∂μsμ=η2T(∂iUj+∂jUi−23δij∇⋅U)2+ζT(∇⋅U)2. (33)

Non-decrease of entropy requires that both viscosities be non-negative.

We are assuming that the interactions are well localized in space and time. We are assuming that they are, for practical purposes, point or contact interactions. Then the energy-momentum tensor is written as a sum of independent contributions.

 Tμν(x)=∑a∫d3p(2π)3pμapνaEafa(x,p) (34)

Allow the system to be slightly out of equilibrium. This means that is not constant in space and time, but that departures from local equilibrium are small. Then we can write

 fa(x,p)=feqa(Uαpα/T)[1+ϕa(x,p)] (35)

and so

 ΔTμν=∑a∫d3p(2π)3pμapνaEafeqa(Uαpα/T)ϕa(x,p) (36)

where . There is a constraint on in that the Landau-Lifshitz condition (26) must be satisfied. It is customary and natural to use the same tensorial decomposition of (28) when expressing as a function of space-time and momentum.

 ϕa=−Aa∂ρUρ+Caμν(DμUν+DνUμ+23Δμν∂ρUρ) (37)

Here in general will depend on the scalar . The tensor could in principle be a linear combination of and . However, the former gives zero contribution. Therefore we write where will in general depend on the scalar .

It is worthwhile emphasizing that the expansion of and in terms of the first order derivatives of the flow velocity is only an approximation. It is referred to as the first order dissipative fluid dynamics. Inclusion of second order derivatives goes under the names of Müller and Israel and Stewart. The second order theory is under intense investigation due to its usefulness in describing high energy nuclear collisions where space-time gradients are not necessarily small. It is also worth emphasizing that these same quantities are zero for an equilibrated system in uniform flow.

It is now a straightforward matter to equate the two expressions for the dissipative part of the energy-momentum tensor. It is advantageous to work in the local rest frame of the fluid.

 ζ=13∑a∫d3p(2π)3|p|2Eafeqa(Ea/T)Aa(Ea) (38)
 η=215∑a∫d3p(2π)3|p|4Eafeqa(Ea/T)Ca(Ea) (39)

How do we determine the and ? The idea is to use the Boltzmann equation where the term is evaluated using the local equilibrium distribution . This is nonzero whenever the flow velocity is changing in space-time. It will act as a source for the collision term on the other side of the Bolztmann equation. With classical statistics the Boltzmann equation reads as follows.

 E−1apμa∂μfeqa=feqa∑bcd11+δab∫d3pb(2π)3d3pc(2π)3d3pd(2π)3feqbW(a,b|c,d)
 +∑bc∫d3pb(2π)3d3pc(2π)3feqcW(c|a,b){ϕc−ϕa−ϕb} (40)

The first task is to compute the left-hand side of the Boltzmann equation. With we have

 ∂μfeqa=−1Tfeqapν(∂μUν−1TUν∂μT). (41)

Using the conservation equations for energy and momentum, , and entropy (since the viscous terms are neglected at this order), , we may deduce that

 Δμν1T∂νT=−Uα∂αUμ (42)

which has solution

 1T∂μT=Uα∂αUμ+ξUμ∂αUα (43)

where is a function of which is undetermined by Eq. (42). It can be determined by substituting the above expression into the conservation equations. This gives , where is the heat capacity per unit volume and is the square of the sound velocity. Thus

 pμ∂μfeqa=−1Tfeqapμpν[∂μUν−Uν(Uα∂α)Uμ]
 =−12Tfeqapμpν[(DμUν+DνUμ+23Δμν∂ρUρ)−23Δμν∂ρUρ
 +2v2sUμUν∂ρUρ]. (44)

After substituting in the structure of the ’s and grouping terms we get

 Aa(∂ρUρ)−Caμν(DμUν+DνUμ+23Δμν∂ρUρ)=0 (45)

where

 Aa=13EaT[(pαaUα)2(1−3v2s)−m2a]
 +∑bc∫d3pb(2π)3d3pc(2π)3feqbW(c|a,b){Ac−Aa−Ab} (46)

and

 Cμνa=pμapνa2EaT+∑cd∫d3pc(2π)3d3pd(2π)3W(a|c,d){Ccpμcpνc+Cdpμdpνd−Capμapνa}
 +∑bc∫d3pb(2π)3d3pc(2π)3feqbW(c|a,b){Capμapνa+Cbpμbpνb−Ccpμcpνc}
 +∑bcd11+δab∫d3pb(2π)3d3pc(2π)3d3pd(2π)3feqbW(a,b|c,d)
 ×{Ccpμcpνc+Cdpμdpνd−Capμapνa−Cbpμbpνb}. (47)

The bulk viscosity is very small or zero in several limits. The first is the conformal limit, which means that the theory has no dimensional parameters, such as mass or intrinsic energy scale. Then and , and so the first term on the right hand side of Eq. (46), the source term, vanishes and so do the . The second is the nonrelativistic limit of a single species of particle. Then and (plus corrections of higher order in ). Once again the first term on the right hand side of Eq. (46) vanishes (to lowest order in ) and so the bulk viscosity should be very small. These arguments do not apply to the shear viscosity since the source term does not involve the equation of state.

## 4 Landau-Lifshitz Condition

The equation (46) does not have a unique solution as it stands. For example, consider elastic scattering for just one type of particle. Starting with one solution we can generate an infinite number of other solutions by making the shift , where and are arbitrary constants. These constants are associated with particle conservation () and energy conservation (). This has been noted in the literature before. It may be restated in more physical terms. To return a system to kinetic and chemical equilibrium after a change in volume, one might either change the number of particles while keeping the average energy per particle fixed, or one might change the average energy per particle while keeping the total number of particles fixed. Now it is apparent that this ambiguity is associated with the Landau Lifshitz condition (26), which is also sometimes called the condition of fit when solving (46).

Consider an arbitrary set of particle species and all possible reactions allowed by the symmetries. Make the shift . The constant must be the same for all species of particle. The constants are just like chemical potentials; they satisfy the same relationships among themselves. Since we are restricting our considerations to systems with zero net quantum numbers, such as electric charge and baryon number, it is obvious that the are all zero, just as all chemical potentials are zero. The constant acts like an inverse temperature and is as yet undetermined.

Suppose that we have a particular solution to (46); does it satisfy the Landau Lifshitz condition (26)? The general solution would be . Using Eq. (36) the Landau Lifshitz condition for the term is

 ∑a∫d3p(2π)3feqa(Ea/T)Ea[Apara(Ea)−bEa]=0. (48)

Here it is useful to know the contributions to the pressure, energy density, entropy density and heat capacity from a single species of particle.

 Pa = T∫d3p(2π)3feqa(Ea/T) ϵa = ∫d3p(2π)3Eafeqa(Ea/T) sa = 13T2∫d3p(2π)3|p|2feqa(Ea/T) cVa = 1T2∫d3p(2π)3E2afeqa(Ea/T). (49)

Now the coefficient is determined in terms of integrals of the particular solutions.

 b=1T2cV∑a∫d3p(2π)3feqa(Ea/T)EaApara(Ea). (50)

If the particular solutions already happen to satisfy the Landau Lifshitz condition, then . Substitution of into Eq. (38), with as determined above, gives an expression for the bulk viscosity.

 ζ=13∑a∫d3p(2π)3Eafeqa(Ea/T)Apara(Ea)(|p|2−3v2sE2a) (51)

Notice that if the particular solutions happen to satisfy the Landau Lifshitz condition then Eq. (51) reduces to Eq. (38).

There is no ambiguity with the in the shear viscosity because of the tensorial structure of the integrand in Eq. (47). Physically the reason has to do with the fact that shear viscosity is associated with the response to changes in shape at fixed volume whereas bulk viscosity is associated with the response to changes in volume at fixed shape.

## 5 Relaxation Time Approximation

Consider the Boltzmann equation (21). Let us suppose that all species of particles for all values of momentum are in equilibrium except for species with momentum . Replace all phase space distributions with their equilibrium values except for , which we allow to be out of equilibrium by a small amount. Thus we write . This is the momentum-dependent relaxation time approximation. We approximate the Boltzmann equation by

 ∂fa(x,t,pa)∂t+va⋅∇fa(x,t,pa)=−ωa(Ea)δfa(x,t,pa) (52)

where

 ωa(Ea)=∑bcd11+δab∫d3pb(2π)3d3pc(2π)3d3pd(2π)3W(a,b|c,d)feqb
 +∑cd∫d3pc(2π)3d3pd(2π)3W(a|c,d)+∑bc∫d3pb(2π)3d3pc(2π)3W(c|a,b)feqb (53)

is the frequency of interaction. The equilibration time is defined as

 τa(E)=ω−1a(E). (54)

The deviation is related to the function defined in Eq. (35) by

 δfa(x,p)=feqa(x,p)ϕ(x,p). (55)

Therefore we can substitute Eqs. (37) and (44) into Eq. (52) to solve for the functions and , where .

 Apara(Ea)=τa(Ea)3TEa[(1−3v2s)E