Perturbative renormalization of the electric field correlator

The momentum diffusion coefficient of a heavy quark in a hot QCD plasma can be extracted as a transport coefficient related to the correlator of two colour-electric fields dressing a Polyakov loop. We determine the perturbative renormalization factor for a particular lattice discretization of this correlator within Wilson's SU(3) gauge theory, finding a ~12% NLO correction for values of the bare coupling used in the current generation of simulations. The impact of this result on existing lattice determinations is commented upon, and a possibility for non-perturbative renormalization through the gradient flow is pointed out.


Introduction
If a system in thermodynamic equilibrium is displaced slightly by means of some external perturbation, it tends to relax back to equilibrium. The rate at which this happens depends on the nature of the perturbation. Among the simplest perturbations is that the current J µ associated with some conserved particle species is not aligned with the flow velocity u µ of the heat bath: (η µν − u µ u ν )J ν = 0. Then the relaxation rate is called the kinetic equilibration rate; more generally it can be a function of the wave vector associated with the perturbation of J µ , with the rate-proper defined through the long-wavelength limit.
In a hot QCD plasma produced in heavy ion collision experiments, with a lifetime of ∼ 10 fm/c and a temperature of ∼ 300 MeV, heavy (charm and bottom) quarks provide examples of conserved particle species. Indeed their weak decays are too slow by many orders of magnitude to play a role. Moreover the particle and antiparticle number densities are to a good approximation conserved separately, given that pair creations and annihilations are also too slow to take place in a typical event, unless T > ∼ 600 MeV [1].
Heavy quarks are naturally displaced from kinetic equilibrium, given that they are generated in an initial hard process which has no knowledge of the thermal state and the flow that it develops during ∼ 1 fm/c. It is empirically observed, however, that the process of kinetic equilibration takes place during the fireball expansion: heavy quark jets get quenched, and eventually heavy quarks participate in hydrodynamic flow almost as efficiently as light quarks do (cf., e.g., refs. [2][3][4]). This calls for a theoretical computation of their kinetic equilibration rate [5], which can subsequently be inserted into phenomenological models permitting for comparisons with data (cf. ref. [6] for a review).
A large body of theoretical and model-based work has already been carried out on heavy quark kinetic equilibration. Focussing here on the deconfined phase, it has been observed that next-to-leading order (NLO) corrections are large and positive [7,8], and that the AdS/CFT setup suggests rapid equilibration [9][10][11]. However ultimately the problem should be studied with methods of lattice QCD. It turns out that in the heavy-quark limit, applicable at least for the bottom quarks, the lattice challenge appears to be relatively manageable, because the observable in question reduces to a purely gluonic "electric field correlator" (cf. eq. (2.2)) whose associated spectral function is expected to be smoother than for almost any other transport observable [12,13]. Indeed first lattice measurements making use of this observation have been carried out within quenched QCD [14][15][16][17].
There are several important issues which the existing simulations have not solved conclusively. One is taking the continuum limit: for a systematic extrapolation the electric field correlator requires a finite renormalization factor, which should ultimately be determined nonperturbatively. Another is analytic continuation: after the continuum limit has been taken, a short-distance singularity can in principle be subtracted [18] and the result subsequently subjected to a well-defined analytic continuation algorithm [19], however implementing this in practice requires exquisite statistical precision [20], so that in the existing studies theoretically motivated ansätze have been employed for extracting the transport coefficient [15][16][17].
The issue addressed in the present paper is that of renormalization. Specifically we compute the renormalization factor for the electric field correlator at 1-loop order in lattice perturbation theory [21,22], and briefly comment on the possibilities for non-perturbative renormalization. With such ingredients and improved statistical precision, the current orderof-magnitude estimate [17] can hopefully be promoted towards a quantitative level.
The plan of this paper is the following. After defining the basic observable in sec. 2, the technical implementation of the computation is outlined in sec. 3. Results for individual Feynman diagrams are documented in sec. 4. The results are put together in sec. 5 for our final expression for the perturbative renormalization factor. We offer an outlook and comments on non-perturbative renormalization in sec. 6.

Basic definitions
The physical observable underlying our considerations is the 2-point correlator of the time derivatives of the spatial components of the conserved vector current, evaluated in an ensemble at a finite temperature T . After taking the heavy quark-mass limit [12], the correlator reduces to the 2-point function of coloured Lorentz forces [11].
We first define the correlator in continuum notation. Letting U (τ b , τ a ) be a temporal Wilson line at the spatial position r ≡ 0 and defining colour-electric fields through where D µ = ∂ µ + ig B A µ is a covariant derivative and g B is the bare gauge coupling of a dimensionally regularized SU(N c ) gauge theory, the electric field correlator reads [12] G E,cont (τ ) ≡ − 1 3 .

(2.2)
Here β ≡ 1/T denotes the inverse temperature. The global Z(N c ) symmetry of the SU(N c ) gauge theory is assumed to be spontaneously or explicitly broken, so that the denominator of eq. (2.2) is non-zero. The lattice discretization of G E is not unique; we have in mind the proposal from ref. [12] which was argued to have small discretization effects, viz.
where lines indicate parallel transporters, and a is the lattice spacing. We focus on pure SU(N c ) gauge theory in the following, setting N f = 0. In this case the Z(N c ) symmetry is spontaneously broken in the high-temperature phase. Even though the correlator is finite after coupling constant renormalization [12,23], the continuum and lattice-regularized expressions do differ by a finite factor in every order of perturbation theory. In the following we determine this factor at NLO, i.e. at O(α s ). Numerical measurement should be multiplied by this factor (cf. eq. (5.4)) in order to obtain the continuum result. If the spectral function (ρ E,cont ) corresponding to the continuum result can subsequently be extracted, then the corresponding transport coefficient is known as the "momentum diffusion coefficient": The kinetic equilibration rate (also known as the drag coefficient η D ) is given by where M kin is a so-called kinetic mass of the heavy quarks. The single factor κ yields two different kinetic equilibration rates of phenomenological interest, one for the charm and another for bottom quarks. In the charm case 1/Γ kin could conceivably be as small as a few fm/c [17], which could then offer for a partial qualitative explanation for the experimentally observed hydrodynamic flow observed in D meson spectra after hadronization.

Technical outline
Denoting by U µ (X) a link matrix which changes in gauge transformations as whereμ is a unit vector in the µ-direction, the observable of eq. (2.3) can be expressed as 1 For a perturbative analysis, the link matrices are expressed in terms of gauge potentials as where g 0 is the bare lattice gauge coupling and T a are Hermitean generators of SU(N c ), normalized as Tr (T a T b ) = δ ab /2. The gauge potentials can be Fourier-represented as where the sum-integration measure is as defined in eq. (A.2). We use the standard Wilson action, with well-known momentum-space Feynman rules (cf. refs. [21,22]). As mentioned before, eq. (2.3) is unambiguous only when the Z(N c ) center symmetry is explicitly or spontaneously broken. We have carried out the renormalization computation at a high temperature, when the latter is the case. Even though formally justifying the use of the weak-coupling expansion, the finite temperature has the side effect that Matsubara zero modes play a prominent role; indeed, they lead to effects of O(1/a) or O(a) in many individual terms. We have verified the cancellation of these structures, which then also serves as a partial crosscheck of our computation.

Results for individual diagrams
As an illustration of the computational procedure, we give in this section results for the individual Feynman diagrams, both in dimensional regularization and in lattice regularization. In each case only those terms which differ in the two schemes are shown; the full results for dimensional regularization, including the scheme-independent parts, can be found in ref. [23]. For simplicity the calculations were carried out in Feynman gauge; their sum is gauge-independent.
In the expressions below, the common sum-integration "measure" is suppressed. Here i enumerates the spatial directions, D is the space-time dimension, and K = (k n , k) is an imaginary-time four-momentum, with k n denoting a Matsubara frequency. The momentum variables have been so chosen that the components of K can be assumed small compared with the ultraviolet cutoff, i.e. |K| ≪ 1 a , where a denotes the lattice spacing. The renormalized gauge coupling of the MS scheme is denoted by g 2 , and C F ≡ (N 2 c − 1)/(2N c ). The contributions of the different Feynman diagrams read as follows. A renormalization contribution is obtained by expressing the bare continuum and lattice couplings in terms of the renormalized MS coupling. 2 For the lattice coupling this relation was determined in ref. [24] and has been extended up to 2-loop level in ref. [25]. For our purposes, eq. (1.1) of ref. [25] is best re-expressed as andμ is the scale parameter of the MS scheme. The coefficients P 1 , P 2 read where unspecified notation is explained in appendix A.1. In practice P 2 emerges when considering the continuum limit, aK ≪ 1, of the following expression [25]: (4.6) We obtain (in the graphical notation used below, the circular solid line denotes the Polyakov loop to which gauge fields can attach; blobs are the colour-electric field operators; and wiggly lines represent gauge field propagators): As a second set of contributions, we collect together effects originating from expanding the Wilson lines between the electric fields, or in the denominator, to quadratic order. Here two 2 We find it convenient to first present the results in terms of the renormalized MS coupling g 2 , even though at the end of the computation g 2 will be re-expressed through the bare lattice coupling g 2 0 .
new constants, denoted by P 3 and p 1 and defined in eqs. (A.5) and (A.11), appear: The unspecified "finite" terms correspond to structures which are integrable even in the absence of an ultraviolet regulator and therefore agree in the two schemes; they can be extracted from the continuum results presented in ref. [23], but are not shown here because of their complicated appearance.
A simple contribution, only appearing on the lattice side, concerns a "tadpole" correction to the electric fields: There are also two graphs which are finite or vanish, and therefore yield no renormalization contribution: . (4.10) The self-energy contribution, which is in agreement with the classic self-energy computation in ref. [26], reads: The remaining graphs yield: (4.12) cont : finite , cont : finite + k 2 The last contribution turns out to be quite tedious to extract, but given its simple final appearance we refrain from dwelling on more details here.

Determination of the renormalization factor
Summing together eqs. (4.7)-(4.14) all appearances of the "three-dimensional" lattice constant p 1 , defined in eq. (A.11), cancel. We are left with cont : k 2 n 1 + Subtracting the two and noting that in the remaining terms the sum-integration measure in eq. (4.1) implies that we can substitute k 2 with the price of inessential contact terms, we obtain (here also D → 4) This difference needs to be cancelled by a renormalization factor Z E in order for the lattice and continuum results to agree: Eq. (5.3) leads to the next-to-leading order perturbative expression for Z E : In the last step we made use of the non-trivial relation determined in eq. (A.10). This simple expression, obtained through a substantial effort, constitutes our final result. In ref. [15], a provisional result for Z E was given, based on a lattice HQET computation of the type in ref. [27], concerning the renormalization of the spatial components of the Noether current. Our result in eq. (5.5) agrees with the part proportional to C F P 1 in the estimate of ref. [15], however that result has a large additional term proportional to −C F p 1 ; this structure is the same that appears in connection with the mass renormalization of a static quark. We are not surprised by the difference, given that G E does not directly correlate the spatial components of the Noether current but their time derivatives. Taking time derivatives involves a specific discretization, which interferes with the movement of the static quark of mass correction ∼ C F p 1 /a by a distance ∼ a in the time direction. It is non-trivial to account for such effects, and conceivable that something got lost in the HQET estimate.
Inserting a numerical value for the coefficient P 1 (cf. eq. (4.4)) and N c = 3 for the number of colours, and using eq. (4.2) to re-express g 2 through g 2 0 , we obtain Recalling g 2 0 = 6/β 0 , where β 0 is the bare lattice coupling, and noting that values in the range β 0 ∼ 7...8 have been used in simulations attempting to make a continuum extrapolation [17], a rather modest ∼ 12% renormalization effect is found. This can be compared, for instance, with the renormalization factor associated with a local version of the vector current: in that case the O(g 2 0 ) correction [28][29][30] is slightly larger than ours, 3 but nevertheless it dominates the full result including all O(g 4 0 ) effects [31,32]. Therefore we may expect a ∼ 5% theoretical uncertainty from effects of O(g 4 0 ) in the present case. Finally we note that in ref. [17] a provisional version of eq. (5.6) was used in connection with a lattice analysis. Unfortunately that version contained an algebraic error, whereby the numerical value was estimated as Z E ≃ 1 + g 2 0 × 0.079. For β 0 ∼ 7...8 the difference with respect to eq. (5.6) is of a similar magnitude as our estimate for O(g 4 0 ) corrections. However, given that systematic errors were estimated generously in ref. [17] and that they are predominantly associated with analytic continuation from the imaginary-time correlator to a spectral function, the final results should not be significantly affected within their 50% uncertainties. (Note in particular that one of the ansätze considered in ref. [17], denoted by 3a, left the overall normalization of the ultraviolet part of the spectral function open, which roughly corresponds to leaving Z E open.)

Conclusions and outlook
The purpose of this technical contribution has been to report on the result of a 1-loop computation in lattice perturbation theory for the renormalization factor defined in eq. (5.4). The result for the case that the action is discretized according to Wilson's classic prescription and the observable is discretized as specified in eq. (2.3) is given in eq. (5.6). The numerical magnitude of this correction is rather modest, and suggests that future updates of simulations such as those in ref. [17] can attempt approaching the continuum limit, once the results are multiplied by the perturbative renormalization factor that we have determined.
For a definitive lattice result, renormalization should be promoted to a non-perturbative level. One possibility for this is to make use of Yang-Mills gradient flow [35]: following ref. [36], we expect that if the electric fields and Wilson lines between them are computed from gauge fields at a positive flow time t > 0, then the continuum limit can be taken at each t and leads to a finite regularization-independent result. Subsequently the limit t → 0 + needs to be taken, but this should not lead to any divergences since G E,cont is finite [12,23]. A program of a similar type has been suggested for the energy-momentum tensor some time ago [37,38], and numerical tests have been carried out in the meanwhile [39]. Recently gradient flow has also been used for defining renormalized Polyakov loop expectation values [40,41].
Ultimately, for physical QCD, dynamical quarks need to be included in the analysis. The discretization choices made in this context are of a great variety, so computations within lattice perturbation theory become cumbersome and non-universal. Therefore we have not embarked on the inclusion of dynamical quarks; this is probably sensible anyways only in combination with non-perturbative renormalization. with N τ denoting the temporal extent in lattice units. In the zero-temperature limit we denote q n → q 0 , and the measure becomes In the renormalization factors we can always use zero-temperature measures, given that thermal effects are finite and independent of the regularization scheme.

A.2. Lattice constants
Apart from the "four-dimensional" constants P 1 and P 2 defined in eqs. (4.4) and (4.5) and from an analogous "three-dimensional" constant p 1 defined in eq. (A.11), a "mixed" constant also makes an appearance in the analysis: i andQ 2 ≡q 2 0 +q 2 . This plays a role in the sum-integral which in dimensional regularization reads where n B (q) ≡ 1/[exp(βq) − 1] is the Bose distribution. On the lattice it evaluates to This sum-integral appears (only) in the 1st and 2nd diagrams of eq. (4.8).
Remarkably, P 3 as defined in eq. (A.5) is related to the constants P 1 and P 2 as defined in eqs. (4.4) and (4.5). One way to see this is to represent the propagators as 1/∆ = ∞ 0 dt e −t∆ and then to carry out the angular integrals as π −π dθ 2π e t cos θ = I 0 (t), where I 0 is a modified Bessel function. For P 2 and P 3 two variables t 1 , t 2 may be introduced for the two propagators, but we can subsequently substitute variables as x ≡ t 1 + t 2 and y ≡ t 1 . The expression for P 3 contains the integral x 0 dy e −y I 0 (y) = xe −x I 0 (x) + I 1 (x) . (A.9) Inserting I 1 (x) = I ′ 0 (x) yields Apart from the coefficients P 1 , P 2 and P 3 defined in eqs. (4.4), (4.5) and (A.5), the threedimensional lattice integral [33,34]  also appears in the results for the individual graphs. It cancels in the sum, cf. eq. (5.2).

A.3. Frequently needed relations
There are many identities between various lattice integrals that are helpful for simplifying the computation; a collection sufficient for our purposes can be found in appendix B of ref. [25]. Let us here just recall two examples that appear particularly often:

A.4. Thermal sums
A basic thermal sum appearing is T qn e iqnτ q 2 n + ω 2 = a e kx + e (Nτ −k)x 2 e Nτ x − 1 sinh x , x = 2 asinh aω 2 , k N τ = τ β . (A.13) Methods for carrying out sums of this type were discussed in ref. [42]. From this sum others can be obtained through the limits ω → 0 and/or τ → 0 and through the omission of the Matsubara zero mode, in which case we denote the variable by q ′ n . In particular,