Charm contribution to bulk viscosity

In the range of temperatures reached in future heavy ion collision experiments, hadronic pair annihilations and creations of charm quarks may take place within the lifetime of the plasma. As a result, charm quarks may increase the bulk viscosity affecting the early stages of hydrodynamic expansion. Assuming thermalization, we estimate the charm contribution to bulk viscosity within the same effective kinetic theory framework in which the light parton contribution has been computed previously. The time scale at which this physics becomes relevant is related to the width of the transport peak associated with the trace anomaly correlator, and is found to be<20 fm/c for T>600 MeV.


Introduction
Recently concrete thoughts about a possible successor to the LHC have been aired [1]. One of the issues under discussion is whether there is a case for a Heavy Ion Collision program at energies much higher than those achievable at the LHC [2]. An example of a possible new physics observable is that if temperatures up to ∼1 GeV could be reached, even charm quarks might become a chemically equilibrated part of the Quark-Gluon Plasma. As the system cools down, their interactions slow down much faster than those of the light partons. Therefore, we might expect a chemical freeze-out process within an otherwise thermal medium, akin to that associated with some dark matter scenarios in cosmology.
The purpose of the present paper is to investigate the role of charm quarks within the weak-coupling expansion. More precisely, we consider a physical observable, the bulk viscosity, which is intimately related to chemical equilibration.
Determining transport coefficients such as the bulk viscosity is a notoriously difficult task even within perturbation theory. Starting from massive quarks (of mass M) and very a e-mail: laine@itp.unibe.ch b e-mail: sohrabi@itp.unibe.ch low temperatures (T M/π ), such that thermal resummations might be expected not to play a role, it turns out that a field-theoretic determination of the rate of pair annihilation and creation necessitates a 3-loop computation [3]. The situation becomes worse at higher temperatures T M/g, where g 2 ≡ 4πα s , when all-order resummations are needed for obtaining the correct leading-order result [4]. It has been demonstrated, however, that the same results can be obtained relatively economically by making use of effective kinetic theory [5,6], which describes quasi-particles having temperature-dependent properties. 1 We adopt this framework in the present paper.
Our goal is to estimate the charm quark contribution to bulk viscosity, as well as their chemical equilibration rate, at temperatures between about 200 and 700 MeV. In this regime the charm quarks, with a pole mass M ∼ 1.5 GeV, can be considered dilute and perhaps also kinetically equilibrated, as is suggested by experiment [9,10]. Ideally, it would be nice to also consider temperatures in the ultrarelativistic regime T M/g in order to be able to compare with a previous computation [11], but this is not achieved in the present paper. Note that if chemical equilibration were to take place, charm quarks should be included in the equation of state, where they would have a substantial influence already in the temperature range considered [12][13][14][15] (they could be relevant even if not fully equilibrated [16]).
The plan of this paper is the following. After outlining the general kinetic theory approach in Sect. 2, the specific setup relevant for our problem is detailed in Sect. 3. Thermally averaged annihilation rates are evaluated in Sect. 4. A numerical solution of the basic equations is presented in Sect. 5, whereas Sect. 6 offers some conclusions and an outlook.

General approach
We start by outlining the general strategy of the analysis in the language of kinetic theory and linear algebra. The precise implementation will be given in Sect. 3. Our notation and basic philosophy follow closely those presented in Ref. [11], to which we refer for further details.
Suppressing all indices, let S be a source vector enforcing a desired "adiabatic" deviation from equilibrium. 2 The system responds, as dictated by a linearised collision matrix C (C T = C), by deviating the phase space distributions by an amount χ around their equilibrium values: (2.1) The bulk viscosity ζ is given by the projection of χ in the direction of S: Let us assume that we know the normalised eigenfunctions v n of C (v T m v n = δ mn ) and the corresponding eigenvalues λ n ∈ R.
In general, the collision matrix has zero modes (Cξ = 0) and is therefore not invertible. However, if the system can equilibrate, as is expected, then there should be a non-trivial solution χ = 0 to S = Cχ . In this case the zero modes are necessarily orthogonal to S: S T ξ = χ T C T ξ = χ T Cξ = 0. Therefore they do not contribute to ζ and can be omitted. In other words, we can restrict to a subspace orthogonal to the zero modes and then invert C in a spectral representation: We observe that the dominant contribution to ζ is given by the smallest eigenvalues, provided that the corresponding v n have a non-vanishing projection in the direction of S. The physical realisation of this picture is the following. For bulk viscosity, the source S breaks conformal invariance. If QCD is approximated as a weakly coupled theory with N f = 3 massless flavours and the massive charm quark, then S is non-zero only because of the presence of the charm quark, and it is proportional to e −M/T . Elastic scatterings as well as 1 → 2 splittings of the light partons produced in charm decays have rates much faster than the charm pair annihilations and creations. It is the latter processes which lead to the smallest eigenvalues and are therefore considered in the following.

Specific setup
Bulk viscosity can be defined as a transport coefficient corresponding to the operator representing (minus) the trace anomaly, where T μν is the energy-momentum tensor. The part T 00 does not contribute to the bulk viscosity, given that its space average is exactly conserved. It can therefore be left out, or subtracted with a suitable coefficient. Within effective kinetic theory [6], there is a structure related to Eq. 3.1 [5,11], representing a shift in pressure ( p = T ii ) when subtracting the contribution of pressure perturbations carrying energy where v 2 s is the sound speed squared). For a momentum mode k and a particle of species a, the relevant subtracted perturbation is of the form [11] Here v a k = ∇ k E a k is a group velocity; k ≡ |k|; β ≡ 1/T ; and E a k is an on-shell energy. The associated source vector will be denoted by where f a is the equilibrium (Bose or Fermi) distribution of a particle of type a; and + and − correspond to bosons and fermions, respectively.
For v 2 s = 1/3 and massless particles with E a k = k, S a vanishes identically. More generally, as pointed out in Ref. [11], all S a vanish identically in a conformal theory. However, with massive quarks S a is non-zero even at leading order, as we now show.
In free SU (N c ) gauge theory with N f massless flavours and one massive quark, v 2 s evaluates to 3 . Both for M π T and M π T this reduces to 1 3 . For M π T , the corrections are exponentially suppressed: . (3.5) In this case the source vectors of Eq. 3.3 read, for charm quarks (S c ), gluons (S g ) and light quarks (S q ), respectively, Let us elaborate on the nature of the approximation that led to Eqs. 3.6-3.8. It is well known that when considering massive particles at finite temperature, two different types of thermal effects appear: exponential corrections and power corrections. As an example, consider the integral Given that e −E k /T < 1, this series is absolutely convergent. The individual terms in the series lead to where K is a modified Bessel function. The modified Bessel function is of the form However, the coefficients c i grow factorially; the series is asymptotic, with a zero radius of convergence. The purpose of the present paper is to compute the bulk viscosity to leading order in e −M/T but to all orders in T /M, avoiding asymptotic series. This may be called a dilute approximation, in contrast to a non-relativistic approximation which would also truncate the power corrections to a finite order in T /M. It can easily be verified that the source of Eqs. 3.6-3.8 carries no energy, . This is related to the discussion above Eq. 2.3 concerning zero modes; E a k plays the role of ξ . Indeed the full collision term of a Boltzmann equation vanishes at any temperature if equilibrium distributions are used; the temperature derivative of this statement asserts that χ a (k) ∝ E a k is a zero-mode of the linearised collision matrix C. So, Eq. 3.11 represents a crosscheck that S a (k) sources a deviation from equilibrium in a subspace orthogonal to the zero mode.
It may be noted that for T M, all the terms in Eqs. 3.6-3.8 are parametrically of a similar magnitude ∼ e −M/T . However, the rates at which different particles respond to this perturbation are different: number changing rates for charm quarks are ∼ α 2 s e −2M/T , whereas for light partons they are α 2 s for k M, or α 2 s e −M/T for k ∼ M. Putting the sources S a on the left-hand side of the equation (cf. Eq. 2.1) and the rates on the right-hand side, we may estimate the magnitudes of the perturbations χ a . We get χ c ∼ e M/T / α 2 s χ g , χ q . For this reason χ g , χ q give an exponentially suppressed contribution in Eq. 2.2. In other words, we can assume light partons to be completely thermalised (χ g , χ q ∼ 0) and consider in essence only the block C cc of the linearised collision matrix.
An important further point concerns kinetic equilibration. The elastic scatterings of charm quarks on light partons are ∼ α 2 s e −M/T . These reactions are again faster than those corresponding to chemical equilibration, so we may assume full kinetic equilibrium to be reached. Full kinetic equilibrium implies that χ c is a k-independent constant. Even though we consider this situation for our results, we do show the k-dependence in our equations for the moment, in order to take the limit of a constant χ c in a correct way through a consideration to be carried out in Sect. 5.
For later reference, we add an additional oscillatory time dependence to the problem, characterised by a frequency ω = 0. Then Eq. 2.1 becomes [17,18] (3.12) and the bulk viscosity is obtained as (3.13)

Scattering matrix elements
In order to solve Eq. 3.12, we need to determine the linearised collision matrix C cc . The relevant processes are shown in Fig. 1. Summing over all indices, the well-known expressions for the matrix elements squared read (cf. e.g. Refs. [19,20]) where s, t, u are the standard kinematic variables. The bulk viscosity is given by Eq. 3.13, where χ c (k 1 ; ω) is determined (up to a discussion concerning kinetic equilibration, cf. Sect. 5) from Eq. 3.12, viz.
Recalling Eq. 3.6 and noting that, because of energy conservation, f B ( p i ) and f F ( p i ) are exponentially suppressed, the basic equation can be reduced to where we have defined an average over the angles between k 1 and k 2 as The integrals in Eq. 4.5 can be carried out in closed form. k 2 ), where Φ b originates from the bosonic final states (the amplitude M 1 in Fig. 1) and Φ f from the fermionic ones (the amplitude M 2 ), the fermionic part reads For In the bosonic case we get In the extreme limit T M/π , the integrals in Eqs. 3.13, 4.4 are saturated by k 1 , k 2 ∼ √ 2T M M. Then we may invoke an approximation in which the energies appearing in the exponentials are set to their non-relativistic forms, , and the energies in the prefactors are expanded to first non-trivial order in k 1 , k 2 . This leads to power-suppressed thermal corrections. As mentioned in the paragraph following Eq. 3.8, in our main results (Sect. 5) we avoid making such an approximation because of its questionable convergence. However, for future reference, let us work out the limiting values in the remainder of the present section.
Setting k 1 , k 2 → 0 where their effects are subleading; making use of the expressions for Φ f (0, 0) and Φ b (0, 0) mentioned above; and carrying out the integral over k 2 , Eq. 4.4 can be solved as where we have defined (following Ref. [3]) Inserting Eqs. 3.6 and 4.9 into Eq. 3.13 yields where the charm quark susceptibility was defined as 3/2 e −M/T . Equation 4.11 agrees with an expression given in Ref. [21] and is used for normalisation below.

Numerical solution
The purpose of this section is to present a numerical evaluation of the final integrals that are needed for determining the bulk viscosity to leading order in e −M/T but to all orders in T /M. As has been discussed in Sect. 3, for theoretically consistent results this is to be done by assuming full kinetic equilibrium, whereby the function χ c is constant. However, a constant value needs to be imposed in a correct way. Allowing first for general k-dependence, we define a quadratic form as The bulk viscosity is formally given as the extremal value: At this point we can restrict the function space to that of constant functions. Then Eq. 4.4 yields, for the observable defined in Eq. 3.13, (5.5) and the weight Φ is from Eqs. 4.6 and 4.8.
In Fig. 2, we plot the results for ζ obtained with the above procedure. The physics conclusion is that when the temperature increases, ζ /ζ 0 , where ζ 0 is from Eq. 4.11, decreases below unity. Our approximation, which assumed e −M/T  Fig. 3 The inverse of chem in physical units, for M = 1.5 GeV and an effective thermal gauge coupling from Ref. [22] (we set MS 360 MeV from Ref. [23]; this leads to α s 0.3 at T /M 0.2 and α s 0.2 at T /M 0.5). For comparison we also plot the result from Ref. [3] where partial lattice input [15,24] has been included. The insert shows the ω-dependence of the transport peak from Eq. 5.2. The axes are normalised to Eqs. 4. 10 and 4.11 small, is reliable at most for T M/π . Nevertheless, it seems conceivable that the solution eventually extrapolates to the result obtained for T M α s T in Ref. [11] (shown with a grey band).
The frequency dependence of the transport peak is plotted in the insert of Fig. 3. Frequency has been normalised to chem from Eq. 4.10; chem is plotted in physical units in Fig. 3. The inverse of the width of the Lorentzian shape serves as an estimate for the microscopic time scale above which the charm contribution to the bulk viscosity plays a role within an effective hydrodynamic description. If the actual time scale characterizing the hydrodynamic expansion of the system is shorter than the microscopic one, a different approach, with the charm density appearing as an additional dynamical variable on par with the temperature and flow velocity, would be needed.

Conclusions
We have shown that in the limit M π T , the bulk viscosity of deconfined QCD plasma grows as ζ 0.04M 4 /(α 2 s T ), cf. Eq. 4.11 and Fig. 2 which shows that Eq. 4.11 provides for a reasonable estimate in this temperature range. Compared with the result ζ 0.25α 2 s T 3 for M α s T [11], we note that massive quarks dominate the QCD bulk viscosity if M 1.6α s T . If the system lived for an infinitely long time, this would be true for charm quarks at all temperatures that can conceivably be reached in current or future heavy ion collision experiments. However, it is very important to keep in mind that dynamical reactions need to take place in order for this contribution to the bulk viscosity to be physically relevant for the hydrodynamical evolution of a system with a finite lifetime. 3 The time scale needed for the bulk viscosity to "build up" can be estimated from the inverse of the width of the transport peak associated with the trace anomaly and is parametrically of the form τ 0.45M 1/2 e M/T /(α 2 s T 3/2 ), cf. Eq. 4.10 and Fig. 3 which suggests that Eq. 4.10 is a conservative estimate in the temperature range considered. The associated processes take place within the lifetime of the plasma generated in heavy ion collision experiments only if a temperature T 600 MeV can be reached for an extended period. For T 400 MeV charm and anticharm quarks constitute separate conserved charges like in a statistical model [27]. In between, a partial chemical equilibration may take place, and a dynamical simulation, perhaps as a further ingredient in an existing one for kinetic equilibration as reviewed in Ref. [28], would be needed for judging the practical importance of charm quarks. For a recent study of the practical effects of a bulk viscosity, see e.g. Ref. [29].
Apart from their effect on bulk viscosity, chemically equilibrated charm quarks would also be relevant for basic thermodynamic quantities such as the equation of state [12][13][14][15]. In addition they could change the shear viscosity by a noticeable amount. 3 This requirement is much easier to satisfy in cosmology than in heavy ion collision experiments, because the expansion rate of the Universe is inversely proportional to the Planck mass and thus much smaller. Correspondingly effects from particle mass thresholds could have a favourable influence e.g. on scalar field friction coefficients directly proportional to the bulk viscosity [25,26].
We end by stressing that obtaining complete leading-order results for T ∼ M necessitates implementing more complicated resummations than have been considered here [11]. The same is also true if we decrease the temperature down to T α 2 s M [30], which is in principle a regime relevant for bottom quarks. Beyond resummed computations, it would also be interesting to address the corresponding physics on the lattice. Comparisons of NLO perturbation theory [21,31] and lattice simulations [24] for imaginary-time correlators suggest that lattice data are already close to the continuum limit, so that future progress may be expected.