# Viscous Coefficients of a Hot Pion Gas

Authors: Sourav Sarkar
Journal: Advances in High Energy Physics
 (2013)

---

## Abstract

The steps essentially involved in the evaluation of transport coefficients in linear response theory using Kubo formulas are to relate the definingretarded correlation function to the corresponding time-ordered one and to evaluate the latter in the conventional perturbation expansion. Here we evaluate the viscosities of a pion gas carrying out both the steps in the real-time formulation. We also obtain the viscous coefficients by solving the relativistic transport equation in the Chapman-Enskog approximation to leading order. An in-medium ππ cross-section is used in which spectral modifications are introduced in the propagator of the exchanged ρ.


---

## Body

## 1. Introduction

One of the most interesting results from experiments at the Relativistic Heavy Ion Collider (RHIC) is the surprisingly large magnitude of the elliptic flow of the emitted hadrons. Viscous hydrodynamic simulations of heavy ion collisions require a rather small value ofη/s, η being the coefficient of shear viscosity and s the entropy density, for the theoretical interpretation of this large collective flow. The value being close to 1/4π, the quantum lower bound for this quantity [1], matter produced in these collisions is believed to be almost a perfect fluid [2].This finding has led to widespread interest in the study of nonequilibrium dynamics, especially in the microscopic evaluation of the transport coefficients of both partonic as well as hadronic forms of strongly interacting matter. In the literature one comes across basically two approaches that have been used to determine these quantities. One is the kinetic theory method in which the nonequilibrium distribution function which appears in the transport equation is expanded in terms of the gradients of the flow velocity field. The coefficients of this expansion which are related to the transport coefficients are then perturbatively determined using various approximation methods. The other approach is based on response theory in which the nonequilibrium transport coefficients are related by Kubo formulas to equilibrium correlation functions. They are then perturbatively evaluated using the techniques of thermal field theory. Alternatively, the Kubo formulas can be directly evaluated on the lattice [3] or in transport cascade simulations [4] to obtain the transport coefficients.Thermal quantum field theory has been formulated in the imaginary as well as real-time [5–9]. For time independent quantities such as the partition function, the imaginary time formulation is well suited and stands as the only simple method of calculation. However, for time dependent quantities like two-point correlation functions, the use of this formulation requires a continuation to imaginary time and possibly back to real-time at the end. On the other hand, the real-time formulation provides a convenient framework to calculate such quantities, without requiring any such continuation at all.A difficulty with the real-time formulation is, however, that all two-point functions take the form of2×2 matrices. But this difficulty is only apparent. Such matrices are always diagonalisable, and it is the 11  component of the diagonalised matrix that play the role of the single function in the imaginary time formulation. It is only in the calculation of this 11  component to higher order in perturbation that the matrix structure appears in a nontrivial way.In the literature transport coefficients are evaluated using the imaginary time formulation [10–12]. Such a coefficient is defined by the retarded correlation function of the components of the energy-momentum tensor. As the conventional perturbation theory applies only to time-ordered correlation functions, it is first necessary to relate the two types of correlation functions using the Källen-Lehmann spectral representation [13–16]. We find this relation directly in real-time formulation. The time-ordered correlation function is then calculated also in the covariant real-time perturbative framework to finally obtain the viscosity coefficients of a pion gas.We also calculate the viscous coefficients in a kinetic theory framework by solving the transport equation in the Chapman-Enskog approximation to the leading order. This approach being computationally more efficient [12] has been mostly used in the literature to obtain the viscous coefficients. The ππ cross-section is a crucial dynamical input in these calculations. Scattering amplitudes evaluated using chiral perturbation theory [17, 18] to lowest order have been used in [19, 20], and unitarization improved estimates of the amplitudes were used in [21] to evaluate the shear viscosity. Phenomenological scattering cross-section using experimental phase shifts has been used in [20, 22–24] in view of the fact that the ππ cross-section estimated from lowest order chiral perturbation theory is known to deviate from the experimental data beyond centre of mass energy of 500 MeV primarily due to the ρ pole which dominates the cross-section in the energy region between 500–1000 MeV. All these approaches have used a vacuum cross section. To construct an in-medium cross-section we employ an effective Lagrangian approach which incorporates ρ and σ meson exchange in ππ scattering. Medium effects are then introduced in the ρ propagator through one-loop self-energy diagrams [25].In Section2 we derive the spectral representations for the retarded and time-ordered correlation functions in the real-time version of thermal field theory. We also review the formulation of the nonequilibrium density operator and obtain the expressions for the viscosities in terms of equilibrium (retarded) two-point functions. The time-ordered function is then calculated to lowest order with complete propagators in the equilibrium theory. In Section 3 we briefly recapitulate the expressions for the viscosities obtained by solving the Uehling-Uhlenbeck transport equation in the kinetic theory framework. We then evaluate the ππ cross-section in the medium briefly discussing the one-loop ρ self-energy due to πh(h=π,ω,h1,a1) loops evaluated in the real-time formulation discussed previously. We end with a summary in Section 4.


## 2. Viscous Coefficients in the Linear Response Theory

### 2.1. Real-Time Formulation

In this section we review the real-time formulation of equilibrium thermal field theory leading to the spectral representations of bosonic two-point functions [16]. This formulation begins with a comparison between the time evolution operator e-iH(t2-t1) of quantum theory and the Boltzmann weight e-βH=e-iH(τ-iβ-τ) of statistical physics, where we introduce τ as a complex variable. Thus while for the time evolution operator, the times t1 and t2(t2>t1) are any two points on the real line, the Boltzmann weight involves a path from τ to τ-iβ in the complex time plane. Setting this τ=-T, where T is real, positive, and large, we can get the contour C shown in Figure 1, lying within the region of analyticity in this plane and accommodating real-time correlation functions [6, 8].Figure 1
The contourC in the complex time plane used here for the real-time formulation.Let a general bosonic interacting field in the Heisenberg representation be denoted byΦl(x), whose subscript l collects the index (or indices) denoting the field component and derivatives acting on it. Although we will call its two-point function as propagator, Φl(x) can be an elementary field or a composite local operator. (If Φl(x) denotes the pion field, it will, of course, not have any index.)The thermal expectation value of the productΦl(x)Φl′†(x′) may be expressed as
(1)〈Φl(x)Φl′†(x′)〉=1Z∑m,n‍e-βEm×〈m|Φl(x)|n〉〈n|Φl′†(x′)|m〉,
where 〈O〉 for any operator O denotes equilibrium ensemble average:
(2)〈O〉=Tr(e-βHO)Z,Z=Tre-βH.
Note that we have two sums in (1), one to evaluate the trace and the other to separate the field operators. They run over a complete set of states, which we choose as eigenstates |m〉 of four-momentum Pμ. Using translational invariance of the field operator:
(3)Φl(x)=eiP·xΦl(0)e-iP·x,
we get
(4)〈Φl(x)Φl′†(x′)〉=1Z∑m,n‍e-βEmei(km-kn)·(x-x′)×〈m|Φl(0)|n〉〈n|Φl′†(0)|m〉.
Its spatial Fourier transform is
(5)∫‍d3xe-ik·(x-x′)〈Φl(x)Φl′†(x′)〉=(2π)3Z∑m,n‍e-βEmei(Em-En)(τ-τ′)δ3(km-kn+k)×〈m|Φl(0)|n〉〈n|Φl′†(0)|m〉,
where the times τ, τ′ are on the contour C. We now insert unity on the left of (5) in the form
(6)1=∫-∞∞‍dk0′δ(Em-En+k0′).
(We reserve k0 for the variable conjugate to the real-time.) Then it may be written as
(7)∫d3xe-ik·(x-x′)〈Φl(x)Φl′†(x′)〉=∫‍dk0′2πe-ik0′(τ-τ′)Mll′+(k0′,k),
where the spectral function M+ is given by [kμ′=(k0′,k)] as follows:
(8)Mll′+(k′)=(2π)4Z∑m,n‍e-βEmδ4(km-kn+k′)×〈m|Φl(0)|n〉〈n|Φl′†(0)|m〉.In just the same way, we can work out the Fourier transform of〈Φl′†(x′)Φl(x)〉 as follows:
(9)∫d3xe-ik·(x-x′)〈Φl′†(x′)Φl(x)〉=∫‍dk0′2πe-ik0′(τ-τ′)Mll′-(k0′,k),
with a second spectral function M- given by
(10)Mll′-(k′)=(2π)4Z∑m,n‍e-βEmδ4(kn-km+k′)×〈m|Φl′†(0)|n〉〈n|Φl(0)|m〉.
The two spectral functions are related by the KMS relation [26, 27]
(11)Mll′+(k)=eβk0Mll′-(k),
in momentum space, which may be obtained simply by interchanging the dummy indices m, n in one of Mll′±(k) and using the energy conserving δ-function.We next introduce thedifference of the two spectral functions:
(12)ρll′(k)≡Mll′+(k)-Mll′-(k)
and solve this identity and the KMS relation (11) for Mll′±(k) as follows:
(13)Mll′+(k)={1+f(k0)}ρll′(k),Mll′-(k)=f(k0)ρll′(k),
where f(k0) is the distribution-like function
(14)f(k0)=1eβk0-1,-∞<k0<∞.
In terms of the true distribution function
(15)n(|k0|)=1eβ|k0|-1,
it may be expressed as
(16)f(k0)=f(k0){θ(k0)+θ(-k0)}=nϵ(k0)-θ(-k0).With the previous ingredients, we can build the spectral representations for the two types of thermal propagators. First consider thetime-ordered one:
(17)-iDll′(x,x′)=〈TcΦl(x)Φl′†(x′)〉=θc(τ-τ′)〈Φl(x)Φl′†(x′)〉+θc(τ′-τ)〈Φl′†(x′)Φl(x)〉.
Using (7), (9), and (13), we see that its spatial Fourier transform is given by [6]
(18)Dll′(τ-τ′,k)=i∫-∞∞‍dk0′2πρll′(k0′,k)×e-ik0′(τ-τ′){θc(τ-τ′)+f(k0′)}.AsT→∞, the contour of Figure 1 simplifies, reducing essentially to two parallel lines, one the real axis and the other shifted by -iβ/2, points which will be denoted, respectively, by subscripts 1 and 2, so that τ1=t,τ2=t-iβ/2 [8]. The propagator then consists of four pieces, which may be put in the form of a 2×2 matrix. The contour ordered θ′s may now be converted to the usual time ordered ones. If τ, τ′ are both on line 1 (the real axis), the τ and t orderings coincide, θc(τ1-τ1′)=θ(t-t′). If they are on two different lines, the τ ordering is definite, θc(τ1-τ2′)=0,θc(τ2-τ1′)=1. Finally if they are both on line 2, the two orderings are opposite, θc(τ2-τ2′)=θ(t′-t).Back to real-time, we can work out the usual temporal Fourier transform of the components of the matrix to get(19)Dll′(k0,k)=∫-∞∞‍dk0′2πρll′(k0′,k)Λ(k0′,k0),
where the elements of the matrix Λ are given by [16]
(20)Λ11=-Λ22*=1k0′-k0-iη+2πif(k0′)δ(k0′-k0),Λ12=Λ21=2πieβk0′/2f(k0′)δ(k0′-k0).
Using relation (16), we may rewrite (20) in terms of n as follows:
(21)Λ11=-Λ22*=1k0′-k0-iηϵ(k0)+2πinϵ(k0)δ(k0′-k0),Λ12=Λ21=2πin(1+n)ϵ(k0)δ(k0′-k0).The matrixΛ and hence the propagator Dll′ can be diagonalised to give
(22)Dll′(k0,k)=U(D¯ll′00-D¯ll′*)U,
where D¯ll′ and U are given by
(23)D¯ll′(k0,k)=∫-∞∞‍dk0′2πρll′(k0′,k)k0′-k0-iηϵ(k0),U=(1+nnn1+n).
Equation (22) shows that D¯ can be obtained from any of the elements of the matrix D, say D11. Omitting the indices ll′, we get
(24)ReD¯=ReD11,ImD¯=tanh(β|k0|2)ImD11.Looking back at the spectral functionsMll′± defined by (8) and (10), we can express them as usual four-dimensional Fourier transforms of ensemble average of the operator products, so that ρll′ is the Fourier transform of that of the commutator:
(25)ρll′(k0,k)=∫‍d4yeik·(y-y′)〈[Φl(y),Φl′(y′)]〉,
where the time components of y and y′ are on the real axis in the τ-plane. Taking the spectral function for the free scalar field
(26)ρ0=2πϵ(k0)δ(k2-m2),
we see that D¯ becomes the free propagator, D¯(k0,k)=-1/(k2-m2).We next consider theretarded thermal propagator:
(27)Dll′R(x,x′)=iθc(τ-τ′)〈[Φl(x,τ),Φl′(x′,τ′)]〉,
where again τ, τ′ are on the contour C (Figure 1). Noting equations (7), (9), and (12) the three-dimensional Fourier transform may immediately be written as(28)Dll′R(τ-τ′,k)=iθc(τ-τ′)∫-∞∞‍dk0′2πe-ik0′(τ-τ′)ρll′(k0′,k).
As before we isolate the different components with real-times and take the Fourier transform with respect to real-time. Thus for the 11 component we simply have
(29)Dll′R(t-t′,k)11=iθ(t-t′)∫-∞∞‍dk0′2πe-ik0′(t-t′)ρll′(k0′,k),
whose temporal Fourier transform gives
(30)Dll′R(k0,k)11=∫-∞∞‍dk0′2πρll′(k0′,k)k0′-k0-iη.
This 11 component suffices for us, but we also display the complete matrix as follows:(31)Dll′R(k0,k)=(Dll′R(k)110ρll′(k){nn+1θ(k0)+n+1nθ(-k0)}-Dll′R*(k)11).Though we deal with matrices in real-time formulation, it is the 11 component that is physical. Equations (23) and (30) then show that we can continue the time-ordered two-point function into the retarded one by simply changing the iϵ prescription as follows:
(32)Dll′R(k0+iη,k)11=D¯ll′(k0+iηϵ(q0)⟶k0+iη,k).
The point to note here is that for the time-ordered propagator, it is the diagonalised matrix and not the matrix itself, whose 11 component can be continued in a simple way.


### 2.2. Transport Coefficients

We now use the linear response approach to arrive at expressions of the transport coefficients as integrals of retarded Green’s functions over space. We follow the method proposed by Zubarev [28], which is excellently reviewed in [10]. Here the system is supposed to be in the hydrodynamical stage where the mean free time of the constituent particles is much shorter than the relaxation time of the whole system under consideration. Thus local equilibrium will be attained quickly, while global equilibrium will be approaching gradually. Since the system is assumed to be not far from equilibrium, we may retain only linear terms in space-time gradients of thermodynamical parameters, like temperature and velocity fields. We assume the energy-momentum of the system to be conserved as follows:
(33)∂μTμν(x,t)=0.
The nonequilibrium density matrix operator is constructed in the Heisenberg picture, where it is independent of time as follows:
(34)dρdt=0.
Following Zubarev, we construct the operator B(x,t) as follows:
(35)B(x,t)=ϵ∫-∞t‍dt1eϵ(t1-t)Fν(x,t1)T0ν(x,t1),(ϵ⟶0+),
where Fν(x,t)=β(x,t)uν(x,t). Here β(x,t) is a Lorentz invariant quantity defining the local temperature, and uν(x,t) is the four-velocity field of the fluid
(36)uν(x,t)uν(x,t)=1.
The construction (35), which smooths out the oscillating terms, resembles the one used in the formal theory of scattering [28, 29] and selects out the retarded solution.The expression (35) is actually independent of t; the time derivative is
(37)ddtB(x,t)=ϵFν(x,t)T0ν(x,t)-ϵ2∫-∞t‍dt1eϵ(t1-t)Fν(x,t1)T0ν(x,t1).
As T0ν and Fν are finite, the right hand side of (37) goes to zero as ϵ→0. Also integrating (35) by parts, we get
(38)B(x,t)=Fν(x,t)T0ν(x,t)-∫-∞t‍dt1eϵ(t1-t)(FνdT0νdt+dFνdtT0ν).
We now consider the space integral of (38). Using the energy-momentum conservation rule (33), we integrate the second term in (38) by parts and neglect the surface integrals to get
(39)∫‍d3xB(x,t)=∫‍d3xFν(x,t)T0ν(x,t)-∫‍d3x∫-∞t‍dt1eϵ(t1-t)Tμν(x,t1)∂μFν(x,t1)≡A-B,
where we have abbreviated the first and second terms by A and -B, respectively. Then the nonequilibrium statistical density matrix is given by
(40)ρ=e-A+B/Tre-A+B.
The first term A in (39) characterises local equilibrium:
(41)ρ0=e-A/Tre-A
where the second term B including the thermodynamical force ∂μFν describes deviation from equilibrium.In order to expandρ in a series in B we define the function
(42)Q(τ)=e-(1-τ)Aeτ(-A+B)
such that the boundary conditions at τ=0 and τ=1 correspond to the equilibrium and nonequilibrium density matrices as follows:
(43)Q(τ=0)=e-A=ρ0,Q(τ=1)=e-A+B=ρ.
We then differentiate Q(τ) with respect to τ to get
(44)dQ(τ)dτ=e-(1-τ)ABe(1-τ)AQ(τ)
which can be integrated to give
(45)Q(τ)=Q(0)+∫0τ‍dτ′e-(1-τ′)ABe(1-τ′)AQ(τ′).
It can be solved iteratively. Keeping up to the first order term (linear response) and setting τ=1, we get the required result
(46)ρ=ρ0[1+∫01‍e-τABeτA].Applying this formula to the energy-momentum tensor, we get its response to the thermodynamical forces as [10]
(47)〈Tμν(x,t)〉=〈Tμν(x,t)〉0+∫‍d3x′∫-∞t‍dt′eϵ(t1-t)×(Tμν(x,t),Tρσ(x′,t′))∂ρFσ(x′,t′),
where
(48)(Tμν(x,t),Tρσ(x′,t′))=∫01‍dτ{〈Tμν(x,t)e-τATρσ(x′,t′)eτA〉0-〈Tμν(x,t)〉0〈Tρσ(x′,t′)〉0}
is the correlation function to be evaluated. As the correlation is assumed to vanish as t′→-∞, it can be put in terms of the conventional retarded Green’s function. Omitting indices, it is
(49)(T(x,t),T(x′,t′))=1β∫-∞t‍dt′〈T(x,t),T(x′,t′)〉ret
with
(50)〈T(x,t),T(x′,t′)〉ret≡iθ(t-t′)〈[T(x,t),T(x′,t′)]〉0.We now use (47) to obtain the expectation value of the viscous-shear stress part of the nonequilibrium energy momentum tensor which is given by
(51)Tμν=Tμν(0)+πμν+(Pμuν+Pνuμ),
where Tμν(0)=(ϵ+p)uμuν-gμνp is the equilibrium part, πμν is the viscous-shear stress tensor, and Pμ is the heat current. Also, with a view to separate scalar, vector and tensor processes, the quantity Tρσ∂ρFσ in (39) is expanded as
(52)Tρσ∂ρFσ=βπρσ∂ρuσ+βPρ(β∂ρβ+u·∂uρ)-βp~∂·u
with p~=p-cs2ϵ, cs being the sound velocity. Using now the fact that the correlation function between operators of different ranks vanishes in an isotropic medium, one can write from (47)
(53)〈πμν(x,t)〉=〈πμν(x,t)〉0+∫‍d3x′∫-∞t‍dt′eϵ(t1-t)×(πμν(x,t),πρσ(x′,t′))β(x′,t′)∂ρuσ(x′,t′)
with 〈πμν(x,t)〉0=0. Following Hosoya et al. [10], we write the correlation function as
(54)(πμν,πρσ)=110[ΔμρΔνσ+ΔμσΔνρ-23ΔμνΔρσ]×(παβ,παβ),
where Δμν=gμν-uμuν. Assuming now that changes in the thermodynamic forces are small over the correlation length of the two-point function, the factor β∂ρuσ can be taken out of the integral giving finally:
(55)〈πμν(x,t)〉=η[ΔρμΔσν(∂ρuσ+∂σuρ)-23ΔμνΔρσ∂ρuσ],
where
(56)η=110∫-∞0‍dt1eϵ(t1-t)∫-∞t1‍dt′×∫‍d3x′〈παβ(x,t),παβ(x′,t′)〉ret=110∫-∞0‍dt1eϵt1∫-∞t1‍dt′i∫‍d3x′θ(-t′)×〈[παβ(0,0),παβ(x′,t′)]〉0;(x⟶0,t⟶0).Again, starting with the pressurep(x,t) on the l.h.s of (47) and following the steps as described previously, we obtain
(57)〈p(x,t)〉=〈p(x,t)〉0-ζ∂ρuρ(x,t),
where the bulk viscosity ζ is given in terms of a retarded correlation function by
(58)ζ=∫-∞0‍dt1eϵ(t1-t)∫-∞t1‍dt′∫‍d3x′〈p~(x,t),p~(x′,t′)〉ret.
Here p~(x,t)=p(x,t)-cs2ϵ(x,t) with ϵ(x,t) the energy density and cs2=〈∂p/∂ϵ〉0.Recall that〈⋯〉0 denotes equilibrium ensemble average. From now on we will drop the subscript “0” on the correlation functions.


### 2.3. Perturbative Evaluation

Clearly the spectral forms and their interrelations derived in Section2.1 hold also for the two-point function appearing in (56) for the shear viscosity. We begin with four-dimensional Fourier transforms. To calculate the 11 element of the retarded two-point function
(59)Π11R(q)=i∫‍d4xeiq(x-x′)θ(t-t′)×〈[παβ(x,t),παβ(x′,t′)]〉,
we consider the corresponding time-ordered one
(60)Π11(q)=i∫‍d4xeiq(x-x′)〈Tπαβ(x,t)παβ(x′,t′)〉,
which can be calculated perturbatively. The viscous stress tensor can be extracted from the energy momentum tensor using the formula
(61)παβ(x)=(ΔαρΔβσ-13ΔαβΔρσ)Tρσ(x),
where Tμν=-gμνℒ+(∂ℒ/∂(∂μϕ→))·∂νϕ→ in which ϕ→(x) denotes the pion triplet. We take the lowest-order chiral Lagrangian given by [18]
(62)ℒ=12∂μϕ→·∂μϕ→-12mπ2ϕ→·ϕ→+ℒint,ℒint=-16fπ2[∂μϕ→·∂μϕ→ϕ→·ϕ→-ϕ→·∂μϕ→ϕ→·∂μϕ→]+mπ24!fπ2(ϕ→·ϕ→)2.
The time-ordered correlator, to leading order, is then given by Wick contractions of pion fields in παβ which is obtained as
(63)παβ(x)=(ΔαρΔβσ-13ΔαβΔρσ)∂ρϕ→(x)·∂σϕ→(x).
In the so-called skeleton expansion, these contractions are expressed in terms of complete propagators (see Figure 2) to get
(64)Π11(q)=i∫‍d4l(2π)4N(l,q)D11(l)D11(l-q),
where D11 is given by (19) and N(l,q) is determined by the derivatives acting on the pion fields as follows:
(65)N(l,q)=-gπ[l2(l-q)2+13{l·(l-q)}2],
where the pion isospin degeneracy factor gπ=3.Figure 2
The first term in the so-called skeleton expansion of the two-point function. Heavy lines denote full propagators.To work out thel0 integral in (64), it is more convenient to use Λ11 given by (20) than by (21). Closing the contour in the upper or lower half l0-plane, we get
(66)Π11(q)=∫‍d3l(2π)3N(l,q)∫‍dk0′2πρ(k0′,l)×dk0′′2πρ(k0′′,l-q)K(q0,k0′,k0′′),
where
(67)K={1+f(k0′)}f(k0′′)q0-(k0′-k0′′)+iη-f(k0′){1+f(k0′′)}q0-(k0′-k0′′)-iη.
The imaginary part of Π11 arises from the factor K as follows:
(68)ImK=-π[{1+f(k0′)}f(k0′′)+f(k0′){1+f(k0′′)}]×δ(q0-(k0′-k0′′))=-πcoth(βq02){f(k0′′)-f(k0′)}×δ(q0-(k0′-k0′′)),
where its real part is given by the principal value integrals.Having obtained the real and imaginary parts ofΠ11(q), we use relations similar to (24) to build the 11 element of the diagonalised Π matrix:
(69)Π¯=∫‍d3l(2π)3N(l,q)∫‍dk0′2πρ(k0′,l)∫‍dk0′′2πρ(k0′′,l-q)×{1+f(k0′)}f(k0′′)-f(k0′){1+f(k0′′)}q0-(k0′-k0′′)+iηϵ(q0).
Finally Π¯ can be continued to Π11R by a relation similar to (31) as follows:
(70)Π11R=∫‍d3l(2π)3N(l,q)×∫‍dk0′2πρ(k0′,l)dk0′′2πρ(k0′′,l-q)×{1+f(k0′)}f(k0′′)-f(k0′){1+f(k0′′)}q0-(k0′-k0′′)+iη.
Note that in (69), (70), we retain the f(k0′)f(k0′′) terms in the numerator to put it in a more convenient form. Change the signs of k0′ and k0′′ in the first and second term, respectively. Noting relations like 1+f(-k0)=-f(k0) and ρ(-k0)=-ρ(k0), we get
(71)Π11R(q)=∫‍d3l(2π)3N(l,q)×∫‍dk0′2πdk0′′2πρ(k0′,l)ρ(k0′′,l-q)×f(k0′)f(k0′′)W(q0,k0′+k0′′),
where
(72)W=1q0+k0′+k0′′+iη-1q0-(k0′+k0′′)+iη.Returning to the expression (56) for η, we now get the three-dimensional spatial integral of the retarded correlation function by setting q=0 in (59) and Fourier inverting with respect to q0 as follows:
(73)i∫‍d3x′θ(-t′)〈[παβ(0→,0),παβ(x′,t′)]〉=-∫‍dq0eiq0t′Π11R(q0,q=0).
This completes our use of the real-time formulation to get the required result. The integrals appearing in the expression for η have been evaluated in [10, 11], which we describe in the following for completeness.As shown in [10], the integral over t1, t′, and q0 in (56) and (73) may be carried out trivially to give
(74)η=i10ddq0Π11R(q0)|q0=0.
The q0 dependence of Π11R is contained entirely in W as follows:
(75)ddq0W(q0)|q0=0=-1(k0′+k0′′-iη)2+1(k0′+k0′′-iη)2=2πiδ′(k0′+k0′′).
Changing the integration variables in (71) from k0′, k0′′ to k¯0=k0′+k0′′ and k0=(1/2)(k0′-k0′′), we get
(76)η=∫‍d3l(2π)3N(l)∫‍dk0(2π)2F(k0,l),
where
(77)F(k0,l)=ddk¯0{ρ(k¯02+k0,l)ρ(k¯02-k0,l)×f(k¯02+k0)f(k¯02-k0)}|k¯0=0.It turns out that the integral overk0 becomes undefined, if we try to evaluate F(k0) with the free spectral function ρ0(k) given by (26). As pointed out in [10], we have to take the spectral function for the complete propagator that includes the self-energy of the pion, leading to its finite width Γ in the medium
(78)ρ(k0,l)=1i[1(k0-iΓ)2-ω2-1(k0+iΓ)2-ω2],ω=l2+mπ2.
Note that this form of the spectral function trivially follows on replacing iη (where η→0+) with iΓ in the free spectral function (26) which can be written as
(79)ρ0(k0,l)=1i[1(k0-iη)2-ω2-1(k0+iη)2-ω2].ThenF(k0,l) becomes
(80)F=-8k02eβk0(eβk0-1)2βΓ2{(k0-iΓ)2-ω2}2{(k0+iΓ)2-ω2}2,
having double poles at k0=2πin/β for n=±1,±2,… and also at k0=±ω±iΓ. The integral over k0 may now be evaluated by closing the contour in the upper/lower half plane to get
(81)∫-∞+∞‍dk0(2π)2F(k0,l)=-18πβω2Γn(ω){1+n(ω)},
where we retain only the leading (singular) term for small Γ. In this approximation (76) gives
(82)η=gπβ30π2∫0∞‍dll6ω2Γn(ω){1+n(ω)}.Proceeding analogously as mentioned above, the lowest order contribution to the bulk viscosity can be obtained as [10]
(83)ζ=gπβ4π2∫0∞‍dll2(l2/3-cs2ω2)2ω2Γn(ω){1+n(ω)}.The widthΓ(l) at different temperatures is known [30] from chiral perturbation theory. The quantity Γ can also be interpreted as the collision frequency, and the inverse of which is the relaxation time τ. For collisions of the form π(l)+π(k)→π(l′)+π(k′) this is given by (see, e.g., [25])(84)Γ(l)=τ-1(l)=∫‍k2dk2π2ωks(s-4mπ2)2ωln(ωk)×(1+n(ωl′))(1+n(ωk′))12∫‍dΩdσdΩ,
where dσ/dΩ is the ππ cross-section. Note that the lowest order formulae for the shear and bulk viscosities obtained previously in the linear response approach coincide with the expressions which result from solving the transport equation in the relaxation-time approximation.


## 2.1. Real-Time Formulation

In this section we review the real-time formulation of equilibrium thermal field theory leading to the spectral representations of bosonic two-point functions [16]. This formulation begins with a comparison between the time evolution operator e-iH(t2-t1) of quantum theory and the Boltzmann weight e-βH=e-iH(τ-iβ-τ) of statistical physics, where we introduce τ as a complex variable. Thus while for the time evolution operator, the times t1 and t2(t2>t1) are any two points on the real line, the Boltzmann weight involves a path from τ to τ-iβ in the complex time plane. Setting this τ=-T, where T is real, positive, and large, we can get the contour C shown in Figure 1, lying within the region of analyticity in this plane and accommodating real-time correlation functions [6, 8].Figure 1
The contourC in the complex time plane used here for the real-time formulation.Let a general bosonic interacting field in the Heisenberg representation be denoted byΦl(x), whose subscript l collects the index (or indices) denoting the field component and derivatives acting on it. Although we will call its two-point function as propagator, Φl(x) can be an elementary field or a composite local operator. (If Φl(x) denotes the pion field, it will, of course, not have any index.)The thermal expectation value of the productΦl(x)Φl′†(x′) may be expressed as
(1)〈Φl(x)Φl′†(x′)〉=1Z∑m,n‍e-βEm×〈m|Φl(x)|n〉〈n|Φl′†(x′)|m〉,
where 〈O〉 for any operator O denotes equilibrium ensemble average:
(2)〈O〉=Tr(e-βHO)Z,Z=Tre-βH.
Note that we have two sums in (1), one to evaluate the trace and the other to separate the field operators. They run over a complete set of states, which we choose as eigenstates |m〉 of four-momentum Pμ. Using translational invariance of the field operator:
(3)Φl(x)=eiP·xΦl(0)e-iP·x,
we get
(4)〈Φl(x)Φl′†(x′)〉=1Z∑m,n‍e-βEmei(km-kn)·(x-x′)×〈m|Φl(0)|n〉〈n|Φl′†(0)|m〉.
Its spatial Fourier transform is
(5)∫‍d3xe-ik·(x-x′)〈Φl(x)Φl′†(x′)〉=(2π)3Z∑m,n‍e-βEmei(Em-En)(τ-τ′)δ3(km-kn+k)×〈m|Φl(0)|n〉〈n|Φl′†(0)|m〉,
where the times τ, τ′ are on the contour C. We now insert unity on the left of (5) in the form
(6)1=∫-∞∞‍dk0′δ(Em-En+k0′).
(We reserve k0 for the variable conjugate to the real-time.) Then it may be written as
(7)∫d3xe-ik·(x-x′)〈Φl(x)Φl′†(x′)〉=∫‍dk0′2πe-ik0′(τ-τ′)Mll′+(k0′,k),
where the spectral function M+ is given by [kμ′=(k0′,k)] as follows:
(8)Mll′+(k′)=(2π)4Z∑m,n‍e-βEmδ4(km-kn+k′)×〈m|Φl(0)|n〉〈n|Φl′†(0)|m〉.In just the same way, we can work out the Fourier transform of〈Φl′†(x′)Φl(x)〉 as follows:
(9)∫d3xe-ik·(x-x′)〈Φl′†(x′)Φl(x)〉=∫‍dk0′2πe-ik0′(τ-τ′)Mll′-(k0′,k),
with a second spectral function M- given by
(10)Mll′-(k′)=(2π)4Z∑m,n‍e-βEmδ4(kn-km+k′)×〈m|Φl′†(0)|n〉〈n|Φl(0)|m〉.
The two spectral functions are related by the KMS relation [26, 27]
(11)Mll′+(k)=eβk0Mll′-(k),
in momentum space, which may be obtained simply by interchanging the dummy indices m, n in one of Mll′±(k) and using the energy conserving δ-function.We next introduce thedifference of the two spectral functions:
(12)ρll′(k)≡Mll′+(k)-Mll′-(k)
and solve this identity and the KMS relation (11) for Mll′±(k) as follows:
(13)Mll′+(k)={1+f(k0)}ρll′(k),Mll′-(k)=f(k0)ρll′(k),
where f(k0) is the distribution-like function
(14)f(k0)=1eβk0-1,-∞<k0<∞.
In terms of the true distribution function
(15)n(|k0|)=1eβ|k0|-1,
it may be expressed as
(16)f(k0)=f(k0){θ(k0)+θ(-k0)}=nϵ(k0)-θ(-k0).With the previous ingredients, we can build the spectral representations for the two types of thermal propagators. First consider thetime-ordered one:
(17)-iDll′(x,x′)=〈TcΦl(x)Φl′†(x′)〉=θc(τ-τ′)〈Φl(x)Φl′†(x′)〉+θc(τ′-τ)〈Φl′†(x′)Φl(x)〉.
Using (7), (9), and (13), we see that its spatial Fourier transform is given by [6]
(18)Dll′(τ-τ′,k)=i∫-∞∞‍dk0′2πρll′(k0′,k)×e-ik0′(τ-τ′){θc(τ-τ′)+f(k0′)}.AsT→∞, the contour of Figure 1 simplifies, reducing essentially to two parallel lines, one the real axis and the other shifted by -iβ/2, points which will be denoted, respectively, by subscripts 1 and 2, so that τ1=t,τ2=t-iβ/2 [8]. The propagator then consists of four pieces, which may be put in the form of a 2×2 matrix. The contour ordered θ′s may now be converted to the usual time ordered ones. If τ, τ′ are both on line 1 (the real axis), the τ and t orderings coincide, θc(τ1-τ1′)=θ(t-t′). If they are on two different lines, the τ ordering is definite, θc(τ1-τ2′)=0,θc(τ2-τ1′)=1. Finally if they are both on line 2, the two orderings are opposite, θc(τ2-τ2′)=θ(t′-t).Back to real-time, we can work out the usual temporal Fourier transform of the components of the matrix to get(19)Dll′(k0,k)=∫-∞∞‍dk0′2πρll′(k0′,k)Λ(k0′,k0),
where the elements of the matrix Λ are given by [16]
(20)Λ11=-Λ22*=1k0′-k0-iη+2πif(k0′)δ(k0′-k0),Λ12=Λ21=2πieβk0′/2f(k0′)δ(k0′-k0).
Using relation (16), we may rewrite (20) in terms of n as follows:
(21)Λ11=-Λ22*=1k0′-k0-iηϵ(k0)+2πinϵ(k0)δ(k0′-k0),Λ12=Λ21=2πin(1+n)ϵ(k0)δ(k0′-k0).The matrixΛ and hence the propagator Dll′ can be diagonalised to give
(22)Dll′(k0,k)=U(D¯ll′00-D¯ll′*)U,
where D¯ll′ and U are given by
(23)D¯ll′(k0,k)=∫-∞∞‍dk0′2πρll′(k0′,k)k0′-k0-iηϵ(k0),U=(1+nnn1+n).
Equation (22) shows that D¯ can be obtained from any of the elements of the matrix D, say D11. Omitting the indices ll′, we get
(24)ReD¯=ReD11,ImD¯=tanh(β|k0|2)ImD11.Looking back at the spectral functionsMll′± defined by (8) and (10), we can express them as usual four-dimensional Fourier transforms of ensemble average of the operator products, so that ρll′ is the Fourier transform of that of the commutator:
(25)ρll′(k0,k)=∫‍d4yeik·(y-y′)〈[Φl(y),Φl′(y′)]〉,
where the time components of y and y′ are on the real axis in the τ-plane. Taking the spectral function for the free scalar field
(26)ρ0=2πϵ(k0)δ(k2-m2),
we see that D¯ becomes the free propagator, D¯(k0,k)=-1/(k2-m2).We next consider theretarded thermal propagator:
(27)Dll′R(x,x′)=iθc(τ-τ′)〈[Φl(x,τ),Φl′(x′,τ′)]〉,
where again τ, τ′ are on the contour C (Figure 1). Noting equations (7), (9), and (12) the three-dimensional Fourier transform may immediately be written as(28)Dll′R(τ-τ′,k)=iθc(τ-τ′)∫-∞∞‍dk0′2πe-ik0′(τ-τ′)ρll′(k0′,k).
As before we isolate the different components with real-times and take the Fourier transform with respect to real-time. Thus for the 11 component we simply have
(29)Dll′R(t-t′,k)11=iθ(t-t′)∫-∞∞‍dk0′2πe-ik0′(t-t′)ρll′(k0′,k),
whose temporal Fourier transform gives
(30)Dll′R(k0,k)11=∫-∞∞‍dk0′2πρll′(k0′,k)k0′-k0-iη.
This 11 component suffices for us, but we also display the complete matrix as follows:(31)Dll′R(k0,k)=(Dll′R(k)110ρll′(k){nn+1θ(k0)+n+1nθ(-k0)}-Dll′R*(k)11).Though we deal with matrices in real-time formulation, it is the 11 component that is physical. Equations (23) and (30) then show that we can continue the time-ordered two-point function into the retarded one by simply changing the iϵ prescription as follows:
(32)Dll′R(k0+iη,k)11=D¯ll′(k0+iηϵ(q0)⟶k0+iη,k).
The point to note here is that for the time-ordered propagator, it is the diagonalised matrix and not the matrix itself, whose 11 component can be continued in a simple way.


## 2.2. Transport Coefficients

We now use the linear response approach to arrive at expressions of the transport coefficients as integrals of retarded Green’s functions over space. We follow the method proposed by Zubarev [28], which is excellently reviewed in [10]. Here the system is supposed to be in the hydrodynamical stage where the mean free time of the constituent particles is much shorter than the relaxation time of the whole system under consideration. Thus local equilibrium will be attained quickly, while global equilibrium will be approaching gradually. Since the system is assumed to be not far from equilibrium, we may retain only linear terms in space-time gradients of thermodynamical parameters, like temperature and velocity fields. We assume the energy-momentum of the system to be conserved as follows:
(33)∂μTμν(x,t)=0.
The nonequilibrium density matrix operator is constructed in the Heisenberg picture, where it is independent of time as follows:
(34)dρdt=0.
Following Zubarev, we construct the operator B(x,t) as follows:
(35)B(x,t)=ϵ∫-∞t‍dt1eϵ(t1-t)Fν(x,t1)T0ν(x,t1),(ϵ⟶0+),
where Fν(x,t)=β(x,t)uν(x,t). Here β(x,t) is a Lorentz invariant quantity defining the local temperature, and uν(x,t) is the four-velocity field of the fluid
(36)uν(x,t)uν(x,t)=1.
The construction (35), which smooths out the oscillating terms, resembles the one used in the formal theory of scattering [28, 29] and selects out the retarded solution.The expression (35) is actually independent of t; the time derivative is
(37)ddtB(x,t)=ϵFν(x,t)T0ν(x,t)-ϵ2∫-∞t‍dt1eϵ(t1-t)Fν(x,t1)T0ν(x,t1).
As T0ν and Fν are finite, the right hand side of (37) goes to zero as ϵ→0. Also integrating (35) by parts, we get
(38)B(x,t)=Fν(x,t)T0ν(x,t)-∫-∞t‍dt1eϵ(t1-t)(FνdT0νdt+dFνdtT0ν).
We now consider the space integral of (38). Using the energy-momentum conservation rule (33), we integrate the second term in (38) by parts and neglect the surface integrals to get
(39)∫‍d3xB(x,t)=∫‍d3xFν(x,t)T0ν(x,t)-∫‍d3x∫-∞t‍dt1eϵ(t1-t)Tμν(x,t1)∂μFν(x,t1)≡A-B,
where we have abbreviated the first and second terms by A and -B, respectively. Then the nonequilibrium statistical density matrix is given by
(40)ρ=e-A+B/Tre-A+B.
The first term A in (39) characterises local equilibrium:
(41)ρ0=e-A/Tre-A
where the second term B including the thermodynamical force ∂μFν describes deviation from equilibrium.In order to expandρ in a series in B we define the function
(42)Q(τ)=e-(1-τ)Aeτ(-A+B)
such that the boundary conditions at τ=0 and τ=1 correspond to the equilibrium and nonequilibrium density matrices as follows:
(43)Q(τ=0)=e-A=ρ0,Q(τ=1)=e-A+B=ρ.
We then differentiate Q(τ) with respect to τ to get
(44)dQ(τ)dτ=e-(1-τ)ABe(1-τ)AQ(τ)
which can be integrated to give
(45)Q(τ)=Q(0)+∫0τ‍dτ′e-(1-τ′)ABe(1-τ′)AQ(τ′).
It can be solved iteratively. Keeping up to the first order term (linear response) and setting τ=1, we get the required result
(46)ρ=ρ0[1+∫01‍e-τABeτA].Applying this formula to the energy-momentum tensor, we get its response to the thermodynamical forces as [10]
(47)〈Tμν(x,t)〉=〈Tμν(x,t)〉0+∫‍d3x′∫-∞t‍dt′eϵ(t1-t)×(Tμν(x,t),Tρσ(x′,t′))∂ρFσ(x′,t′),
where
(48)(Tμν(x,t),Tρσ(x′,t′))=∫01‍dτ{〈Tμν(x,t)e-τATρσ(x′,t′)eτA〉0-〈Tμν(x,t)〉0〈Tρσ(x′,t′)〉0}
is the correlation function to be evaluated. As the correlation is assumed to vanish as t′→-∞, it can be put in terms of the conventional retarded Green’s function. Omitting indices, it is
(49)(T(x,t),T(x′,t′))=1β∫-∞t‍dt′〈T(x,t),T(x′,t′)〉ret
with
(50)〈T(x,t),T(x′,t′)〉ret≡iθ(t-t′)〈[T(x,t),T(x′,t′)]〉0.We now use (47) to obtain the expectation value of the viscous-shear stress part of the nonequilibrium energy momentum tensor which is given by
(51)Tμν=Tμν(0)+πμν+(Pμuν+Pνuμ),
where Tμν(0)=(ϵ+p)uμuν-gμνp is the equilibrium part, πμν is the viscous-shear stress tensor, and Pμ is the heat current. Also, with a view to separate scalar, vector and tensor processes, the quantity Tρσ∂ρFσ in (39) is expanded as
(52)Tρσ∂ρFσ=βπρσ∂ρuσ+βPρ(β∂ρβ+u·∂uρ)-βp~∂·u
with p~=p-cs2ϵ, cs being the sound velocity. Using now the fact that the correlation function between operators of different ranks vanishes in an isotropic medium, one can write from (47)
(53)〈πμν(x,t)〉=〈πμν(x,t)〉0+∫‍d3x′∫-∞t‍dt′eϵ(t1-t)×(πμν(x,t),πρσ(x′,t′))β(x′,t′)∂ρuσ(x′,t′)
with 〈πμν(x,t)〉0=0. Following Hosoya et al. [10], we write the correlation function as
(54)(πμν,πρσ)=110[ΔμρΔνσ+ΔμσΔνρ-23ΔμνΔρσ]×(παβ,παβ),
where Δμν=gμν-uμuν. Assuming now that changes in the thermodynamic forces are small over the correlation length of the two-point function, the factor β∂ρuσ can be taken out of the integral giving finally:
(55)〈πμν(x,t)〉=η[ΔρμΔσν(∂ρuσ+∂σuρ)-23ΔμνΔρσ∂ρuσ],
where
(56)η=110∫-∞0‍dt1eϵ(t1-t)∫-∞t1‍dt′×∫‍d3x′〈παβ(x,t),παβ(x′,t′)〉ret=110∫-∞0‍dt1eϵt1∫-∞t1‍dt′i∫‍d3x′θ(-t′)×〈[παβ(0,0),παβ(x′,t′)]〉0;(x⟶0,t⟶0).Again, starting with the pressurep(x,t) on the l.h.s of (47) and following the steps as described previously, we obtain
(57)〈p(x,t)〉=〈p(x,t)〉0-ζ∂ρuρ(x,t),
where the bulk viscosity ζ is given in terms of a retarded correlation function by
(58)ζ=∫-∞0‍dt1eϵ(t1-t)∫-∞t1‍dt′∫‍d3x′〈p~(x,t),p~(x′,t′)〉ret.
Here p~(x,t)=p(x,t)-cs2ϵ(x,t) with ϵ(x,t) the energy density and cs2=〈∂p/∂ϵ〉0.Recall that〈⋯〉0 denotes equilibrium ensemble average. From now on we will drop the subscript “0” on the correlation functions.


## 2.3. Perturbative Evaluation

Clearly the spectral forms and their interrelations derived in Section2.1 hold also for the two-point function appearing in (56) for the shear viscosity. We begin with four-dimensional Fourier transforms. To calculate the 11 element of the retarded two-point function
(59)Π11R(q)=i∫‍d4xeiq(x-x′)θ(t-t′)×〈[παβ(x,t),παβ(x′,t′)]〉,
we consider the corresponding time-ordered one
(60)Π11(q)=i∫‍d4xeiq(x-x′)〈Tπαβ(x,t)παβ(x′,t′)〉,
which can be calculated perturbatively. The viscous stress tensor can be extracted from the energy momentum tensor using the formula
(61)παβ(x)=(ΔαρΔβσ-13ΔαβΔρσ)Tρσ(x),
where Tμν=-gμνℒ+(∂ℒ/∂(∂μϕ→))·∂νϕ→ in which ϕ→(x) denotes the pion triplet. We take the lowest-order chiral Lagrangian given by [18]
(62)ℒ=12∂μϕ→·∂μϕ→-12mπ2ϕ→·ϕ→+ℒint,ℒint=-16fπ2[∂μϕ→·∂μϕ→ϕ→·ϕ→-ϕ→·∂μϕ→ϕ→·∂μϕ→]+mπ24!fπ2(ϕ→·ϕ→)2.
The time-ordered correlator, to leading order, is then given by Wick contractions of pion fields in παβ which is obtained as
(63)παβ(x)=(ΔαρΔβσ-13ΔαβΔρσ)∂ρϕ→(x)·∂σϕ→(x).
In the so-called skeleton expansion, these contractions are expressed in terms of complete propagators (see Figure 2) to get
(64)Π11(q)=i∫‍d4l(2π)4N(l,q)D11(l)D11(l-q),
where D11 is given by (19) and N(l,q) is determined by the derivatives acting on the pion fields as follows:
(65)N(l,q)=-gπ[l2(l-q)2+13{l·(l-q)}2],
where the pion isospin degeneracy factor gπ=3.Figure 2
The first term in the so-called skeleton expansion of the two-point function. Heavy lines denote full propagators.To work out thel0 integral in (64), it is more convenient to use Λ11 given by (20) than by (21). Closing the contour in the upper or lower half l0-plane, we get
(66)Π11(q)=∫‍d3l(2π)3N(l,q)∫‍dk0′2πρ(k0′,l)×dk0′′2πρ(k0′′,l-q)K(q0,k0′,k0′′),
where
(67)K={1+f(k0′)}f(k0′′)q0-(k0′-k0′′)+iη-f(k0′){1+f(k0′′)}q0-(k0′-k0′′)-iη.
The imaginary part of Π11 arises from the factor K as follows:
(68)ImK=-π[{1+f(k0′)}f(k0′′)+f(k0′){1+f(k0′′)}]×δ(q0-(k0′-k0′′))=-πcoth(βq02){f(k0′′)-f(k0′)}×δ(q0-(k0′-k0′′)),
where its real part is given by the principal value integrals.Having obtained the real and imaginary parts ofΠ11(q), we use relations similar to (24) to build the 11 element of the diagonalised Π matrix:
(69)Π¯=∫‍d3l(2π)3N(l,q)∫‍dk0′2πρ(k0′,l)∫‍dk0′′2πρ(k0′′,l-q)×{1+f(k0′)}f(k0′′)-f(k0′){1+f(k0′′)}q0-(k0′-k0′′)+iηϵ(q0).
Finally Π¯ can be continued to Π11R by a relation similar to (31) as follows:
(70)Π11R=∫‍d3l(2π)3N(l,q)×∫‍dk0′2πρ(k0′,l)dk0′′2πρ(k0′′,l-q)×{1+f(k0′)}f(k0′′)-f(k0′){1+f(k0′′)}q0-(k0′-k0′′)+iη.
Note that in (69), (70), we retain the f(k0′)f(k0′′) terms in the numerator to put it in a more convenient form. Change the signs of k0′ and k0′′ in the first and second term, respectively. Noting relations like 1+f(-k0)=-f(k0) and ρ(-k0)=-ρ(k0), we get
(71)Π11R(q)=∫‍d3l(2π)3N(l,q)×∫‍dk0′2πdk0′′2πρ(k0′,l)ρ(k0′′,l-q)×f(k0′)f(k0′′)W(q0,k0′+k0′′),
where
(72)W=1q0+k0′+k0′′+iη-1q0-(k0′+k0′′)+iη.Returning to the expression (56) for η, we now get the three-dimensional spatial integral of the retarded correlation function by setting q=0 in (59) and Fourier inverting with respect to q0 as follows:
(73)i∫‍d3x′θ(-t′)〈[παβ(0→,0),παβ(x′,t′)]〉=-∫‍dq0eiq0t′Π11R(q0,q=0).
This completes our use of the real-time formulation to get the required result. The integrals appearing in the expression for η have been evaluated in [10, 11], which we describe in the following for completeness.As shown in [10], the integral over t1, t′, and q0 in (56) and (73) may be carried out trivially to give
(74)η=i10ddq0Π11R(q0)|q0=0.
The q0 dependence of Π11R is contained entirely in W as follows:
(75)ddq0W(q0)|q0=0=-1(k0′+k0′′-iη)2+1(k0′+k0′′-iη)2=2πiδ′(k0′+k0′′).
Changing the integration variables in (71) from k0′, k0′′ to k¯0=k0′+k0′′ and k0=(1/2)(k0′-k0′′), we get
(76)η=∫‍d3l(2π)3N(l)∫‍dk0(2π)2F(k0,l),
where
(77)F(k0,l)=ddk¯0{ρ(k¯02+k0,l)ρ(k¯02-k0,l)×f(k¯02+k0)f(k¯02-k0)}|k¯0=0.It turns out that the integral overk0 becomes undefined, if we try to evaluate F(k0) with the free spectral function ρ0(k) given by (26). As pointed out in [10], we have to take the spectral function for the complete propagator that includes the self-energy of the pion, leading to its finite width Γ in the medium
(78)ρ(k0,l)=1i[1(k0-iΓ)2-ω2-1(k0+iΓ)2-ω2],ω=l2+mπ2.
Note that this form of the spectral function trivially follows on replacing iη (where η→0+) with iΓ in the free spectral function (26) which can be written as
(79)ρ0(k0,l)=1i[1(k0-iη)2-ω2-1(k0+iη)2-ω2].ThenF(k0,l) becomes
(80)F=-8k02eβk0(eβk0-1)2βΓ2{(k0-iΓ)2-ω2}2{(k0+iΓ)2-ω2}2,
having double poles at k0=2πin/β for n=±1,±2,… and also at k0=±ω±iΓ. The integral over k0 may now be evaluated by closing the contour in the upper/lower half plane to get
(81)∫-∞+∞‍dk0(2π)2F(k0,l)=-18πβω2Γn(ω){1+n(ω)},
where we retain only the leading (singular) term for small Γ. In this approximation (76) gives
(82)η=gπβ30π2∫0∞‍dll6ω2Γn(ω){1+n(ω)}.Proceeding analogously as mentioned above, the lowest order contribution to the bulk viscosity can be obtained as [10]
(83)ζ=gπβ4π2∫0∞‍dll2(l2/3-cs2ω2)2ω2Γn(ω){1+n(ω)}.The widthΓ(l) at different temperatures is known [30] from chiral perturbation theory. The quantity Γ can also be interpreted as the collision frequency, and the inverse of which is the relaxation time τ. For collisions of the form π(l)+π(k)→π(l′)+π(k′) this is given by (see, e.g., [25])(84)Γ(l)=τ-1(l)=∫‍k2dk2π2ωks(s-4mπ2)2ωln(ωk)×(1+n(ωl′))(1+n(ωk′))12∫‍dΩdσdΩ,
where dσ/dΩ is the ππ cross-section. Note that the lowest order formulae for the shear and bulk viscosities obtained previously in the linear response approach coincide with the expressions which result from solving the transport equation in the relaxation-time approximation.


## 3. Viscous Coefficients in the Kinetic Theory Approach

The kinetic theory approach is suitable for studying transport properties of dilute systems. Here one assumes that the system is characterized by a distribution function which gives the phase space probability density of the particles making up the fluid. Except during collisions, these (on-shell) particles are assumed to propagate classically with well-defined position, momenta, and energy. It is possible to obtain the nonequilibrium distribution function by solving the transport equation in the hydrodynamic regime by expanding the distribution function in a local equilibrium part along with nonequilibrium corrections. This expansion in terms of gradients of the velocity field is used to linearize the transport equation. The coefficients of expansion which are related to the transport coefficients satisfy linear integral equations. The standard method of solution involves the use of polynomial functions to reduce these integral equations to algebraic ones.


### 3.1. Transport Coefficients at First Chapman-Enskog Order

The evolution of the phase space distribution of the pions is governed by the (transport) equation:(85)pμ∂μf(x,p)=C[f],
where C[f] is the collision integral. For binary elastic collisions p+k→p′+k′ which we consider, this is given by [22]
(86)C[f]=∫‍dΓkdΓp′dΓk′×[f(x,p′)f(x,k′){1+f(x,p)}×{1+f(x,k)}-f(x,p)f(x,k)×{1+f(x,p′)}{1+f(x,k′)}]W,
where the interaction rate is as follows:
(87)W=s2dσdΩ(2π)6δ4(p+k-p′-k′),
and dΓq=d3q/(2π)3q0. The 1/2 factor comes from the indistinguishability of the initial state pions.For small deviation from local equilibrium, we write, in the first Chapman-Enskog approximation,(88)f(x,p)=f(0)(x,p)+δf(x,p),δf(x,p)=f(0)(x,p)[1+f(0)(x,p)]ϕ(x,p),
where the equilibrium distribution function (in the new notation) is given by
(89)f(0)(x,p)=[e(pμuμ(x)-μπ(x))/T(x)-1]-1,
with T(x), uμ(x), and μπ(x) representing the local temperature, flow velocity, and pion chemical potential, respectively. Putting (88) in (85), the deviation function ϕ(x,p) is seen to satisfy
(90)pμ∂μf(0)(x,p)=-ℒ[ϕ],
where the linearized collision term
(91)ℒ[ϕ]=f(0)(x,p)∫‍dΓkdΓp′dΓk′f(0)(x,k)×{1+f(0)(x,p′)}{1+f(0)(x,k′)}×[ϕ(x,p)+ϕ(x,k)-ϕ(x,p′)-ϕ(x,k′)]W.
Using the form of f(0)(x,p) as given in (89) on the left side of (90) and eliminating time derivatives with the help of equilibrium thermodynamic laws, we arrive at [31]
(92)[Q∂νuν+pμΔμν(p·u-h)(T-1∂νT-Duν)-〈pμpν〉〈∂μuν〉]f(0)(1+f(0))=-Tℒ[ϕ],
where D=uμ∂μ, ∇μ=Δμν∂ν, Δμν=gμν-uμuν, and 〈·〉 indicates a space-like symmetric and traceless combination. In this equation,
(93)Q=-13mπ2+(p·u)2{43-γ′}+p·u{(γ′′-1)h-γ′′′T},
where(94)γ′=(S20/S21)2-(S30/S21)2+4z-1S20S31/(S21)2+z-1S30/S21(S20/S21)2-(S30/S21)2+3z-1S20S31/(S21)2+2z-1S30/S21-z-2,γ′′=1+z-2(S20/S21)2-(S30/S21)2+3z-1S20S31/(S21)2+2z-1S30/S21-z-2,γ′′′=S20/S21+5z-1S31/S21-S30S31/(S21)2(S20/S21)2-(S30/S21)2+3z-1S20S31/(S21)2+2z-1S30/S21-z-2,with z=mπ/T and h=mπS31/S21. The functions Snα(z) are integrals over Bose functions [31] and are defined as Snα(z)=∑k=1∞‍ekμπ/Tk-αKn(kz), Kn(x) denoting the modified Bessel function of order n. The left hand side of (90) is thus expressed in terms of thermodynamic forces with different tensorial ranks. In order to be a solution of this equation ϕ must also be a linear combination of the corresponding thermodynamic forces. It is typical to take ϕ as
(95)ϕ=A∂·u+Bμ∇μν(T-1∂νT-Duν)-Cμν〈∂μuν〉
which on substitution into (92) and comparing coefficients of the (independent) thermodynamic forces on both sides yields the set of equations:
(96)ℒ[A]=-Qf(0)(p){1+f(0)(p)}T,ℒ[Cμν]=-〈pμpν〉f(0)(p){1+f(0)(p)}T
ignoring the equation for Bμ which is related to thermal conductivity. These integral equations are to be solved to get the coefficients A and Cμν. It now remains to link these to the viscous coefficients ζ and η. This is achieved by means of the dissipative part of the energy-momentum tensor resulting from the use of the nonequilibrium distribution function (88) in(97)Tμν=∫‍dΓppμpνf(p)=Tμν(0)+ΔTμν,
where
(98)ΔTμν=∫‍dΓpf(0)(1+f(0))Cαβ〈pαpβ〉〈∂μuν〉+∫‍dΓpf(0)(1+f(0))QAΔμν∂σuσ.
Again, for a small deviation ϕ(x,p), close to equilibrium, so that only first-order derivatives contribute, and the dissipative tensor can be generally expressed in the form [32, 33]
(99)ΔTμν=-2η〈∂μuν〉-ζΔμν∂σuσ.Comparing, we obtain the expressions of shear and bulk viscosity as follows:(100)η=-110∫‍dΓpCμν〈pμpν〉f(0)(p){1+f(0)(p)},ζ=-∫‍dΓpQAf(0)(p){1+f(0)(p)}.
The coefficients A and Cμν are perturbatively obtained from (96) by expanding in terms of orthogonal polynomials which reduces the integral equations to algebraic ones. After a tedious calculation using the Laguerre polynomial of 1/2 integral order, the first approximation to the shear and bulk viscosity comes out as
(101)η=T10γ02c00,ζ=Tα22a22,
where
(102)γ0=-10S32(z)S21(z),c00=16{I1(z)+I2(z)+13I3(z)},α2=z32[13(S30S21-z-1)+(S20S21+3zS31S21){(1-γ′′)S31S21+γ′′′z-1}-(43-γ′){S30S21+15z-2S32S21+2z-1}],a22=2z2I3(z).
The integrals Iα(z) are given by
(103)Iα(z)=z4[S21(z)]2e(-2μπ/T)×∫0∞‍dψcosh3ψsinh7ψ×∫0π‍dΘsinΘ12dσdΩ(ψ,Θ)∫02π‍dϕ×∫0∞‍dχsinh2αχ∫0π‍dθsinθ×e2zcoshψcoshχ(eE-1)(eF-1)(eG-1)(eH-1)Mα(θ,Θ),
where μπ is the chemical potential of pions. The exponents in the Bose functions are given by
(104)E=z(coshψcoshχ-sinhψsinhχcosθ)-μπT,F=z(coshψcoshχ-sinhψsinhχcosθ′)-μπT,G=E+2zsinhψsinhχcosθ,H=F+2zsinhψsinhχcosθ′,and the functionsMα(θ,Θ) represent
(105)M1(θ,Θ)=1-cos2Θ,M2(θ,Θ)=cos2θ+cos2θ′-2cosθcosθ′cosΘ,M3(θ,Θ)=[cos2θ-cos2θ′]2.
The relative angle θ′ is defined by cosθ′=cosθcosΘ-sinθsinΘcosϕ.Note that the differential cross-section which appears in the denominator is the dynamical input in the expressions forη and ζ. It is this quantity that we turn to in the next section.


### 3.2. Theππ Cross-Section with Medium Effects

The strong interaction dynamics of the pions enters the collision integrals through the cross-section. In Figure3, we show the ππ cross-section as a function of the centre of mass energy of scattering. The different curves are explained in the following. The filled squares referred to as experiment is a widely used resonance saturation parametrization [23, 34] of isoscalar and isovector phase shifts obtained from various empirical data involving the ππ system. The isospin averaged differential cross-section is given by
(106)dσ(s)dΩ=4qcm2[19sin2δ00+59sin2δ02+13·9sin2δ11cos2θ],
where
(107)δ00=π2+arctan(E-mσΓσ/2),δ11=π2+arctan(E-mρΓρ/2),δ02=-0.12pmπ.
The widths are given by Γσ=2.06p and Γρ=0.095p((p/mπ)/(1+(p/mρ)2))2 with mσ=5.8mπ and mρ=5.53mπ.Figure 3
Theππ cross-section as a function of centre of mass energy. The dotted line indicates the cross-section obtained using (109) which agrees well with the experimental values (106) shown by filled squares. The dashed and solid lines depict the in-medium cross-section for ππ and πh(h=π,ω,h1,a1) loops, respectively, in the ρ self-energy evaluated at T= 160 MeV.To get a handle on the dynamics, we now evaluate theππ cross-section involving ρ and σ meson exchange processes using the interaction Lagrangian as follows:
(108)ℒ=gρρ→μ·π→×∂μπ→+12gσmσπ→·π→σ,
where gρ=6.05 and gσ=2.5. In the matrix elements corresponding to s-channel ρ and σ exchange diagrams which appear for total isospin I=1 and 0, respectively, we introduce a decay width in the corresponding propagator. We get [25]
(109)ℳI=0=gσ2mσ2[3s-mσ2+imσΓσ+1t-mσ2+1u-mσ2]+2gρ2[s-ut-mρ2+s-tu-mρ2],ℳI=1=gσ2mσ2[1t-mσ2-1u-mσ2]+gρ2[2(t-u)s-mρ2+imρΓρ+t-su-mρ2-u-st-mρ2],ℳI=2=gσ2mσ2[1t-mσ2+1u-mσ2]+gρ2[u-st-mρ2+t-su-mρ2].
The differential cross-section is then obtained from dσ/dΩ=|ℳ|2¯/64π2s, where the isospin averaged amplitude is given by |ℳ|2¯=(1/9)∑‍(2I+1)|ℳI|2¯.The integrated cross-section, after ignoring theI=2 contribution is shown by the dotted line (indicated by “vacuum”) in Figure 3 and is seen to agree reasonably well with the experimental cross-section up to a centre of mass energy of about 1 GeV beyond which the theoretical estimate gives higher values. We hence use the experimental cross-section beyond this energy.After this normalisation to data, we now turn to the in-medium cross-section by introducing the effective propagator for theρ in the previous expressions for the matrix elements. This is obtained in terms of the self-energy by solving the Dyson equation and is given by
(110)Dμν=Dμν(0)+Dμσ(0)ΠσλDλν,
where Dμν(0) is the vacuum propagator for the ρ meson and Πσλ is the self-energy function obtained from one-loop diagrams shown in Figure 4. The standard procedure [35] to solve this equation in the medium is to decompose the self-energy into transverse and longitudinal components. For the case at hand the difference between these components is found to be small and is hence ignored. We work with the polarization averaged self-energy function defined as
(111)Π=13(2ΠT+q2ΠL),
where
(112)ΠT=-12(Πμμ+q2q-2Π00),ΠL=1q-2Π00,Π00≡uμuνΠμν.
The in-medium propagator is then written as
(113)Dμν(q0,q→)=-gμν+qμqν/q2q2-mρ2-ReΠ(q0,q→)+iImΠ(q0,q→).
The scattering, decay, and regeneration processes which cause a gain or loss of ρ mesons in the medium are responsible for the imaginary part of its self-energy. The real part on the other hand modifies the position of the pole of the spectral function.Figure 4
π
h self-energy diagrams where h stands for π, ω, h1, and a1 mesons.As discussed in Section2.1, in the real-time formulation of thermal field theory the self-energy assumes a 2×2 matrix structure of which the 11 component is given by
(114)Πμν11(q)=i∫d4k(2π)4Nμν(q,k)Dπ11(k)Dh11(q-k),
where D11 is the 11 component of the scalar propagator given by D11(k)=Δ(k)+2πif(0)(k)δ(k2-m2). It turns out that the self-energy function mentioned above can be obtained in terms of the 11 component through the relations [35, 36]
(115)ReΠμν=ReΠμν11,ImΠμν=ϵ(q0)tanh(βq02)ImΠμν11.
Tensor structures associated with the two vertices and the vector propagator are included in Nμν and are available in [37] where the interactions were taken from chiral perturbation theory. It is easy to perform the integral over k0 using suitable contours to obtain
(116)Πμν(q0,q)=∫d3k(2π)314ωπωh×[(1+f(0)(ωπ))N1μν+f(0)(ωh)N3μνq0-ωπ-ωh+iηϵ(q0)+-f(0)(ωπ)N1μν+f(0)(ωh)N4μνq0-ωπ+ωh+iηϵ(q0)+f(0)(ωπ)N2μν-f(0)(ωh)N3μνq0+ωπ-ωh+iηϵ(q0)+-f(0)(ωπ)N2μν-(1+f(0)(ωh))N4μνq0+ωπ+ωh+iηϵ(q0)(1+f(0)(ωπ))N1μν+f(0)(ωh)N3μνq0-ωπ-ωh+iηε(q0)],
where f(0)(ω)=1/(e(ω-μπ)/T-1) is the Bose distribution function with arguments ωπ=k2+mπ2 and ωh=(q-k)2+mh2. Note that this expression is a generalized form for the in-medium self-energy obtained by Weldon [38]. The subscript i(=1,…,4) on Nμν in (116) corresponds to its values for k0=ωπ, -ωπ, q0-ωh, and q0+ωh, respectively. It is easy to read off the real and imaginary parts from (116). The angular integration can be carried out using the δ-functions in each of the four terms in the imaginary part which define the kinematically allowed regions in q0 and q where scattering, decay, and regeneration processes occur in the medium leading to the loss or gain of ρ mesons [37]. The vector mesons ω, h1, and a1 which appear in the loop have negative G-parity and have substantial 3π and ρπ decay widths [39]. The (polarization averaged) self-energies containing these unstable particles in the loop graphs have thus been folded with their spectral functions as follows:
(117)Π(q,mh)=1Nh∫(mh-2Γh)2(mh+2Γh)2‍dM2×1πIm[1M2-mh2+iMΓh(M)]Π(q,M)
with Nh=∫(mh-2Γh)2(mh+2Γh)2‍dM2(1/π)Im[1/(M2-mh2+iMΓh(M))]. The contributions from the loops with heavy mesons may then be considered as a multipion contribution to the ρ self-energy.The in-medium cross-section is now obtained by using the fullρ-propagator (113) in place of the usual vacuum propagator Dμν(0) in the scattering amplitudes. The long dashed line in Figure 3 shows a suppression of the peak when only the ππ loop is considered. This effect is magnified when the πh loops (solid line indicated by multipion) are taken into account and is also accompanied by a small shift in the peak position. Extension to the case of finite baryon density can be done using the spectral function computed in [40] where an extensive list of baryon (and antibaryon) loops are considered along with the mesons. A similar modification of the ππ cross-section for a hot and dense system was seen also in [41].We plotη versus T in Figure 5 obtained in the Chapman-Enskog approximation showing the effect of the in-medium ρ propagation in the pion gas [25]. We observe ~10% change at T=150MeV due to medium effects compared to the vacuum when all the loops in the ρ self-energy are considered. The effect reduces with temperature to less than 5% at 100 MeV.Figure 5
The shear viscosity as a function of temperature in the Chapman-Enskog approximation. The dash-dotted line indicates use of the vacuum cross-section, and the dashed and solid lines correspond to in-medium cross-section for theππ and multipion cases respectively.We noted in Section2 that the lowest order result for η in the response theory framework coincides with that obtained in the relaxation time approximation which is in fact the simplest way to linearize the transport equation. Here one assumes that f(x,p) goes over to the equilibrium distribution f(0)(x,p) as a result of collisions, and this takes place over a relaxation time τ(p) which is the inverse of the collision frequency defined in (84). The right hand side of (85) is then given by -Ep[f(x,p)-f(0)(x,p)]/τ(p) which subsequently leads to the expressions (82) and (83) for the shear and bulk viscosities [42]. In Figure 6 we show the temperature dependence of η in the relaxation time approximation. The values in this case are lower than that obtained in the Chapman-Enskog method, though the effect of the medium is larger. In addition to the fact that the expressions for the viscosities are quite different in two approaches, the difference in the numerical values obtained in the two cases also depends significantly on the energy dependence of the ππ cross-section [43].Figure 6
The shear viscosity as a function of temperature in the relaxation time approximation. The dash-dotted and solid lines correspond to the use of in-medium cross-sections in (84) for ππ and multipion loops, respectively. The dashed line represents the vacuum case.In Figure7 we show the numerical results for the bulk viscosity of a pion gas as function of T. It is seen from an analysis of the left hand side of the transport equation that, while the shear viscosity depends on elastic processes, bulk viscosity is sensitive to number changing processes. However, in heavy ion collision experiments matter is known to undergo early chemical freezeout. Number changing (inelastic) processes having much larger relaxation times go out of equilibrium at this point and a temperature dependent chemical potential results for each species so as to conserve the number corresponding to the measured particle ratios. We hence use a temperature dependent pion chemical potential taken from [44] in this case. It is interesting to observe that ζ decreases with T in contrast to η which increases. The trend followed by ζ is similar to the findings of [45]. Additional discussions concerning the temperature dependence of viscosities for a chemically frozen pion gas are available in [31].Figure 7
The bulk viscosity as a function of temperature for a chemically frozen pion gas. The dashed and solid lines correspond to the use of in-medium cross-sections in (84) for ππ and multipion loops, respectively. The dotted line represents the vacuum case.


## 3.1. Transport Coefficients at First Chapman-Enskog Order

The evolution of the phase space distribution of the pions is governed by the (transport) equation:(85)pμ∂μf(x,p)=C[f],
where C[f] is the collision integral. For binary elastic collisions p+k→p′+k′ which we consider, this is given by [22]
(86)C[f]=∫‍dΓkdΓp′dΓk′×[f(x,p′)f(x,k′){1+f(x,p)}×{1+f(x,k)}-f(x,p)f(x,k)×{1+f(x,p′)}{1+f(x,k′)}]W,
where the interaction rate is as follows:
(87)W=s2dσdΩ(2π)6δ4(p+k-p′-k′),
and dΓq=d3q/(2π)3q0. The 1/2 factor comes from the indistinguishability of the initial state pions.For small deviation from local equilibrium, we write, in the first Chapman-Enskog approximation,(88)f(x,p)=f(0)(x,p)+δf(x,p),δf(x,p)=f(0)(x,p)[1+f(0)(x,p)]ϕ(x,p),
where the equilibrium distribution function (in the new notation) is given by
(89)f(0)(x,p)=[e(pμuμ(x)-μπ(x))/T(x)-1]-1,
with T(x), uμ(x), and μπ(x) representing the local temperature, flow velocity, and pion chemical potential, respectively. Putting (88) in (85), the deviation function ϕ(x,p) is seen to satisfy
(90)pμ∂μf(0)(x,p)=-ℒ[ϕ],
where the linearized collision term
(91)ℒ[ϕ]=f(0)(x,p)∫‍dΓkdΓp′dΓk′f(0)(x,k)×{1+f(0)(x,p′)}{1+f(0)(x,k′)}×[ϕ(x,p)+ϕ(x,k)-ϕ(x,p′)-ϕ(x,k′)]W.
Using the form of f(0)(x,p) as given in (89) on the left side of (90) and eliminating time derivatives with the help of equilibrium thermodynamic laws, we arrive at [31]
(92)[Q∂νuν+pμΔμν(p·u-h)(T-1∂νT-Duν)-〈pμpν〉〈∂μuν〉]f(0)(1+f(0))=-Tℒ[ϕ],
where D=uμ∂μ, ∇μ=Δμν∂ν, Δμν=gμν-uμuν, and 〈·〉 indicates a space-like symmetric and traceless combination. In this equation,
(93)Q=-13mπ2+(p·u)2{43-γ′}+p·u{(γ′′-1)h-γ′′′T},
where(94)γ′=(S20/S21)2-(S30/S21)2+4z-1S20S31/(S21)2+z-1S30/S21(S20/S21)2-(S30/S21)2+3z-1S20S31/(S21)2+2z-1S30/S21-z-2,γ′′=1+z-2(S20/S21)2-(S30/S21)2+3z-1S20S31/(S21)2+2z-1S30/S21-z-2,γ′′′=S20/S21+5z-1S31/S21-S30S31/(S21)2(S20/S21)2-(S30/S21)2+3z-1S20S31/(S21)2+2z-1S30/S21-z-2,with z=mπ/T and h=mπS31/S21. The functions Snα(z) are integrals over Bose functions [31] and are defined as Snα(z)=∑k=1∞‍ekμπ/Tk-αKn(kz), Kn(x) denoting the modified Bessel function of order n. The left hand side of (90) is thus expressed in terms of thermodynamic forces with different tensorial ranks. In order to be a solution of this equation ϕ must also be a linear combination of the corresponding thermodynamic forces. It is typical to take ϕ as
(95)ϕ=A∂·u+Bμ∇μν(T-1∂νT-Duν)-Cμν〈∂μuν〉
which on substitution into (92) and comparing coefficients of the (independent) thermodynamic forces on both sides yields the set of equations:
(96)ℒ[A]=-Qf(0)(p){1+f(0)(p)}T,ℒ[Cμν]=-〈pμpν〉f(0)(p){1+f(0)(p)}T
ignoring the equation for Bμ which is related to thermal conductivity. These integral equations are to be solved to get the coefficients A and Cμν. It now remains to link these to the viscous coefficients ζ and η. This is achieved by means of the dissipative part of the energy-momentum tensor resulting from the use of the nonequilibrium distribution function (88) in(97)Tμν=∫‍dΓppμpνf(p)=Tμν(0)+ΔTμν,
where
(98)ΔTμν=∫‍dΓpf(0)(1+f(0))Cαβ〈pαpβ〉〈∂μuν〉+∫‍dΓpf(0)(1+f(0))QAΔμν∂σuσ.
Again, for a small deviation ϕ(x,p), close to equilibrium, so that only first-order derivatives contribute, and the dissipative tensor can be generally expressed in the form [32, 33]
(99)ΔTμν=-2η〈∂μuν〉-ζΔμν∂σuσ.Comparing, we obtain the expressions of shear and bulk viscosity as follows:(100)η=-110∫‍dΓpCμν〈pμpν〉f(0)(p){1+f(0)(p)},ζ=-∫‍dΓpQAf(0)(p){1+f(0)(p)}.
The coefficients A and Cμν are perturbatively obtained from (96) by expanding in terms of orthogonal polynomials which reduces the integral equations to algebraic ones. After a tedious calculation using the Laguerre polynomial of 1/2 integral order, the first approximation to the shear and bulk viscosity comes out as
(101)η=T10γ02c00,ζ=Tα22a22,
where
(102)γ0=-10S32(z)S21(z),c00=16{I1(z)+I2(z)+13I3(z)},α2=z32[13(S30S21-z-1)+(S20S21+3zS31S21){(1-γ′′)S31S21+γ′′′z-1}-(43-γ′){S30S21+15z-2S32S21+2z-1}],a22=2z2I3(z).
The integrals Iα(z) are given by
(103)Iα(z)=z4[S21(z)]2e(-2μπ/T)×∫0∞‍dψcosh3ψsinh7ψ×∫0π‍dΘsinΘ12dσdΩ(ψ,Θ)∫02π‍dϕ×∫0∞‍dχsinh2αχ∫0π‍dθsinθ×e2zcoshψcoshχ(eE-1)(eF-1)(eG-1)(eH-1)Mα(θ,Θ),
where μπ is the chemical potential of pions. The exponents in the Bose functions are given by
(104)E=z(coshψcoshχ-sinhψsinhχcosθ)-μπT,F=z(coshψcoshχ-sinhψsinhχcosθ′)-μπT,G=E+2zsinhψsinhχcosθ,H=F+2zsinhψsinhχcosθ′,and the functionsMα(θ,Θ) represent
(105)M1(θ,Θ)=1-cos2Θ,M2(θ,Θ)=cos2θ+cos2θ′-2cosθcosθ′cosΘ,M3(θ,Θ)=[cos2θ-cos2θ′]2.
The relative angle θ′ is defined by cosθ′=cosθcosΘ-sinθsinΘcosϕ.Note that the differential cross-section which appears in the denominator is the dynamical input in the expressions forη and ζ. It is this quantity that we turn to in the next section.


## 3.2. Theππ Cross-Section with Medium Effects

The strong interaction dynamics of the pions enters the collision integrals through the cross-section. In Figure3, we show the ππ cross-section as a function of the centre of mass energy of scattering. The different curves are explained in the following. The filled squares referred to as experiment is a widely used resonance saturation parametrization [23, 34] of isoscalar and isovector phase shifts obtained from various empirical data involving the ππ system. The isospin averaged differential cross-section is given by
(106)dσ(s)dΩ=4qcm2[19sin2δ00+59sin2δ02+13·9sin2δ11cos2θ],
where
(107)δ00=π2+arctan(E-mσΓσ/2),δ11=π2+arctan(E-mρΓρ/2),δ02=-0.12pmπ.
The widths are given by Γσ=2.06p and Γρ=0.095p((p/mπ)/(1+(p/mρ)2))2 with mσ=5.8mπ and mρ=5.53mπ.Figure 3
Theππ cross-section as a function of centre of mass energy. The dotted line indicates the cross-section obtained using (109) which agrees well with the experimental values (106) shown by filled squares. The dashed and solid lines depict the in-medium cross-section for ππ and πh(h=π,ω,h1,a1) loops, respectively, in the ρ self-energy evaluated at T= 160 MeV.To get a handle on the dynamics, we now evaluate theππ cross-section involving ρ and σ meson exchange processes using the interaction Lagrangian as follows:
(108)ℒ=gρρ→μ·π→×∂μπ→+12gσmσπ→·π→σ,
where gρ=6.05 and gσ=2.5. In the matrix elements corresponding to s-channel ρ and σ exchange diagrams which appear for total isospin I=1 and 0, respectively, we introduce a decay width in the corresponding propagator. We get [25]
(109)ℳI=0=gσ2mσ2[3s-mσ2+imσΓσ+1t-mσ2+1u-mσ2]+2gρ2[s-ut-mρ2+s-tu-mρ2],ℳI=1=gσ2mσ2[1t-mσ2-1u-mσ2]+gρ2[2(t-u)s-mρ2+imρΓρ+t-su-mρ2-u-st-mρ2],ℳI=2=gσ2mσ2[1t-mσ2+1u-mσ2]+gρ2[u-st-mρ2+t-su-mρ2].
The differential cross-section is then obtained from dσ/dΩ=|ℳ|2¯/64π2s, where the isospin averaged amplitude is given by |ℳ|2¯=(1/9)∑‍(2I+1)|ℳI|2¯.The integrated cross-section, after ignoring theI=2 contribution is shown by the dotted line (indicated by “vacuum”) in Figure 3 and is seen to agree reasonably well with the experimental cross-section up to a centre of mass energy of about 1 GeV beyond which the theoretical estimate gives higher values. We hence use the experimental cross-section beyond this energy.After this normalisation to data, we now turn to the in-medium cross-section by introducing the effective propagator for theρ in the previous expressions for the matrix elements. This is obtained in terms of the self-energy by solving the Dyson equation and is given by
(110)Dμν=Dμν(0)+Dμσ(0)ΠσλDλν,
where Dμν(0) is the vacuum propagator for the ρ meson and Πσλ is the self-energy function obtained from one-loop diagrams shown in Figure 4. The standard procedure [35] to solve this equation in the medium is to decompose the self-energy into transverse and longitudinal components. For the case at hand the difference between these components is found to be small and is hence ignored. We work with the polarization averaged self-energy function defined as
(111)Π=13(2ΠT+q2ΠL),
where
(112)ΠT=-12(Πμμ+q2q-2Π00),ΠL=1q-2Π00,Π00≡uμuνΠμν.
The in-medium propagator is then written as
(113)Dμν(q0,q→)=-gμν+qμqν/q2q2-mρ2-ReΠ(q0,q→)+iImΠ(q0,q→).
The scattering, decay, and regeneration processes which cause a gain or loss of ρ mesons in the medium are responsible for the imaginary part of its self-energy. The real part on the other hand modifies the position of the pole of the spectral function.Figure 4
π
h self-energy diagrams where h stands for π, ω, h1, and a1 mesons.As discussed in Section2.1, in the real-time formulation of thermal field theory the self-energy assumes a 2×2 matrix structure of which the 11 component is given by
(114)Πμν11(q)=i∫d4k(2π)4Nμν(q,k)Dπ11(k)Dh11(q-k),
where D11 is the 11 component of the scalar propagator given by D11(k)=Δ(k)+2πif(0)(k)δ(k2-m2). It turns out that the self-energy function mentioned above can be obtained in terms of the 11 component through the relations [35, 36]
(115)ReΠμν=ReΠμν11,ImΠμν=ϵ(q0)tanh(βq02)ImΠμν11.
Tensor structures associated with the two vertices and the vector propagator are included in Nμν and are available in [37] where the interactions were taken from chiral perturbation theory. It is easy to perform the integral over k0 using suitable contours to obtain
(116)Πμν(q0,q)=∫d3k(2π)314ωπωh×[(1+f(0)(ωπ))N1μν+f(0)(ωh)N3μνq0-ωπ-ωh+iηϵ(q0)+-f(0)(ωπ)N1μν+f(0)(ωh)N4μνq0-ωπ+ωh+iηϵ(q0)+f(0)(ωπ)N2μν-f(0)(ωh)N3μνq0+ωπ-ωh+iηϵ(q0)+-f(0)(ωπ)N2μν-(1+f(0)(ωh))N4μνq0+ωπ+ωh+iηϵ(q0)(1+f(0)(ωπ))N1μν+f(0)(ωh)N3μνq0-ωπ-ωh+iηε(q0)],
where f(0)(ω)=1/(e(ω-μπ)/T-1) is the Bose distribution function with arguments ωπ=k2+mπ2 and ωh=(q-k)2+mh2. Note that this expression is a generalized form for the in-medium self-energy obtained by Weldon [38]. The subscript i(=1,…,4) on Nμν in (116) corresponds to its values for k0=ωπ, -ωπ, q0-ωh, and q0+ωh, respectively. It is easy to read off the real and imaginary parts from (116). The angular integration can be carried out using the δ-functions in each of the four terms in the imaginary part which define the kinematically allowed regions in q0 and q where scattering, decay, and regeneration processes occur in the medium leading to the loss or gain of ρ mesons [37]. The vector mesons ω, h1, and a1 which appear in the loop have negative G-parity and have substantial 3π and ρπ decay widths [39]. The (polarization averaged) self-energies containing these unstable particles in the loop graphs have thus been folded with their spectral functions as follows:
(117)Π(q,mh)=1Nh∫(mh-2Γh)2(mh+2Γh)2‍dM2×1πIm[1M2-mh2+iMΓh(M)]Π(q,M)
with Nh=∫(mh-2Γh)2(mh+2Γh)2‍dM2(1/π)Im[1/(M2-mh2+iMΓh(M))]. The contributions from the loops with heavy mesons may then be considered as a multipion contribution to the ρ self-energy.The in-medium cross-section is now obtained by using the fullρ-propagator (113) in place of the usual vacuum propagator Dμν(0) in the scattering amplitudes. The long dashed line in Figure 3 shows a suppression of the peak when only the ππ loop is considered. This effect is magnified when the πh loops (solid line indicated by multipion) are taken into account and is also accompanied by a small shift in the peak position. Extension to the case of finite baryon density can be done using the spectral function computed in [40] where an extensive list of baryon (and antibaryon) loops are considered along with the mesons. A similar modification of the ππ cross-section for a hot and dense system was seen also in [41].We plotη versus T in Figure 5 obtained in the Chapman-Enskog approximation showing the effect of the in-medium ρ propagation in the pion gas [25]. We observe ~10% change at T=150MeV due to medium effects compared to the vacuum when all the loops in the ρ self-energy are considered. The effect reduces with temperature to less than 5% at 100 MeV.Figure 5
The shear viscosity as a function of temperature in the Chapman-Enskog approximation. The dash-dotted line indicates use of the vacuum cross-section, and the dashed and solid lines correspond to in-medium cross-section for theππ and multipion cases respectively.We noted in Section2 that the lowest order result for η in the response theory framework coincides with that obtained in the relaxation time approximation which is in fact the simplest way to linearize the transport equation. Here one assumes that f(x,p) goes over to the equilibrium distribution f(0)(x,p) as a result of collisions, and this takes place over a relaxation time τ(p) which is the inverse of the collision frequency defined in (84). The right hand side of (85) is then given by -Ep[f(x,p)-f(0)(x,p)]/τ(p) which subsequently leads to the expressions (82) and (83) for the shear and bulk viscosities [42]. In Figure 6 we show the temperature dependence of η in the relaxation time approximation. The values in this case are lower than that obtained in the Chapman-Enskog method, though the effect of the medium is larger. In addition to the fact that the expressions for the viscosities are quite different in two approaches, the difference in the numerical values obtained in the two cases also depends significantly on the energy dependence of the ππ cross-section [43].Figure 6
The shear viscosity as a function of temperature in the relaxation time approximation. The dash-dotted and solid lines correspond to the use of in-medium cross-sections in (84) for ππ and multipion loops, respectively. The dashed line represents the vacuum case.In Figure7 we show the numerical results for the bulk viscosity of a pion gas as function of T. It is seen from an analysis of the left hand side of the transport equation that, while the shear viscosity depends on elastic processes, bulk viscosity is sensitive to number changing processes. However, in heavy ion collision experiments matter is known to undergo early chemical freezeout. Number changing (inelastic) processes having much larger relaxation times go out of equilibrium at this point and a temperature dependent chemical potential results for each species so as to conserve the number corresponding to the measured particle ratios. We hence use a temperature dependent pion chemical potential taken from [44] in this case. It is interesting to observe that ζ decreases with T in contrast to η which increases. The trend followed by ζ is similar to the findings of [45]. Additional discussions concerning the temperature dependence of viscosities for a chemically frozen pion gas are available in [31].Figure 7
The bulk viscosity as a function of temperature for a chemically frozen pion gas. The dashed and solid lines correspond to the use of in-medium cross-sections in (84) for ππ and multipion loops, respectively. The dotted line represents the vacuum case.


## 4. Summary and Conclusion

To summarize, we have calculated the shear viscosity coefficient of a pion gas in the real-time version of thermal field theory. It is simpler to the imaginary version in that we do not have to continue to imaginary time at any stage of the calculation. As an element in the theory of linear response, a transport coefficient is defined in terms of a retarded thermal two-point function of the components of the energy-momentum tensor. We derive Källen-Lehmann representation for any (bosonic) two-point function of both time-ordered and retarded types to get the relation between them. Once this relation is obtained, we can calculate the retarded function in the Feynman-Dyson framework of the perturbation theory.Clearly the method is not restricted to transport coefficients. Any linear response leads to a retarded two-point function, which can be calculated in this way. Also quadratic response formulae have been derived in the real-time formulation [46].We have also evaluated the viscous coefficients in the kinetic theory approach to leading order in the Chapman-Enskog expansion. Here we have incorporated an in-mediumππ cross-section and found a significant effect in the temperature dependence of the shear viscosity.The viscous coefficients and their temperature dependence could affect the quantitative estimates of signals of heavy ion collisions particularly where hydrodynamic simulations are involved. For example, it has been argued in [47] that corrections to the freeze-out distribution due to bulk viscosity can be significant. As a result the hydrodynamic description of the pT spectra and elliptic flow of hadrons could be improved by including a realistic temperature dependence of the viscous coefficients. Such an evaluation essentially requires the consideration of a multicomponent gas preferably containing nucleonic degrees of freedom, so that extensions to finite baryon chemical potential can be made. Work in this direction is in progress.



---

*Source: 627137-2013-06-24.xml*