Bulk Viscosity of Hot Quark Plasma from Non-Equilibrium Statistical Operator

We provide a discussion of the bulk viscosity of two-flavor quark plasma, described by the Nambu--Jona-Lasinio model, within the framework of Kubo-Zubarev formalism. This discussion, which is complementary to our earlier study, contains a new, detailed derivation of the bulk viscosity in the case of multiple conserved charges. We also provide some numerical details of the computation of the bulk viscosity close to the Mott transition line, where the dissipation is dominated by decays of mesons into quarks and their inverse processes. We close with a summary of our current understanding of this quantity, which stresses the importance of loop resummation for obtaining the qualitatively correct result near the Mott line.


Introduction
Transport coefficients of hot and dense quark plasma are key inputs in the hydrodynamical description of the heavy-ion experiments at Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC). The matter created in these experiments exhibits a very small ratio of the shear viscosity to the entropy density, which is close to the lower bound placed by the uncertainty principle [1] and conjectured on the basis of AdS/CFT duality [2].
The bulk viscosity describes the dissipation in cases where pressure falls out of its equilibrium value on uniform expansion or contraction of fluid. It vanishes in several cases, e.g., for an ultrarelativistic or nonrelativistic gas interacting weakly with local forces via binary collisions, as well as in strongly coupled systems with conformal symmetry. Because at high energies QCD is almost conformally symmetric, the bulk viscosity of quark-gluon plasma is small in the perturbative regime [3][4][5][6]. At low energies, the conformal symmetry is broken by the quark mass and/or by dimensionful regularization of the ultraviolet divergences, in which case the bulk viscous effects can become important. Large values of bulk viscosity were found, in particular, close to the chiral phase transition line [7,8]. Computations of the QCD bulk viscosity in the strongly coupled regime where carried out using various methods including lattice simulations [5,9,10], quasiparticle Boltzmann transport [11][12][13][14][15] and the Kubo formalism [16][17][18].
The focus of this contribution, which is complementary to our earlier study [16], is the bulk viscosity of quark matter in the non-perturbative regime as it is realized close to the chiral phase transition line. The case of bulk viscosity is special because it requires a resummation of an infinite series of loop diagrams, whereas the remaining coefficient (shear viscosity, thermal and electrical conductivities) are given by the one-loop result only, see Ref. [19]. Specifically, our aim here is to provide details of the computation of this quantity which complement our earlier publication [16]. First, we provide a formal derivation of the bulk viscosity coefficient from the Zubarev formalism of

Bulk Viscosity Formula from the Non-Equilibrium Statistical Operator
The coefficient of the bulk viscosity was computed within the NESO method for the case of a system without conserved charges in the seminal paper by Hosoya et al. [23]. Our purpose here is to extend that derivation to the case of systems with multiple conserved charges.
Hydrodynamics of relativistic quantum fluids is described by the energy-momentum tensor T µν (x) and currents of conserved chargesN µ a (x). We consider the general case of multiple conserved charge flavors (e.g., baryonic, electric, etc.) which are labeled by the index a. In this case, these conservation laws take the form For systems in the hydrodynamic regime, one can introduce local thermodynamic variables, such as temperature T(x) ≡ β −1 (x), chemical potentials µ a (x) and fluid 4-velocity u ν (x) as smooth functions of the space-time coordinates x ≡ (x, t). Below we will specify the (matching) conditions which are necessary to identify these quantities. Our next step is to define a NESO aŝ with operators defined asÂ where the covariant quantities (c-numbers) are defined as Note that the limit ε → +0 in Equation (4) should be taken only after the thermodynamic limit. The proper order of taking the thermodynamic and ε → +0 limits guarantees that the NESO in Equation (2) satisfies the Liouville equation with an infinitesimal source term, which breaks the time-reversibility of that equation and chooses its retarded solution for positive values of ε [20,21,23]. Equations (2)-(6) generalize the analogous expressions of Refs. [23,24] to the case of a system with multiple conserved charges.
The first term in the exponent in Equation (2) corresponds to the local equilibrium part of the statistical operator, defined asρ The second term in the exponent of Equation (2) is the non-equilibrium part of the statistical operator, which can be interpreted as a thermodynamic "force". For small deviations from equilibrium, it can be treated as a small perturbation. Expanding the NESO around the local equilibrium distribution and keeping linear in the operatorB terms we find The statistical average of any operatorX(x) can be written now as where X (x) l is the local equilibrium average and a two-point correlation function has been defined as [23,24] X (x),Ĉ( The final point of our general discussion of the NESO method is the procedure by which the quantities β ν and α a are matched with the relevant thermodynamic variables in an arbitrary non-equilibrium state. This can be achieved by the following matching conditions [20,21,24] where the operators of the energy and charge densities are defined asˆ = u µ u νT µν andn a = u µN µ a . In these expressions, the fluid 4-velocity u µ , which is normalized to unity u µ u µ = 1, should be "tied" to a physical current. This could be either the energy flow, which specifies the Landau-Lifshitz frame [25] or the charge flow, which specifies the Eckart frame [26]. In the Landau frame u µ T µν = ˆ u ν , whereas in the Eckart frame N µ a = n a u µ . The conditions (11) define the temperature and the chemical potentials of components as non-local functionals of ˆ (x) and n a (x) [27]. However, the hydrodynamic description requires thermodynamic parameters as local functions of the energy and charge densities. In practice, this difficulty is circumvented by dividing the fluid into elements which are in local statistical equilibrium and each of which is independent of the other [28]. In practice, the local equilibrium values ˆ l and n a l in Equation (11) are then evaluated assuming formally constant values of β and µ a , which are identified by matching ˆ l and n a l to the real values of these quantities ˆ and n a at any given point x. In this way one can construct a fictitious local equilibrium state, characterized by the thermodynamic parameters β(x) and µ a (x), such that it reproduces the local values of the energy and charge densities at every point of the space and time.

Decomposition into Different Dissipative Processes
To identify the different dissipative processes, we now exploit the common decompositions of the energy-momentum tensor and the charge currents into the ideal and dissipative partŝ wherep is the operator of pressure; ∆ µν = g µν − u µ u ν is the projection operator onto the 3-space orthogonal to u µ and has the properties The dissipative quantitiesπ µν ,q µ andĵ µ a are the operators of the shear stress tensor, energy diffusion flux and charge diffusion fluxes, respectively, and they satisfy the following conditions The operators on the right-hand sides of Equations (12) and (13) can be obtained via the projections ofT µν andN which in the local rest frame [u µ = (1, 0, 0, 0)] read In Equation (17) we introduced a fourth-rank traceless projector orthogonal to u µ The hydrodynamic quantities π µν , q µ and j µ a are obtained as the statistical averages of the corresponding operators over the NESO according to Equations (9) and (10). In local equilibrium the averages of the dissipative operators vanish: To compute the non-equilibrium averages of these operators it is convenient to write the operator C given by Equation (5) as a sum of contributions of different dissipative processes according to Equations (12) and (13). Similar decompositions were performed in Refs. [23,24]. Recalling the properties (14) and (15) we obtain where we introduced the covariant time-derivative D = u ρ ∂ ρ , the covariant spatial derivative ∇ ρ = ∆ ρλ ∂ λ , the fluid expansion rate θ = ∂ ρ u ρ and the velocity stress tensor via σ λρ = ∆ λραβ ∂ α u β .
As a first approximation, we can eliminate the terms Dβ, Dα a and Du λ using the equations of ideal hydrodynamics where p = p( , n a ) is the equilibrium pressure fixed by an equation of state, and h = + p is the enthalpy density. Choosing and n a as independent thermodynamic variables and using the first two equations in (23) we can write In the next step we exploit the thermodynamic relations to obtain Substituting these relations back into Equations (24) and (25) we obtain Now the first three terms in Equation (22) can be combined as followŝ

Kubo Formula for the Bulk Viscosity
By definition, bulk viscous pressure Π measures the deviation of the thermodynamic pressure from its equilibrium value, which results from the expansion or compression of the fluid. Therefore, it might appear at a first glance that the bulk viscous pressure should be identified as p − p l . However, it is easy to see that such a definition would be erroneous. To understand the problem, we go back to the matching conditions (11), which define the local equilibrium state. As explained above, these conditions are satisfied only if the local equilibrium distribution function is evaluated formally assuming uniform background values of the thermodynamic parameters, i.e., as if these were constant in space and time with the given values β(x) and µ a (x). Because the local equilibrium distribution (7) is actually a functional of non-uniform thermodynamic parameters, the average values ˆ l and n a l in the full computation are shifted from the actual values of = ˆ and n a = n a by additional gradient terms ∆ ≡ ˆ − ˆ l and ∆n a ≡ n a − n a l , which were neglected in Equation (11). These shifts bring, in their turn, an additional shift in the equilibrium part of the pressure p l , which should not be included in the bulk viscous pressure [4,29,30]. Thus, the bulk viscous pressure should be defined as the difference between the actual non-equilibrium pressure p and the equilibrium pressure p( , n a ), which is not equal to p l ≡ p( ˆ l , n a l ): where we kept only the linear terms. Then, the bulk viscous pressure is given by where we used the definition ofp * given by Equation (31). From Equations (9), (22) and (30) we then obtain where we dropped the correlators between operators of different rank, because they vanish in isotropic medium according to Curie's theorem [31]. Introducing the bulk viscosity as we rewrite Equation (34) as The correlator (35) can be evaluated using uniform background values of thermodynamic parameters, i.e., as if the system is in global thermal equilibrium. Finally, the bulk viscosity can be cast in the form of a Kubo formula [23,24] where the two-point retarded equilibrium Green's function in the zero-wavenumber limit is given by where the square brackets denote a commutator.

Bulk Viscosity within the Two-Flavor NJL Model
In this section, we illustrate the computation of the bulk viscosity following Ref. [16]; in doing so we will provide some numerical details not exposed earlier. The Lagrangian of the two-flavor NJL model contains scalar-isoscalar and pseudoscalar-isovector channels of interactions among quarks and is given by where ψ = (u, d) T is the Dirac field for u and d quarks, m 0 = 5.5 MeV is the current-quark mass, G = 10.1 GeV −2 is the effective coupling and τ is the vector of Pauli matrices in the space of isospin. The NJL model is regularized with a three-momentum cut-off Λ = 0.65 GeV. Assuming isospin-symmetry, the only conserved current is the net particle current given bŷ The energy-momentum tensor readŝ The relevant operatorp * which enters the Kubo formula (37) with the correlator given by Equation (38) in the local rest frame reads (see Equations (18) and (31)) Inserting Equation (42) where we omitted the arguments of the operators. Substituting here the explicit expressions forT µν andN µ and switching to the imaginary-time (Matsubara) formalism via the substitutions t → −iτ, with two-point correlation functions defined as where ω n = 2πnT, n = 0, ±1, . . . are the bosonic Matsubara frequencies; T τ is the imaginary time-ordering operator; i/ ∂ τ ≡ −γ 0 ∂ τ − iγ j ∂ j , andâ andb are either constants or γ-matrices contracted with partial derivatives. (Note that the correlators which arise from the interaction part of Equation (39) vanish for ω = 0 because of the energy conservation, see Ref. [16] for details). Figure 1 illustrates diagrammatically the series of the loop diagrams which contribute to the correlation function given by Equation (45).

Resummation of the Feynman Diagrams
The class of leading-order diagrams which contributes to the correlation function (45) is identified according to the O(1/N c ) power-counting scheme. In this scheme each diagram is selected according to its power with respect to the color number N c , which is determined by the following rules [22]: (a) each quark loop contributes a factor of N c , which arises from the trace over the color space; (b) each coupling G contributes a factor of 1/N c . It is easy to see that the leading-order diagrams in the correlation function (45) are of the order of O(N 1 c ) and involve loop diagrams without vertex corrections, i.e., those of the type shown in the first and the second lines in Figure 1. Indeed, the factor N c associated with each additional loop is compensated by the factor 1/N c from an interaction insertion. Therefore, we conclude that we need to resume an infinite chain of loop diagrams without vertex corrections. To carry out the resummation, define the single-loop diagram in the momentum space as where D(p, iω l ) is the full (i.e., dressed) quark propagator defined in the imaginary time, and the summation is over the fermionic Matsubara frequencies ω l = π(2l + 1)T − iµ, l = 0, ±1, . . . , where µ is the quark chemical potential. The traces should be taken in Dirac, color, and flavor spaces. The resummation then leads to where Γ = 1 (the correlators with the pseudoscalar vertex ∝ γ 5 vanish because the relevant traces vanish in the Dirac space). The frequency arguments in Equation (47) were omitted for the sake of brevity. The effective coupling in Equation (47) is related to the bare coupling G viã The diagrams involving vertex corrections, such as the one shown in the third line of Figure 1, are of higher order in the O(1/N c ) power-counting scheme. Thus, the computation of the leading-order contribution to the bulk viscosity reduces to the calculation of the series of loop diagrams defined by Equation (47), which in turn requires the evaluation of the single-loop diagram given by Equation (46). To carry out the sum over the Matsubara frequencies in Equation (46) one needs the frequency dependence of the operatorsâ andb which arises whenâ,b ∝ ∂ τ [see Equation (44)]. Indeed, in the frequency space, such dependence translates as ∂ τ → −iω l ≡ −i(ω l + ω n /2). For such cases, we separate the iω l -dependence by formally factorizing the frequency dependence into a function f (iω l ), i.e., we writeâ...b... = f (iω l )â 0 ...b 0 ..., whereâ 0 andb 0 do not depend on iω l . After summation over the Matsubara frequencies and subsequent analytical continuation (i.e., iω n = ω + iδ) we obtain (see Appendix A for details) whereñ(ε) = n(ε) − 1/2 with n(ε) = [e β(ε−µ) + 1] −1 being the Fermi distribution for quarks. Finally, we separate the real and imaginary parts in Equation (49) by exploiting the Dirac identity From Equations (49) and (50) we find and Using Equation (52) we can compute the imaginary part of Equation (47) where To compute the traces in Equations (49) and (51) one needs to exploit the Dirac decomposition of the spectral function where m is the quark mass. The coefficients A s , A 0 and A v can be expressed in terms of the relevant components of the quark self-energy according to the relations [19,32,33] where i = ImΣ i , r i = ReΣ i , i = s, 0, v, and From now on we will neglect the irrelevant real parts of the self-energy, which lead to momentum-dependent corrections to the constituent quark mass in next-to-leading order O(N −1 c ) and will keep only the imaginary parts which were computed in Refs. [19,32,33]. The three components of the quark self-energy are identified via With this input, we can now calculate the relevant correlators entering Equation (44). Computing the traces and performing the angular integrations in Equations (49) and (51), we find, e.g., where A i ≡ A i (p, ε ); N c = 3 and N f = 2 are the color and flavor numbers, respectively. The remaining correlation functions can be computed in analogy to Equations (64) and (65); the explicit expressions are given in Ref. [16].
Inserting the relevant correlators in Equations (37), (44), (53)-(58), we can write the bulk viscosity as where each of the three terms arises from the corresponding terms in Equation (53). The first (single-loop) term is given by The multiloop contributions are given by where the effective couplingḠ is given byḠ with the polarization loop and Here the functions a , b , c , y are obtained from a, b, c, y defined above by substitution ε → ε . Equations (66)-(73) express the bulk viscosity of quark plasma in terms of the components of its spectral function. It is remarkable that the multiloop contributions to the bulk viscosity vanish trivially in the chirally symmetric case, where m 0 = 0. Indeed, for massless quarks and temperatures T > T c , where T c is the critical temperature of the chiral phase transition, we find x = 0 and a, a ∝ m = 0 in Equation (70). As a result, we have ζ 1,2 = 0 and ζ = ζ 0 according to Equations (68) and (66).

Numerical Results
We use the Lorentz components of the quark spectral function, obtained previously in Refs. [19,32,33], to evaluate numerically the bulk viscosity. We concentrate on the region of the phase diagram which is located above the Mott transition temperature T M , which is defined as the threshold temperature at any given chemical potential above which the meson decay into two on-mass-shell quarks is kinematically allowed. It is identical to the chiral phase transition temperature T c in the chiral limit m 0 = 0.
To gain insight into the numerical results for the bulk viscosity, it is useful first to analyze the integrands of Equations (67) and (70)-(73). The integrands of ζ 0 , I 1 and I 2 are shown in Figure 2. Each of these integrands develops a peak structure at p |ε|, whereby the height of the peak rapidly increases with |ε|. The momentum integrals are increasing functions of |ε| as long as |ε| ≤ Λ, and decreasing for larger energies (because of the momentum cut-off p ≤ Λ, see Figure 2). Note that for µ = 0 the integrands are even functions of ε, which reflects the quark-antiquark symmetry. The factor ∂n(ε)/∂ε breaks this symmetry for non-vanishing chemical potentials by increasing the contribution of quarks. Thus, we see that the dominant contribution to the bulk viscosity comes from the modes with p |ε|, whereby the quark contribution dominates the antiquark contribution at non-zero µ.
Now we turn to the three-dimensional integrals R 0 andR given by Equations (70) and (73). Their integrands are strongly peaked at p |ε| |ε |, and have two smaller maxima located at p |ε| and p |ε | in the cases where |ε| = |ε |, see Figure 3. As a result, the momentum integrals of these expressions obtain the main contribution form energies ε ±ε. We also observe that the height of each peak increases with |ε| for |ε| ≤ Λ and sharply drops beyond the cut-off. The peak structures seen above reflect the quasiparticle-like nature of the excitations, which however have non-zero width because of the meson decay and recombination processes included in our consideration.    The results for the bulk viscosity are shown in Figure 5. The multiloop contributions ζ 1 and ζ 2 dominate over the one-loop contribution ζ 0 in the regime close to the Mott transition line, where all three components rapidly decrease with the temperature and density. At high enough temperatures, the one-loop contribution scales as T 3 and dominates over the multiloop contributions. The net bulk viscosity which is the sum of the one-loop and multiloop contributions then exhibits a shallow minimum as a function of temperature. From the analysis above, we thus conclude that the single-loop approximation is justified only at sufficiently high temperatures where multiloop diagrams are suppressed. We now comment briefly on the case where the chiral symmetry is intact, i.e., m 0 = 0. In this case, the multiloop contributions vanish automatically, as explained in Section 3. The bulk viscosity is then given by the single-loop contribution ζ 0 taken in the limit m → 0, which is shown in Figure 5 by the dotted lines for zero and finite chemical potentials. Contrary to the case where m 0 = 0, here ζ 0 is smooth at the Mott temperature and increases with the temperature ∝ T 3 in the whole temperature-density range. We thus conclude that the explicit chiral symmetry breaking is essential for the correct description of the bulk viscosity in the low-temperature region of the phase diagram, especially in the region close to the chiral phase transition line.
For completeness, we also compare our results with the shear viscosity η, which was computed previously in Ref. [19] (see also Refs. [32,33]) by employing the same formalism and approximations. Figure 6 shows the dependence of the ratios ζ/s and η/s, where s is the entropy density, on the temperature for several values of the chemical potential [16,19]. For comparison, we also show the AdS/CFT lower bound 1/4π on the η/s ratio [2]. We see that both ratios decrease rapidly with the temperature, but the slope of this decrease is larger in the case of the bulk viscosity in the region which is close to the Mott transition line. In this regime, the bulk viscosity exceeds the shear viscosity by factors ζ/η 5 ÷ 20. Thus, in the low-temperature regime close to the Mott transition line the bulk viscosity is the dominant source of dissipation. It is worth stressing that had we kept only the one-loop contribution to the bulk viscosity, it would have been negligible compared to the shear viscosity. As ζ drops much faster than η with the temperature, the shear viscosity is the dominant dissipation channel at high temperatures. Our numerical results can be fitted using the formula where y = µ/µ 0 , and µ 0 = 0.345 GeV is the value of the chemical potential at which the Mott line terminates. The coefficients a, b, c, d depend only on the chemical potential and are given by The relative error of the fit formula (74) is ≤ 10% for chemical potentials µ ≤ 0.2 GeV. The bulk viscosity according to our fit formula is shown in Figure 5 by empty circles. The fit formula above should be complemented by a fit to the Mott temperature as a function of the chemical potential, which is given by the formula where T 0 = T M (µ = 0) = 0.213 GeV, and γ = 2.7. The relative accuracy of the formula (79) is ≤ 3% for chemical potentials µ ≤ 0.32 GeV.

Conclusions
In this contribution, we provided a general derivation of the Kubo formula for the bulk viscosity from Zubarev's formalism of NESO generalized to systems with multiple conserved charges. The method was then illustrated on the example of computation of the bulk viscosity of quark matter in the framework of the two-flavor NJL model. The previous discussion of Ref. [16] has been supplemented by further details.
The key finding of our work is that at low temperatures and close to the Mott transition line the overall multiloop contribution to the bulk viscosity is larger than the one-loop contribution. This is in contrast to the results found for the shear viscosity and the thermal and electrical conductivities, for which the single-loop approximation gives the leading-order result [19]. We have shown that the bulk viscosity decreases with the temperature and the chemical potential in this regime, attains a minimum and then increases again at higher temperatures where the one-loop contribution becomes dominant.
Phenomenologically interesting is the fact that the bulk viscosity provides the main source of dissipation of stresses close to the Mott line as it exceeds the shear viscosity in this regime by factors of 5 ÷ 20. The bulk viscosity drops faster than the shear viscosity as the temperature increases and it becomes negligible above a certain value of the temperature. Finally, we observed that in the chiral symmetric case, where the quark masses vanish above the (critical) Mott temperature, the picture is different. In this case, the multiloop contributions to the bulk viscosity vanish, and consequently the bulk viscosity becomes negligible compared to the shear viscosity in the entire temperature-density plane. We show now that the contribution of the large circle to the integral (A3) vanishes. Because of the sum rule ∞ −∞ dεA(p, ε) = const, the quark propagator for large |z| has the scaling D ∝ z −1 . Therefore, for large |z| the integrand in Equation (A3) scales as ∝ñ(z)z k−2 (recall that f ∝ z k ). For the Fermi distribution function we have the asymptotics n(z) → Rez→∞ 0 and n(z) → Rez→−∞ 1. Substituting z = Re iφ , dz = iRe iφ dφ, and performing the limit R → ∞, we can write for the integral along the circle z k−2ñ (z) = R k−1 4π π/2 −π/2 dφe i(k−1)φ − 3π/2 π/2 dφe i(k−1)φ = R k−1 π sin(k−1)π/2 k−1 , k = 1, 0, k = 1.
As seen from Equation (A4), the integral vanishes in the limit R → ∞ if k = 0, 1. If k = 2, we have S R ∝ R → ∞, which is purely real and can be dropped since in this case we need only the imaginary parts of S[ f ] functions. As a result, we need to keep only the integrals along the two branch cuts shown in Figure A1. Then, using also Equation (A2), we obtain for Equation (A3) where we took into account that n(ε − iω n ) = n(ε). Substituting the spectral representation (A1) into Equation (A5), changing the variables ε ↔ ε in the first term and performing analytical continuation via iω n → ω + iδ, we find Substituting this expression into Equation (46), we obtain Equation (49) of the main text.