Critical behavior of the bulk viscosity in QCD

We study the behavior of the bulk viscosity $\zeta$ in QCD near a possible critical endpoint in the phase diagram. We verify the expectation that $(\zeta/s)\sim a(\xi/\xi_0)^{x_\zeta}$, where $s$ is the entropy density, $\xi$ is the correlation length, $\xi_0$ is the non-critical correlation length, $a$ is a constant and $x_\zeta\simeq 3$. Using a recently developed equation of state that includes a critical point in the universality class of the Ising model we estimate the constant of proportionality $a$. We find that $a$ is typically quite small, $a\sim O(10^{-4})$. We observe, however, that the result is sensitive to the commonly made assumption that the Ising temperature axis is approximately aligned with the QCD baryon chemical potential axis. If this is not the case, then the critical $\zeta/s$ can approach the non-critical value of $\eta/s$, where $\eta$ is the shear viscosity, even if the enhancement of the correlation length is modest, $\xi/\xi_0\sim 2$.


I. INTRODUCTION
There are several active programs dedicated to exploring the phase diagram of Quantumchromodynamics (QCD) using collisions of heavy ions [1]. A central feature of the phase diagram is the possible presence of a critical endpoint of a first order phase transition between the hadronic phase and the quark gluon plasma (QGP) phase. Experimentally, a critical point is expected to manifest itself in terms of a non-trivial beam energy, rapidity, or system size dependence of fluctuation observables [2]. In a static system fluctuation observables are controlled by the critical equation of state [3][4][5]. The equation of state near the QCD critical point is expected to be in the liquid-gas (Ising) universality class. The equation of state of the Ising model is known from lattice simulations [6], and accurate parametrizations are available in the literature [7]. More recently, there have been efforts to map the Ising equation of state onto the QCD phase diagram, taking into account information from lattice QCD about the equation of state and the susceptibilities at zero baryon chemical potential [8,9].
The dynamic behavior of fluctuations is expected to be governed by model H in the classification of dynamical critical phenomena by Hohenberg and Halperin [10][11][12]. Model H is a hydrodynamic theory that describes the coupling of the order parameter field to a conserved momentum density. It predicts the dynamic critical exponent for the relaxation of the order parameter, and the critical behavior of the transport coefficients, the shear viscosity η, the bulk viscosity ζ, and the thermal conductivity κ. In model H, fluctuations of the order parameter with wave number q ∼ ξ −1 , where ξ is the correlation length, relax on a time scale τ ∼ ξ z , where z 3. This behavior is intermediate between ordinary diffusion (z 2), and the relaxation of a conserved charge not coupled to fluctuations of the fluid velocity (z 4). There is a very mild divergence in the shear viscosity, and more pronounced critical behavior in the thermal conductivity and bulk viscosity [10,[13][14][15][16][17], Physical effects related to critical transport phenomena have been observed in ordinary fluids. For example, the critical behavior of the bulk viscosity manifests itself in sound attenuation near the critical endpoint [17].
Recently, a number of authors have investigated the dynamic evolution of fluctuations in an expanding QCD medium. This includes studies of non-critical correlation functions [18][19][20], simulations of critical stochastic diffusion in an expanding medium [21], and deterministic frameworks for the evolution of two-point [22][23][24] or higher n-point functions [25] near a critical point.
There is a general expectation that the large critical exponent in equ. (1), combined with the strong deviation of the QCD equation of state from scale invariance, will lead to a substantial enhancement of the bulk viscosity, and to large effects on the evolution of a heavy ion collision near a critical end point [26,27]. Our goal in the present work is to study this problem more quantitatively, based on the equation of state constructed by Parotto et al. [9]. We will verify the expected scaling behavior, and estimate the overall coefficient.
We will also study the relaxation of the bulk pressure near a QCD critical point. These results complement earlier studies of non-critical contributions to the bulk viscosity in QCD [28][29][30][31].

II. FLUCTUATIONS OF THE ENERGY DENSITY AND PRESSURE
Fluctuations of the order parameter are governed by an entropy functional S = d 3 x s.
The entropy functional of the Ising model is a function of the densities where and ψ are the Ising energy density and order parameter. The corresponding intensive variables are the reduced temperature r and the magnetic field h, We are following the notation of Landau and Lifshitz [32] as well as Akamatsu et al. [23].
The analogous canonical pair in QCD is where (e, n) are the energy and baryon density, β = 1/T is the inverse temperature and µ is baryon chemical potential. We will assume that there is a map between the intensive variables in QCD and the corresponding Ising variables, see Fig. 1 and Sect. III. Fluctuations of the QCD pressure are given by Thermodynamic variables are defined in the text. The mapsR b A and R A b are defined in equ. (25,26).
The map between (r, h) and (−β, βµ) induces a map between the QCD densities (e, n) and the Ising densities ( , ψ). We assume that the singular part of the QCD entropy density is proportional to the Ising entropy density, s sing (e, n) = As Is ( (e, n), ψ(e, n)). A common assumption is that the images of the r and h axes in the QCD phase diagram are approximately orthogonal, and that the Ising temperature axis is almost aligned with the QCD chemical potential axis 1 , ∂ /∂e 0 [8,9,33]. This implies that Onuki observed that the Ising entropy functional contains a tri-linear coupling between and ψ 2 , and as a result fluctuations of the pressure couple to ψ 2 [16,17]. This means that correlation functions of δP are controlled by the order parameter relaxation rate, and the slow relaxation of ψ leads to a large enhancement in the bulk viscosity. 1 We will reconsider this assumption in Sect. VI.
The presence of a tri-linear coupling between and ψ 2 in the Ising entropy S[ψ, ] can be understood from the relation between entropy and Gibbs free energy 2 [34] exp where E is the total energy, 0 is the background energy density, and T is the temperature.
For small and ψ we can expand the entropy functional as where k, v, u and γ are constants. C 0 is the heat capacity, and S 0 and T 0 are the background entropy and temperature. The free energy functional is , so the Legendre transform generates the T dependence of the correlation length. Note that in this Gaussian approximation γ is a constant, independent of T . We will see in the following Section that, in general, γ scales with T − T 0 according a non-trivial critical exponent.

III. CRITICAL ENTROPY FUNCTIONAL
We determine the critical entropy density using the Ising equation of state constructed by Zinn-Justin [7]. We write the Gibbs free energy as [7,35]  where r = (T − T c )/T c . The constants h 0 and M 0 will be specified below, α 0.11 is the specific heat exponent, and g(θ) is a function that is given in App. A. We use the (R, θ) where β 0.33 is the order parameter exponent.
We can expand the Gibbs free energy for small r and ψ. We find where Note that in the mean field approximation we have α = 0 and β = 1/2, so that G ∼ 1 2 rψ 2 , as expected. Also note that 2 − α − 2β =γ is the susceptibility exponent. The Ising energy density is The entropy functional determines the tri-linear coupling defined in the previous section.
We get There are two differences compared to the mean field result in the previous section. The first is that the tri-linear coupling vanishes near the transition point, with a critical scaling controlled by the exponent 3 (1 − 2β) 0.34. We will see that this is more than compensated by the divergence in the order parameter relaxation rate. The second is that there is an amplitude ratio, that means the coupling is different on the first order (r < 0) and crossover (r > 0) of the transition. We find Note that the scaling with r can be converted to a scaling relation involving the correlation length using |r| ∼ ξ −1/ν , with ν 0.63. Indeed, Zinn-Justin provides a scaling form of the correlation length where ξ 0 is an overall scale and g ξ (θ) (1 − 5θ 2 /18) [7].
Parotto et al. construct a map X A (X b ) from intensive QCD variables X a = (−β, βµ) to Ising variables X A = (r, h). The specific map considered in [9] is a simple linear relation wherew,ρ and α 1,2 are parameters. Most of the work in the existing literature [8,9,33] assumes that α 1 0 and α 2 π/2. We will discuss this assumption in more detail below, but for now we will assume that α 1 = 0 and α 2 = π/2. Parotto et al. choosew = 1,ρ = 2 and A = T 3 c . They also use M 0 = 0.605 and h 0 = 0.394. In order to determine the map between the conserved densities we use the fact that the transformation between the intensive variables can be specified in terms of the matrix (see 3 Note that in the mean field approximation β = 1/2 and γ ± is a constant. Also note that the scaling exponent of γ ± agrees with the value determined using renormalization group arguments [16,34]. Halperin et al. also provide a diagrammatic argument that explains why the critical scaling of γ is related to the specific heat and susceptibility exponents. This map induces a relationship between the densities x A (x b ), described by Since the densities are conjugate to the intensive variables we must have In the simple case that α 1 = 0 and α 2 = π/2 we obtain where we have used that the normalization of the Ising and QCD entropy functional differ by a factor A.

IV. CRITICAL BULK VISCOSITY
The bulk viscosity is determined by the Kubo relation where G ij,kl R (ω, k) is the retarded correlation function of the (spatial components) of the stress tensor The critical behavior arises from the contribution of the non-equilibrium pressures to the trace of the stress tensor, 1 3 Π ii = δP , with δP given in equ. (11). In order to compute the right hand side of the Kubo relation equ. (29) we consider the symmetric correlation where we have defined c = (nT Aa n γ) 2 . In statistical field theory we can factorize the correlation function and determine the retarded correlator using the fluctuation-dissipation relation. We get where ∆ S,R are the symmetric and retarded correlation function of the order parameter ψ.
Here, λ 0 is a transport coefficient and θ is a stochastic force. The noise correlator is The diffusion constant is D = λ 0 /χ 0 where the susceptibility is given by The retarded and symmetric order parameter correlation functions are where Γ k = λ 0 k 2 /χ k . The bulk viscosity is We note that the order parameter field ψ is dimensionless. We can absorb all dimensionful parameters in equ. (33) into a length scale, the non-critical correlation length ξ 0 , and a non-critical relaxation time t 0 . These are two free parameters in the model H description of the QCD critical point, in addition to the parameters that appear in the equation of state.
As a rough estimate of ξ 0 we will use the entropy density of the system, and assume that sξ 3 0 = 1. The relaxation time is determined by the diffusion constant, D 0 = t 0 /ξ 2 0 . The integral in equ. (38) is easily performed in the ω → 0 limit. We find where we have taken into account the scaling of γ ± ≡ γ R ± r 1−2β , see equ. (21). The critical exponent x ζ = 3.92 is larger than the one found by Onuki, see equ. (1). We will show in Sect. VI that this is due to the fact that in equ. (33) we have neglected the coupling of the order parameter to the momentum density of the fluid. The coupling γ R ± a n was determined in equ. (21) and (28). For typical values of the diffusion constant t 0 T is of order one; we have used t 0 1.8 fm from [23]. The coefficient 3/(32π) ∼ 3 · 10 −2 is a loop factor. The main uncertainty in equ. (39) is the factor n/s. For typical equations of state that put the critical endpoint within reach of the RHIC beam energy scan program this coefficient is quite small.
For example, the equation of state shown in Fig. 6 of Parotto et al. [9] has (n/s) crit 1/25.
We then get We conclude that the critical enhancement in the bulk viscosity is not large, unless the correlation length becomes quite large. We will consider the growth in the effective correlation length in the next Section, and we will consider a possible scenario that overcomes the suppression by n/s in Sect. VI.

V. BULK VISCOSITY IN AN EXPANDING SYSTEM
In a heavy ion collision the growth of the correlation length is limited by the failure of the evolution to go precisely through the critical point, and by the finite time available for the correlation length to reach its equilibrium value [23,36]. The recent study [23] concludes that the second effect dominates, and that ξ/ξ 0 is limited by ξ KZ /ξ 0 where ξ KZ is the Kibble-Zurek scale. Akamatsu et al. [23] estimate that ξ KZ /ξ 0 ∼ 1.33, so that the critical bulk viscosity in equ. (40) does not become large.
There is a second effect which appears even if the correlation length does become large.
Critical slowing down also implies that the bulk pressure only relaxes slowly to the equilibrium value ∆P ζ( ∇ · u), where u is the fluid velocity [22]. The response ∆P (t) to a time-dependent bulk stress ( ∇ · u)(t ) is given by where G(t) is the Fourier transform of equ. (38) The time integral of G(t) is given by the bulk viscosity ζ, and the asymptotic behavior of G(t) for large t is governed by hydrodynamic tails G(t) ∼ 1/t 3/2 [18,20,[37][38][39][40][41][42]. The tail of the critical contribution is which has an even stronger dependence on the correlation length than the static value of the bulk viscosity. We note, however, that the asymptotic behavior only sets in at parametrically We may also consider a bulk stress that is turned on at t = 0 and then study the evolution of ∆P (t) towards its asymptotic value ∆P (∞) = ζ( ∇· u). The relaxation is not exponential.
Instead, we find very slow relaxation where, again, the asymptotic behavior requires t ξ 4 . For small t we find that ∆P grows linearly in t, but the linear growth involves a much smaller power of the correlation length. This behavior is shown in Fig. 2. We observe that even for a modest enhancement of the correlation length, axes of the plot are the baryon chemical potential µ and the temperature T in units of the critical temperature in the chiral limit. The model is described in [46], and the parameters were chosen as explained in the text. The contour lines show the magnitude of the chiral susceptibility.

VI. REFINEMENTS
In this section we will discuss two refinements to the calculation described in the previous sections. The first involves the treatment of the order parameter relaxation rate. In model H it has been shown that the relaxation rate Γ k is governed by the coupling of ψ to the momentum density of the fluid [10]. Consider, for example, the choice ψ = s/n. Then the conservation laws generate a coupling of the form (∂ψ)/(∂t) = π · ∇ψ + . . ., where π is the momentum density of the fluid. The momentum relaxation rate has a very weak divergence, so that it makes sense to approximate that rate by the non-critical shear viscosity η 0 of the fluid. This is known as the Kawasaki approximation. Within this approximation [10,17,43,44] where We will also use a refined parameterization of the order parameter susceptibility where η 0.036 is the correlation function exponent and χ 0 = ξ 2 0 (ξ/ξ 0 ) 2−η . In this context we will continue to define ξ 0 by (sξ 3 0 ) = 1, but we will take t 0 = c ξ 0 , where c is the speed of light, and the relaxation time scale is set by the non-critical value of η 0 /s. We can now recompute the critical contribution to the bulk viscosity in equ. (38). We find where is the result of a numerical integral, and we have used the hyperscaling relation νd = 2−α = 2β + γ, where d = 3 is the number of spatial dimensions. For the Kawasaki function we have z = 3, which is close to the result in the -expansion, z = 3.05. Using the values of α and ν quoted above we find z − α/ν 2.8, which agrees with Onuki's result [16]. For We note that the smallness of the pre-exponent in equ. (51) is mainly due to the small factor (n/s). The appearance of this factor is due to the assumption about the orientation of the Ising axes in the QCD phase diagram made above equ. (7). It was recently observed that close to the chiral limit this assumption is not correct, and that in models of QCD this observation also applies to realistic values of the quark masses [45]. To illustrate this fact we show in Fig. 3 the critical endpoint and the orientation of the Ising axes in the random matrix model introduced in [46]. The overall scale was fixed by the requirement that T χ , the phase transition temperature in the chiral limit, and T pc , the the pseudo-critical temperature for the physical value of the pion and kaon mass, is the same as in lattice QCD, T χ /T pc (132 +3 −6 MeV )/(156.5±1.5 MeV ) [47]. In the random matrix model this correspond to a quark mass m q 10 MeV. We find α 1 105 • and α 2 165 • . In this case the dominant contribution to δP arises from with a e = (∂ )/(∂(δe)). Using the mapping relation in equ. (23,24) together with equ. (27) we find a e = βρw sin(α 1 ) and With sin 2 (α 1 ) 1/4 from Fig. 3 we observe that, at least on the first order side of the transition, the critical bulk viscosity can reach ζ/s ∼ 1/(4π), comparable to the non-critical shear viscosity 4 .

VII. CONCLUSIONS AND OUTLOOK
In this work we have studied the critical bulk viscosity in QCD. We find that the result is very sensitive to the orientation of the Ising axes in the QCD phase diagram. It is often assumed that the Ising temperature axis is aligned with the QCD chemical potential axis, and that the Ising magnetic field axis is orthogonal to that. In that case, the QCD bulk viscosity is suppressed by a factor (n/s) 2 . If we include a possible misalignment of the Ising axes our final result is given in equ. (53). We observe that there is significant dependence on what side of the phase transition we are considering; the critical contribution to the bulk viscosity is about an order of magnitude larger on the first order side of the transition.
In this regime the critical enhancement of the bulk viscosity may become comparable to the non-critical shear viscosity, even for a modest enhancement of the correlation length, ξ/ξ 0 ∼ 2.
We have also studied the response function for the bulk pressure in the critical regime.
As expected, relaxation is very slow. The bulk pressure approaches its equilibrium value as 1/ √ t, and the initial rise of the bulk pressure involves a much smaller power of the correlation length. Note that in the present work we have only considered the Fourier transform of the equilibrium response. A more complete treatment in an expanding system would be based on computing the response from the real time correlation functions of the order parameter in an expanding medium [20,23]. Also note that we focused on the leading critical enhancement in the bulk viscosity. In a real heavy ion collision it may well be the case that the bulk viscosity is dominated by contributions that are more weakly singular near the critical endpoint.
We would like to thank M. Stephanov for useful discussions. This work was supported in part by the US Department of Energy grant DE-FG02-03ER41260 and by the BEST (Beam Energy Scan Theory) DOE Topical Collaboration.