Viscous coefficients of a hot pion gas

The steps essentially involved in the evaluation of transport coefficients in linear response theory using Kubo formulas are to relate the defining {\em retarded} correlation function to the corresponding {\em 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 {\em 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 $\pi\pi$ cross-section in used in which spectral modifications are introduced in the propagator of the exchanged $\rho$.


I. 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 non-equilibrium 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 non-equilibrium 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 non-equilibrium 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][6][7][8][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 of 2 × 2 matrices. But this difficulty is only apparent: Such matrices are always diagonalisable and it is the 11-component of the diagonalised matrix that plays 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 non-trivial way.
In the literature transport coefficients are evaluated using the imaginary time formulation [10][11][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][14][15][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 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 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 have been used in [20,[22][23][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 Sec. II 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 non-equilibrium 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 Sec. III 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 = π, ω, h 1 , a 1 ) loops evaluated in the real-time formulation discussed above. We end with a summary in Sec. IV.

A. 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 t 1 and t 2 (t 2 > t 1 ) 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 Fig. 1, lying within the region of analyticity in this plane and accommodating real time correlation functions [6,8].
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 shall 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 where O for any operator O denotes equilibrium ensemble average; Note that we have two sums in (2.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, Its spatial Fourier transform is where the times τ, τ ′ are on the contour C. We now insert unity on the left of eq. (2.5) in the form (We reserve k 0 for the variable conjugate to the real time.) Then it may be written as where the spectral function M + is given by In just the same way, we can work out the Fourier transform of with a second spectral function M − is given by The two spectral functions are related by the KMS relation [26,27] in momentum space, which may be obtained simply by interchanging the dummy indices m, n in one of M ± ll ′ (k) and using the energy conserving δ-function.
We next introduce the difference of the two spectral functions, 11) and solve this identity and the KMS relation (2.10) for M ± ll ′ (k), In terms of the true distribution function it may be expressed as With the above ingredients, we can build the spectral representations for the two types of thermal propagators. First consider the time-ordered one, As T → ∞, the contour of Fig. 1 simplifies, reducing essentially to two parallel lines, one the real axis and the other shifted by −iβ/2, points on 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, 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 where the elements of the matrix Λ are given by [16] Using relation (2.15), we may rewrite (2.19) in terms of n, The matrix Λ and hence the propagator D ll ′ can be diagonalised to give where D ll ′ and U are given by Eq. (2.21) shows that D can be obtained from any of the elements of the matrix D, say D 11 . Omitting the indices ll ′ , we get ReD = ReD 11 , ImD = tanh(β|k 0 |/2)ImD 11 . (2.23) Looking back at the spectral functions M ± ll ′ defined by (2.7, 2.9), 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, 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, we see that D becomes the free propagator, D(k 0 , k) = −1/(k 2 − m 2 ). We next consider the retarded thermal propagator where again τ , τ ′ are on the contour C (Fig. 1). Noting eqs. (2.6,2.8,2.11) the three dimensional Fourier transform may immediately be written as 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 whose temporal Fourier transform gives This 11-component suffices for us, but we also display the complete matrix, . (2.30) Though we deal with matrices in real time formulation, it is the 11-component that is physical. Eqs. (2.22) and (2.29) then show that we can continue the time-ordered two-point function into the retarded one by simply changing the iǫ prescription, 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.

B. 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, (2.32) The non-equilibrium density matrix operator is constructed in the Heisenberg picture, where it is independent of time, Following Zubarev, we construct the operator B(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, The construction (2.34), which smooths out the oscillating terms resemble the one used in the formal theory of scattering [28,29] and selects out the retarded solution.
The expression (2.34) is actually independent of t; the time derivative is As T 0ν and F ν are finite, the right hand side of (2.36) goes to zero as ǫ → 0. Also integrating (2.34) by parts, we get We now consider the space integral of (2.37). Using the energy-momentum conservation rule (2.32), we integrate the second term in (2.37) by parts and neglect the surface integrals to get where we have abbreviated the first and second terms by A and −B respectively. Then the non-equilibrium statistical density matrix is given by The first term A in eq. (2.38) characterises local equilibrium, while 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 such that the boundary conditions at τ = 0 and τ = 1 correspond to the equilibrium and non-equilibrium density matrices, which can be integrated to give It can be solved iteratively. Keeping up to the first order term (linear response) and setting τ = 1 we get the required result Applying this formula to the energy-momentum tensor, we get its response to the thermodynamical forces as [10] T µν ( 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 We now use eq. (2.46) to obtain the expectation value of the viscous-shear stress part of the non-equilibrium energy momentum tensor which is given by 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 (2.38) is expanded as with p = p − c 2 s ǫ, c s being the sound velocity. Using now the fact that the correlation function between operators of different ranks vanish in an isotropic medium, one can write from (2.46) with π µν (x, t) 0 = 0. Following Hosoya [10] we write the correlation function as (π µν , π ρσ ) = 1 10 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 Again, starting with the pressure p(x, t) on the l.h.s of eq. (2.46) and following the steps as described above we obtain (2.56) where the bulk viscosity ζ is given in terms of a retarded correlation function by Here p(x, t) = p(x, t) − c 2 s ǫ(x, t) with ǫ(x, t) the energy density and c 2 s = ∂p/∂ǫ 0 . Recall that · · · 0 denotes equilibrium ensemble average. From now on we shall drop the subscript '0' on the correlation functions.

C. Perturbative evaluation
Clearly the spectral forms and their inter-relations derived in Sec. IIA hold also for the two-point function appearing in eq. (2.55) for the shear viscosity. We begin with four-dimensional Fourier transforms. To calculate the 11-element of the the retarded two-point function we consider the corresponding time-ordered one, which can be calculated perturbatively. The viscous stress tensor can be extracted from the energy momentum tensor using the formula denotes the pion triplet. We take the lowest order chiral Lagrangian given by [18] The time-ordered correlator, to leading order, is then given by Wick contractions of pion fields in π αβ which is obtained as In the so-called skeleton expansion, these contractions are expressed in terms of complete propagators (see Fig. 2) to get, where D 11 is given by eq. (2.18) and N (l, q) is determined by the derivatives acting on the pion fields, where the pion isospin degeneracy factor g π = 3. To work out the l 0 integral in eq. (2.63), it is more convenient to use Λ 11 given by eq. (2.19) than by eq. (2.20). Closing the contour in the upper or lower half l 0 -plane we get (2.66) The imaginary part of Π 11 arises from the factor K, while 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 eq. (2.23) to build the 11-element of the diagonalised Π matrix, .

(2.68)
Finally Π can be continued to Π R 11 by a relation similar to eq. (2.30), Note that in eqs. (2.68,2.69) we retain the f (k ′ 0 )f (k ′′ 0 ) terms in the numerator to put it in a more convenient form. Change the signs of k ′ 0 and k ′′ 0 in the first and second term respectively. Noting relations Returning to the expression (2.55) for η, we now get the three-dimensional spatial integral of the retarded correlation function by setting q = 0 in eq. (2.58) and Fourier inverting with respect to 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 Refs. [10,11], which we describe below for completeness. As shown in Ref. [10], the integral over t 1 , t ′ and q 0 in eqs. (2.55) and (2.72) may be carried out trivially to give The q 0 dependence of Π R 11 is contained entirely in W , Changing the integration variables in eq. (2.70) from It turns out that the integral over k 0 becomes undefined, if we try to evaluate F (k 0 ) with the free spectral function ρ 0 (k) given by eq. (2.25). As pointed out in Ref. [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, (2.77) Note that this form of the spectral function trivially follows on replacing iη (where η → 0 + ) with iΓ in the free spectral function (2.25) which can be written as (2.78) Then F (k 0 , l) becomes having double poles at k 0 = 2πin/β for n = ±1, ±2, · · · and also at k 0 = ±ω ± iΓ. The integral over k 0 may now be evaluated by closing the contour in the upper/lower half-plane to get where we retain only the leading (singular) term for small Γ. In this approximation eq. (2.75) gives Proceeding analogously as above, the lowest order contribution to the bulk viscosity can be obtained as [10] The width Γ(l) at different temperatures is known [30] from chiral perturbation theory. The quantity Γ can also be interpreted as the collision frequency, 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]) where dσ dΩ is the ππ cross-section. Note that the lowest order formulae for the shear and bulk viscosities obtained above in the linear response approach coincide with the expressions which result from solving the transport equation in the relaxation-time approximation.

III. 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 non-equilibrium distribution function by solving the transport equation in the hydrodynamic regime by expanding the distribution function in a local equilibrium part along with non-equilibrium 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.

A. Transport coefficients at first Chapman-Enskog order
The evolution of the phase space distribution of the pions is governed by the (transport) equation where C[f ] is the collision integral. For binary elastic collisions p + k → p ′ + k ′ which we consider, this is given by [24] C where the interaction rate, and dΓ q = d 3 q (2π) 3 q0 . 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 where the equilibrium distribution function (in the new notation) is given by with T (x), u µ (x) and µ π (x) representing the local temperature, flow velocity and pion chemical potential respectively. Putting (3.3) in (3.1) the deviation function φ(x, p) is seen to satisfy where the linearized collision term Using the form of f (0) (x, p) as given in (3.4) on the left side of (3.5) and eliminating time derivatives with the help of equilibrium thermodynamic laws we arrive at [31] [ where D = u µ ∂ µ , ∇ µ = ∆ µν ∂ ν , ∆ µν = g µν − u µ u ν and · indicates a space-like symmetric and traceless combination. In this equation with z = m π /T and h = m π S 1 3 /S 1 2 . The functions S α n (z) are integrals over Bose functions [31] and are defined as S α n (z) = ∞ k=1 e kµπ /T k −α K n (kz), K n (x) denoting the modified Bessel function of order n. The left hand side of (3.5) 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 which on substitution into (3.7) and comparing coefficients of the (independent) thermodynamic forces on both sides, yields the set of equations 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 non-equilibrium distribution function (3.3) in Again, for a small deviation φ(x, p), close to equilibrium, so that only first order derivatives contribute, the dissipative tensor can be generally expressed in the form [32,33] Comparing, we obtain the expressions of shear and bulk viscosity, The coefficients A and C µν are perturbatively obtained from (3.13) and (3.14) 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 come out as , c 00 = 16{I 1 (z) + I 2 (z) + 1 3 and The integrals I α (z) are given by where µ π is the chemical potential of pions. The exponents in the Bose functions are given by E = z(cosh ψ cosh χ − sinh ψ sinh χ cos θ) − µ π /T F = z(cosh ψ cosh χ − sinh ψ sinh χ cos θ ′ ) − µ π /T G = E + 2z sinh ψ sinh χ cos θ H = F + 2z sinh ψ sinh χ cos θ ′ , (3.25) and the functions M α (θ, Θ) represent 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 we turn to in the next section.

B. The ππ cross-section with medium effects
The strong interaction dynamics of the pions enters the collision integrals through the cross-section. In Fig. 3 we show the ππ cross-section as a function of the centre of mass energy of scattering. The different curves are explained below. 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 sin 2 δ 0 0 + 5 9 sin 2 δ 2 0 + 1 3 · 9 sin 2 δ 1 1 cos 2 θ (3.27) where 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 π . To get a handle on the dynamics we now evaluate the ππ cross-section involving ρ and σ meson exchange processes using the interaction Lagrangian 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], (3.30) The differential cross-section is then obtained from dσ dΩ = |M| 2 /64π 2 s where the isospin averaged amplitude is given by |M| 2 = 1 9 (2I + 1)|M I | 2 . The integrated cross-section, after ignoring the I = 2 contribution is shown by the dotted line (indicated by 'vacuum') in Fig. 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 above expressions for the matrix elements. This is obtained in terms of the self-energy by solving the Dyson equation and is given by µν is the vacuum propagator for the ρ meson and Π σλ is the self energy function obtained from one-loop diagrams shown in Fig. 4. The standard procedure [35] to solve this equation in the medium is to decompose the selfenergy 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 The in-medium propagator is then written as D µν (q 0 , q) = −g µν + q µ q ν /q 2 q 2 − m 2 ρ − ReΠ(q 0 , q) + iImΠ(q 0 , 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.
As discussed in Sec. IIA, 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 where D 11 is the 11-component of the scalar propagator given by D 11 (k) = ∆(k) + 2πif (0) (k)δ(k 2 − m 2 ). It turns out that the self-energy function mentioned above can be obtained in terms of the 11-component through the relations [35,36] Re Π µν = Re Π 11 µν ImΠ µν = ǫ(q 0 ) tanh(βq 0 /2)Im Π 11 µν . (3.36) 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 k 0 using suitable contours to obtain where f (0) (ω) = 1 e (ω−µπ )/T −1 is the Bose distribution function with arguments ω π = k 2 + m 2 π and ω h = (q − k) 2 + m 2 h . 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 (3.37) correspond to its values for k 0 = ω π , −ω π , q 0 − ω h , q 0 + ω h respectively. It is easy to read off the real and imaginary parts from (3.37). 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 q 0 and q where scattering, decay and regeneration processes occur in the medium leading to the loss or gain of ρ mesons [37]. The vector mesons ω, h 1 and a 1 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, . The contributions from the loops with heavy mesons may then be considered as a multi-pion contribution to the ρ self-energy. The in-medium cross-section is now obtained by using the full ρ-propagator (3.34) in place of the usual vacuum propagator D (0) µν in the scattering amplitudes. The long dashed line in Fig. 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 multi-pion) 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 anti-baryon) 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 Fig. 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 = 150 MeV 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.
We noted in Sec. II 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 (2.83). The right hand side of eq. (3.1) is then given by −E p [f (x, p) − f (0) (x, p)]/τ (p) which subsequently leads to the expressions (2.81) and (2.82) for the shear and bulk viscosities [42]. In Fig. 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].
In Fig. 7 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 freeze-out. 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].

IV. 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 p T 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 multi-component 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.