Critical behavior of the bulk viscosity in QCD

TL;DR

This paper investigates the critical behavior of bulk viscosity in QCD near a potential critical point, estimating its divergence and sensitivity to assumptions about the phase diagram's alignment.

Contribution

It provides a quantitative estimate of the bulk viscosity's critical divergence in QCD and highlights the importance of the phase diagram's orientation assumptions.

Findings

Bulk viscosity diverges as $(\xi/\xi_0)^{x_\xi}$ with $x_\xi \,\approx 3$.

The proportionality constant $a$ is approximately $10^{-4}$.

Critical bulk viscosity can approach shear viscosity if the phase diagram assumptions are altered.

Abstract

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$.

Equations

$$\eta \sim {\xi }^{0.05},\kappa \sim {\xi }^{0.9},\zeta \sim {\xi }^{2.8}.$$

$$x^A=(ϵ,\psi ),$$

$$X_A=-\frac{\partial s}{\partial x^A}=(r,h).$$

$$x^a=(e,n),X_a=(-\beta ,\beta \mu ),$$

$$\delta P={\frac{\partial P}{\partial \beta }|}_{\beta \mu }\delta \beta +{\frac{\partial P}{\partial \beta \mu }|}_{\beta }\delta (\beta \mu )=-\frac{e+P}{\beta }\delta \beta +\frac{n}{\beta }\delta (\beta \mu ).$$

$$\delta P=\frac{e+P}{\beta }\frac{\partial s}{\partial (\delta e)}-\frac{n}{\beta }\frac{\partial s}{\partial (\delta n)}.$$

$$\delta P=\frac{e+P}{\beta }\frac{\partial \psi }{\partial (\delta e)}\frac{\partial s^{Is}}{\partial \psi }-\frac{n}{\beta }\frac{\partial ϵ}{\partial (\delta n)}\frac{\partial s^{Is}}{\partial ϵ}.$$

$$\exp (-\beta G[\psi ,r])=\int \mathcal{D}ϵ\exp (S[\psi ,ϵ]-E/T),E=\int d^{3}x(ϵ+{ϵ}_0),$$

$$S[\psi ,ϵ]=-\int d^{3}x\left\{k(\nabla \psi {)}^2+\frac{v}{2}{\psi }^2+\frac{u}{4}{\psi }^4+\gamma ϵ{\psi }^2+\frac{1}{2C_0}{ϵ}^2\right\}+S_0+\frac{E}{T_0},$$

$$\beta G[\psi ,r]=\int d^{3}x\left\{k(\nabla \psi {)}^2+\frac{\tilde{v}}{2}{\psi }^2+\frac{\tilde{u}}{4}{\psi }^4\right\}$$

$$\delta P=nTA{a}_{nϵ}\gamma {\psi }^2$$

$$G[\psi ,r]=h_0{M}_0{R}^{2-\alpha }g(\theta ),$$

$$\psi$$

$$r$$

$$h=h_0{R}^{\beta \delta }\tilde{h}(\theta ),$$

$$G[\psi ,r]\simeq h_0{M}_0\left\{g_{\pm }{r}^{2-\alpha }+\frac{1}{2}{m}_{\pm }^2{r}^{2-\alpha -2\beta }(\psi -{\psi }_0{)}^2+\dots \right\},$$

$$g_{\pm}=\left\{\begin{array}[]{cc}-0.84&r>0\\ -1.58&r<0\end{array}\right.\hskip 28.45274ptm^{2}_{\pm}=\left\{\begin{array}[]{cc}1.00&r>0\\ 4.77&r<0\end{array}\right.\hskip 28.45274pt\psi_{0}=\left\{\begin{array}[]{cc}0&r>0\\ \psi_{\pm}&r<0\end{array}\right.$$

$$ϵ=\frac{\partial G}{\partial r}\simeq h_0{M}_0\left\{g_{\pm }(2-\alpha )r^{1-\alpha }+\frac{1}{2}{m}_{\pm }^2(2-\alpha -2\beta )r^{1-\alpha -2\beta }(\psi -{\psi }_0{)}^2+\dots \right\},$$

$$s[\psi ,ϵ]$$

$${\gamma }_{\pm }=\frac{m_{\pm }^2}{2(-g_{\pm })}\frac{2-\alpha -2\beta }{(2-\alpha )(1-\alpha )}|r{|}^{1-2\beta }.$$

$$\gamma_{\pm}=\left\{\begin{array}[]{cc}0.43\,r^{1-2\beta}&r>0\\ 1.10\,|r|^{1-2\beta}&r<0\end{array}\right.$$

$$\xi ={\xi }_0{R}^{-\nu }{g}_{\xi }^{1/2}(\theta ),$$

$$\frac{T-T_c}{T_c}$$

$$\frac{\mu -{\mu }_c}{T_c}$$

$${\bar{R}}_A^b=\frac{\partial X_A}{\partial X_b}.$$

$$R_b^A=\frac{\partial x^A}{\partial x^b}.$$

$$R_b^A{\bar{R}}_C^b={\delta }_C^A,$$

$$a_{nϵ}=\frac{\partial ϵ}{\partial (\delta n)}=\frac{\bar{\rho }\bar{w}}{A}$$

$$\zeta (\omega )=-\frac{1}{9\omega }\mathrm{Im}{G}_R^{ii,jj}(\omega ,0)$$

$$G_R^{ij,kl}(\omega ,k)=-i\int dt{d}^{3}x{e}^{-i(\omega t-kx)}⟨[{\Pi }^{ij}(0,0),{\Pi }^{kl}(x,t)]⟩\theta (t).$$

$$G_S(\omega ,k)=c\int dt{d}^{3}x{e}^{-i(\omega t-kx)}⟨{\psi }^2(0){\psi }^2(t,x)⟩$$

Full text

I Introduction

There are several programs dedicated to exploring the phase diagram of Quantumchromodynamics (QCD) at heavy ion accelerator laboratories around the world Braun-Munzinger:2015hba . A central feature of the phase diagram is a possible 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 Stephanov:1998dy .

In a static system fluctuation observables are controlled by the critical equation of state Stephanov:2008qz ; Asakawa:2009aj ; Stephanov:2011pb . 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 Engels:2002fi , and accurate parametrizations are available in the literature Zinn-Justin:2002ru . 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 Nonaka:2004pg ; Parotto:2018pwx .

The dynamic behavior of fluctuations is expected to be governed by model H in the classification of dynamical critical phenomena by Hohenberg and Halperin Son:2004iv ; Hohenberg:1977ym ; Folk:2006ve . Model H is a hydrodynamic theory that describes to 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 $\eta$, the bulk viscosity $\zeta$, and the thermal conductivity $\kappa$. In model H fluctuations of the order parameter with wave number $q\sim\xi^{-1}$, where $\xi$ is the correlation length, relax on a time scale $\tau\sim\xi^{z}$, where $z\simeq 3$. This behavior is intermediate between ordinary diffusion ($z\simeq 2$), and critical relaxation of a conserved charge not coupled to fluctuations of the fluid velocity ($z\simeq 4$). There is a very mild divergence in the shear viscosity, and more pronounced critical behavior in the thermal conductivity and bulk viscosity Kadanoff:1968 ; Hohenberg:1977ym ; Kroll:1981 ; Ferrell:1985 ; Onuki:1997 ; Onuki:2002 ,

[TABLE]

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 Onuki:2002 .

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 Akamatsu:2016llw ; Akamatsu:2017rdu ; Martinez:2018wia , simulations of critical stochastic diffusion in an expanding medium Nahrgang:2018afz , and deterministic frameworks for the evolution of two-point Stephanov:2017ghc ; Akamatsu:2018vjr ; An:2019osr or higher n-point functions Mukherjee:2015swa 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 Karsch:2007jc ; Monnai:2016kud . Our goal in the present work is to study this problem more quantitatively, based on the equation of state constructed by Parotto et al. Parotto:2018pwx . 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 Arnold:2006fz ; Moore:2008ws ; Lu:2011df ; Dusling:2011fd .

II Fluctuations of the energy density and pressure

Fluctuations of the order parameter are governed by an entropy functional $S=\int d^{3}x\,s$. The entropy functional of the Ising model is a function of the densities

[TABLE]

where $\epsilon$ and $\psi$ are the Ising energy density and order parameter. The corresponding intensive variables are the reduced temperature $r$ and the magnetic field $h$,

[TABLE]

We are following the notation of Landau and Lifshitz Landau:smI as well as Akamatsu et al. Akamatsu:2018vjr . The analogous canonical pair in QCD is

[TABLE]

where $(e,n)$ are the energy and baryon density, $\beta=1/T$ is the inverse temperature and $\mu$ 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

[TABLE]

Using the fact that $(\beta,\beta\mu)$ is conjugate to $(e,n)$ we obtain

[TABLE]

The map between $(r,h)$ and $(\beta,\beta\mu)$ induces a map between the QCD densities $(e,n)$ on the Ising densities $(\epsilon,\psi)$. We assume that the singular part of the QCD entropy density is proportional to the Ising entropy density, $s^{\it sing}(e,n)=As^{\it Is}(\epsilon(e,n),\psi(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 axis111We will reconsider this assumption in Sect. V., $\partial\epsilon/\partial e\simeq 0$ Nonaka:2004pg ; Bluhm:2016byc ; Parotto:2018pwx . This implies that

[TABLE]

Onuki observed that the Ising entropy functional contains a tri-linear coupling between $\epsilon$ and $\psi^{2}$, and as a result fluctuations of the pressure couple to $\psi^{2}$ Onuki:1997 ; Onuki:2002 . This means that correlation functions of $\delta P$ are controlled by the order parameter relaxation rate, and the slow relaxation of $\psi$ leads to a large enhancement in the bulk viscosity.

The presence of a tri-linear coupling between $\epsilon$ and $\psi^{2}$ in the Ising entropy $S[\psi,\epsilon]$ can be understood from the relation between entropy and Gibbs free energy222Note that there is some disagreement in the literature on whether $G[\psi,r]$ should be called free energy or Gibbs free energy. Zinn-Justin refers to this functional as a free energy and equ. (3) of Asakawa et al. Nonaka:2004pg as well as Parotto et al. Parotto:2018pwx follows that notation, but Akamatsu et al. Akamatsu:2018vjr uses the term Gibbs free energy, and reserves the term free energy for $F[h,r]$. Halperin:1974

[TABLE]

where $E$ is the total energy, $\epsilon_{0}$ is the background energy density, and $T$ is the temperature. For small $\epsilon$ and $\psi$ we can expand the entropy functional as

[TABLE]

where $k,v,u$ and $\gamma$ 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

[TABLE]

and $\tilde{v}=v+2\gamma C_{0}(T-T_{0})/T_{0}^{2}$, so the Legendre transform generates the $T$ dependence of the correlation length. Note that in this Gaussian approximation $\gamma$ is a constant, independent of $T$. We will see in the following Section that, in general, $\gamma$ scales with $T-T_{0}$ according a non-trivial critical exponent.

Finally, combining equ. (7) and equ. (9) we obtain the slow mode contribution to fluctuations of the pressure

[TABLE]

where we have defined $a_{n\epsilon}=(\partial\epsilon)/(\partial(\delta n))$.

III Critical entropy functional

We determine the critical entropy density using the Ising equation of state constructed by Zinn-Justin Zinn-Justin:2002ru . We write the Gibbs free energy as Schofield:1969 ; Zinn-Justin:2002ru

[TABLE]

where $r=(T-T_{c})/T_{c}$. The constants $h_{0}$ and $M_{0}$ will be specified below, $\alpha=0.11$ is the specific heat exponent, and $g(\theta)$ is a function that is given in App. A. We use the $(R,\theta)$ coordinates

[TABLE]

where $\beta=0.33$ is the order parameter exponent. Here, $R\in[0,\infty)$ and $\theta\in[-\theta_{0},\theta_{0}]$, where $\theta_{0}$ is determined by the condition $\tilde{h}(\theta_{0})=0$. Furthermore, $\tilde{h}(\theta)$ is a function that appears in the magnetic equation of state

[TABLE]

$\delta=4.78$ is the external field exponent, and $\tilde{h}(\theta)$ is also specified in App. A. For this parameterization we find $\theta_{0}=1.154$. The value of $\theta_{0}$ determines the boundaries of the first order region for $r<0$, $\psi_{\pm}=\pm M_{0}R^{\beta}\theta_{0}$.

We can expand the Gibbs free energy for small $r$ and $\psi$. We find

[TABLE]

where

[TABLE]

Note that in the mean field approximation we have $\alpha=0$ and $\beta=1/2$, so that $G\sim\frac{1}{2}r\psi^{2}$, as expected. Also note that $2-\alpha-2\beta=\tilde{\gamma}$ is the susceptibility exponent. The Ising energy density is

[TABLE]

We can now determine an entropy functional that describes fluctuations of the energy density at constant $\psi$ by performing a Legendre transformation, $s[\psi,\epsilon]=G[\psi,r]-r\epsilon$. Expanding $s$ for small values of the arguments gives

[TABLE]

We can now read off the tri-linear coupling defined in the previous section. We get

[TABLE]

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 exponent333Note that this exponent agrees with the exponent determined using renormalization group arguments Halperin:1974 ; Onuki:1997 . Halperin et al. provide a diagrammatic argument that explains why the critical scaling of $\gamma$ is related to the specific heat and susceptibility exponents. $(1-2\beta)\simeq 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

[TABLE]

Note that the scaling with $r$ can be converted to a scaling relation involving the correlation length using $|r|\sim\xi^{-1/\nu}$, with $\nu\simeq 0.63$. Indeed, Zinn-Justin provides a scaling form of the correlation length

[TABLE]

where $\xi_{0}$ is an overall scale and $g_{\xi}(\theta)\simeq(1-5\theta^{2}/18)$ Zinn-Justin:2002ru .

Parotto et al. construct a map $X_{A}(X_{b})$ from intensive QCD variables $X_{a}=(-\beta,\beta\mu)$ to Ising variables $X_{A}=(r,h)$. The specific map considered in Parotto:2018pwx is a simple linear relation

[TABLE]

where $\bar{w},\bar{\rho}$ and $\alpha_{1,2}$ are parameters. Most of the work in the existing literature Nonaka:2004pg ; Parotto:2018pwx ; Bluhm:2016byc assumes that $\alpha_{1}\simeq 0$ and $\alpha_{2}\simeq\pi/2$. We will discuss this assumption in more detail below, but for now we will assume that $\alpha_{1}=0$ and $\alpha_{2}=\pi/2$. Parotto et al. choose $\bar{w}=1,\bar{\rho}=2$ and $A=T_{c}^{3}$. 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 Fig. 1)

[TABLE]

This map induces a relationship between the densities $x_{A}(x_{b})$, described by

[TABLE]

Since the densities are conjugate to the intensive variables we must have

[TABLE]

which implies that the matrix $R$ is the inverse of $\bar{R}$. In the simple case that $\alpha_{1}=0$ and $\alpha_{2}=\pi/2$ we obtain

[TABLE]

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

[TABLE]

where $G_{R}^{ij,kl}(\omega,k)$ is the retarded correlation function of the (spatial components) of the stress tensor

[TABLE]

The critical behavior arises from the contribution of the non-equilibrium pressures to the trace of the stress tensor, $\frac{1}{3}\Pi^{ii}=\delta P$, with $\delta P$ given in equ. (11). In order to compute the right hand side of the Kubo relation equ. (29) we consider the symmetric correlation function

[TABLE]

where we have defined $c=(nTAa_{n\epsilon}\gamma)^{2}$. In statistical field theory we can factorize the correlation function and determine the retarded correlator using the fluctuation-dissipation relation. We get

[TABLE]

where $\Delta_{S,R}$ are the symmetric and retarded correlation function of the order parameter $\psi$. The slow dynamics of $\psi$ is governed by a diffusion equation Hohenberg:1977ym

[TABLE]

Here, $\lambda_{0}$ is a transport coefficient and $\theta$ is a stochastic force. The noise correlator is

[TABLE]

The diffusion constant is $D=\lambda_{0}/\chi_{0}$ where the susceptibility is given by

[TABLE]

The retarded and symmetric order parameter correlation functions are

[TABLE]

where $\Gamma_{k}=\lambda_{0}k^{2}/\chi_{k}$. The bulk viscosity is

[TABLE]

We note that the order parameter field $\psi$ is dimensionless. We can absorb all dimensionful parameters in equ. (33) into a length scale, the non-critical correlation length $\xi_{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 $\xi_{0}$ we will use the entropy density of the system, and assume that $s\xi_{0}^{3}=1$. We will comment on this choice in Sect. VII. The relaxation time is determined by the diffusion constant, $D_{0}=t_{0}/\xi_{0}^{2}$.

The integral in equ. (38) is easily performed in the $\omega\to 0$ limit. We find

[TABLE]

where we have taken into account the scaling of $\gamma_{\pm}\equiv\gamma^{R}_{\pm}r^{1-2\beta}$, see equ. (21). The critical exponent $x_{\zeta}=3.92$ is larger than the one found by Onuki, see equ. (1). We will show in Sect. V 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 $\gamma^{R}_{\pm}a_{n\epsilon}$ 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}\simeq 1.8$ fm from Akamatsu:2018vjr . The coefficient $3/(32\pi)\sim 3\cdot 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. Parotto:2018pwx has $n/s\simeq 1/25$. We then get

[TABLE]

We conclude that the critical enhancement in the bulk viscosity is not large, unless the correlation length becomes quite large. We will study the growth in the effective correlation length in the Sect. VI, and we will consider a possible scenario that overcomes the suppression by $n/s$ in the following Section.

V 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 $\Gamma_{k}$ is governed by the coupling of $\psi$ to the momentum density of the fluid Hohenberg:1977ym . Consider, for example, the choice $\psi=s/n$. Then the conservation laws generate a coupling of the form $(\partial\psi)/(\partial t)=\vec{\pi}\cdot\vec{\nabla}\psi+\ldots$, where $\vec{\pi}$ 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 $\eta_{0}$ of the fluid. This is known as the Kawasaki approximation. Within this approximation Kawasaki:1970 ; Hohenberg:1977ym ; Onuki:2002

[TABLE]

where

[TABLE]

We will also use a refined parameterization of the order parameter susceptibility

[TABLE]

where $\eta\simeq 0.036$ is the correlation function exponent and $\chi_{0}=\xi_{0}^{2}(\xi/\xi_{0})^{2-\eta}$. In this context we will continue to define $\xi_{0}$ by $(s\xi_{0}^{3})=1$, but we will take $t_{0}=c\,\xi_{0}$, where $c$ is the speed of light, and the relaxation time scale is set by the non-critical value of $\eta_{0}/s$. We can now recompute the critical contribution to the bulk viscosity in equ. (38). We find

[TABLE]

where

[TABLE]

is the result of a numerical integral, and we have used the hyperscaling relation $\nu d=2-\alpha=2\beta+\gamma$, 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 $\epsilon$-expansion, $z=3.05$. Using the values of $\alpha$ and $\nu$ quoted above we find $z-\alpha/\nu\simeq 2.8$, which agrees with Onuki’s result Onuki:1997 . For $\eta_{0}/s\simeq 1/(4\pi)$ the estimate in equ. (44) is numerically close to the result given above. We find

[TABLE]

We note that the smallness of the pre-exponent in equ. (46) 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 Pradeep:2019ccv . To illustrate this fact we show in Fig. 2 the critical endpoint and the orientation of the Ising axes in the random matrix model introduced in Halasz:1998qr . The overall scale was fixed by the requirement that $T_{\chi}$, 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_{\chi}/T_{pc}\simeq(132^{+3}_{-6}\,{\it MeV})/(156.5\pm 1.5\,{\it MeV})$ Ding:2019prx . In the random matrix model this correspond to a quark mass $m_{q}\simeq 10$ MeV. We find $\alpha_{1}\simeq 105^{\circ}$ and $\alpha_{2}\simeq 165^{\circ}$. In this case the dominant contribution to $\delta P$ arises from

[TABLE]

with $a_{e\epsilon}=(\partial\epsilon)/(\partial(\delta e))$. Using the mapping relation in equ. (23,24) together with equ. (27) we find $Aa_{e\epsilon}=\beta\bar{\rho}\bar{w}\sin(\alpha_{1})$ and

[TABLE]

With $\sin^{2}(\alpha_{1})\simeq 1/4$ from Fig. 2 we observe that, at least on the first order side of the transition, the critical bulk viscosity can reach $\zeta/s\sim 1/(4\pi)$, comparable to the non-critical shear viscosity444Note that Parotto et al. Parotto:2018pwx use $\alpha_{1}=3.85^{o}$, based on the curvature of the freezeout curve. This corresponds to $\sin(\alpha_{1})^{2}\simeq 4\cdot 10^{-3}$.. Note that equ. (46,48) are limiting cases, corresponding to $\sin(\alpha_{1})^{2}\sim 0$ and \sin(\alpha_{1})^{2}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\sim}\hss}\raise 1.0pt\hbox{>}}(n/s)^{2}, of a more general relation that involves both $a_{n\epsilon}$ and $a_{e\epsilon}$. The more general result follows from the coupling $\delta P=\bar{\rho}\bar{w}[s\sin(\alpha_{1})+n(\cos(\alpha_{1})-\beta\mu\sin(\alpha_{1}))]TA\gamma\psi^{2}$.

VI 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 Berdnikov:1999ph ; Akamatsu:2018vjr . The recent study Akamatsu:2018vjr concludes that the second effect dominates, and that $\xi/\xi_{0}$ is limited by $\xi_{KZ}/\xi_{0}$ where $\xi_{KZ}$ is the Kibble-Zurek scale. Akamatsu et al. Akamatsu:2018vjr estimate that $\xi_{KZ}/\xi_{0}\sim 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 $\Delta P\simeq\zeta(\vec{\nabla}\cdot\vec{u})$, where $\vec{u}$ is the fluid velocity Stephanov:2017ghc . The response $\Delta P(t)$ to a time-dependent bulk stress $(\vec{\nabla}\cdot\vec{u})(t^{\prime})$ is given by

[TABLE]

where $G(t)$ is the Fourier transform of equ. (38)

[TABLE]

In the following, we will use the simple form $\Gamma_{k}=\lambda_{0}k^{2}\chi_{k}^{-1}$. The result is easily generalized to the Kawasaki approximation in equ. (41). The time integral of $G(t)$ is given by the bulk viscosity $\zeta$, and the asymptotic behavior of $G(t)$ for large $t$ is governed by hydrodynamic tails $G(t)\sim 1/t^{3/2}$ Alder:1967 ; Ernst:1971 ; Pomeau:1975 ; Kovtun:2003vj ; Kovtun:2012rj ; Martinez:2017jjf ; Akamatsu:2016llw ; Martinez:2018wia . The tail of the critical contribution is

[TABLE]

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 late time, $(t/t_{0})\gg(\xi/\xi_{0})^{4}$.

We may also consider a bulk stress that is turned on at $t=0$ and then study the evolution of $\Delta P(t)$ towards its asymptotic value $\Delta P(\infty)=\zeta(\vec{\nabla}\cdot\vec{u})$. The relaxation is not exponential. Instead, we find very slow relaxation

[TABLE]

where, again, the asymptotic behavior requires $t\gg\xi^{4}$. For small $t$ we find that $\Delta P$ grows linearly in $t$,

[TABLE]

but the linear growth involves a much smaller power of the correlation length. This behavior is shown in Fig. 3. We observe that even for a modest enhancement of the correlation length, $\xi/\xi_{0}=1.5$, the critical bulk pressure relaxes on a time scale much slower than the natural equilibration time $t_{0}$.

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. (48). 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, $\xi/\xi_{0}\sim 2$.

We note that the overall scale for $\zeta/s$ is set by the two parameters $\eta_{0}/s$ and $s\xi_{0}^{3}$. While $\eta_{0}/s$ is fairly well constrained by data, this is not the case for $s\xi_{0}^{3}$. We have assumed that $s\xi_{0}^{3}\simeq 1$, based on the idea that in a hadron gas the correlation length cannot be shorter than the distance between particles. However, the exact value of $\xi_{0}$ will depend on the composition of the gas, and the precise observable used to define the correlation length. For example, the recent work of Akamatsu et al. Akamatsu:2018vjr uses $\xi_{0}\simeq 1.2$ fm. Combined with a freezeout value $sT^{3}\simeq 7$ Cleymans:2005xv this number corresponds to $s\xi_{0}^{3}\simeq 5$.

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/\sqrt{t}$, and the initial rise of the bulk pressure involves a much smaller power of the correlation length than the scaling of $\zeta$. 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 Akamatsu:2018vjr ; Martinez:2018wia .

Finally, we may compare our results to simple scaling relations that have been discussed in the literature. Weinberg proposed the relation $\zeta\sim(c_{s}^{2}-1/3)^{2}\eta$, where $c_{s}$ is the speed of sound Weinberg:1972 . This relation is satisfied in the weak coupling limit of QCD Arnold:2006fz , and it holds in simple relaxation time models of kinetic theory Dusling:2011fd . Weinberg’s relation predicts an enhancement of the bulk viscosity near the critical point, but no critical divergence of the type studied in this work. Based on an analysis of spectral sum rules, Karsch et al. suggested a different scaling relation, $\zeta/s\sim(1/c_{s}^{2}-3)$ Karsch:2007jc . This formula predicts a mild divergence of the bulk viscosity, $\zeta/s\sim\xi^{\alpha/\nu}$, much weaker than our result $\zeta/s\sim\xi^{z-\alpha/\nu}$.

Implementations of second order non-conformal relativistic fluid dynamics take into account a bulk viscous relaxation time $\tau_{\zeta}$. A typical assumption, motivated by kinetic theory, is that $\tau_{\zeta}\sim\zeta/P$ Dusling:2011fd ; Monnai:2016kud . This relation captures in an approximate way the critical slowing down of the response discussed in Sect. VI, but a simple relaxation time approximation cannot capture the hydrodynamic tail given in equ. (51).

We would like to thank M. Stephanov and Yi Yin 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.

Appendix A Ising equation of state

We use a parameterization of the Ising equation of state constructed by Zinn-Justin Zinn-Justin:2002ru . We write

[TABLE]

with

[TABLE]

and $g(1)=0.0424455$ as well as

[TABLE]

The magnetic equation of state is $h=h_{0}R^{\beta\delta}\tilde{h}(\theta)$ with

[TABLE]

and

[TABLE]