Quark Mass Dependence of Heavy Quark Diffusion Coefficient from Lattice QCD

We present the first study of the quark mass dependence of the heavy quark momentum and spatial diffusion coefficients using lattice QCD with light dynamical quarks corresponding to a pion mass of 320 MeV. We find that, for the temperature range 195 MeV $<T<$ 293 MeV, the spatial diffusion coefficients of the charm and bottom quarks are smaller than those obtained in phenomenological models that describe the $p_T$ spectra and elliptic flow of open heavy flavor hadrons.

Introduction.-Heavy-ionexperiments at high energies hint toward a rapid thermalization of heavy (charm and bottom) quarks.These observations are surprising since the relaxation time of a heavy quark immersed in a quark-gluon plasma (QGP) is expected to be M/T times larger than the relaxation time of the light bulk degrees of freedom constituting the QGP, where M is the heavy quark mass [1,2] and T is the temperature of the QGP.These experimental observations corroborate the picture that the QGP created in such high-energy heavy-ion collisions is an almost perfect fluid, see Refs.[3][4][5].This makes the heavy quark diffusion coefficient one of the fundamental transport properties of the QGP, along with other transport coefficients such as shear and bulk viscosities.
The heavy quark momentum diffusion coefficient κ is defined as the average momentum transfer squared to a heavy quark from the medium per unit time, and thus characterizes the kinetic relaxation of heavy quarks toward thermal equilibrium.One can also define the spatial heavy quark diffusion coefficient D s in terms of the conserved net heavy flavor number through the usual Kubo formula.This spatial heavy quark diffusion coefficient will depend on M , and also has a well defined limit when M goes to zero.Close to equilibrium, in the limit M ≫ T , the momentum and spatial heavy quark diffusion coefficients are related [2,6], where ⟨p 2 ⟩ is thermal averaged momentum squared of the heavy quark.
Experimentally measured p T spectra and the elliptic flows of open charm and open bottom hadrons provide information on the degree of thermalization of the heavy quarks in the QGP, see Refs.[3][4][5] for reviews.D s can be estimated by fitting these experimental measurements using phenomenological transport models, see Refs.[3][4][5].
A key ingredient of these transport models is an effective momentum-dependent heavy-quark diffusion coefficient, which in turn may depend on some effective in-medium cross sections.The D s appearing in the Kubo formula can be obtained as the zero-momentum limit of this effective diffusion coefficient.
Very recently, D s has been calculated in 2+1 flavor QCD for infinitely heavy quarks [7].This QCD result for the infinite-mass limit turns out to be smaller than the phenomenological estimate of D s for the charm and bottom quarks.This begs the question of whether or not these discrepancies arise solely due to the mass dependence of D s .In this Letter we address this question by presenting first lattice QCD calculations of the massdependent D s with light dynamical quarks.
Theoretical framework.-ForM ≫ T , κ can be calculated using heavy quark effective theory [6,8,9].In this framework κ is expressed in terms of the correlation function of chromo-electric (E) and chromo-magnetic (B) fields connected by fundamental Wilson lines [6,8,9].
In this effective theory where κ E,B (T ) = 2T lim ω→0 [ρ E,B (ω, T )/ω] [6,9], ρ E,B are the spectral functions corresponding to the E and B field correlation functions, and ⟨v 2 ⟩ is the mean-squared thermal velocity of the heavy quark [6].The quark mass dependence of κ enters through ⟨v 2 ⟩.At the leading order in 1/M , ⟨v 2 ⟩ = 3T /M .In this way, κ B controls the quark mass dependence of κ.
Lattice QCD calculations of κ B rely on accessing ρ B (ω, T ) from the correlator [6] where β = 1/T is the inverse temperature, τ is the Euclidean time separation of the B operators, and U (τ 1 , τ 2 ) is a thermal Wilson line connecting the B-fields located at Euclidean times τ 1 and τ 2 .ρ B is related to G B (τ, T ) via the integral equation Lattice QCD setup.-In the present calculation we use the same lattice QCD setup and ensembles that were used for the calculation of κ E [7], specifically, 2+1 flavors of quarks in the Highly Improved Staggered Quark fermionic action [10] and the tree-level improved Lüscher-Weisz gauge action [11,12] with physical values of the kaon and 320 MeV pion masses and at T =195, 220, 251, 293 and 352 MeV.At each temperature (except 352 MeV) we use three lattice spacings (a) to carry out continuum extrapolations (a → 0) of G B .Further details are provided in Supplemental Material [13].
The B-fields are discretized on the lattice as For measurements of G B we use a Symanzik-improved version [14] of gradient flow [15][16][17][18].Guided by our experience [19,20] and perturbative QCD (pQCD) [21] we limit the gradient-flow time (τ F ) within the range √ 8τ F < τ /3.We find that gradient flow improves the signal-to-noise ratio of G B .
At 1-loop level [22] in pQCD G B has a nontrivial anomalous dimension, and gradient flow serves as a nonperturbative renormalization scheme for G B .The continuum-extrapolated G B are renormalized in the gradient-flow scheme at the scale µ F = 1/ √ 8τ F .The renormalization-group invariant physical correlator, G phys.

B
, is obtained via the one-loop pQCD matching [23] G phys.
(5) The power-law corrections arising from mixing with highdimension operators are removed through the τ F → 0 extrapolations at each T .The matching function, Z match , involves three components: matching from the gradientflow to the MS scheme at a scale μτF , matching between MS-renormalized thermal QCD to the static quark effective theory at a scale μT , and running of the anomalous dimension of the operator from μT to μτF .If Z match were known up to all orders in pQCD, the dependence of G phys.

B
on the scales μT and μτF would exactly cancel.Since Z match is known only up to one loop (NLO), we estimate the uncertainty from unknown higher-order effects in the matching by varying the values of each scale; for μT we consider the two choices μT = 2πT or 19.18T , and for μτF we consider μτF = µ F or 1.4986µ F .Further details and the expression for Z match are given in Supplemental Material [13].
Data analysis.-TheB-field correlator G B scales with a strong negative power of τ since ρ B ∝ ω 3 at large ω.To mitigate this as well as the lattice artifacts and distortion due to gradient flow [7] we normalize G B with , where G LO B is the tree-level pQCD results for G B at nonzero lattice spacing and gradient-flow time, C F is the Casimir factor, and g 2 is the strong coupling.For brevity we suppress the arguments of the normalization correlator.
Our lattice data are not G B (τ, T, τ F ) directly, but G B (τ, T, N τ , τ F ) with N τ the number of points across the lattice time direction, which is related to the lattice spacing via N τ a = 1/T .Therefore, before performing the τ F → 0 limit in Eq. ( 5), we must first take the N τ → ∞ limit to find G B (τ, T, τ F ).Because the available τ values also depend on the lattice spacing, we first perform a τ -interpolation of G B (τ, T, N τ , τ F )/G norm on the two coarsest lattices, separately at each (T, τ F ) pair, to establish the value at the τ values available on the finest lattice.This then allows the N τ → ∞ extrapolation, which is performed independently at each (τ, T, τ F ) triple.In extrapolating to N τ → ∞ we assume that the discretization errors scale as 1/N 2 τ = (aT ) 2 .For each τ T , to take the τ F → 0 limit we multiply the continuum-extrapolated results for G B /G norm by Z match and do linear extrapolations in τ F , as suggested by NLO pQCD [30].To avoid potentially large discretization effects and to keep the gradient-flow scale smaller than all relevant physical scales in the problem, we restrict to the range of flow times 0.25 ≤ √ 8τ F /τ ≤ 0.30 [7].At T = 352 MeV we also perform such flow extrapolation based on the single lattice.
As an example, in Fig. 1 (left panel) we show the resulting G phys.

B
for different choices of μτF and μT .In B can be found in Supplemental Material [13].
G phys.

B
are fitted to Eq. ( 4) to obtain κ B .κ B is encoded in the infrared region of the spectral function, ρ ir B (ω) = ωκ B /(2T ).For the ultraviolet region of spectral function, ρ uv B (ω, µ), we use leading order and NLO vacuum pQCD results.To convert the NLO pQCD results [23,31] from the MS scheme, at the scale µ, to the renormalization-group invariant physical scheme we use the matching NLO Wilson coefficient [32], c 2 B (µ, μT ), to the static quark effective theory.To account for possible higher order effects we also introduce another (in addition to κ B ) ω-independent fit parameter, K.The resulting choices for the UV part of the physical spectral function are ρ uv,phys.
Lattice QCD (LQCD) results for the spatial diffusion coefficients, Ds, for the charm, bottom and infinitely-heavy quarks compared with those from the quasi-particle model (QPM) [24] and the T-matrix approach [25,26].Also shown are the infinitely-heavy quark diffusion coefficients from the ALICE collaboration's phenomenological estimate [27,28], NLO perturbative calculation [29] and AdS/CFT estimate at a certain value of λ [8].The width of the NLO perturbative QCD band corresponds to the variation of the renormalization scale from µ = 2πT (upper boundary of the orange band) and µ = 4πT (lower boundary).
over between the two regimes.(2) The smooth maximum

B
(ω, µ) and µ, and (3) three interpolating models ρ max , ρ smax , and ρ plaw , we carry out 24 different fits for K and κ B on each bootstrap sample of gauge configurations at each T .The final result for κ B at each T is obtained from the median and 68% confidence limit of the distribution of all the bootstrap samples over gauge configurations and fit forms and, thus, include both statistical as well as systematic errors arising from different model and scales choices.−2.41 T 3 , which are of similar magnitude to κ E obtained for the same ensembles [7].The κ B for 2+1 flavor QCD turns out be much larger than those from quenched QCD [20,31] at the same values of T /T c .For further details on the fits and results see Supplemental Material [13].
The ⟨v 2 ⟩ appearing in Eq. ( 2) can be obtained either from the low-frequency part of the spectral function corresponding to the net-flavor current [9,36] or in the quasi-particle model with a temperature dependent quark mass [37].In quenched QCD it has been shown that a quasi-particle model with a temperature dependent heavy quark mass fitted to the heavy quark number susceptibility gives a ⟨v 2 ⟩ that agrees with the one obtained from the low-frequency part of the net heavy quark current spectral function [37].Therefore, in this work we adopt the quasi-particle model to calculate ⟨v 2 ⟩.For the temperature dependent charm quark is obtained from the continuum-extrapolated lattice QCD results for the charm susceptibility [38].No lattice QCD results for bottom quark susceptibility are available presently.Therefore, we simply fix the effective bottom quark mass to 4.8 GeV.The ⟨p 2 ⟩ needed to obtain D s , cf.Eq. ( 1), are estimated from the quasi-particle model in the same way as for the ⟨v 2 ⟩.For further details, see Supplemental Material [13].There we also show that the ⟨v 2 ⟩, ⟨p 2 ⟩ and the values of D s are not too sensitive to the precise choice of the bottom quark mass.
Final results for D s are summarized in Fig. 2. We find a slight increase in D s with decreasing heavy quark mass.We compare our results with those obtained from a phenomenological quasi-particle model (QMP) [24] and the T-matrix approach [25,26].We find that D s in QCD is smaller than the results of these calculations and shows smaller dependence on the heavy quark mass, with the exception of the T-matrix result at the highest temperature considered by us.Our result for D s is also smaller than other phenomenological estimates [27,28], which do not take into account the quark mass dependence.Finally, for completeness, we show the AdS/CFT estimate [8] of D s with a certain value of λ and the result of the NLO perturbative calculation [29] in the limit M → ∞.
Conclusion.-We presented the 2+1 flavor lattice QCD calculations of the quark mass dependence of the momentum and spatial heavy quark diffusion coefficients.In the temperature range 195 MeV ≤ T ≤ 352 MeV the quark mass dependence turns out to be quite small.In conjunction with the previous results for an infinitelyheavy quark [7], the present calculations provide the first non-perturbative QCD results for the charm and bottom quark diffusion coefficients in the QGP.These non-perturbative QCD results will serve as critical inputs to and benchmarks for various dynamical models to study thermalization of charm and bottom quarks in the strongly-coupled medium created in heavy-ion collisions at the Large Hadron Collider and the Relativistic Heavy-Ion Collider.
All data from our calculations, presented in the figures of this paper, can be found in [42].

Supplemental Materials
In these supporting materials we provide the details mentioned in the Letter for the calculation of the finite mass correction to the heavy quark momentum diffusion coefficient via B-field correlators.We start with Sec.I by giving additional details to the lattice setup.In Sec.II we describe the calculation of the matching factor Z match .Sec. III is devoted to the discussions on the continuum and flow time extrapolations of the B-field correlators.In Sec.IV we detail the modeling of the spectral function.The obtained κ B is compared to other estimates and κ E in Sec.V. We close by checking the dependence of 2πT D s on the thermal charm and bottom quark mass in Sec.VI.The gauge configurations in this study were generated using a Rational Hybrid Monte Carlo (RHMC) algorithm [43,44].The configurations are saved after every 10 trajectories with acceptance rate ≈ 80%.The dynamical strange quark mass m s is at its physical value and the light quarks (up/down) are carrying mass m l = m s /5, corresponding to the pion mass of 320 MeV.The masses are obtained from the parametrization of the lines of constant physics from [45].Here we consider four temperatures T =195, 220, 251 and 293 MeV.At each temperature three lattices are used to perform the continuum extrapolation.The finest lattices are of size 96 3 × N τ , where N τ = 24, 28 and 36.All the finest lattices have the same β = 8.249, corresponding to 1/a = 7.036 GeV determined using the r 1 -scale [46], with r 1 = 0.3106 fm taken from [47].The two coarser lattices are of size 64 3 × N τ , where N τ =20, 22/24, generated with a β value that gives the desired temperature, again determined through the r 1 scale.The lattice parameters of the ensembles used in this calculation, including the number of gauge configurations are summarized in Tab.I.
To account for the possible auto-correlation residing in the correlators we estimate the auto-correlation time by analyzing the Polyakov loop at the maximum flow time allowed in this study √ 8τ F T = 0.15.The Polyakov loop, as part of the correlator, exhibits stronger auto-correlation than the correlator itself.Thus its auto-correlation can serve as a conservative estimate for the auto-correlation of the correlator.The auto-correlation also becomes stronger at larger flow time.Therefore our strategy of calculating the auto-correlation is robust.The obtained integrated auto-correlation time is ∼ 10 − 30 in terms of the number of configurations.Then we group the configurations of that size to make independent bins that are passed to the bootstrap analysis afterwards.
When presenting the temperature dependence of the obtained heavy quark momentum diffusion coefficient, it is conventional in the literature to express it in terms of T /T c .T c in the quenched case means the confinement/deconfinement transition temperature.In full QCD it makes more sense to associate T c with the chiral crossover temperature, which can be determined by locating the peak location of the disconnected chiral susceptibility.A modelling of the peak location of the disconnected chiral susceptibility suggests that T c = 180(2) MeV [48] for the quark masses used in this study.

II. DETERMINATION OF THE MATCHING FACTOR
The matching factor from the B-field correlator in the gradient flow scheme to the physical correlator is obtained using a three-step matching procedure.We first match the gradient flow scheme to the MS scheme at some UV scale μτF [23].Then evolve the B-field operators from scale μτF to some thermal scale μT , and finally match to the physical correlator using the result of Ref. [32].The matching factor reads [23] ln where γ 0 = 3/8π 2 is the leading-order anomalous dimension of a B-field.In the above formula we use MS coupling constant g 2 MS that is calculated using the RunDec package [49,50] with N f = 3, Λ MS =0.339 GeV [51] at 5-loop.This is justified because at this order the difference between the MS coupling constant and gradient flow coupling constant can be neglected.Different choices of the scales lead to four different Z match (μ T , μτF , µ F ), which are shown as a function of √ 8τ F /r 0 in Fig. 3. r 0 = 0.469 fm is taken from [52].Changing μT at fixed μτF amounts to an overall multiplicative shift, while changing μτF changes the shape of the curve.However, the small-τ F limiting behavior is not sensitive to μτF .This indicates that the obtained κ B will show a small difference among different μτF choices.

III. CONTINUUM AND FLOW TIME EXTRAPOLATIONS
As discussed in the main text, for the analysis of the B-field correlator and the continuum extrapolations we need the normalization correlator, The normalization correlator in the continuum at zero flow time reads [9]: Its lattice version at finite flow time is calculated in the same way as for the chromo-electric correlator in [7] with the difference that the observable matrix M AB µν here reads [53] M AB µν (p, τ ) = After tree level improvement and normalization, which can be done effectively by dividing the the unmatched correlator by the lattice version of the normalized correlator at finite flow time, our results are shown in Fig. 4 at two selected flow times √ 8τ F /τ = 0.25, 0.3 for different temperatures.We show the results for the finest and for the coarsest lattice.We see from the figure that the discretization effects are larger than for the chromo-electric correlator [7].
To perform the continuum extrapolations the correlators at different lattice spacings are first interpolated using cubic spline to the same separations τ T that are present on the finest lattice.In this way the lattice spacing effects can be seen for each τ T .In Fig. 5 they are shown for some selected τ T s at T = 195, 293 MeV.One is the lowest and the other is the highest temperature available in this calculation.As an example we pick a fixed relative flow time √ 8τ F /τ = 0.3.As in Ref. [7] the interpolated correlators are extrapolated to the continuum limit by fitting to the Ansatz f in a combined way, assuming that the slope decreases with increasing τ T [7].Again, it can be seen that, the lattice spacing effects are very minor, similar as in the case of the E-field correlators [7].
The obtained continuum correlators multiplied by the matching factors are independently extrapolated to the τ F → 0 limit linearly in the window 0.25 ≤ √ 8τ F /τ ≤ 0.3.Some of the instances are shown in Fig. 6 for T = 195 MeV with different μT /T and μτF /µ F .We see from the figure that the change of μτF changes significantly the slope of the zero flow time extrapolation, while changing μT /T effectively amounts to a multiplicative factor in the extrapolated value.In Fig. 6 we also show the zero flow time extrapolation of the normalized B-field correlator at T = 293 MeV.We see that for the same choices of μT /T and μτF /µ F the slope of the zero flow time extrapolation is somewhat smaller at high temperature.

IV. SPECTRAL RECONSTRUCTION
As discussed in the main text, to obtain κ B from the B-field correlators we need to model the spectral function.At high energy the spectral functions can be reliably calculated in perturbation theory.Therefore, we could use the perturbative result to model the high energy part of the spectral function.The spectral function was calculated at NLO in Ref. [31].Evaluating the NLO contribution numerically we found that the thermal effects are quite small, except at very small ω/T values, where perturbation theory is not valid.Therefore, we use the zero temperature limit of the perturbative spectral function in MS scheme at large ω: Here N c = 3, N f = 3 and C F = (N 2 c − 1)/(2N c ).We note that Eq. ( 10) contains −2π 2 /3 rather than −8π 2 /3 as appears in Ref. [31] due to a mistake in the original calculation.To obtain the spectral function for the physical scheme we need to multiply the above expressions by c B (µ, μT ) for which we use the 1-loop renormaliziation group inspired resummed expression After multiplying ρ uv,NLO B by c 2 B (µ, μT ) the MS scale choice µ cancels at leading order, but there will be a residual µ-dependence from higher order contributions.So in practice the choice of µ has an effect.We consider two choices of µ.The first choice is µ 2 = 4ω 2 + µ 2 DR , where µ DR is the thermal scale of dimensionally reduced effective field theory [35].Our second choice is µ 2 = (0.13306ω) 2 + µ 2 DR , which minimizes the NLO contribution for very large ω.Both choices correspond to a smooth variation of the scale from the natural thermal scale to the frequency in the UV, thus avoiding large logs for large ω.We also use the above choices of µ for ρ uv,LO B .g 2 MS (µ) is calculated in the same way as it in Eq. (6).To account for the missing higher order effects we multiply the UV part of the spectral function by a constant factor, K.
With these choices and definitions we fit the models to our lattice data.All models contain two fit parameters: κ B /T 3 and K.In this section we show fit results using µ 2 = (0.13306ω) 2 + µ 2 DR , μT = 19.18Tand μτF /µ F = 1 as an example.In Fig. 7 we show the fitted spectral functions for the lowest temperature T = 195 MeV on the left and the highest temperature T = 293 MeV on the right.It can be seen that using LO   function makes only a small difference, much smaller than the overall model dependence of the spectral function.This is because beyond the IR scale LO and NLO perturbative spectral function roughly differs by a multiplicative factor, which can be compensated by the adjustment of the K-factor in the fitting process, see Fig. 8.The other two temperatures are similar.The ratio of the model correlators and the lattice correlator is shown in Fig. 9, again for T = 195 MeV and T = 293 MeV.We can see that the model spectral function describes our lattice data well within errors.The obtained κ B from different models and scales forms a distribution, see Fig. 10.We calculate the weighted average and use it as our final estimate for κ B .

V. COMPARISON OF κE AND κB
In Fig. 11 we compare κ E and κ B obtained in quenched lattice QCD and 2 1-flavor lattice QCD.The quenched κ E are taken from [19,34,54] calculated using gradient flow or multilevel algorithm.The full QCD κ E is taken from Ref. [7].For quenched κ B we select two estimates, one from multilevel algorithm [31] and the other from gradient flow at finite τ F [20].We present the figures in terms of both T /T c and T .One can see that κ E and κ B are of similar size, in both quenched limit and full QCD, when shown in T /T c .The full QCD results are much larger at small T /T c .When presented in T , both κ E /T 3 and κ B /T 3 follow a decreasing pattern, connecting smoothly to the quenched results at high temperature.To obtain κ and D s we need to evaluate the averaged velocity squared, ⟨v 2 ⟩, and the averaged momentum squared, ⟨p 2 ⟩ of the heavy quark.This can be done assuming that the heavy quarks are well defined quasi particles in the hot medium with a potentially temperature dependent mass, M [37].Thus we need to define the masses of the charm and bottom quarks.
To obtain the charm quark mass, we fit the continuum-extrapolated lattice QCD results of the charm quark number susceptibility [38] to the quasi particle model [37] where E 2 p = M 2 + p 2 .Here we use Boltzmann approximation because the heavy quark masses are significantly larger than the temperature.The lattice data of the susceptibility and the fit results are shown in the left panel of Fig. 12 and the extracted charm quark mass is shown in the right panel.It can be seen that the quasi particle model describes the lattice data rather well.We also see that the temperature variation of the charm quark mass between T = 190 MeV and T = 300 MeV is 309 MeV.Then the mass is interpolated to the desired temperatures, listed in Table II.In the table we indicate the uncertainty of the charm masses due to the uncertainties in χ charm mass, ⟨v 2 ⟩ and ⟨p 2 ⟩ can be calculated again in quasi-particle model [37]: In  to the uncertainty M c are small.
We consider two values 4.5 GeV and 4.8 GeV for the bottom quark.The difference between these two values of the bottom quark mass should cover possible thermal effects, see above.In Table II we also show ⟨v 2 ⟩, ⟨p 2 ⟩/(3M T ) and D s for bottom quarks.Again the effects of the choice of M b are small.[19,34,54] and the 2+1 flavor κE is taken from [7].The quenched κB are taken from [20,31].The κE points are slightly shift horizontally for better visibility.12. Left: lattice results of the charm quark susceptibility [38] and the fit to quasi particle model.Right: the extracted charm quark mass from the fits.

Fig. 1 (
right panel) we show the T dependence of G phys.B for μT /T = 19.18 and μτF /µ F = 1.0.Further details on and examples of G phys.

μτFFIG. 3 .
FIG. 3. Z match (μ T , μτ F , µF) calculated for different scales at T = 195 MeV as an example.The grey band covers the flow time range valid for the τF → 0 extrapolation, see Sec.III.

g 2 0 20 FIG. 4 .
FIG. 4. The B-field correlators in the gradient flow scheme for different temperatures calculated on the finest (open symbols) and coarsest lattices (filled symbols) at relative flow time √ 8τF/τ = 0.25 (left) and 0.30 (right).The lattice results have been normalized by G norm obtained at non-zero flow time and lattice spacing.

FIG. 5 .
FIG. 5.The lattice spacing dependence of the B-field correlator at T = 195 MeV (left) and T = 293 MeV (right) at √ 8τF/τ = 0.3.The dashed lines show the central values of the continuum extrapolations.

FIG. 9 .
FIG. 9.The ratio of the fit correlators to the lattice correlators at T = 195 MeV (left) and = 293 MeV (right).
FIG.12.Left: lattice results of the charm quark susceptibility[38] and the fit to quasi particle model.Right: the extracted charm quark mass from the fits.

TABLE I .
The temperatures used in this calculation and the β values, bare quark masses, lattice sizes and the number of configurations at each temperature.
The fitted spectral functions using different models for T = 195 MeV (left) and T = 293 MeV (right).
Table II we show ⟨v 2 ⟩, ⟨p 2 ⟩/(3M T ) and D s for charm quarks.Note that the uncertainties in these quantities due

TABLE II .
The dependence of ⟨v 2 ⟩, ⟨p 2 ⟩/(3M T ) and 2πT Ds on the quark mass.The four columns in the middle are for the bottom quark.The uncertainties of 2πT Ds are from κ solely.The first row at each temperature is for m b =4.5 GeV while the second row for 4.8 GeV.The four columns on the right are for the charm quark.Note that the uncertainties in the first brackets of the last column inherit from κ and the second brackets from the uncertainties in the charm quark mass.