Charmonium dissociation and heavy quark transport in hot quenched lattice QCD

We study the properties of charmonium states at finite temperature in quenched lattice QCD on large and fine isotropic lattices. We perform a detailed analysis of charmonium correlation and spectral functions both below and above Tc. Our analysis suggests that the S wave states disappear at about 1.5 Tc. The charm diffusion coefficient is estimated and found to be approximately 1/{\pi}T at 1.5Tc {\leq} T {\leq} 3Tc.


Introduction
In the hot and deconfined QCD medium it has been conjectured that, unlike light mesons, the heavy mesons might survive and dissociate only at some higher temperature due to Debye screening [1]. Therefore the suppression of their yield in nucleus-nucleus compared to proton-nucleus or protonproton collisions may serve as probe for the properties of the medium. The experiments carried out at the SPS and LHC at CERN and the RHIC at BNL have indeed observed J/ψ suppression [2]. The interpretation of experimental data, however, is not as straightforward and in order to disentangle the cold nuclear matter effects on the one hand and hot medium effects on the other, it is crucial to have a good understanding of the behavior of heavy quarks and quarkonia in the hot medium. In addition, a substantial elliptic flow of heavy quarks has been observed [3]. As the heavy quark diffusion coefficient D can be related to the ratio of shear viscosity to entropy density η/s [4,5], a non-perturbative estimate of D is highly desirable.
From the theoretical point of view, the meson spectral function at finite temperature [6] contains all the information on the hadron properties in the thermal medium. This explicitly includes the information on bound states and transport properties. As a consequence the non-perturbative computation of spectral functions is of key interest.
In the following we report on a recent success in this direction in quenched lattice QCD. We present the spectral functions in the pseudoscalar and vector channels and estimate the heavy quark diffusion constant. See [7] for a recent overview and [8] for a more detailed and extended discussion.

Euclidean correlation functions and the charmonium spectral function
At vanishing momentum charmonium correlation functions are defined via current operators by: Here J H is e.g. a local mesonic current operatorq(τ, x)Γ H q(τ, x), where Γ H = γ i , γ 5 for vector (V ii ) and pseudo-scalar (PS ) channels, respectively. Notice ... T denotes a thermal average. On the lattice the temperature T is related to Euclidean temporal extent aN τ by T = 1/(aN τ ), where a is the lattice spacing. The above correlation function is exactly related to the spectral function via: where the kernel K is given by K(τ, T, ω) = cosh(ω(τ − 1 2T ))/sinh( ω 2T ). The spectral function ρ(ω) contains all the information on the in medium properties of hadrons and its calculation therefore is our primary goal. The dissociation temperature can be read off from the deformation of the spectral function, while the heavy quark diffusion constant D is related to the vector spectral function via the Kubo formula: where χ 00 is the quark number susceptibility.
To study and estimate the charmonium spectral functions in the pseudoscalar and vector channels as temperature is increased, we first examine the behavior of the correlator itself in a number of ratios and differences. It should be clear that in this fashion we can only obtain a qualitative idea on the properties of the spectral function. The estimates are then used to construct the necessary default models for an analysis based on the maximum entropy method (MEM) [9]. Using MEM we invert Eq. 2 and compute the desired spectral function. A detailed analysis of our systematics and statistical uncertainties can be found in [8].

Numerical results
For our numerical investigation the standard Wilson plaquette action was used in the gauge sector, while for the (valence) charm quarks we choose the non-perturbatively O(a) improved Wilson clover fermion action. Note, here the mass of vector meson is tuned to be almost the physical J/ψ mass. Employing a fixed scale approach we measure correlation functions on very fine quenched lattices sized 128 3 × 96, 128 3 × 48, 128 3 × 32 and 128 3 × 24. These correspond to the temperatures 0.73 T c , 1.46 T c , 2.20 T c and 2.93 T c , whereby T c 270MeV in quenched QCD. Note, the required lattice parameters are given in Tab. 1.

Charmonium correlation functions
To begin our analysis we show the correlation functions in the pseudoscalar and vector channels in physical units in Fig. 1. Clearly all correlators agree throughout most of the available Euclidean time distance except when they near their respective midpoints. This small effect in the large distance region is the essential piece for understanding the thermal modification of the bound state and transport regions of the spectral function, while the correlation function at small to intermediate distances is dominated by the continuum contribution. As a consequence a more sophisticated analysis highlighting the thermal modifications by taking e.g. appropriate ratios or differences is necessary.
To analyze the properties of the underlying spectral functions on the correlator level more thoroughly we turn to studying it in conjuncture with the so called reconstructed correlator. In a reconstructed correlator the temperature dependence of the kernel in Eq. 2 is taken into account, while the spectral function remains unchanged. It is particularly interesting to compare a zero or low temperature spectral function with a thermal one and as it turns out in this case it is possible to calculate the reconstructed correlator analytically as follows, first one constructs [10]: At this point one would generally need the spectral function ρ(ω, T ) at a reference temperature T and consequently the evaluation of Eq. (4) would suffer from the uncertainty of the determination of that spectral function. However there is a useful and exact relation, which is a generalization of  Ref. [11]: where and N τ are the number of time slices in the temporal directions at temperature T and T , respectively. τ denotes the time slice of the correlation function at temperature T while τ denotes the time slice of the correlation function at temperature T . The sum of τ on the right hand side of Eq. (5) starts from τ = τ with a step length N τ to the upper limit N τ − N τ + τ. After putting ρ(ω, T ) into both sides of the above relation and performing the integration over ω, one immediately arrives at [8]: with this relation we are therefore capable of computing the reconstructed correlator directly from our low temperature results at N τ = 96 without having to extract the spectral function, thereby bypassing a systematic source of error in this quantity. Given this insight we use Eq. 6 to compute the ratio of the thermal and reconstructed correlators G(τ)/G rec (τ) at 1.46T c , 2.20T c and 2.93T c in the pseudoscalar and vector channels, where the results are shown in Fig. 2. Concentrating first on the pseudoscalar case in Fig. 2(left) we observe that the ratios decrease monotonically at the available temperatures as the Euclidean time separation increases. Note, at the largest distances, the ratios deviate from unity by about 5%, 8% and 12% at 1.46, 2.20 and 2.93 T c , respectively. Turning to the right of Fig. 2 and therefore the V ii channel we see an entirely different behavior. Here the ratios increase monotonically with increasing distance, reaching a deviation of about 16% at the large distances. These results are interesting as there is no emerging transport contribution in the PS channel, therefore our results hint at the decrease in this channel being due to bound state modification, while a large part of the modification in the V ii channels is due to the emerging transport contribution.
diff sub Figure 3. The ratio G diff (τ, T )/G diff rec (τ, T ) (G sub (τ, T )/G sub rec (τ, T )) for the PS (left) and V ii (right) channels as a function of the Euclidean distance τ at T = 1.46, 2.20 and 2.93 T c . The superscript "x" denotes either "diff" or "sub".
To study this more closely we define two ratios, which are approximately independent of a zero mode contribution: The first equals the ratio of the time derivative of the correlators to the time derivative of the reconstructed correlators at τ + 1/2 [12], while the second is the ratio of midpoint subtracted correlators [13]. As noted above both ratios are approximately independent of a zero mode contribution, consequently they should be almost independent of transport phenomena. The corresponding results in the PS channel are given in Fig. 3(left) and we find that both G di f f (τ)/G di f f rec (τ) and G sub (τ)/G sub rec (τ) give very similar results. Throughout the available temperature range we retain the decreasing trend with increasing distance from the analysis of G(τ)/G rec (τ), although the magnitude of the deviation has decreased, the values at the largest distance being shifted up by about 3% at both 1.46 T c and 2.20 T c and about 6% at 2.93 T c . As there is no transport contribution in this channel we would have expected fairly small changes from the behavior of G(τ)/G rec (τ) to G di f f (τ)/G di f f rec (τ) or G sub (τ)/G sub rec (τ) and indeed we see this confirmed. Turning to the V ii channel on the right of Fig. 3 we first note that G di f f (τ)/G di f f rec (τ) and G sub (τ)/G sub rec (τ) give similar results as in the PS case. The behavior compared to G(τ)/G rec (τ) on the other hand has been greatly modified. As such the data for 1.46T c instead of increasing now is almost flat, while 2.20T c and 2.93T c show a decreasing trend similar to the PS case. From this we infer that most of the temperature dependence in G(τ)/G rec (τ) in the V ii channel is in fact due to the transport contribution.
Another way to look at the data is to take the difference G(τ) − G rec (τ), implying: consequently this difference between the measured correlator and the reconstructed correlator provides information on the difference of spectral functions below and above T c and the corresponding results are given in Fig. 4. Note in both V ii and PS channels we observe that G(τ) − G rec (τ) increases with τT at the available temperatures. In the PS case on the left of Fig. 4 however all values are negative, while for the V ii on the right they are positive. As G(τ) − G rec (τ) probes the difference of the two spectral functions Δρ(ω) the results in the PS channel require Δρ(ω) < 0 in a certain ω region. With no transport contribution in this channel and given the positivity condition of the spectral functions this implies the low temperature bound state peak overcompensates any remnants of the peak at finite temperature. At the same time in the V ii channel the positive values indicate a more intricate situation, as the emergence of transport leads to a positive contribution, while an overcompensation of the bound state would lead to a negative part. Note, the weakest τ-dependence is seen in the V ii channel at 1.46T c . Here we observe an almost flat behavior with positive values throughout. Assuming that there is a complete cancellation of bound state peaks at this temperature in Δρ(ω) we may try to fit the transport contribution to the Ansatz defined in Eq. 3. In this case there is only one free parameter and the uncertainty in the heavy quark mass, for M = 1.0GeV and M = 1.8GeV we obtain: Even though it is not convincing to use a one parameter fit here to obtain a robust value for the heavy quark diffusion coefficient, these results are in the right ballpark and provide useful input for our subsequent MEM analysis.

Charmonium spectral functions via MEM
Clearly the comparison of the measured correlator with the reconstructed correlator can only give a rough idea of the magnitude of any medium effects at a certain temperature. To really explore the properties of the thermal charmonium states, one has to go to the level of spectral functions. Here we do this by using the maximum entropy method (MEM) to extract the spectral functions directly from the correlator data. In the MEM analysis presented here, we use as number of frequencies N ω = 8000, the minimum energy aω min = 0.000001 and implement the modified kernel EPJ Web of Conferences 00061-p.6 K * = tanh(ω/2) · K [14,15] in order to explore the low energy behavior of spectral function [16]. The default model in the PS channel is a normalized free lattice spectral function, while in the V ii channel it is a normalized free lattice spectral function with a transport peak based on Eq. 3 and the diffusion coefficient estimated above.
The results of this analysis are given in Fig. 5, where we show the PS (left) and V ii (right) spectral functions divided by ω 2 at the available temperatures above T c and also below T c . Note, throughout we show the spectral functions in PS and V ii channels with error bands given by the statistical uncertainties, while the error of the peak positions, see [8], are given as horizontal error bars. In Fig. 5(left) we observe that at 0.73 T c the spectral function in the PS channel has large uncertainties in the amplitude at the point which corresponds to the ground state peak location in the mean spectral function. However, even at the lower end of the error band, the amplitude is still larger than the peak amplitudes at the higher temperatures within the errors. Unlike the large uncertainties shown in the amplitude of the peak height, the peak location of the ground state peak at 0.73 T c is well determined. A Jackknife analysis yields m η c = 3.31(4) GeV. At 1.46 T c this peak is shifted by about 0.8 GeV to around 4.1 GeV. At 2.23 T c there is hardly a peak structure that can be identified within the statistical uncertainties. At 2.93 T c the spectral function flattens further. Thus this picture suggests that η c is melted already at 1.46 T c .
Turning to the right of Fig. 5 we analyze the resonance part of the spectral function in the V ii channel. We observe that the peak location of the spectral function at 0.73 T c does not have an overlap with the peak location of the spectral function at 1.46 T c and the amplitudes between these two differ a lot (see horizontal error bars). At both 2.20 T c and 2.93 T c there is hardly any peak structure. As such this picture indicates that also J/ψ is already dissociated at 1.46 T c .
Focusing on the very low frequency part of the spectral function given in the vector channel, we show the spectral function divided by ωT in Fig. 6(left). Notice here the statistical uncertainties on the amplitude of the peak are relatively small. The charm diffusion coefficient is related to the amplitude ICFP 2012 00061-p.7  Figure 6. Left: statistical uncertainties of transport peaks at T > T c . Right: the resulting charm diffusion coefficients. The boxes stand for statistical error estimated from Jackknife method while the bars stand for systematic uncertainties from MEM analyses. To compare also the results from perturbation theory [17], Langevin dynamics [5], Ads/CFT [18] and HQET [21,22] are given.
of the transport peak at vanishing frequency through the Kubo formula (Eq. 3). The corresponding estimates for the heavy diffusion coefficient D are summarized in the right plot of Fig. 6. The boxes denote the statistical uncertainties and the error bars reflect systematic uncertainties obtained from the analyses discussed in [8]. Note the charm diffusion coefficient obtained at 1.46 T c is the most reliable one among the three temperatures above T c since more prior information is known at this temperature. At 2.20 T c and 2.93 T c , due to the lack of precise prior information and also less number of data points that can be used in the MEM analyses, the uncertainties on the charm diffusion coefficient thus might be underestimated.
For comparison we highlight results from alternative approaches to the calculation of the heavy quark diffusion coefficient in Fig. 6(right), showing the perturbative next-to-leading order results of [17], those from a computation in AdS/CFT [18], as well as those from Langevin dynamics derived in [5]. Additionally we show results from a recently proposed method to compute the momentum diffusion coefficient κ from a purely gluonic operator derived from heavy quark effective theory (HQET) [19,20]. Hereby κ is related to the heavy quark diffusion coefficient by κ = 2T 2 /D. First calculations from two groups on this operator agree fairly well [21,22] and compare favourably with our results.

Summary
In summary our analysis of charmonium correlation functions on large isotropic quenched gauge configurations using Wilson clover fermions implies that both J/ψ and η c have dissociated by a temperature of 1.46T c . At the same time we estimated the charmonium diffusion coefficient and found a value 2πT D 1.8 − 2.3 in the available temperature range.