Continuum extrapolation of quarkonium correlators at non-zero temperature

In the investigation of in-medium modifications of quarkonia and for determining heavy quark diffusion coefficients, correlation functions play a crucial role. For the first time we perform a continuum extrapolation of charmonium and bottomonium correlators in the vector channel based on non-perturbatively clover-improved Wilson fermions in quenched lattice QCD. Calculations were done at 4 different lattice spacings, spatial extents between 96 and 192, aspect ratio from 1/6 to 1/2, for 5 temperatures between $T/T_c = 0.75$ and $T/T_c = 2.2$. We interpolate between different quark masses to match to the same vector meson mass over different lattice setups. Afterwards we extrapolate the renormalized correlators to the continuum. While we find a strong temperature dependence for charmonium, bottomonium states are only slightly affected at higher temperatures.


Introduction
In the investigation of properties of the quark gluon plasma (QGP), which is created at heavy ion collisions, quarkonium bound states are of great interest. Due to the color Debye screening, these states are expected to melt at a certain dissolution temperature and, therefore, suppression of J/ψ yields at heavy ion collisions compared to p-p-collisions can be a signal for the existence of a quark gluon plasma [1].
Indeed, J/ψ-suppression has been found at all three large colliders SPS, RHIC and LHC and also Υ suppression has been reported at the LHC [2][3][4][5][6][7]. Nevertheless, for a better theoretical interpretation of the experiments, one needs more information about the in-medium modification of quarkonium bound states.
In addition, measurements at RHIC and LHC revealed a non-zero elliptic flow for heavy quarks [8][9][10]. This indicates a relaxation towards thermal equilibrium and gives evidence to a collective motion of heavy quarks. For a better understanding of these hydrodynamic effects, knowledge about the transport properties of heavy quarks within the medium is needed. Especially, the heavy quark diffusion coefficient D is of interest, as it is related to the energy loss of a heavy quark in the medium. So far leading order perturbative calculations are not in agreement with the experimental data [11] and next-to-leading order calculations show that one can expect large corrections by higher orders [12]. Therefore, first principle studies, such as lattice QCD, are necessary to understand the nonperturbative effects in these quantities.
Both, in medium modification as well as transport properties are encoded in the spectral function of quarkonia. Unfortunately, the spectral function and the corresponding Euclidean correlator are related by an integral and a lot of details of the spectral function are averaged away. Moreover, due to the limitation of the τ-direction, only few lattice points are available. Therefore, very fine lattices as well as good statistics are needed to extract structures of the spectral function. So far, there were several attempts using very fine quenched lattices, though no continuum extrapolation has been realized [13][14][15]. In [16] a continuum extrapolation for the heavy quark momentum diffusion coefficient based on a color electric correlator has been carried out. However, this approach is limited to the heavy quark mass limit and as this analysis is not based on meson correlators, in-medium properties of bound states cannot be extracted by this approach.
For the first time, we have carried out a continuum extrapolation for charmonium and bottomonium correlation functions in quenched lattice QCD. In this paper we will show the details and steps that have been been used to perform the continuum extrapolation. A similar analysis for light quarks can be found in [17,18].

Meson correlation functions and the spectral function
The main object of interest in this study, the Euclidean mesonic correlation function, is defined as where the mesonic current is given by Here, ψ are the fermion fields and Γ H = γ 5 , γ µ , 1, γ 5 γ µ corresponds to the meson channel we are looking at (pseudo-scalar (P), vector (V), scalar (S), axial-vector (AV), respectively). In this study, we perform the continuum extrapolation for the example of the vector channel, but, except for the renormalization, the process works analogously for the other channels. For results for the pseudoscalar channel, see [19]. The above correlator is projected to finite momentum using the lattice Fourier transformation, with n = τ,x = (x, y, z) andp = (p x , p y , p z ) for temporal correlators or n = z,x = (x, y, τ) and p = (p x , p y , p τ ) for spatial correlators, respectively. In the following we restrict top = 0 and, except for the mass determination, to temporal correlators. The temporal correlator is connected to the spectral function via The spectral function contains all information on the medium properties of mesons. Whereas many information, like dissociation temperatures, can be extracted from the deformation of the spectral functions at higher temperature, the diffusion constant D can be directly extracted from the vector spectral function: Thereby, the index at the spectral function stems from the µ-index from Γ V = γ µ and χ q is the quark number susceptibility given by While this is in principle constant in τ due to current conservation, lattice cut-off effects are noticeable for short τ-ranges. Therefore we extract χ q from G V 00 (τT = 0.5).
Since the quark number susceptibility is affected by the same renormalization as G V , it is used to circumvent the uncertainties of the renormalization constants using renormalization independent ratios by dividing all correlators by the quark number susceptibility at T = 2.25T c .
In order to remove the exponential drop of the correlator we normalize all correlators by the free correlator which is given by The choice of m q is arbitrary in this case, as this normalization is only temporary and will be removed after the final continuum extrapolation. In the following we use m q = 1.5GeV for charmonium and m q = 5GeV for bottomonium.

Lattice setup
The correlators were measured on quenched configurations generated with the standard Wilson gauge action on large isotropic lattices. All configurations are separated by 500 sweeps which consist of  [20] with updated coefficients to include a wider β-range [19]. The correlators have been calculated using non-perturbatively clover-improved Wilson fermions [21]. We computed the correlators for six different κ-values for all lattices except for the N σ = 120 lattice, where we used four different κ-values. The κ-values have been roughly distributed between bottomonium and charmonium using the vector screening masses from the spatial correlator in the deeply confined phase at 0.75T c for tuning. The screening masses have been determined by averaging over the masses from correlated two state fits for different fit intervals. As a fit Ansatz we use G(n τ ) = A 1 cosh(m 1 (n τ − N τ /2)) + A 2 cosh(m 2 (n τ − N τ /2)).
For the fit intervals, we fixed the upper limit to N σ /2 − N σ /24 and gradually increased the lower limit. Afterwards we averaged over all fit results that lied within a plateau of the screening mass and in the corresponding χ 2 /d.o. f. of the fit (See figure 1). Due to the high costs of calculating the correlators, the quark mass tuning is not optimal. This will be corrected by a quark mass interpolation later on.
For the finest lattice, we used four additional quark sources for the two κ-values closest to bottomonium and charmonium to increase the statistics at these important data points. κvalues closest to charmonium. We find a very strong mass dependence: Even for small differences in the corresponding screening mass, we get a rather large deviation of the correlators. In the right plot we show the correlators whose measured screening mass lies closest to the J/ψ mass. As it can be seen, the ordering of the correlators does not match to the ordering of the lattice sizes. This makes it impossible to perform a continuum extrapolation based directly on the lattice correlators. As this effect stems from the slight uncertainties of the quark mass tuning, we have to make sure that all correlators are matched to exactly the same quark mass. This is realized by interpolating between the correlators corresponding to different screening masses. To do so, we plot the correlators at each τT separately against the measured screening mass from the T = 0.75T c lattice (see right plot of figure 3). Note that this means that also for higher temperatures we use the screening mass from the lowest temperature lattice.

Interpolating to the same vector meson mass
The actual interpolation is done by fitting the data with an exponential of a quadratic function: Then, the new interpolated correlators are computed by inserting the same reference mass m ref /T for each τT and each lattice size into the above equation (vertical line in the right plot of figure  3). As also the temperature differs slightly between the lattices, we choose m/T from the finest lattice as the reference mass for bottomonium and charmonium.. Because of the large differences between the correlators of different mass, we did not use all six available different masses. Instead we only used the closest four correlators for charmonium and the closest three for bottomonium.
In principle this method can be performed for each configuration individually. However, due to fluctuations between the individual measurements, this gives larger errors. Instead we used bootstrap samples with a size of 1000 samples for the mass interpolation. When the choice of the samples is the same for each τT , the bootstrap method preserves the correlations within the correlator.    For the three coarsest lattices, the data points arise from a B-spline interpolations to match to the data points of the finest lattice at each τT . Right: The mass interpolated lattice correlators and the corresponding continuum extrapolation. We also show the data point for T 2 /χ q,cont. , at zero distance. In both plots, the correlators have been normalized with the free correlator defined in equation 7 with a quark mass m q = 1.5GeV and with the quark number susceptibility χ q at T = 2.25T c .
With this method we are now able to calculate the correlator for any arbitrary meson mass between charmonium and bottomonium and can retune the correlators for all lattice sizes after the actual computation of the correlator and, therefore, can compensate uncertainties that stem from slight errors in the tuning of the quark masses.

Extrapolation to the continuum
Having performed the mass interpolation on the three coarsest lattices (as we used m ref /T from the finest lattice, we do not need to interpolate the finest lattice), we end up with well ordered correlators (right plot in figure 4) and we are ready to perform the continuum extrapolation. This is realized using a linear extrapolation in 1/N τ 2 on each bootstrap sample (see left plot of figure 4). The three coarsest lattices have been B-spline interpolated in order to estimate all correlators at distances available on the finest lattice. Due to the cut off effects at short distances, we restricted the extrapolation to values larger than τT = 0.08.
The final results for all temperatures can be found in figure 5. We find that the correlators agree for short distances of τT at almost all temperatures. Only at T = 2.25T c for charmonium we find deviations. Nevertheless, for larger distances we see clear temperature effects for charmonium. On the contrary, for bottomonium we find temperature effects only for T = 2.25T c , while for lower temperatures, the bound state seems to be independent of the temperature.

Conclusion
For the first time, we were able to perform a continuum extrapolation of temporal Euclidean vector meson correlation functions for charmonium and bottomonium in quenched lattice QCD. To do so, we interpolated the correlators to exactly the same vector meson mass and then performed a continuum extrapolation linear in 1/N 2 τ . Within the temperature range T = 0.75T c − 2.25T c , we observe strong temperature effects for charmonium, but only slight effects for bottomonium at T = 2.25T c .
The continuum extrapolated correlators are now ready for further investigations: Using qualitative motivated fit Ansätze one can find which parts of the spectral function contribute to the temperature dependence and which not. Furthermore, statistical Bayesian methods, such as Stochastical Analytical Interference (SAI) [22,23], the Stochastic Optimization Method (SOM) [24] or the Maximum Entropy Method (MEM) [25] can be used to extract details of the spectral function. These methods have already been tested with the lattice data of the finest lattice [26] and will be applied to the actual continuum correlators.