The errant life of a heavy quark in the quark-gluon plasma

In the high-temperature phase of QCD, the heavy quark momentum diffusion constant determines, via a fluctuation-dissipation relation, how fast a heavy quark kinetically equilibrates. This transport coefficient can be extracted from thermal correlators via a Kubo formula. We present a lattice calculation of the relevant Euclidean correlators in the gluon plasma, based on a recent formulation of the problem in heavy-quark effective field theory (HQET). We find a $\approx20%$ enhancement of the Euclidean correlator at maximal time separation as the temperature is lowered from $6T_c$ to $2T_c$, pointing to stronger interactions at lower temperatures. At the same time, the correlator becomes flatter from $6T_c$ down to $2T_c$, indicating a relative shift of the spectral weight to lower frequencies. A recent next-to-leading order perturbative calculation of the correlator agrees with the time dependence of the lattice data at the few-percent level. We estimate how much additional contribution from the $\omega\lesssim T$ region of the perturbative spectral function would be required to bring it in agreement with the lattice data at $3.1T_c$.


Introduction
Imagine a heavy quark moving through the plasma of light quarks and gluons with more than its fair share of thermal kinetic energy 3 2 k B T . Interactions with the medium will gradually cause it to slow down, at a rate determined by the drag coefficient η. It also experiences stochastic interactions that balance the drag force in such a way that the heavy quark's kinetic energy tends to 3 2 k B T . In this paper, we study from first principles the strength of these stochastic interactions, which is characterized by a parameter called the momentum diffusion coefficient and denoted by κ. This is a quantity that is of phenomenological interest. Measurements at the Relativistic Heavy Ion Collider (RHIC) have shown that heavy quarks display substantial elliptic flow [1,2], implying stronger medium interactions than extrapolated weak-coupling calculations would suggest. It has been estimated [3] that a diffusion coefficient D 1/T , or equivalently κ/T 3 2, is required in order to accommodate the RHIC data. Theoretically, we will be concerned with collisional energy loss, which is expected to dominate for small P/M; however, one should keep in mind that in heavy-ion collisions, the quarks are not produced at rest; hence, radiative energy loss can play a role too (see [4] for a recent review). A separate motivation is that the dynamics of the heavy-quark probe are somewhat simpler to study than the diffusion of transverse momentum carried by the constituents of the plasma, which is characterized by shear viscosity. There are reasons for expecting that, once the heavy-quark diffusion constant has been determined, it can be used to calibrate other transport coefficients, such as shear viscosity, because the ratios of transport coefficients are more stable predictions of the weak-coupling expansion [5].
We briefly review the known analytic results for the momentum diffusion coefficient. To strict leading order, κ is given by the following equation [5] Unfortunately, for realistic values of the Debye screening mass m D , this expression is negative. Carrying out the integrals without expanding in m D (m 2 D = g 2 T 2 (N c /3 + N f /6) in leading order), positive values are obtained, for instance at g 2 ≈ 3, κ/T 3 ≈ 0.4 for N f = 3 but 0.23 for N f = 0 [5]. Remarkably, a next-to-leading order (NLO) computation has been carried out [6], with the result κ/T 3 ≈ 2.7 for g 2 ≈ 3.0 for N f = 3. It is unclear at this point how useful the expansion is.
The N = 4 super Yang-Mills (SYM) theory provides an interesting testing ground for analytic methods. The result obtained by holographic methods in the strong-coupling limit of 3,7,8], where λ ≡ g 2 N c is the 't Hooft coupling. On the other hand, the NLO weak-coupling result has also been worked out [9], Although the apparent convergence is again poor, it is encouraging that the strong-coupling result crosses the weak-coupling result at an intermediate coupling of λ ≈ 10.2. This suggests that the values obtained from NLO expressions such as (2) are in the right ball-park 1 .

3
In the next section, we review the steps that lead to a Kubo formula for the momentum diffusion coefficient κ in the static limit of the probe-quark. In section 3, we present our numerical results for the Euclidean correlators obtained by Monte-Carlo simulations. We conclude the paper in section 4.

The spectral function of the heavy-quark current
This section is mainly a review of the recent literature on the subject of heavy-quark diffusion, with some comments pertinent to calculations performed in Euclidean space. The basic definitions of thermal correlators and the relations among them are gathered in the appendix; for unexplained notation, see the appendix. The main equations of the linear response framework are also reviewed there.
To study the diffusion of a conserved charge such as the heavy-quark number N = dx n(t, x), the Hamiltonian can be perturbed by In the hydrodynamic treatment of the problem, Fick's law j = −D∇n and the conservation equation ∂ t n + ∇ · n = 0 lead to the diffusion equation, whose solution in Fourier space takes the formñ Here, χ s (k) = β dx e −ikx n(t, x)n(0) and χ s ≡ χ s (0), is the particle number susceptibility. Via (A.20), this determines the retarded correlator G nn R (ω) for small ω and k, Before proceeding further, it is also instructive to write down the correlator in the real-time domain, Using equation (A.10), the contribution of this exponential tail to the Euclidean correlator (say) at t = β/2 is χ s (k), independently of the diffusion coefficient D. This is a manifestation of the difficulty of extracting transport information from Euclidean correlators. The longitudinal part of the current correlator (i.e. j z j z if k = (0, 0, k)) is related to the density correlator by the current conservation equation, k). Thus the current spectral function reads implying in particular the Kubo formula, The diffusion of a heavy quark (one has in mind the charm or preferably, from a theoretical point of view, the bottom quark) in the quark-gluon plasma is characterized by a time scale 4 M/T 2 , which is long compared with the thermal time scale of 1/T . For this reason, it is expected that a classical Langevin equation should appropriately describe the thermalization of heavy quarks [5]. See [11] and references therein for a derivation. The heavy quark's classical equations of motion are For a given ξ(t), the equation is easily solved to give The equipartition of energy requires p 2 2M to be 3 2 T in equilibrium; the drag and fluctuation coefficients are thus related by the fluctuation-dissipation relation The mean square distance covered by the particle is also easily worked out. For a thermal initial distribution of momenta, This equation describes both the early-time directed motion, 1 3 x 2 (t) = 1 3 v 2 t 2 , 1 3 v 2 = T M , and the late-time diffusive motion, 1 3 x 2 (t) = 2Dt [12]. Let now P(t, x) be the probability that a heavy quark starts at the origin at t = 0 and moves a distance x over a time t. If the distribution of heavy quarks at time zero is n(0, x), at time t it will be given by the convolution or equivalently If one assumes the noise to be Gaussian distributed, then the probability distribution P(t, x) is Gaussian [12], with a width given by (14), and therefore so is P(t, k). Applying the general rule (A.20), one then obtains the retarded correlator The resulting spectral functions were obtained numerically in [12]. In particular, the k = 0 spectral function takes the form The spectral structure that is obtained from the Langevin equation is expected to arise for a sufficiently heavy diffusing particle. Conversely, the presence of a transport peak allows one to define the quantities appearing in the Langevin equation directly from the spectral function [13]. The effective mean square velocity is then given by where is a cutoff that separates the scale η from the correlation time of the medium (which is typically of order T , or gT at weak coupling). The 'kinetic mass' M kin is further defined so as to satisfy the equipartition theorem, M kin v 2 = 3T . Finally, the momentum diffusion coefficient κ(M) can be defined as A weak-coupling calculation shows that while κ(M) and D are only weakly dependent on M, the drag coefficient η ∼ T 2 /M×a power of the coupling constant is parametrically small compared to the medium time scale. By re-expressing the right-hand side of (20) in terms of the correlator of two heavy-quark currents, the authors of [3,13] were able to formulate the task of computing κ in the static limit of heavy-quark effective theory (HQET). In that limit, the leading contribution is the 'force-force' correlator, where the force is given by the Lorentz expression g E. The relevant Euclidean correlator reads, after evaluating the fermion line contractions, where the color parallel transporters U (t 2 , t 1 ) in the fundamental representation are propagators of static quarks. In particular, the Polyakov loop appears in the denominator of (21). The momentum diffusion coefficient is given by the low-frequency limit of the corresponding spectral function via equation (A.12), The obvious advantage of this formulation is that the large scale M has disappeared from the problem. The spectral function ρ HQET has been studied in detail at NLO in perturbation theory [14]. Remarkably, even in the weak-coupling limit, the function is smooth as small frequencies. This is in contrast with the narrow transport peaks that are found at weak coupling in e.g. the shear channel. This property represents a clear advantage for numerical studies of the spectral function: the form of the kernel in equation (A.12) makes the Euclidean correlator very insensitive to the functional form of the spectral function at ω T . On the downside, while the spectral function of the current-current correlator grows as N c 12π 2 ω 2 at large frequencies, the spectral function of the E-field correlator grows even faster, like g 2 C F 6π 2 ω 3 [13]. This implies that the low-frequency part makes a comparatively small contribution to the Euclidean correlator studied in the next section.

Force-force correlators from lattice gauge theory
In this section, we describe a first calculation of the HQET 'force-force' correlator (21) by means of Monte-Carlo simulations in Euclidean space. The calculation is performed in the 6 deconfined phase of SU(3) gauge theory. We employ the isotropic Wilson action [15], where the 'plaquette' P µν is the product of four link variables U µ (x) ∈ SU(3) around an elementary cell in the (µ, ν) plane. The size of the lattice is N τ × N 3 σ , with periodic boundary conditions in all directions. As a local update algorithm, we use the standard combination of heatbath and over-relaxation [16]- [19] sweeps in a ratio increasing from 5 to 7 as the lattice spacing is decreased. No multi-level algorithm [20,21] was used here, although we expect that for sufficiently small lattice spacing, such an algorithm will be beneficial, since the fluctuations of the E fields will be UV dominated. To set the scale, we use the parameterization of the Sommer scale r 0 ≈ 0.5 fm given in [22] based on the data [23], and convert between r 0 and T c by using r 0 T c = 0.746(7) ( [24] and references therein).
There is a lot of freedom in discretizing the HQET Lagrangian (see [25] for an introduction). As color parallel transporters, we choose the lattice gauge links obtained after one iteration of HYP smearing [26]. This action has been used extensively at zero temperature [27], and our choice of smearing parameters is known in the literature as the HYP1 action [28]. The original motivation for choosing this action is that it reduces the size of the UV-divergent selfenergy, thus leading to a reduction in statistical errors. At the same time, the cut-off effects were found to be controllable. In [29], the chromo-magnetic field B was also discretized with HYP-link variables, yielding also a benefit in statistical error reduction. Following this example, we use HYP-links for the chromo-electric field as well.
While the renormalization factor for the B field was computed in [29] (it is scale dependent), the corresponding factor for the chromo electric field remains to be computed. The same methods apply to the renormalization of this factor, but in the meantime we use a preliminary way of normalizing the chromo electric operator. We expect that at short time separations t, perturbation theory provides an accurate prediction for G E (t). Therefore the absolute normalization of the chromo electric field can be obtained by requiring that the lattice correlator match the NLO perturbative prediction at some reference time t ref . The difficulty of this procedure is that t ref must still be large enough in lattice units for the discretization errors to be under control. In this work, we made the compromise to choose t ref = 1/4T at a temperature of about 3.1T c on a 16 × 64 3 lattice. In physical units, this represents a separation of about 0.05 fm, and four lattice spacings in lattice units. At one value of the bare coupling g 2 0 , this normalization is then valid for other values of N τ . In this way, we have varied the temperature by changing N τ from 8 to 22 at fixed 6/g 2 0 = 7.483. The chromo electric field correlators are displayed in figure 1. Only those points with t/a 4 are displayed. At shorter separations, the lattice correlators exhibit non-monotonicity (and are even negative at separation 0 and 1 lattice spacing), a not unexpected lattice artefact. It is clear that with the procedure adopted here, we cannot disentangle discretization errors from the bare-coupling dependence of the renormalization factor of E. Our results should accordingly be regarded as preliminary. Eventually, the renormalization factor should be computed along the lines of [29], where the B operator was treated instead.
To repeat, the data at 3.1T c have been calibrated to match the O(g 4 ) result of [14] at t = 1/4T . In figure 1, different temperatures are accessed by varying N τ . On the logarithmic scale of the plot, the temperature dependence of the Euclidean correlator normalized by T 4 is weak. Furthermore, the NLO perturbative prediction provides a rather good description of the t dependence at 3.1T c . The temperature dependence of the mid-lattice data points (normalized as  [14]. This result has been used to determine the normalization of the E field by matching the lattice data to it at t T = 1 4 and 3.1T c . The statistical error bars are smaller than the data symbols. in figure 1) is shown in much greater detail in figure 2. The temperature variation between 2T c and 6T c is of the order of 20%. The sign of the variation is the same as predicted by the perturbative expression: the magnitude of the force-force correlator in the gluon plasma increases as the temperature decreases. Although we do not yet have non-perturbative control over the absolute normalization of the correlators, the latter cancels out in the relative fall-off of the correlator. An observable that  (24). The '+' denotes the O(g 4 ) prediction at 3T c [14]. The '×' denotes the function (t) that results from adding the low-frequency contribution ρ(ω) = 1 π κ tanh( ω 2T )θ ( − |ω|) to the O(g 4 ) spectral function, with parameters κ and given in the caption. The statistical error bars are smaller than the data symbols, but discretization errors are probably non-negligible, particularly at smaller t values.
measures this is the quantity (t) 0 defined by We remark that (t) has a continuum limit, in which tanh (β/2 − t) = − d dt log G E (t). One can interpret it as the location of a delta function in the spectral function, which by itself reproduces the local fall-off of the Euclidean correlator. The function (t) is displayed in figure 3 for three different temperatures, where the temperature is varied this time by changing the bare coupling g 2 0 at fixed N τ = 16. We note that because the statistical samples of the numerator and the denominator in equation (24) are highly correlated, the numerical results for these ratios have uncertainties at the few-per mille level. Figure 3 also shows the perturbative prediction at 3.1T c corresponding to the curve appearing in figure 1; it is based directly on equation (24) rather than on the continuum version of this equation, to allow for a more direct comparison with the lattice data. The lattice (t) differs from the perturbative one only by a few per cent, namely the latter falls off slightly more steeply. Since the difference is very small, it could partly be due to discretization effects. To reduce their influence, we concentrate on the largest values of t. One may ask, nonetheless, how large a difference in the transport coefficient κ could this discrepancy possibly correspond to. An estimate is obtained by adding a low-frequency correction to the perturbative spectral function, Other functional forms than (25), such as a Breit-Wigner curve, would perhaps be more realistic, but would not change our conclusions in any significant way. Adding such a term to the spectral function has the effect of making the Euclidean correlator flatter, and indeed, 9 by adjusting κ, one can obtain good agreement for the largest two t values between the perturbative prediction modified by equation (25) and the lattice data. At T = 3.1T c , we find that, for = T and κ/T 3 = 0.352(38), agreement is obtained with the lattice data at the two largest t values. This represents a substantial enhancement of κ over the leading-order perturbative value mentioned in the introduction. An equally good agreement is obtained if one chooses = 2T , which leads to κ/T 3 ≈ 0.204 (22).

Conclusion
To summarize, we have found that the Euclidean force-force correlator evaluated at t = β/2 admits a ≈ 20% increase as the temperature is lowered from 6 to 2T c . At 3.1T c its t dependence is described at the few-per cent level by the recent NLO perturbative result [14]. The somewhat flatter behavior seen in the lattice data can be explained by an enhancement of the spectral function at low frequencies, and adopting this explanation leads to a substantial increase in κ over the leading-order result obtained in [5]. While it is too early to draw phenomenological conclusions, at present the increase appears to be not quite sufficient to explain the experimentally observed elliptic flow of heavy quarks, as discussed in the introduction. It will be interesting to see the results brought by the current lead-lead collisions at the LHC. The ALICE experiment has very recently reported [30] an integrated elliptic flow (i.e. of light constituents) about 30% larger than at RHIC. From a technical point of view, it is encouraging that precision at the few-per-mille level can be achieved on the force-force correlator on a 16 × 64 3 lattice. Given this statistical precision, it is essential, as a next step, to control the discretization errors at a comparable level. One could then study more systematically and quantitatively the implications of the lattice data for the spectral function. A mandatory step in controlling the discretization errors is to compute the renormalization factor of the E field non-perturbatively. The motivation to perform these calculations is now quite strong. play a particularly important role in finite-temperature physics. The integral transform over the positive half-axis is analytic in the half-complex plane Im(ω) > 0. We will refer to it as the retarded correlator. The spectral function, which is really a distribution, is defined as For B = A † , the spectral function is identically related to the imaginary part of the retarded correlator, The Euclidean correlator is defined as It obeys the Kubo-Martin-Schwinger relation, and therefore admits representation as a Fourier series, where ω = 2π T and we have dropped the label specifying the operators A and B.
In frequency space, the Euclidean and retarded correlators are related by The = 0 case has to be treated somewhat more carefully (see [31] for details). In coordinate space, the relation between the Euclidean correlator and the real-time correlator is .
Finally, the most commonly used relation between Euclidean and real-time correlators is the mixed coordinate-frequency space relation, The retarded correlator G AB R is important because it is related to the response of operator A to a perturbation of the system by operator B. One then finds that to linear order in f , the expectation value of A in the perturbed system is Equation (A.16) is the master formula of the linear response theory. It shows that the retarded correlator G AB R determines the response of an observable A to a time-dependent external field that couples to B. A source term of the form is often adopted to study how system relaxes back to equilibrium after having been perturbed adiabatically. The static susceptibility is defined as the expectation value of A at t = 0, This formula shows that the relaxation of the observable A back to its equilibrium value and the retarded correlator G AB R (ω) are in one-to-one correspondence. Since the late-time relaxation is expected to be described by hydrodynamic evolution, this equation can be exploited to obtain a prediction of the small-ω functional form of G R .