Thermal quarkonium physics in the pseudoscalar channel

The pseudoscalar correlator is an ideal lattice probe for thermal modifications to quarkonium spectra, given that it is not compromised by a contribution from a large transport peak. We construct a perturbative spectral function incorporating resummed thermal effects around the threshold and vacuum asymptotics above the threshold, and compare the corresponding imaginary-time correlators with continuum-extrapolated lattice data for quenched SU(3) at several temperatures. Modest differences are observed, which may originate from non-perturbative mass shifts or renormalization factors, however no resonance peaks are needed for describing the quenched lattice data for charmonium at and above T ~ 1.1Tc ~ 350 MeV. For comparison, in the bottomonium case a good description of the lattice data is obtained with a spectral function containing a single thermally broadened resonance peak.


Introduction
Even though potential models and sum rules have suggested that charmonia resonances (J/ψ, η c ) have dissolved by the time that the temperature of a QCD plasma reaches 300 MeV or so (cf. e.g. refs. [1][2][3][4] for overviews), it has been difficult to consolidate this picture with direct lattice measurements of thermal 2-point correlators corresponding to vector and pseudoscalar operators. A part of the problem is that spectral information is well hidden in an imaginary-time measurement. For instance, in the vector channel (J/ψ), the imaginary-time correlator gets a substantial contribution from a "transport peak" located at small frequencies |ω| ≪ T 2 /M [5,6], where T is the temperature and M is a heavy-quark mass. This reflects interesting physics of heavy quark diffusion (cf. ref. [7,8] and references therein), but has little to do with the fate of quarkonium resonances at ω ∼ 2M , even though it severely hampers the corresponding imaginary-time measurement [6]. 1 In contrast to the vector channel, no transport peak is expected to be present in the pseudoscalar channel (η c ) [14,15]. The inertness of the pseudoscalar correlator then led ref. [16] to conclude that there is no sign of melting of the corresponding charmonium resonances (another recent analysis can be found in ref. [17]). However, another lattice study [18] did find small changes in the pseudoscalar case (in contrast to ref. [16], ref. [18] made use of the quenched approximation; we return to this in sec. 7). Moreover, for the vector channel, it was subsequently shown [19] that within statistical and systematic uncertainties the results of ref. [18] are compatible with perturbative predictions involving no charmonium resonances (only a modest threshold enhancement). This demonstrates that J/ψ and η c dissociation is a possibility and that a very fine resolution is necessary for constraining this physics.
The purpose of the current study is to carry out considerations such as those in refs. [18,19] but for the pseudoscalar correlator and with improved resolution. 2 The pseudoscalar correlator reflects the physics of the η c but couples to the same non-relativistic operators as the vector channel as far as spin-independent terms are concerned.
Our presentation is organized as follows. After defining the basic setup in sec. 2, the vacuum asymptotics of the pseudoscalar spectral function is reviewed in sec. 3, whereas sec. 4 contains a discussion of thermal effects around the quark-antiquark threshold. The lattice setup, quark-mass interpolation, renormalization, and continuum extrapolation are summarized in sec. 5. The perturbative and lattice results are compared in sec. 6, and conclusions and an outlook are offered in sec. 7. Two appendices detail the thermal correlator at strict next-toleading order (NLO), as a function of imaginary time and real frequency, respectively.

Basic setup
We consider QCD with one heavy quark, of bare mass M B and MS mass m(μ), and N f essentially massless dynamical quarks (later on we set N f = 0). The pseudoscalar imaginarytime correlator is defined as where ψ denotes the heavy quark Dirac spinor, and ... c indicates that only the connected (i.e. flavour non-singlet) contraction is included in the thermal average. Defined this way, the correlator is believed to be finite after mass and gauge coupling renormalization, even though care needs to be taken with the regularization of γ 5 (see below).
In the presence of a spatial ultraviolet regulator, the correlator can be Fourier-transformed. Denoting by ω n a bosonic Matsubara frequency, the Fourier transform is given byG P (ω n ) ≡ 1/T 0 dτ e iω n τ G P (τ ). The corresponding spectral function is formally given by the cut ρ P (ω) = ImG P (−i[ω + i0 + ]) , (2.2) and the correlator of eq. (2.1) is reproduced through (2.3) If present, a narrow "transport peak", ρ P /ω ∝ δ(ω), would yield a constant (τ -independent) contribution to G P (τ ). However, we have checked by an explicit computation (cf. appendix A) that G P contains no τ -independent contribution at NLO. This can be partly understood from the partially conserved axial current (PCAC) relation where T a is a traceless matrix in flavour space. According to eq. (2.4) the pseudoscalar correlator can be obtained from a corresponding axial charge correlator through derivatives, ∂ 2 τ G A 0 = 4G P , so that constant parts appearing in G A 0 get deleted. However, at finite temperature the argument is not rigorous, since in principle G A 0 could contain ∼ τ ( 1 T − τ ) which would yield a constant G P after the second derivative. As mentioned, up to NLO no such term is found. A physical reason for the absence of a constant part is that they only appear in cases in which the operators are related to a conserved current in some limit (for example, the scalar density equals the fermionic part of the "trace anomaly" T µ µ ). In order to estimate ρ P , it is useful to have a physical picture in mind. As a cut, the spectral function describes a hadronic decay width of a pseudoscalar meson. The spatial average in eq. (2.1) implies that the meson is at rest with respect to the heat bath. The energies of the decay products are of the order ω/2. For ω ≫ 2M ≫ πT , where M denotes a pole mass, all thermal effects are small (exponentially suppressed at leading order, powersuppressed in general [22]), because thermal motion represents a minor correction to the kinematics of the decay products. In this regime, the spectral function can be extracted from vacuum computations. Decreasing the energy to ω ∼ 2M , we are entering the threshold region. In this situation even the vacuum computations become less precise, because the decay products move slowly and have much time for final-state interactions, leading to the so-called Sommerfeld effect (cf. ref. [23] and references therein). Once the momenta of the decay products are of the same order as the Debye mass, M v ∼ m D ∼ gT where g ≡ √ 4πα s , thermal corrections are of order unity [24]. Below the threshold, thermal corrections represent the dominant physics: the real-time static potential develops an imaginary part [25][26][27] which leads to the dissociation of quarkonium resonances if T ≫ α s M [28][29][30].
This rich physics implies that a number of different techniques are needed for a reasonably precise estimate of the pseudoscalar spectral function. The vacuum computations are discussed in sec. 3 and the thermal ones in sec. 4.

Vacuum contribution above the threshold
We start with a brief non-expert review of vacuum results for the pseudoscalar spectral function, relevant for M v ≫ m D . In this regime thermal effects are power-suppressed [22]. Many computations have been carried out in the "on-shell" (i.e. "pole-mass") scheme, keeping the gauge coupling in the MS scheme. Like above, we may then define the pseudoscalar operator as P = M Bψ iγ 5 ψ, 3 and re-expand then M B in terms of the pole mass M [33][34][35].
The vacuum ρ P has been estimated for a general ω/M up to O(α 2 s ) [36,37], and for the asymptotics at ω ≫ M and for certain moments of the spectral function up to O(α 3 s ) [38][39][40]. The relation of the QCD spectral function to the corresponding NRQCD one, needed for near-threshold resummation, is known up to O(α 2 s ) [41]. Unfortunately these results show poor convergence, even for the bottom quark case [42].
To illustrate the problem, we write the spectral function in the form The subscript c refers to a connected quark contraction. The R-function is expanded as where C F ≡ (N 2 c − 1)/(2N c ) and T F ≡ 1/2. The asymptotics of the different contributions read (cf. e.g. refs. [43,44]; ζ n ≡ ζ(n)) Order by order these describe the full ω-dependence [36,37] well already at modest ω > ∼ 4M , however as a whole the expression does not converge. Indeed, even if we choosē for evaluating α s (μ), so that it is small at ω ≫ M , its decrease ∼ 1/ ln(ω 2 /Λ 2 MS ) is not fast enough to kill the growing ln(ω 2 /M 2 ) appearing in the coefficients in eqs. (3.4)-(3.7).
However, it is not necessary to stick to the on-shell scheme. The reason for using a pole mass is that it is formally necessary for defining a perturbative series around the threshold (cf. appendix B). However, the pole mass may be re-expanded as an MS mass [33][34][35]:   Changing the normalization from eq. (3.1) to and writingR p c as in eq. (3.2), mass dependence drops out from the asymptotics: We now see that the choice in eq. (3.8) removes large logarithms and leads to increasingly small corrections as ω grows. Therefore eq. (3.10) yields a reliable prediction for large ω.
We note in passing that the asymptotics in eqs. (3.11)-(3.15) is identical to what is obtained for the connected part of the scalar density correlator [43,44]. This demonstrates that the scalar and pseudoscalar densities have the same anomalous dimension, which in turn is related to the quark mass anomalous dimension, so that m 2 (μ)R p c (ω,μ) is independent ofμ. Given the poor convergence of eq. (3.9), the question arises how we should relate M and m(μ) to each other. Following old tradition, we parametrize quark masses through their MS values at the scaleμ ref ≡ 2 GeV. The corresponding running mass is obtained from 5-loop running [45,46]. One way to fix the value of M is to require that eqs. (3.1) and (3.10) agree at a value ω = xM where we assume both representations to be reliable, say x = 4...8: The results obtained at O(α 2 s ) are shown in table 1. (There are many other pole mass definitions, however none of them are quite satisfactory for our purposes.)

Thermal contributions around the threshold
Just above the threshold (i.e. for ω−2M ≪ M ), the frequency ω can be parametrized through a relative velocity as M v 2 ≡ ω − 2M . If M v < ∼ gT , thermal effects are of order unity [24]. Because the heat bath breaks Lorentz invariance, their technical treatment is much harder than that of vacuum contributions: within strict perturbation theory, only the NLO level has been reached (cf. appendix B). Unfortunately the NLO expression turns out to be rather useless in practice; it indicates that thermal corrections are power-suppressed as ∼ α s T 2 /M 2 , whereas in reality there are also thermal corrections only suppressed by (non-integer) powers of α s . These originate from soft-gluon mediated effects (thermal heavy quark mass correction, Debye screening, real 2 → 2 scatterings of heavy quarks off plasma particles), and require a resummed treatment. We follow here the implementation of ref. [28].
Before proceeding let us stress that, as is characteristic of resummed treatments, the approach can be justified theoretically only in a certain parametric range (see below). In practice it is also used somewhat outside of this range, and for this purpose it is "enhanced" by a number of phenomenological ingredients (see below), in order to avoid a drastic breakdown in the latter situations.
Within the non-relativistic regime v ≪ 1, the pseudoscalar spectral function is related to the vector channel one: For the vector channel, the spectral function is obtained from a Wightman function C > as where C > is solved from For t > 0 the potential is of the form [25][26][27] where the function represents the effects of real 2 → 2 scatterings of the quark and antiquark off medium particles. For t < 0 the sign of Im V T is reversed. For numerical estimates we insert 2-loop values of m D and a thermal α s from ref. [52] (for m D the 3-loop level has also been reached [53]). At short separations, r ≪ 1/m D , we replace the thermal potential by a vacuum one. At 2-loop level, where γ is the Euler constant and, for N f = 0 [54], The 3-loop potential is also known in analytic form by now [55]. At rm D ≪ 1 the potential of eq. (4.8) is less binding than that of eq. (4.5) because the coupling runs towards zero; at rm D ≫ 1 eq. (4.5) is less binding because of Debye screening. For our numerics, we employ max{V 0 , Re[V T (r)−V T (∞)]} as the r-dependent real part of the potential, whereby V 0 applies at short separations and Re V T at large ones. This regulates the infrared sensitivity of the potential that has been widely discussed in the context of heavy-quark mass definitions. In order to combine the threshold behaviour from eqs. (4.1)-(4.5) with the asymptotics from eq. (3.10), we normalize ρ NRQCD P by a multiplicative factor. Given that the perturbative value of this factor is poorly determined, we rather choose a free coefficient, denoted by A: (4.10) The value of A is chosen so that ρ QCD P attaches to the MS asymptotics from eq. (3.10) continuously and with a continuous first derivative at some 2M < ω < 3M , cf. fig. 1.
Once we are sufficiently below the threshold, 2M − ω ≫ α 2 s M , the description of the spectral function through eqs. (4.1)-(4.5) is no longer applicable. In fact the Schrödinger description overestimates the spectral function in this regime: the correct result is non-zero  at finite temperature, but exponentially suppressed, as discussed below eq. (B.6). In order to model this suppression, we multiply φ of eq. (4.5) by θ(2M −ω)e −|ω−2M |/T . In the parametric range |ω − 2M | ∼ α 2 s M ≪ gT, α s M, πT ≪ M for which an approach based on eqs. (4.1)-(4.5) can be justified [27][28][29], this amounts to a higher-order effect. The range ω ≤ 1.0M , which has very little effect after adopting this recipe, is cut off from the numerical evaluation of eq. (2.3). These choices are purely phenomenological; formally they amount to effects which are of higher order in α s than our computation.
For small frequencies, the prefactor in eq. (4.2) results in an additional suppression of the spectral function. However other exponentially small (when ω ∼ 2M ≫ πT ) effects have been omitted in the non-relativistic solution of C > . We "reconstruct" the exponential factor by taking it over from the tree-level computation in which, as indicated by eq. (5.2), it is given by tanh GeV, when the non-relativistic expansion is at most marginally viable, this leads to a ∼ 10% reduction of the threshold enhancement obtained from the Schrödinger equation. Of course, this reduction is partly compensated for by the matching in eq. (4.10). In any case, our "best estimates" for the thermal spectral function at different m(μ ref ), obtained as outlined above, are shown in fig. 1.
We note in passing that, in accordance with the original findings of ref. [28], a small resonance peak can be observed at the largest quark masses at this temperature. We return to a discussion of this point in sec. 6.

Ensemble and statistics
We have carried out lattice simulations in quenched SU(3) gauge theory with O(a) improved Wilson quarks as valence quarks, in order to measure the connected pseudoscalar correlator and take the continuum limit. At each lattice spacing, this involves an interpolation to a bare quark mass that reproduces the physical J/ψ or Υ mass at a low temperature. The basic techniques (action, simulation algorithm), as well as tests suggesting that finite-volume effects are reasonably small, are described in ref. [18]. The progress in the past 5 years amounts to simulating several fine lattice spacings, and on each of them several values of the bare quark mass. Details concerning the ensemble are collected in table 2.
There is a well-known problem with simulations close to the continuum limit with periodic boundary conditions, namely the "freezing" of topological degrees of freedom [56]. At low temperatures, this effect is unphysical. Our measurements were separated by 500 sweeps, each consisting of 1 heatbath and 4 overrelaxation updates. The initial thermalization consisted of 2000 sweeps at β = 7.192 and of 5000 sweeps at β = 7.793. With this setup we do observe some evolution between the topological sectors up to β ≃ 6.8, however at the values β > ∼ 7.2 that play a role in our analysis, no evolution takes place. Even though the resulting uncertainties should be minor at high temperatures, where the physical topological susceptibility is small, the freezing does affect our low-temperature runs as well, particularly our scale setting [57] 4 though the parameter r 0 [62]. In some sense, the values of r 0 /a in our study should be interpreted as "r 0 /a at fixed trivial topology". It is not known how much these differ from the corresponding r 0 /a in the physical θ-vacuum, however in principle the differences are suppressed by the inverse volume and therefore perhaps moderate.

Tuning of quark mass
For any given lattice spacing, the first step is to interpolate results to the physical bare quark mass. We have determined quark propagators at six different bare quark masses (or hopping parameters; cf. table 2), roughly speaking between the charm and bottom masses, with a denser grid around the former. As the corresponding observable we measure the screening mass corresponding to the spatial components of the vector current (including only transverse components with respect to the measurement direction). The measurement is carried out deep in the confined phase (T ∼ 0.75T c in table 2), and is interpreted as being a good approximation for a measurement at zero temperature. Thus, the screening mass is identified as the J/ψ or Υ mass (we use J/ψ rather than η c for scale setting because it is a 4 We take the opportunity to note that the values of (t 0 /a 2 ) Clover cited in ref. [57] were marred by a bug, and are too large by a few %. This does not affect any of the conclusions of ref. [57], because the clover results were excluded due to their peculiar volume dependence. We thank Lukas Mazur for locating the bug.  s × N τ ) included in the current analysis. The measurements of r 0 /a at low temperatures are based on our own measurements of large Wilson loop expectation values, and correspond to a sector of fixed trivial topology, due to the freezing of topological degrees of freedom [56]; we refrain from estimating systematic uncertainties. When citing physical units, we make use of r 0 = 0.47(1) fm from ref. [50]. Conversions to units of T c are based on r 0 T c = 0.7457(45) from ref. [57]. Our new data for r 0 /a, together with the older one summarized in table 2 of ref. [57], which include classic results from refs. [58,59] narrower resonance, however in practice this difference is inessential on our resolution). The screening mass is extracted from a two-exponential fit, with an ensemble of fit ranges chosen by hand close to where a minimal χ 2 /d.o.f. ∼ 1 can be found. The result is denoted by m V ; the corresponding values are plotted in fig. 2(left) in lattice units.
Proceeding now to the thermal runs in the temporal direction (τ ), the correlators were measured at the same values of the bare quark masses. Therefore each measurement can be assigned to a specific value of the zero temperature (screening) mass m V . The data is subsequently interpolated quadratically in m V , where G free P is defined in eq. (5.2). The four κ-values closest to the desired point are used for this interpolation (or, in some cases, mild extrapolation). Having determined the coefficients α i (τ T ), we finally set m V → m J/ψ or m V → m Υ , in order to get the physical result at the given τ and T . The procedure is illustrated in fig. 2(right) for four values of τ T . The interpolated correlator for the physical mass can be read off from the dotted line.

Normalization of imaginary-time correlators
With the data interpolated to a physical quark mass, the next task is to extrapolate to the continuum limit. This is facilitated by normalizing the lattice data to a function which captures most of the rapid τ -dependence at small τ ≪ 1/T .
Because of a non-zero anomalous dimension, there is no exact "free result" for the pseudoscalar correlator like, say, for the vector channel correlator in the chiral limit. However, we can define where the expression in curly brackets is the tree-level spectral function, and the remaining factors amount to those in eq. (2.3). Because of the anomalous dimension, a result normalized through eq. (5.2) does not go to a constant value at τ ≪ 1/T , though we expect it to be slowly varying. By definition we set M 1S ≡ 1.5 GeV in eq. (5.2) for the charmonium case, and M 1S ≡ 5.0 GeV for the bottomonium case (the physical value being M 1S ≈ 4.7 GeV).

Renormalization and continuum extrapolation
In the absence of non-perturbative renormalization factors for quenched massive pseudoscalar densities, we have made use of massless perturbative results [63,64], supplemented by a 1loop correction for the mass dependence [65,66]. A tadpole-improved coupling was inserted in these formulae, along the lines discussed in ref. [67]. The difference of the 1-loop and 2-loop expressions from ref. [64] was employed for getting a feeling about the systematic uncertainties of perturbative renormalization. There is an important subtlety related to the pseudoscalar renormalization factor, denoted by Z P . Since Z P is supposed to bring us from the lattice to the MS scheme, it depends on how γ 5 is defined on the latter side. With the 't Hooft-Veltman choice [31], the strict MS pseudoscalar density needs to be multiplied with an additional finite renormalization factor in order for it to have the same anomalous dimension as the scalar density [32], In so-called naive dimensional regularization, no additional factor is needed at low perturbative orders. The 1-loop pseudoscalar renormalization factor computed in ref. [63] is smaller than that in ref. [64], the reason being that it is multiplied by Z 5 in comparison with ref. [64].
We are interested in a pseudoscalar density which corresponds to the perturbative result presented in sec. 3 and has the same anomalous dimension as the scalar density (cf. comments below eq. (3.15)). Then the normalization of ref. [63] is appropriate, i.e. we need Z 5 Z L,MS We have also carried out fits quadratic in 1/N 2 τ , or linear in 1/N 2 τ but restricted to the three finest lattices; the variations are of the same order as the differences from using 1-loop and 2-loop renormalization factors, which we then adopt as our estimate of systematic uncertainties. Right: Continuum-extrapolated correlators. The uncertainties have been obtained by adding together, in quadrature, statistical errors and systematic uncertainties from the difference of using 1-loop and 2-loop renormalization factors. the notation of ref. [64]. For reference we note that numerically Z 5 Z L,MS P ∼ 0.8 in the range of our β's and masses.
The perturbative renormalization factors bring us to the MS scheme at the scaleμ = 1/a. For scale-independent results the pseudoscalar correlator should then be multiplied by a quark mass in the same scheme and at the same renormalization scale, i.e. at which only three lattice spacings are available. Other fit forms (quadratic in 1/N 2 τ , or linear in 1/N 2 τ but restricted to the three finest lattices) have also been successfully tested, cf. caption of fig. 3.
The whole procedure is implemented in the form of a bootstrap analysis. The final continuum-extrapolated results are shown in the right panels of figs. 3 and 4 for charmonium, and in the right panel of fig. 5 for bottomonium.

Modelling and comparison
The purpose of this section is to compare the predictions originating from the spectral functions in fig. 1, shown as the corresponding imaginary-time correlators for the charmonium case in fig. 4(left) and for the bottomonium case in fig. 5(left), with the lattice data shown in figs. 4(right) and 5(right). We restrict ourselves to a rather qualitative discussion here, identifying two main issues which help to explain the differences observed.
Let us start by noting that one striking feature of the lattice data is its apparent inertness: if plotted in units of τ T c , all curves more or less fall on top of each other, cf. fig. 4(right). Such an inertness, however, does not indicate the persistence of bound states; indeed a similar inertness is seen on the perturbative side where no bound states are present, cf. fig. 4(left).
There is one clear difference between the lattice and perturbative sets in fig. 4: the perturbative data display a larger positive "curvature" as a function of τ T c , overshooting the lattice data at τ T c > ∼ 0.2.
This difference may have a relatively simple explanation. In order to demonstrate this, we introduce a "model" spectral function, which is obtained from the perturbative one by adjusting its overall normalization as well as the threshold location: 5 The overall magnitude is adjusted because of uncertainties related to the perturbative renormalization factors on which we rely on the lattice side (cf. sec. 5.4), and the threshold location because pole-type masses are poorly estimated in perturbation theory (cf. sec. 3). Results from simple χ 2 fits to the lattice data, with A and B treated as fit parameters, are shown in fig. 6. An excellent representation of the lattice data can be found in all cases. The same is true, by and large, for the bottomonium case. The fit parameters are collected in table 3. The model spectral functions corresponding to fig. 6 are shown in fig. 7(left), where they are compared with the unmodified perturbative ones. The first essential feature is that some additional spectral weight is needed at very large ω. This accounts for the small difference observed at τ T < ∼ 0.15 between the lattice and perturbative results. This discrepancy may 5 Functional forms which induce a shift only at moderate ω lead to similar results, e.g. A ρ   originate from the perturbative renormalization factors that we have used (cf. sec. 5.4). Nonperturbative renormalization would help to clarify this issue. The second feature observed in fig. 7(left) is that the thresholds shift to larger masses. This is perfectly admissible, given that the procedure we adopted for estimating a vacuum "pole mass" in sec. 3 is subject to large uncertainties, and that thermal mass corrections beyond the perturbative ones included through eq. (4.5) could be substantial. In fact, it is known from lattice studies of renormalized Polyakov loop expectation values that close to T c thermal mass corrections are positive, whereas eq. (4.5) predicts a negative thermal mass correction, which in Polyakov loop measurements is observed only at T > ∼ 3T c [68].
For the bottomonium case, the model spectral functions are shown in 7(right). The perturbative input originates from fig. 1 and contains a resonance peak at T < ∼ 1.5T c . At these larger frequencies non-perturbative mass shifts are barely visible. In contrast the overall normalization is corrected by a larger amount than for the charmonium case, and downwards. This difference is not surprising, given that because of the larger mass, discretization effects, including corrections of O(α 2 s am L ), imply that our renormalization factors and continuum extrapolation are likely to contain larger systematic uncertainties.
We stress that, as fig. 7(left) shows, no resonance peaks are needed for representing the lattice data for the charmonium correlator even at the lowest temperatures in the deconfined phase; a modest threshold enhancement is perfectly sufficient. In contrast a thermally broad-  Table 3: Best fit parameters according to eq. (6.1). The left set corresponds to charmonium, the right to bottomonium. In these fits the errors of the lattice results at different values of τ T , which are dominated by systematic uncertainties, have been treated as independent of each other. Therefore the results are somewhat qualitative in nature, and we refrain from citing errors.
ened resonance peak, as predicted by resummed perturbation theory, may well be present in the bottomonium case at T < ∼ 1.5T c , cf. fig. 7(right). Even though perturbation theory contains inherent uncertainties, such features have been consistently observed in previous computations, however the quantitative properties of the peaks depend on the precise approximation under which the thermal potential and the gauge coupling have been estimated.

Conclusions and outlook
We have presented a resummed perturbative estimate of the thermal quarkonium pseudoscalar spectral function on one hand (secs. 3 and 4), and a continuum-extrapolated quenched lattice measurement of the corresponding imaginary-time correlator on the other (sec. 5). Our main conclusions originate from a comparison of these two computations ( fig. 6). Unambiguous but not overwhelming differences are observed, which call for an interpretation. Two possible culprits have been put forward in sec. 6. First, it is plausible that the perturbative renormalization that we have used for the lattice correlators, and other uncertainties related to the continuum extrapolation at the very large β-values that we have used, could result in inaccuracies on a ∼ 5 − 10% level for the overall normalization of the imaginarytime correlators. This problem is expected to be more severe in the bottomonium case, given that O(α 2 s am L ) corrections could affect the continuum extrapolation. Second, the thermal quark-antiquark threshold location is not accessible to perturbation theory on a quantitative level; treating it rather as a fit parameter significantly improves upon the agreement. The latter problem is, relatively speaking, more severe in the charmonium case, given that the threshold is located at a smaller energy.
In contrast, there is no need to modify the perturbative charmonium spectral function through resonance peaks at any of the deconfined temperatures that we have considered, which is consistent with a rapid dissociation of the η c meson in the deconfined phase of quenched QCD. In the bottomonium case, a thermally broadened η b peak can persist up to ∼ 1.5T c . To rephrase these observations more strictly, our statement is that perturbative spectral functions, which display these features, yield imaginary-time correlators perfectly compatible with lattice data, apart from a modest shift of the threshold location. At the same time the existence of some features beyond the perturbative ones cannot be excluded.
In order to consolidate this tentative picture, several future steps are needed. To remove the uncertainty concerning overall normalization, non-perturbative renormalization is desirable, including the determination of quark mass effects of O(am L ) so that O(a) improvement is complete [69]. Ambiguities originating from the fact that our low-temperature runs, needed for scale setting, are frozen to the trivial topological sector, should be addressed, for instance by taking volume averages in very large volumes [70]. For insight on dissociation patterns, spectral reconstruction techniques could be applied to the continuum-extrapolated data.
We end with an important comment concerning the quenched approximation on which this study was based. Generally, for a comparative value of T /T c , we would expect the quenched deconfined phase to be better described by perturbation theory (i.e. more weakly coupled) than the unquenched one. A physics argument is that in the quenched theory hadrons (glueballs) are heavy, with m 0 ++ ≫ 1 GeV, so that the system needs to be heated up to a high temperature to "dissolve" them. Concretely, T c ≃ 1.24Λ MS > Λ MS , so that α (EQCD) s | T ≃Tc ≃ 0.2 is reasonably "small" [52]. In contrast, the lightest excitations of the unquenched theory are pions, with m π ≪ 1 GeV, and only a modest heating is needed to reach the transition temperature (or crossover). Concretely, T c ≃ 0.45Λ MS < Λ MS , so that α (EQCD) s | T ≃Tc > 0.3 is undoubtedly too large for perturbation theory to apply. To summarize, the "effective" strong coupling, describing soft thermal physics near threshold, is likely to be larger in the unquenched theory, and therefore the physics of η c could differ from that in the quenched case. Indeed, according to an up-to-date potential model study [11], the η c is expected to display a resonance peak at T > T c in the unquenched theory. Therefore our analysis should be repeated for that situation (which has previously been addressed with relativistic lattice techniques in ref. [16], however without a continuum extrapolation or a comparison with resummed perturbation theory). We hope that the present paper helps to trace out a path for tackling this ambitious challenge. In this appendix we compute the imaginary-time correlator G P (τ ) up to NLO in strict perturbation theory. The purpose is to demonstrate that it contains no constant contribution, unlike the corresponding correlators in the vector [19] and scalar [20] channels. Making use of the notation  Here To display the subsequent results, we employ the notation of ref. [19]. The LO result becomes The sign in the denominator is chosen according to whether the particle is a boson (ǫ q ⇒ −) or a fermion (E p ⇒ +). Scheme dependence can be expressed as where δ was defined below eq. (A.4). Inserting the expressions listed in appendix A of ref. [19] for J m 1 m 2 n 1 n 2 and I m 1 m 2 m 3 n 1 n 2 n 3 n 4 n 5 , the remaining NLO contribution reads where È denotes a principal value; n B and n F are Bose and Fermi distributions; ǫ q ≡ q; E p ≡ p 2 + M 2 ; E pq ≡ p 2 + q 2 + 2pqz + M 2 ; E ± pq ≡ E pq | z=± ; ∆ στ ≡ ǫ q + σE p + τ E pq . It is observed from eqs. (A.6)-(A.9) that all terms contain non-trivial τ -dependence, i.e. that there is no constant contribution in G P (τ ) at NLO.
The main goal is to demonstrate how the spectral function behaves below the threshold, i.e. for 2M − ω ≫ α 2 s M , a regime that cannot be addressed with the methods of sec. 4. Physically, the NLO corrections originate from virtual heavy quark self-energy and gluon exchange contributions, as well as from real gluon emissions and absorptions.
The spectral functions corresponding to the master sum-integrals appearing in eq. (A.5) have been worked out in ref. [72]. Making use of the objects S i j (ω) defined there, the spectral function can be written as We have set here ǫ → 0 whenever the master sum-integral that it multiplies is finite. Expressing the result as a sum of a vacuum and thermal part, ρ P (ω)| raw = ρ P (ω)| vac + ρ P (ω)| T , (B.2) where the vacuum part is defined to include also the explicit T 2 visible in eq. (B.1), and inserting the expressions for the functions S i j (ω) from ref. [72], the vacuum part reads The last term in eq. (B.3) is divergent at the threshold ω ∼ 2M . There is a unique choice which avoids this at all T , namely adopting the pole mass scheme (δ ≡ 0) and resumming the explicit thermal correction into an effective mass modifying the LO result [74], We then re-interpret the mass appearing in ρ P (ω)| vac /M 2 as the thermal effective mass. The correction in eq. (B.5) is very small in practice for T ≪ M , M eff −M ∼ α s T 2 /M ; in particular it is smaller than the thermal mass originating from Debye screening in eq. (4.5), which is of order −α 3/2 s T . Consider finally ρ P (ω)| T . Going over to the notation of eq. (3.1), it can be represented as The main importance of eq. (B.6) lies in its first term, which remains non-zero below the threshold (ω < 2M ) and gives the dominant contribution in this regime. Writing ω = 2M + E ′ and assuming |E ′ | ≪ M , the restriction on the integration variable q reads θ(q + E ′ ). Therefore, for E ′ < 0, the result is ∼ α s C F (T /M ) 3/2 n B (|E ′ |) ∼ α s C F (T /M ) 3/2 e −|E ′ |/T .