Charm and beauty in the deconfined plasma from quenched lattice QCD

We present continuum extrapolated results of charmonium and bottomonium correlators in the vector channel at several temperatures below and above $T_c$. The continuum extrapolation jointly performed with the interpolations to have physical values of $J/\psi$ and $\Upsilon$ masses in the confined phase is based on calculations on several large quenched isotropic lattices using clover-improved Wilson valence fermions carrying different quark masses. The extrapolated lattice correlators are confronted with perturbation theory results incorporating resummed thermal effects around the threshold from pNRQCD and vacuum asymptotics above the threshold. An additional transport peak is modelled below the threshold allowing for an estimate of the diffusion coefficients for charm and bottom quarks. We find that charmonium correlators in the vector channel can be well reproduced by perturbative spectral functions above $T_c$ where no resonance peaks for $J/\psi$ are needed at and above 1.1 $T_c$, while for bottomonium correlators a resonance peak for $\Upsilon$ is still needed up to 1.5 $T_c$. By analyzing the transport contribution to the correlators we find that the drag coefficient of a charm quark is larger than that of a bottom quark.


I. INTRODUCTION
Heavy quark-antiquark bound states, quarkonia, have been proposed as a thermometer of quark gluon plasma in heavy ion collisions since they are formed at a very early stage of the collisions and may survive in the deep deconfined phase due to their hierarchically small sizes and large binding energies [1,2]. The suppression of quarkonia yields in AA collisions compared to those in the pp collisions have been observed in RHIC and LHC energies [3][4][5][6][7][8][9][10][11][12], however, its interpretation is still not very clear due to the interplay between the cold and hot nuclear effects [13]. Due to their nonperturbative features it thus is important to understand the fate of quarkonia in the hot medium from lattice QCD computations.
In addition, it was observed that open heavy mesons show an unexpectedly substantial elliptic flow that is comparable to that of light-quark mesons at RHIC [14,15] and LHC [16]. Moreover, heavy quarks are found to lose a significant amount of energy similar to light flavors at sufficiently high transverse momentum, while with decreasing transverse momentum an energy loss hierarchy is expected due to the dead cone effect, in short, heavier quarks suffer less energy loss. Experimentally the nuclear modification factor of open flavor mesons seems to support such a picture [17][18][19]. Phenomenological explanations of these phenomena require a modeling of the heavy quark diffusion in a hot and dense medium. This requires knowledge about the heavy quark diffusion coefficients D [20][21][22][23] which can be determined in lattice QCD calculations as they are encoded in the correlation and spectral functions of quarkonia in the vector channel.
In the heavy quark mass limit recent progress has been made to estimate the heavy quark momentum diffusion coefficient based on continuum extrapolated colorelectric field correlation functions [24][25][26]. The subleading quark mass corrections to this transport coefficient are proportional to a color-magnetic field correlator [27]. In the current study we will utilize full relativistic vector meson correlation functions to estimate the charm and bottom diffusion coefficients.
The spectral functions of quarkonia in the vector channel contain all information about the in-medium hadron properties like the dissociation temperatures of the corresponding bound states and heavy quark diffusion coefficients. However, the spectral function cannot be obtained directly from lattice QCD and is only related to lattice QCD computable Euclidean correlation functions. Investigations on quarkonium spectral functions extracted from two point correlation functions were started about two decades ago [28]. Since the lattice spacing has to be smaller than the inverse of the heavy quark mass to control lattice cutoff effects and the extraction of spectral function requires a large number of data points in the temporal direction of lattices, most studies in lattice QCD focus on charmonium spectral functions and correlation functions, where continuum extrapolated results only exist for those in the pseudoscalar channel [29]. Due to the much larger mass, bottomonium spectral functions in the relativistic formalism have only been studied on highly anisotropic lattices [30,31]. Recently, studies have been carried out using nonrelativistic heavy quark formulations on anisotropic [32] and isotropic lattices [33]. A review of current status of lattice studies on heavy quarkonium in extreme conditions can be found in [34,35].
The main goal of this work is to compare lattice correlators of both charmonium and bottomonium in the vector channel with those integrated from the perturbative arXiv:2108.13693v4 [hep-lat] 23 Dec 2021 spectral functions. The results will be used to investigate the thermal modifications of J/ψ and Υ and the diffusion coefficients of charm and bottom quarks. For this we will start with the construction of perturbative spectral functions, given in Sec. II. In Sec. III we present the lattice setup and describe how we perform the mass interpolation and continuum extrapolation. Section IV is devoted to comparing the lattice and perturbative results in the bound state region of the vector spectral function.
In Sec. V we analyze the transport peak mainly based on a Lorentzian ansatz. In the last section we draw the conclusion. Parts of the study have been presented in various conferences and workshops [36][37][38][39][40][41][42] and in the Ph.D. thesis of Anna-Lena Lorenz [43].

II. SPECTRAL FUNCTIONS IN THE VECTOR CHANNEL
The quarkonium spectral function cannot be obtained directly on the lattice, and it is related to the Euclidean mesonic two point correlation function via the integral equation, where is a temperature (T ) and frequency (ω) dependent integration kernel. Specific to the vector channel that we are considering in this work Γ H = γ µ and thus G H (τ ) = G 00 + G ii (τ ). G 00 is the zeroth component independent of the distance τ [44] and G ii (τ ) the sum of spatial components (we use the Einstein summation convention throughout this paper), The spectral function well above the threshold and around the threshold, on the other hand, can be obtained from the perturbation theory. For frequencies well above the threshold, the spectral function can be described by ultraviolet asymptotics [44]: with R c , a polynomial in α s up to five-loop order R c (ω 2 ) = r 0,0 + r 1,0 α s + (r 2,0 + r 2,1 l) α 2 s + r 3,0 + r 3,1 l + r 3,2 l 2 α 3 s + r 4,0 + r 4,1 l + r 4,2 l 2 + r 4,3 l 3 α 4 where l = ln μ 2 ω 2 andμ is the renormalization scale, whose range can be found in [44]. The coefficients r ij for the vector channel can also be found in [44]. The thermal contributions arising around the threshold can be obtained by applying pNRQCD calculations [45] as The threshold mentioned above is a frequency at which the free quark spectral function switches from vanishing to nonvanishing value [46]. At zero momentum it locates at 2M with M the quark mass. C > is a Wightman function, which is solvable for a real-time static potential from hard thermal loop resummation [47]. The two energy regimes are matched by modifying the pNRQCD result with a factor A match , so that it smoothly connects to the vacuum asymptotics at a certain point ω match . This matching procedure was successfully developed in the pseudoscalar channel in [29]. The resulting spectral function is valid down to frequencies around and above the threshold, and overestimates the regime 2M − ω α 2 s M . An exponential cutoff T is introduced to model the spectral function for the low frequencies. The whole spectral function then reads Note that the above perturbative calculations were carried out in the Minkowski space with metric (+−− −), where ρ V (ω) = ρ ii (ω) − ρ 00 (ω). Considering that ρ 00 (ω) ∼ ωδ(ω) [44], there would be no difference between ρ ii (ω) and ρ V (ω) around or above the threshold region. So in the following analysis we will use Eq. (6) to model ρ ii (ω) in this frequency region. As for the very low frequency region, ρ ii (ω) is supposed to have a transport peak. In the high temperature limit, the transport peak has the following form [46,48,49]: where χ q is the quark number susceptibility and M is the quark mass. In the interacting case, the δ-peak can be smeared into a Lorentzian peak with a finite width of η (drag coefficient) as [49] ρ trans According to [44] this expression overestimates the transport contribution for larger frequencies, so we multiply Eq. (8) with a cutoff function 1/ cosh( ω 2πT ) which becomes unity as ω → 0. Applying the Einstein relation one arrives at the Kubo formula which relates the spectral function and the heavy quark diffusion coefficient D as An estimation for the range of D is possible via its relation to the heavy quark momentum diffusion coefficient κ [50]: κ has been determined from lattice calculations before [24] and hints to a range of 2πT D ∈ [3.71, 6.91] at 1.5T c .
Recently a similar study in a much wider range of temperatures can be found in [51]. Figure 12 of [52], and [28] give an overview of different results for 2πT D for different temperatures.

III. LATTICE SETUP
The spectral function described in the previous section will be compared to continuum extrapolated lattice correlators of both charmonium and bottomonium in the vector channel. To realize the large and fine lattices required for our analysis, we choose the quenched approximation. The configurations are generated with a separation of 500 sweeps each consisting of one heat The aspect ratio is fixed at a certain temperature, and changes from 2 to 6 from the lowest temperature to the highest temperature. The lattice sizes and the number of measured configurations are listed in Table I. The lattice spacing is obtained using r 0 /a and with r 0 = 0.472(5) from [53]. As seen from Table I the lattice spacings used in our simulation are sufficiently small such that both bottom and charm quarks can be accommodated on the lattice. Since the tuning of hopping parameters κ to have the physical masses of J/ψ and Υ is nontrivial, five to six different values of κ at each lattice spacing have been used to compute the correlation functions. Our lattice setup and computations thus make the interpolation of correlators to the continuum limit and the case with physical masses of J/ψ and Υ possible. To compare the correlators from different lattices, they need to be renormalized. In the vector channel, there are different options regarding the renormalization. In addition to perturbative renormalization constants known up to two-loop order [55,56], there are nonperturbatively determined renormalization constants given in [57]. Another possibility is to take the continuum limit of renormalization independent ratios with the quark number susceptibility χ q given by the zeroth component of the vector correlator G 00 . Although the renormalization constants computed from [55,56] are comparable among each other, we decide on using the renormalization independent ratio, e.g. G(τ, T ) divided by the quark number susceptibility for the continuum extrapolation in this work. Here we chose the value of quark number susceptibility at T = 2.25T c (denoted as χ q ) as a normalization, as the susceptibility is more precise at higher temperatures. Using continuum extrapolated results for χ q /χ q at T = 2.25T c and the respective temperature T we obtain the correct normalization in the continuum. We also remark here that the extracted heavy quark diffusion coefficient [cf. Eq. (10)] is also renormalization independent.
In the next step, we need to ensure that the masses and temperatures on the different lattices match. At each lattice spacing we compute the correlation functions for five to six different values of hopping parameters κ related to the bare quark mass. The screening masses obtained at 0.75 T c 2 shows that, due to the nontrivial quark mass tuning, the different lattices do not have the same ground state vector meson mass m V (see Table II and Fig. 1). To overcome this problem, an interpolation in m V /T between the correlators computed at different values of κ is required. We adopt the ansatz where T = 2.25T c and (p, q, r) are fit parameters. We can see that this ansatz describes the data well as shown in the top plot of Fig. 2. Note that for the purpose of guiding the eye, hereafter the correlators are normally shown divided by G f ree (τ T ), a correlator computed from the vector spectral function in the noninteracting case [46,48]. They are calculated at quark mass of 1.5 GeV and 5.0 GeV for charmonium and bottomonium respectively. Then we insert physical J/ψ or Υ mass to the fitted curve and obtain the correlators at physical mass. As an example the interpolated correlators at physical mass of J/ψ on the 144 3 × 48 lattice are shown in the bottom plot of Fig. 2.
2. An example of the interpolation between correlators obtained using different values of hopping parameters κ to obtain the correlator at a physical J/ψ mass on 144 3 × 48 lattices. From the six measured κ values, we chose the four closest ones to the charm quark mass. For every point in τ T , we interpolated with the ansatz Eq. (12) to the value at m V = m J/ψ . Note that the ground state meson mass m V is obtained from two-state fits to the spatial correlators at 0.75 Tc and we use it for all temperatures in the mass interpolation. The top plot shows this interpolation at three example points, while the bottom plot shows the final result of the mass interpolation.
After the mass interpolation, we carry out a combined spline fit via which we are able to interpolate the correlators to the same points in τ T and extrapolate them to the continuum limit at the same time. We choose piecewise polynomials as an ansatz for the spline fit to our correlators, where Here d is the degree of the underlying polynomials, and n is the number of knots which is chosen by hand for different datasets. (τ T ) 0 is a reference point for τ T and it cannot be any of the knots, and t j denotes the position of the knot. a i and c j are spline coefficients which can be determined when fitted to the lattice data. To incorporate lattice cutoff effects one replaces the coefficients a i with certain functions, whose form depends on how the operators concerned are constructed on the lattice. As in this work the O(a)-improved Wilson (clover) fermions are used, one natural choice for the ansatz would be where b 1 and b 2 are fit parameters. To obtain the continuum extrapolated values for the correlators one just needs to take 1/N 2 τ → 0. To estimate the errors, the whole procedure is conducted on bootstrap samples. We show the charmonium and bottomonium correlators on each lattice and the continuum-extrapolated correlators at 1.5 T c in Fig. 3. One can see that lattice cutoff effects are larger at smaller distances. Using the forementioned method we obtain a reliable continuum extrapolation down to τ T = 0.1. To be more certain that no cutoff effects influence our analysis, we decide to start our fits at τ T ≈ 0.2.
The final continuum extrapolated correlators at 0.75, 1.1, 1.3, 1.5 and 2.2 T c for charmonium and bottomonium in the vector channel are summarized in the left and right plots of Fig. 4, respectively. We already can draw some conclusions from the correlators without extracting the spectral functions. We see that correlators at the short distance agree for all temperatures, meaning that the region mostly influenced by the vacuum asymptotic part of the spectral function does not depend on the temperature much. At larger τ , where the threshold region and the transport peak dominate the behavior, the correlators split, indicating that at least one of the two regimes is heavily temperature dependent. We can also see that comparing with bottomonium correlators, charmonium correlators have much stronger temperature dependence, especially for those at long distances which are relevant for the properties of resonance and transport peaks. This is a clear sign that charmonium suffers more thermal modifications than bottomonium. 3 To obtain more quantitative results, in the following section we analyze the lattice data using ansatz con-structed based on the perturbative spectral function described in Sec.II. The fits are conducted on every bootstrap sample to gain a correct error estimate. We then crosscheck the fit results by maximum entropy method (MEM) analyses. For the analyses in bootstrap a covariance matrix C k,l is needed. Since the bootstrap compromises the covariance matrix calculated from the continuum data, we instead use the covariance matrix of the finest lattice and rescale it to the continuum as where δG(τ k ) means the error of the correlators at distance τ k . Here the superscripts "cont" and "lat" stand for the continuum extrapolated and lattice results, respectively.

IV. COMPARISON BETWEEN LATTICE AND PERTURBATIVE RESULTS IN THE BOUND STATE REGION
In [29], it was found that the perturbative spectral function was well suited for describing the lattice correlators in the pseudoscalar channel after introducing corrections for systematic errors. However, the extraction of the information on the fate of J/ψ and Υ from the vector correlators is more complicated since the transport peak lying in the very low frequency region is not described by the perturbative spectral function. In the following we thus divide our analyses in the two different regimes. In this section, we investigate the bound state region, i.e. the intermediate and large frequency part of the spectral function, where the perturbative spectral function is valid. The small frequency part that contains the transport peak is then evaluated in the next section. The complete spectral function and the corresponding correlator are then given by respectively. Here ρ mod ii (ω) is the model spectral function [29] of the form The factors A and B correct for two sources of systematic errors that account for some of the quantitative differences in the comparison between the lattice data and perturbative results. On the lattice side, the renormalization might be off. This is taken care of by the overall normalization factor A. On the perturbative side, the relation between the pole mass and the MS mass is poorly determined which might lead to a slightly smaller or larger threshold location. This is taken care of by the mass shift B.
As the contribution from the transport peak to the correlator, i.e. G trans ii (τ ) is nearly τ independent, one can  thus look into the differences of correlators at neighboring points [58] In G dif f ii (τ /a) the contribution from the transport peak, mostly influencing the correlator at τ T ≈ 0.5, is suppressed. In this way one can directly confront the perturbative results with the lattice data of G dif f ii (τ ) as was done in [29] for correlators in the pseudoscalar channel.  Inserting the model spectral function Eq. (18) into Eq. (2), we obtain an expression for the correlator that is then fitted to the lattice data G dif f ii (τ ). The resulting parameters A and B/T are listed in Table III, and the comparison between the fits and lattice data is shown in the top panel of Fig. 5. One can see that the lattice data is well described by the ansatz Eq. (18). A is close to one and B is small, indicating that the perturbative spectral function is a suitable ansatz. The resulting spectral functions are shown in the top panel of Fig. 6 for charmonium (left) and bottomonium (right). We find that there is no need for a resonance peak to describe the charmonium data in the current temperature window, while for bottomonium one thermally broadened resonance peak represents the data better at T ≤ 1.5T c . The position of this peak almost does not change with temperature. At 2.25 T c the peak structure is gone. As a cross-check we perform MEM analyses using the fit results as default models. The results are shown in the bottom panel of Fig. 6. We can see that for both charmonium and bottomonium at all temperatures, output spectral func-  (17)] obtained from the above fits. The correlators are normalized by G f ree,dif f ii χ q /T 2 in all cases. Here, we observe a difference between the original and model correlators that hints to a transport contribution. For bottomonium this difference is small, while it grows for charmonium at higher temperatures.
tions almost overlap with the inputs, which suggests the perturbative spectral function a good ansatz to describe the lattice correlators. We remark here that the MEM analyses serve only as a consistency check, as the output spectral function from the MEM analyses is known to have large default model dependencies [59].

V. CHARM AND BOTTOM QUARK DIFFUSION COEFFICIENTS
When comparing the fit results of G dif f ii to the original lattice data (see the bottom panel of Fig. 5), we observe a difference more obviously at higher temperatures. This difference is a clear sign of a transport contribution. Qualitatively, we can already draw some conclusions, before analyzing the difference more closely in the following section. As can be seen from the bottom right plot in Fig. 5 the difference is very small for bottomonium. This indicates that the correlator is almost dominated by the bound state region. For charmonium as shown in the bottom left plot of Fig. 5, the difference is much larger compared to the case of bottomonium and it increases with growing temperatures.
In the previous section we obtained an expression for the spectral function in the bound state region which can describe G dif f ii well. We now construct a correlator that describes the contribution from the transport peak by the subtraction The obtained G trans ii (τ T ) for both charmonia and bottomonia at 1.1, 1.3, 1.5 and 2.25 T c are shown in the left and right plot of Fig. 7, respectively.
As expected from the findings in [49], we observe a very weak dependence of G trans ii (τ T ) on τ T at all temperatures, especially for the charm sector. The curvelessness of G trans ii (τ T ) implies a slender hope to reconstruct the transport peak without further information. This is verified when we model the transport peak using a Lorentzian ansatz (see Eq. (8)). We vary the heavy quark diffusion coefficient D in the range 2πT D ∈ [0. 2,4] and make use of the Einstein relation Eq. (9) to obtain the drag coefficient η. For the values of quark masses we use M c = 1.28 GeV and M b = 4.18 GeV [60]. However, it is found that all the choices can describe the lattice data equally well within errors and almost no differences among G trans ii (τ T )/G trans ii (τ T = 0.5) resulting from various values of 2πT D can be seen.

A. Relative magnitude of drag coefficients of charm and bottom quarks
Even though 2πT D cannot be determined by analyzing the curvature of G trans ii (τ T ) in τ T , we can still draw some conclusions on the relative magnitudes of the drag coefficients of charm and bottom quark by comparing charmonium and bottomonium correlators at the midpoint (τ T = 0.5). The procedure to determine the relative magnitude is illustrated as follows. For small ω/T , we expand the kernel and the 1/ cosh(ω/2πT ) cutoff term that is multiplied to Eq. (8) at the midpoint: where c 0 = 2 and c 1 = 3+π 2 12π 2 , for instance. With this and the Lorentzian ansatz we obtain the midpoint correlator: It is clear that the first term 2 arctan( ωcut η ) is the most dominant term and higher orders are negligible for small η/T . With these simplifications, the ratio of the midpoint correlators for charmonium and bottomonium is given by As the ratio of quark masses M b /M c is around 3 [60], and according to the top plot in Fig. 8    should be smaller than 1.
Since arctan(1/x) is a monotonically decreasing function of x for x > 0, we thus have i.e. the drag coefficient of a charm quark is larger than that of a bottom quark in the current temperature window. As one can also observe from the top plot of Fig. 8 the ratio increases with increasing temperature. This could indicate that the difference between η c and η b becomes smaller at higher temperatures.
Since the curvature of G trans ii (τ T ) can hardly provide any information on the heavy quark diffusion coefficient, in the following sections we turn to other two quantities: the midpoint correlator G trans ii (τ T = 0.5) and the thermal moments [58,61,62] through which there is a hope that the transport peak could be reconstructed.

B. Solving transport peak using midpoint correlators
In this section we consider the midpoint correlators which are shown already in Fig. 8. At the midpoint, the integration kernel simplifies to 1/ sinh( ω 2T ). This allows us to compare continuum extrapolated lattice data at the midpoint to the midpoint correlator obtained in the same way as in the above section from the model spectral functions including the Lorentzian ansatz. The difference is that here we use the full kernel function and insert M c = 1.28 GeV and M b = 4.18 GeV [60] for the required masses and assume an error of about 10% that accounts for the uncertainty in the definition of the mass in this approach. As an example we show the estimation of the drag coefficient of a bottom quark at 1.5 T c in Fig. 9 where the middle point correlator is shown as a function of the drag coefficient. The horizontal dashed line and the surrounding green band represent the lattice data of correlator divided by T χ q at the middle point τ T = 0.5, while the solid curve denotes the result obtained using the model spectral function of the Lorentzian form [cf. Bottom: ratio of the transport contribution to the correlator to the complete correlator. As it is seen, the transport contribution only makes up a small fraction of the correlator.
Eq. (8)] and the surrounding purple band denotes the uncertainly arising from the values of quark masses. The lower and upper bound of this overlapping region between the "lattice data" and "model" results shown in Fig. 9 can thus be regarded as the range for the estimated values of η/T . With the estimated range for η/T we are also able to determine D via the Einstein relation Eq. (9). Following this procedure η/T and 2πT D obtained for charm and bottom quarks at different temperatures are listed in Table IV.  The estimate of η/T described above obviously depends on the upper integration limit ω cut . Since the transport contribution described by the Lorentzian ansatz is only valid for small ω, we also investigated the effect of four different upper limits (ω cut = ∞, M, πT and T ) for the integration. For bottomonium, a plateau was reached, where each integration limit gave roughly the same values for η/T . Since the dependence on the integration limit was mild, we choose infinity as the upper bound. As for charmonium, the analysis is more complicated as the transport peak seems to be not well separated from the bound state or continuum region. Charm quark diffusion coefficient 2πT D obtained via the current approach decreases with increasing ω cut /T at all the temperatures considered in the current study, and the second and third largest values of 2πT D obtained using ω cut = M and πT are almost the same, but they are at most about 1.5 times that obtained using ω cut = ∞ at each temperature. This might indicate that the transport peak in the charm sector has a long tail stretching to the large ω region. We thus only show the obtained values of 2πT D and η/T for the charm quark obtained using ω cut /T = ∞ in Table IV. The results obtained for the charm quark thus suffer larger uncertainties than those for the bottom quark. As seen from Table IV the drag coefficient of a charm quark is larger than that of a bottom quark at each temperature (also hold for ω cut πT or M ), which is consistent with our estimate on the relative magnitude of η c and η b in the previous subsection.
C. Solving transport peak using thermal moments From the above section we learn that even using the midpoint correlators it is still difficult to obtain reliable transport coefficient for charmonium. In this section we try to tackle it from the so-called thermal moments which are defined as the Taylor coefficients [58,61,62] when expanding the correlator around the midpoint, where we have neglected the high order contributions in the last line. To get rid of renormalization we further build ratios with which the expansion could be rewritten as To obtain the moments we calculate the curvature from the data The curvature could also be expressed using the ratios Now we can get the ratios R 2n+2,2n H by fits based on Eq. (30). Note that the approximation is valid close to the midpoint, so the fit can only be conducted on points close to τ T = 0.5. At the same time, to have stable fits the fit intervals cannot be too small. For this reason we vary the lower limit τ min T of the fit interval and keep the upper bound at τ T = 0.5. We will see that after a first few values of τ min T , a plateau can be reached and we use the average over the plateau as our final estimate for R 2n+2,2n . Also for stabilities, we choose n = 1 for charmonium and n = 2 for bottomonium in Eq. (30).
We show the first thermal ratios obtained from fits in Fig. 10 as horizontal constant dashed lines. To extract the information on the transport peak, we also calculate the ratios using the spectral function Eq. (17) with varying η/T . For R 2,0 we have with where we use A and B from Table III  In previous subsections we have attempted to estimate 2πT D and η, by either analyzing the midpoint correlators G trans ii or the thermal moments based on Lorentzian ansatz for the transport peak. It is found that when using midpoint correlators we could obtain more reliable results for the bottom quark while when using thermal moments η/T (and also 2πT D) for the charm quark is more accessible. In general all the obtained results support that η c > η b holds true in the current temperature window. In this section we try to combine both results by taking only the most trustworthy ones, namely (i) In Table IV obtained by analyzing the correlator only at the midpoint, results for bottom quark at T ≥ 1.3 T c are chosen.
(ii) In Table V obtained by analyzing the curvature of the correlator via thermal moments, results for bottom quark at 2.25 T c are chosen while those for charm quark at T ≥ 1.5T c are chosen.
We plot the selected results in Fig. 11, as a summary of our analyses for the charm and bottom quark diffusion coefficients. In Fig. 11 2πT D for charm quark are shown at two temperatures, i.e. T = 1.5 and 2.25 T c as red bands, while 2πT D for bottom quark are shown at three highest temperatures, i.e. T = 1.3, 1.5 and 2.25 T c as blue bands. At 2.25 T c we have combined the estimated range for 2πT D of bottom quark obtained in Table IV and Table V. We remark here that the vertical lines denote the possible ranges of the diffusion coefficients arising from the uncertainty of the heavy quark mass used in our analyses, and they do not characterize the size of the statistical error. The charm and bottom quark masses used in our analyses range from the 90% to 110% of their values listed by the Particle Data Group (PDG), and if we amplify the range of the heavy quark masses, the estimated range of the diffusion coefficients would become broader. For example, if we use the heavy quark masses to be 80% to 120% of their PDG values, the estimated range of 2πT D will become [0.14, 0.75] instead of [0.18, 0.51] in Fig.9.
As seen from Fig. 11 there is no significant temperature dependence of 2πT D for both charm bottom quarks. The results of 2πT D are much smaller than those obtained in pQCD with α s ≈ 0.2. Results of 2πT D converted from the static heavy quark momentum diffusion coefficient [cf. Eq. (11)] obtained in [24][25][26] are also shown. These results do not show much temperature dependence as well, and are in general larger than 2πT D of both charm and bottom quarks obtained in the current study.
We also noticed that in Refs. [63,64] a holographic estimate gives 2πT D = 1, but it is for R-charge diffusion. The AdS/CFT calculations for heavy quark diffu- 66]. These estimates from AdS/CFT can be compatible with our results in this study given certain values of g 2 Y M N c .

VI. CONCLUSION
In this work we have computed charmonium and bottomonium correlators in the vector channel at various quark masses on four large and fine isotropic lattices in the quenched approximation at temperatures ranging from 0.75 T c to 2.25 T c . With these data we are able to interpolate the correlators to those with physical J/ψ and Υ mass on the lattice and perform extrapolation to the continuum limit. From our analyses we see a qualitatively good agreement between our continuum ex-  [25] and Altenkort et al. 2021 [26] are for 2πT D of a static quark rather than relativistic charm and bottom quarks considered in the current study. For comparison the results obtained from the (next-leading-order) NLO perturbative QCD with αs ≈ 0.2 (pQCD) are also shown.
trapolated lattice data and the correlator obtained from perturbative spectral functions constructed from matching pNRQCD calculations to vacuum asymptotics. We extended the analysis in [29] to the vector channel, where we divide the spectral function into the bound state region at larger ω and the transport region at small frequencies. To compare perturbation theory results and lattice data in the bound state region, we used the differences of neighboring points in the correlator and fitted a model spectral function accounting for systematic uncertainties. With this, the high frequency part is well described by the perturbative spectral function, as only mild modifications are needed. We find that for charmonium the perturbative spectral function without any resonance peak is sufficient to describe the continuum extrapolated lattice data. For bottomonium on the other hand, a thermally broadened resonance peak is needed to describe the lattice data for temperatures up to 1.5 T c . Our results of spectral functions have been cross-checked using the MEM and were found to be a good description to the correlator.
For the transport contribution we analyzed the midpoint correlators and the thermal moments based on Lorentzian ansatz. We find that the drag coefficient of a charm quark is larger than that of a bottom quark. We also managed to constrain the charm and bottom quark diffusion coefficient D to a possible range. Since different methods have their own (dis)advantages we combine the results by taking only the most reliable results from each method and summarize in Table VI. We find that charm and bottom quark diffusion coefficients, obtained at physical quark mass in this study, are smaller than those converted from lattice calculations of heavy quark momentum diffusion coefficient [24][25][26]. The reason of such discrepancy can be that the heavy quark momentum diffusion coefficient is calculated in the heavy quark mass limit and also the studies [24][25][26] only consider the leading term. Recently the subleading terms of heavy quark momentum diffusion coefficient in T /M have been worked out in [27] and one of them can be estimated from a color-magnetic correlator. This correlator needs to be studied in the future on the lattice and may bring the result closer to our estimates here.

Charm
Bottom All data from our calculations, presented in the figures of this paper, can be found in [67].