1 Introduction

{centering}

Massive vector current correlator in thermal QCD

Y. Burnier and M. Laine

Institute for Theoretical Physics, Albert Einstein Center, University of Bern,
Sidlerstrasse 5, CH-3012 Bern, Switzerland

Abstract

We present an NLO analysis of the massive vector current correlator at temperatures above a few hundred MeV. The physics of this correlator originates from a transport peak, related to heavy quark diffusion, and from the quark-antiquark threshold, related to quarkonium physics. In the bottom case both can be studied with separate effective theories, but for charm these may not be accurate, so a study within the full theory is needed. Working in imaginary time, the NLO correlator can be computed in unresummed perturbation theory; comparing with lattice data, we find good agreement. Subsequently we inspect how non-perturbative modifications of the transport peak would affect the imaginary-time correlator. The massive NLO quark-number susceptibility is also contrasted with numerical measurement.

October 2012

## 1 Introduction

Heavy (charm and bottom) quarks are excellent probes for the properties of the hot QCD plasma generated in heavy ion collision experiments. On the theoretical side, the existence of a mass scale  200 MeV renders the heavy quarks relatively tractable, permitting for an interpolation between the simple static dynamics of the infinite mass limit and the high mobility case manifested by lighter quarks. It is particularly fortunate that two heavy flavours are available, offering for a handle on the functional dependence on . On the experimental side, heavy quarks and quarkonia are readily tagged because of their distinctive leptonic decays. Indeed thermal modifications of the bottomonia spectra were among the first spectacular results produced by the LHC heavy ion program [1].

For the bottom quark case, recent years have seen significant progress in theoretical studies of the main phenomena involved, namely single quark “transport” (diffusion, kinetic equilibration) as well as physics near the quark-antiquark threshold (quarkonium dissociation, chemical equilibration) (cf. ref. [2] for a review and refs. [3][12] for recent contributions). Largely this progress has been achieved through the use of modern effective field theory methods (Heavy Quark Effective Theory, or HQET, for single quarks; Non-Relativistic QCD, or NRQCD, for quarkonium physics). Once properly formulated the effective field theory observables can be measured non-perturbatively with lattice Monte Carlo methods, and indeed first results suggest that these avenues may lead to substantial progress [13][17].

In the charm quark case, however, it is not guaranteed that the heavy quark expansion converges fast enough to yield quantitatively accurate results for all observables of interest. On the other hand, it no longer appears prohibitively expensive to treat charm quarks as “light” degrees of freedom in lattice QCD. Indeed, results have appeared concerning both thermodynamic quantities [18, 19, 20] and imaginary-time correlators relevant for determining dynamical properties of the system [21] (earlier works can be found in refs. [22][25] and references therein). Yet, a modest scale hierarchy between and does exist, and consequently systematic errors, both from lattice artifacts and from the unavoidable analytic continuation [26], are likely to be harder to control than in the light quark case. In fact, given that systematic errors related to analytic continuation remain substantial even for light quarks [27], further crosschecks appear welcome.

The goal of the current study is to derive results within next-to-leading order (NLO) perturbation theory which may help in the analysis of lattice data such as those in ref. [21]. In order to allow for a direct appreciation of Euclidean measurements, we inspect how specific modifications of transport properties and threshold features manifest themselves in imaginary time. Ultimately, in accordance with the philosophy of ref. [28] and practical tests of refs. [27, 29], the goal would be to subtract “trivial” ultraviolet features from continuum-extrapolated lattice data, in order to allow for a model-independent extraction of real-time physics [30].

On the technical side, the current work represents a continuation of our earlier study [31], in which the massive vector current spectral function was computed at NLO in the domain (keeping only those thermal effects which are not exponentially suppressed). Here we keep the full mass dependence, permitting for an extrapolation also to the regime , as well as contributions from the transport peak at which were omitted in ref. [31]. Moreover we work directly in imaginary time, which has the benefit that the usual problems of convergence at small are milder. In particular, at NLO infrared safe results can be obtained without resummations, similarly to what has previously been achieved for gluonic observables [32, 29].

The plan of this paper is the following. After specifying the observables considered and discussing the methods employed (sec. 2), we outline the qualitative structure of our findings in sec. 3. The detailed analytic and numerical results of the strict NLO analysis comprise sec. 4, whereas in secs. 5 and 6 the effects of non-perturbative modifications of the transport peak and quarkonium threshold, respectively, are inspected. Sec. 7 presents our conclusions; appendix A results for all the “master” sum-integrals considered; and appendix B details related to renormalization.

## 2 Observables and methods

### 2.1 Basic definitions

The main quantity considered is the vector current correlator related to a massive flavour. Like in lattice QCD we work in Euclidean signature, with the usual thermal boundary conditions imposed across the time direction. Then the correlator is defined as

 G\tiny\rm{V}(τ) ≡ d∑μ=0∫x⟨(¯ψγμψ)(τ,x)(¯ψγμψ)(0,0)⟩T (2.1) ≡ G00(τ)−Gii(τ),0<τ<1T, (2.2)

where the Dirac matrices are Minkowskian, and a sum over spatial indices is implied. Because of current conservation the charge correlator is independent of , and we denote the corresponding “susceptibility” by

 χ≡∫β0dτG00(τ)=βG00(0),β≡1T. (2.3)

For future reference let us also record the free massless results for these correlators [33]:

 G\scriptsize freeii(τ) ≡ 2NcT3[π(1−2τT)1+cos2(2πτT)sin3(2πτT)+2cos(2πτT)sin2(2πτT)+16], (2.4) χ\scriptsize free ≡ NcT23. (2.5)

Here refers to the number of colours. Later on the group theory factor will also appear. Spacetime dimension is denoted by .

At NLO, we find it convenient to compute the correlators and . The spatial part is then obtained from eq. (2.2). Verifying explicitly the -independence of provides for a nice crosscheck of the computation.

### 2.2 Wick contractions for the vector current correlator

The Wick contractions for are as given in ref. [31]. Denoting

 Q≡(ωn,0),ΔP≡P2+M2, (2.6)

 = 2CA∑∫{P}{(D−2)Q2−4M2ΔPΔP−Q−2(D−2)ΔP}.

Here denotes fermionic Matsubara momenta, and .

At NLO, we have to decide on a meaning of the renormalized mass. Although conceptually subtle, it is technically convenient to employ a pole mass; then the bare mass parameter, , can be expressed as where at NLO

 δM2 = (2.8) = −g2CF∫k12ϵkEpk[(D−2)(Epk−ϵk)+4M2(ϵk+Epk)(ϵk+Epk)2−E2p] (2.9) = −6g2CFM2(4π)2(1ϵ+ln¯μ2M2+43). (2.10)

Here is the scale parameter of the scheme, terms of were omitted, and

 ϵk≡|k|,Ep≡√p2+M2,Epk≡√(p−k)2+M2. (2.11)

Because of Lorentz invariance the vector in eq. (2.8) can be chosen at will, as long as we set after carrying out the -integral; this means that the vector in eq. (2.9) is arbitrary. The corresponding counterterm contribution reads

 = 4CAδM2∑∫{P}{D−2Δ2P−2ΔPΔP−Q+4M2−(D−2)Q2Δ2PΔP−Q}, (2.12)

and it is often convenient to identify of eq. (2.9) as the integration variable of eq. (2.12).

The “genuine” 2-loop graphs amount to

 (D−2)2K2Δ2P−(D−2)2Δ2PΔP−K−2(D−2)K2ΔPΔP−K+4(D−2)M2K2Δ2PΔP−K −2(D−2)K2ΔPΔP−Q+4(D−2)M2−(D−2)2Q2K2Δ2PΔP−Q+2(D−2)K2ΔP−KΔP−Q +4(D−2)ΔPΔP−KΔP−Q−4(D−2)M2−(D−2)2Q2Δ2PΔP−KΔP−Q −16M2+2(D−2)2K⋅Q−4(D−2)Q2K2ΔPΔP−KΔP−Q+16M4−4(D−2)Q2M2K2Δ2PΔP−KΔP−Q −2(D−4)M2+12(D−2)(8−D)Q2+(D−2)K2ΔPΔP−KΔP−QΔP−K−Q +8M4−2(D−4)M2Q2−(D−2)Q4K2ΔPΔP−KΔP−QΔP−K−Q}. (2.13)

### 2.3 Wick contractions for the susceptibility

In the case of the zero components, viz. , the LO correlator reads, in momentum space,

 = 2CA∑∫{P}{2ΔP−Q2+4E2pΔPΔP−Q},

whereas the counterterm graph can be expressed as

 = 4CAδM2∑∫{P}{Q2+4E2pΔ2PΔP−Q−2ΔPΔP−Q−1Δ2P}.

The genuine 2-loop graphs amount to111The appearance of the “energy variables” in the numerators implies that there is a certain redundancy in the basis: (2.16) We have used this to eliminate terms containing in the numerator.

 −D−2K2Δ2P+D−2Δ2PΔP−K+2K2ΔPΔP−K−4M2K2Δ2PΔP−K −2(D−2)K2ΔPΔP−Q+(D−2)(4E2p+Q2)K2Δ2PΔP−Q−2K2ΔP−KΔP−Q +4(D−2)ΔPΔP−KΔP−Q−(D−2)(4E2p+Q2)Δ2PΔP−KΔP−Q −16M2−4K⋅Q+4Q2K2ΔPΔP−KΔP−Q+4M2(4E2p+Q2)K2Δ2PΔP−KΔP−Q −(D−2)(E2p+E2pk−ϵ2k+K2+12Q2)−4M2ΔPΔP−KΔP−QΔP−K−Q +4M2(E2p+E2pk−ϵ2k)+2Q2(E2p+E2pk+M2)+Q4K2ΔPΔP−KΔP−QΔP−K−Q}. (2.17)

### 2.4 Matsubara sums and spatial integrals

The next step is to convert the momentum-space expressions to configuration space. If we denote the above result (eqs. (2.2), (2.12), (2.13)) by , then the conversion is obtained as

 G\tiny\rm{V}(τ)=T∑ωne−iωnτ~G\tiny\rm{V}(ωn), (2.18)

and similarly for . At NLO we are thereby faced with a three-fold Matsubara sum. Making use of standard techniques, reviewed in some detail in ref. [31], these sums can be carried out in a closed form, whereby we are left with integrals over at most two spatial momenta (i.e. two radial directions and one angle). Intermediate results at this stage are displayed for all individual master sum-integrals in appendix A, and for their sums in appendix B, in the latter case with the cancellation of -divergences verified as well.

Two interesting crosschecks are available. First of all, all -dependent terms disappear from , cf. eq. (B.3). Second, individual parts of the expressions contain “contact terms” , where denotes the periodic Dirac-delta. These arise from structures that are independent of , but also from sum-integrals in which appears in the numerator. It can be verified, however, that all contact terms cancel, both at LO and at NLO.

It remains to carry out the spatial integrals. We write

 ∫p,k=∫p,k∫z, (2.19)

where the normalization of the angular variable is chosen so that . Depending on the case, it is convenient to substitute variables as

 ∫z=12pk∫E+pkE−pkdEpkEpk=12pk∫ϵ+pkϵ−pkdϵpkϵpk, (2.20)

where , . The angular integrals are doable in most cases. In addition, it is also possible to carry out partial integrations with respect to the radial directions, which helps to reduce the number of independent terms (for the massless case, see ref. [34] for a recent discussion). Closed massless loop integrals are typically solvable, but in general the integrations remain to be carried out numerically. Our final expressions are given in the next two sections, cf. eqs. (3.2)–(3.5), (4.1), (4.4), (4.5).

## 3 Leading order results and qualitative pattern

In order to illustrate the qualitative structure of the results, we recall in this section the LO expressions for the quantities considered. In general, two types of contributions appear: those that depend on , and those that are constant. To display the -dependence we introduce the periodic dimensionless function

 DEk+1⋯EnE1⋯Ek(τ)≡e(E1+⋯+Ek)(β−τ)+(Ek+1+⋯+En)τ+e(E1+⋯+Ek)τ+(Ek+1+⋯+En)(β−τ)[eβE1±1]⋯[eβEn±1]. (3.1)

In the denominator, the sign is chosen according to whether the particle is a boson or a fermion (at NLO, there is only one boson, with the on-shell energy denoted by , cf. eq. (2.11)). We also denote . With these conventions, the LO result for reads

 G\tiny\rm{LO}\tiny\rm{V}(τ)∣∣∣\tiny\rm{τ-dep.} = −G\tiny\rm{LO}ii(τ)∣∣∣\tiny\rm{% τ-dep.}=−2CA∫p(2+M2E2p)D2Ep(τ), (3.2) G\tiny\rm{LO}\tiny\rm{V}(τ)∣∣∣\tiny\rm{const.} = −4CA∫pM2Tn′\tiny\rm{F{}}(Ep)E2p. (3.3)

Here denotes the Fermi distribution ( denotes the Bose distribution). The susceptibility only contains a -independent part:

 Tχ\tiny\rm{LO}=G\tiny\rm{LO}00=−4CA∫pTn′\tiny\rm{F{}}(Ep). (3.4)

Finally, the constant part of the spatial correlator is obtained through the use of eq. (2.2):

 G\tiny\rm{LO}ii(τ)∣∣∣\tiny\rm{const.}=−4CA∫p(1−M2E2p)Tn′\tiny\rm{F{}}(Ep). (3.5)

To our knowledge none of these leading-order expressions can be integrated in terms of standard elementary functions.

In terms of the spectral function, viz. , the constant contribution in eq. (3.5) arises from (an infinitely narrow) transport peak around , whereas the “fast” -dependence in eq. (3.2) originates from the quark-antiquark continuum at . Around the middle of the Euclidean time interval, . Therefore, both terms contribute in a comparable manner at large Euclidean time separations, even though they originate from completely different types of physics (see also ref. [35]).

At NLO, the same three structures appear as at LO: a constant , as well as a spatial correlator which has both a constant and a -dependent part. It is generally believed that the perturbative series for the constant part of breaks down at some point and that in the full theory, has no projection to the Matsubara zero mode; the reason is that the spatial components of the vector current do not couple to conserved charges. The corresponding “slow” -dependence then reflects the physics of a “smeared” transport peak. Yet it remains true that single quark physics from the transport peak and quark-antiquark physics from the threshold region are expected to contribute in a comparable manner to the Euclidean correlator around . (This is demonstrated explicitly in figs. 3, 4 below.)

## 4 NLO results

### 4.1 Susceptibility

Proceeding to NLO, the result for the susceptibility is obtained from eq. (B.3) after partial integrations:

 Tχ\tiny\rm{NLO}=G\tiny\rm{NLO}00 = 4g2CACF∫pTn′\tiny\rm{F{}}(Ep)p2∫k[n\tiny\rm{B{}}(ϵk)ϵk+n\tiny\rm{F{}}(Ek)Ek(1−M2k2)]. (4.1)

This can be shown to agree with what can be extracted from ref. [36] as . (The substance of the information is already there in ref. [37].) The massless loop evaluates to , and in the limit its effect can be interpreted as an effective mass correction to the LO result of eq. (3.4),  [38]. On the other hand, in the massless limit eq. (4.1) reduces to the well-known correction (cf. [39] and references therein)

 G\tiny\rm{NLO}00\lx@stackrelM≪πT=−g2NcCFT38π2+O(g3). (4.2)

A numerical evaluation of eqs. (3.4), (4.1) is shown in fig. 1. For the gauge coupling we have inserted the 2-loop value for the thermal coupling from ref. [40], and for comparison with quenched lattice data from refs. [41, 21] we assume that , . In the massive case, where no continuum extrapolation is available in ref. [21], we choose to compare with results obtained on a lattice . The quark mass cited,  GeV is phenomenologically on the low side [42]; based on eq. (2.10), viz.

 M≃m\tiny\rm{¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\rm MS}(mc){1+4g2\tiny\rm{¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\rm MS}(mc)CF(4π)2+O(g4)}, (4.3)

we get  GeV but assign a large error to this number. Ref. [21] also cites ,  MeV, so we estimate , but with substantial systematic uncertainties, from lattice artifacts, string tension measurement [43], quenching, as well as perturbative input. With this in mind, the excellent agreement seen in fig. 1 is remarkable, and supports the long-held belief that all quarks, and heavy quarks in particular, are well described by the weak-coupling expansion at surprisingly low temperatures. (In the massless case the good agreement is consistent with the recent unquenched study of ref. [44], in which a similar “dimensionally reduced” resummation scheme was used as here [36], however an almost perfect match already at NLO may be somewhat coincidental.)

### 4.2 Vector current correlator

The vector current correlator is obtained from eqs. (B.1), (B.2) after partial integrations. Its -dependent part reads

 (4.4) + ∫p,k\mathbbmP{∫zDϵkEpEpk(τ)ϵkEpEpkΔ+−Δ−+[−2E2p+M2(Δ−−Δ+++2ϵ2kΔ+−Δ−+)+4ϵ2kM4Δ2++Δ+−Δ−+] +∫z2DEpϵkEpk(τ)ϵkEpEpkΔ++Δ−−[E2p+E2pk−M2(Δ+−Δ−++2ϵ2kΔ++Δ−−)−4ϵ2kM4Δ2−+Δ++Δ−−] + D2Ep(τ)2ϵ3k(2+M2E2p)[−1+E2p(E+pk−E−pk)−pϵk(E+pk+E−pk)2p(E2p−ϵ2k)+ϵ2kM2(E+pk−E−pk)p(E2p−ϵ2k)E+pkE−pk +2E2p−M22pEp(ln∣∣∣(Ep+p)(2p−ϵk)(Ep−p)(2p+ϵk)∣∣∣+ln∣∣∣1−ϵ2k/(Ep+E−pk)21−ϵ2k/(Ep+E+pk)2∣∣∣)+θ(k)] + D2Ep(τ)n\tiny\rm{B{}}(ϵk)ϵk[3M22p2E2p+1pEplnEp+pEp−p +1ϵ2k(2+M2E2p)(−1+2E2p−M22pEplnEp+pEp−p)] + D2Ep(τ)n\tiny\rm{F{}}(Ek)Ek[3M22p2E2p+10p4+16p2M2+3M42(p2−k2)p2E2p+M2pkE2pln∣∣∣p+kp−k∣∣∣ −2E4p+2E2kE2p−M42pk(Ep−Ek)E3p(ln∣∣∣p+kp−k∣∣∣+EkEp+Ekln∣∣∣M2+EpEk+pkM2+EpEk−pk∣∣∣)]},

where denotes a principal value, and etc are defined in eq. (A.2). The constant contribution reads

 G\tiny\rm{NLO}\tiny\rm{V}∣∣∣\tiny\rm{const.}4g2CACF = (4.5) + n\tiny\rm{F{}}(Ek)Ek[1p2−3E2p−M2p2k2−M2E2pE2k+M2(4E2k−M2)2pkE2pE2kln∣∣∣p+kp−k∣∣∣]}.

The spatial correlator subsequently follows from eq. (2.2), with the susceptibility inserted from eq. (4.1).

Although the expression in eq. (4.4) is finite, its numerical evaluation is non-trivial. There are at least three separate challenges: at small various parts of the expression are divergent, and care must be taken in order to avoid significance loss in their cancellation; at there is a pole which is defined in the sense of a principal value; and at large there is a vacuum part which decreases only slowly (although it is integrable). For the reader’s benefit, let us briefly specify how we have dealt with these challenges.

• The small- divergence originates from the terms integrated over in eq. (4.4) and from terms where the integral had already been carried out. For the latter type the divergent part reads

 D2Ep(τ)ϵ3k(2+M2E2p)(−1+2E2p−M22pEplnEp+pEp−p)(12+n\tiny\rm{B{}}(ϵk)), (4.6)

containing both vacuum and Bose-enhanced structures. The integral needs to be carried out precisely enough such that the cancellation takes duly place.

• The principal value integration can be handled for instance by reflecting the range into the range :

 (4.7)

Here has to be evaluated precisely enough for cancellations at to take place.

• A possible way to accelerate the convergence at large is with the help of the function in eq. (4.4). (Note that a power tail only appears in the vacuum part.) The simplest subtraction removes just the leading asymptotic behaviour , e.g.

 θ(k)≡3E2pΘ(k−k\tiny\rm{min})k2+λ2,∫∞0dkθ(k)k=3E2p2λ2ln(1+λ2k2\tiny\rm{% min}), (4.8)

but of course more refined choices can be envisaged. (By replacing by a “gluon mass” it would be possible to take and still carry out exactly, cf. eqs. (B.22), (Appendix B), but the price to pay is that then the small- range has more structure than before, with a would-be divergence only cut off at .)

For a transparent representation of the numerical results, we consider two different normalizations. A simple and theoretically clean possibility is to normalize the results to the free massless expression from eq. (2.4). Another reference point is to make use of the full NLO spectral function in vacuum [45][48]:

 ρ\scriptsize vac\tiny\rm{V}(ω) = −θ(ω−2M)CA(ω2−4M2)12(ω2+2M2)4πω+θ(ω−2M)8g2CACF(4π)3ω2{ (4M4−ω4)L2(ω−√ω2−4M2ω+√ω2−4M2)+(7M4+2M2ω2−3ω4)acosh(ω2M) +ω(ω2−4M2)12[(ω2+2M2)lnω(ω2−4M2)M3−38(ω2+6M2)]}+O(g4),

where the function is defined as

 L2(x)≡4Li2(x)+2Li2(−x)+[2ln(1−x)+ln(1+x)]lnx. (4.10)

There is no transport peak in the vacuum expression, and recalling eq. (2.2), the corresponding spatial correlator can be obtained through

 G\scriptsize recii(τ)≡∫∞0dωπ(−ρ\scriptsize vacV)(ω)cosh(β2−τ)ωsinhβω2. (4.11)

A normalization with respect to a similar “reconstructed” correlator has been used in ref. [21] (see also ref. [49]), and may be useful for phenomenological purposes, although from the theoretical perspective it induces new systematic uncertainties.

In fig. 2 we show our results in both normalizations, compared with lattice data from ref. [21]. (The gauge coupling has been fixed as explained in connection with fig. 1; at very small a different running would appear reasonable but in the absence of an NNLO computation and continuum-extrapolated lattice data, we stick to the simplest choice in the following.) Like in fig. 1, an excellent agreement is found at large time separations, if a value is assumed. (The discrepancy in fig. 2(left) at small is probably due to the missing continuum extrapolation.)

## 5 Modification of the transport peak

As mentioned in sec. 3, the correlator has a constant (-independent) part at LO and at NLO, but within the full dynamics this is expected to turn into a slowly evolving function. The purpose of this section is to estimate how precisely Euclidean data should be measured in order to resolve the slow time dependence.

In order to reach this goal, we model the transport peak through a Lorentzian shape,

 ρ\tiny\rm{(L)}ii(ω)≡3Dχωη2ω2+η21cosh(ω2πT). (5.1)

Here corresponds to the heavy flavour diffusion coefficient. The Lorentzian shape can be correct only at small frequencies, , cf. e.g. ref. [50]; we have chosen to cut it off at large frequencies through the same recipe that has been used in the massless case [27]. The susceptibility is fixed according to ref. [21], . We then vary , in each case tuning the other parameter so as to keep the area under the transport peak fixed at the value predicted by the NLO expression, eqs. (3.5), (4.1), (4.5), (2.2):

 1β∫β0dτG\tiny\rm{(L)}ii(τ)=∫∞0dωπ2ρ\tiny% \rm{(L)}ii(ω)βω≡[G\tiny\rm{LO% }ii+G\tiny\rm{NLO}ii]\tiny\rm{const.}. (5.2)

Subsequently the “correct” is obtained by replacing the part through a -dependent function, , determined by via eq. (4.11). As a guideline, we recall that may be expected according to refs. [14, 15]; that ref. [21] found even ; and that in the massless case values down to appear possible [27].

Results are shown in fig. 3. We observe that if , then there is hope of resolving it with high enough precision of the lattice data (it appears that statistical errors should probably be reduced to 20% of the current ones to be sensitive to the features of the transport peak; obviously improvements in statistical accuracy need to be accompanied by a corresponding decrease in systematic uncertainties). In the case a reliable determination from the massive vector current correlator appears challenging.

## 6 Modification of the threshold region

The second qualitative structure affecting the massive vector current correlator, the quarkonium threshold region inducing a “fast” time dependence from the energy scale , is also expected to undergo drastic changes in the interacting theory. At low temperatures, the spectral function is characterized by quarkonium resonances; at high temperatures, these are expected to move, broaden, and eventually dissolve into a mere threshold enhancement.

Based on our earlier investigations [51, 52], we expect that in the temperature range of interest there is at most one resonance peak in the vector channel spectral function, placed slightly to the left from the free quark-antiquark threshold. Akin to eq. (5.1), we model this by a skewed Breit-Wigner shape,222The approach in this section is schematic; for possible other model shapes see, e.g., ref. [53]. (A general discussion of the sum rule approach can be found in ref. [54].)

 ρ\tiny\rm{(BW)}ii(ω)≡Aω2γ2(ω′)2+γ21cosh(ω′2M),ω′≡ω−2M+ΔM, (6.1)

constructed so as not to contribute to the transport coefficient. To reduce the number of free parameters to two, we set in the following. The contribution from such a peak, determined through eq. (4.11), is added to the thermal NLO result as such, and to the vacuum result of eq. (LABEL:full_vac) with , , keeping the area under the peak roughly invariant. Obviously these choices are arbitrary, but they should nevertheless convey a qualitative impression on the importance of resonance contributions. Based on ref. [31], in which a resonance peak around a threshold was matched to the thermal NLO spectral function above the threshold, we expect that and .

Results are shown in fig. 4. Comparing with fig. 3, it can be seen that a change of the spectral shape around the threshold region affects the Euclidean correlator at smaller than a change of the transport peak. In fact, if the resolution, after continuum extrapolation, were high enough that the lattice and perturbative curves could be subtracted from each other (rather than normalizing them to a function which diverges at short distances), then the two features could be disentangled by inspecting intermediate distances, , in which the threshold region contributes much more prominently than the transport peak. In contrast the two features are difficult to tell apart if reliable data is only available for . (That said, we stress again that the present section is meant as indicative only.)

## 7 Conclusions

The main purpose of this paper has been to compute the massive quark-number susceptibility and vector current correlator at next-to-leading order (NLO) in thermal QCD. The NLO results are shown in eqs. (4.1), (4.4), (4.5), and illustrated numerically in figs. 1, 2.

Our semi-analytic results can be directly compared with numerical lattice Monte Carlo simulations. Although no continuum limit has been reached for the massive vector current correlator [21], the agreement as seen in figs. 1, 2 is quite remarkable. (We have compared with quenched data, because these are “closest” to the continuum limit, but our analytic results are also valid for the unquenched case.)

The good agreement suggests that resummed perturbative NLO computations, through which the Minkowskian spectral function has been determined around the threshold region [51, 52], might be more accurate than sometimes assumed. A similar observation has also been made by comparing NLO computations to lattice data at spatial separations relevant for quarkonium physics [55].

The back side of a remarkable agreement is that the strict NLO result of the present paper reflects the physics of an infinitely narrow transport peak, rather than genuine heavy quark diffusion, and of a threshold singularity, rather than genuine quarkonium resonances — and yet it works well. This implies that even though impressive efforts to extract properties of the transport peak and quarkonium resonances from lattice data are being made, a very fine resolution is needed for getting systematic errors fully under control.

In order to quantify the resolution that lattice simulations should reach for observing substantial deviations from the NLO results, we have finally probed the significance of possible non-perturbative effects. The impact of a smeared transport peak is illustrated in fig. 3; the impact of resonance-like spectral weight below the threshold in fig. 4. Interpreting lattice results such as those in ref. [21] as a deviation from NLO, it appears that a diffusion coefficient could be observed in principle, but that larger values are hard to disentangle. Perhaps, a fruitful approach is to determine the diffusion coefficient by other means [14, 15], subtract the corresponding contribution from lattice data, and hope that any remaining deviations reflect the physics of the threshold region.

Ultimately, going beyond modelling, we suspect that systematic studies are only possible once a continuum extrapolation is available in a broad -range. At this stage vacuum physics can be subtracted, as outlined in ref. [27], and the non-divergent remainder subjected (at least in principle) to a model-independent analytic continuation [30].

## Acknowledgements

We thank H.-T. Ding for helpful discussions and providing us with lattice data from ref. [21], and M. Vepsäläinen for collaboration at initial stages of this work. This work was partly supported by the Swiss National Science Foundation (SNF) under grant 200021-140234.

## Appendix Appendix A Results for individual master sum-integrals

In this appendix we list results for the 2-loop sum-integrals in eqs. (2.13), (2.17), after carrying out the Matsubara sums as well as the Fourier transformation in eq. (2.18). For brevity we refer to the sum-integrals with the notation

 Im1m2m3n1n2n3n4n5(τ)≡T∑ωne−iωnτ∑∫K{P}(M2)m1(Q2)m2(2K⋅Q)m3