Strong coupling constant and heavy quark masses in 2+1 flavor QCD

We present a determination of the strong coupling constant and heavy quark masses in 2+1-flavor QCD using lattice calculations of the moments of the pseudo-scalar quarkonium correlators at several values of the heavy valence quark mass with Highly Improved Staggered Quark (HISQ) action. We determine the strong coupling constant in $\overline{MS}$ scheme at four low energy scales corresponding to $m_c$, $1.5m_c$, $2m_c$ and $3m_c$, with $m_c$ being the charm quark mass and from these we obtain $\Lambda_{\overline{MS}}^{n_f=3}=301 \pm 16$ MeV, which is equivalent to $\alpha_s(\mu=M_Z,n_f=5)=0.1161(12)$. For the charm and bottom quark masses in $\overline{MS}$ scheme we obtain: $m_c(\mu=m_c,n_f=4)=1.2672(84)$ GeV and $m_b(\mu=m_b,n_f=5)=4.188(29)$ GeV.


I. INTRODUCTION
In recent years there was an extensive effort toward accurate determination of QCD parameters since the precise knowledge of these parameters is important for testing the predictions of the Standard Model. Two important examples are the sensitivity of Higgs branching ratios to the heavy quark masses and the strong coupling constant [1,2] and the stability of the Standard Model vacuum [3,4]. Lattice calculations play an increasingly important role in the determination of the QCD parameters as these calculations become more and more precise with advances in computational approaches.
The strong coupling constant α s has been known for a long time. The current Particle Data Group (PDG) values α s (M Z ) = 0.1181 (11) [5] has small errors suggesting that the uncertainties are well under control. However, there are several phenomenological determinations of α s that give much smaller central values and have small errors [6]. Therefore, the errors of α s determinations may not be fully understood. The lattice average of α s provided by Flavor Lattice Averaging Group (FLAG) agrees well with the PDG average [7]. However, there are several lattice determinations reporting smaller values [8][9][10]. The running of the strong coupling constant at lower energy scales is also interesting. For example, for testing the weak coupling approach to QCD thermodynamics through comparison to lattice QCD results [11][12][13][14][15][16][17] one needs to know the coupling constant at a relatively low energy scale ∼ πT , with T being the temperature. The analysis of the τ decay offers the possibility of α s extraction at a low energy scale but there are large systematic uncertainties due to different ways of organizing the perturbative expansion in this method (see Refs. [18][19][20] for a recent work on this topic and references therein). Lattice QCD calculations on the other hand are well suited to map out the running of α s at low energies.
Lattice determination of the charm quark mass significantly improved over the years. Though, the ETMC collaborations reported much larger value for the charm quark mass [21]. Some lattice QCD calculations use 2 or 3 flavors of dynamical quarks [22][23][24], while others use 4 dynamical flavors [21,25,26]. Therefore, additional calculations of the charm quark mass could be useful.
These calculations can also help to better understand the dependence of the charm quark mass on the number of active quark flavors.
Determination of the bottom quark mass is difficult in the lattice simulations due to the large discretization errors caused by powers of m h a, where m h is the bare mass of the heavy quarks. One needs a small lattice spacing to control the corresponding discretization errors. One possibility to deal with this problem is to perform calculations with heavy quark masses smaller than the bottom quark mass and extrapolate to the bottom quark mass guided by heavy quark effective theory [27,28]. Several determinations of the bottom to charm quark mass ratio have been reported and slight inconsistency has been found. The ratio recently obtained by the ETMC collaboration [29] shows smaller value than that previously determined by the HPQCD collaboration [25,27] as well as the recent determination of the Fermilab, MILC and TUMQCD collaborations [28].
In this paper we report on the calculations of α s and the heavy quark masses in 2+1 flavor QCD using the Highly Improved Staggered Quark (HISQ) action and moments of pseudo-scalar quarkonium correlators. We extend the previous work reported in Ref. [10]. This paper is organized as follows. In section II we introduce the details of the lattice setup. In section III we discuss the moments of quarkonium correlators and our main numerical results. The extracted values of the strong coupling constant and heavy quark masses are discussed in section IV and compared to other lattice and non-lattice determinations. The paper is concluded in Sec. V. Some technical details of the calculations are given in the Appendices.

II. LATTICE SETUP AND DETAILS OF ANALYSIS
To determine the heavy quark masses and the strong coupling constant we calculate the pseudoscalar quarkonium correlators in (2 + 1)-flavor lattice QCD. As in our previous study we take advantage of the gauge configurations generated using the tree level improved gauge action [30] and the Highly Improved Staggered Quark (HISQ) action [31] by HotQCD collaboration [32]. The strange quark mass, m s , was fixed to its physical value, while for the light (u and d) quark masses the value m l = m s /20 was used. The latter corresponds to the pion mass m π = 161 MeV in the continuum limit, i.e. the sea quark masses are very close to the the physical value. This is the same set of gauge configurations as used in Ref. [10]. We use additional HISQ gauge configurations with light sea quark masses m l = m s /5, corresponding to the pion mass of 322 MeV at five lattice spacings corresponding to lattice gauge coupling β = 10/g 2 0 = 7.03, 7.825, 8.0, 8.2 and 8.4, generated for the study of the QCD equation of state at high temperatures [12]. This allows to perform calculations at three smaller lattice spacings, namely a = 0.035, a = 0.03 and 0.025 fm, and, also check for sensitivity of the results to the light sea quark masses. As we will see later, the larger than the physical sea quark mass has no effect on the moments of quarkonium correlators.
For the valence charm and bottom quarks we use the HISQ action with the so-called ǫ-term [31], which removes the tree-level discretization effects due to the large quark mass up to O((am) 4 ). The HISQ action with ǫ-term turned out to be very effective for treating the charm quark on the lattice [26, 31,33,34]. The lattice spacing in our calculations has been fixed using r 1 scale defined in terms of the energy of static quark anti-quark pair V(r) as We use the value of r 1 determined in Ref. [35] using the pion decay constant as an input: In the above equation all the sources of errors in Ref. [35] have been added in quadrature. The above value of r 1 corresponds to the value of the scale parameter determined from the Wilson flow w 0 = 0.1749 (14) fm [32]. This agrees very well with determination of Wilson flow parameter by BMW collaboration w 0 = 0.1755(18)(4) fm [36]. It is also consistent with the value r 1 = 0.3133(23)(3) fm reported by HPQCD collaboration within errors [37]. It turns out that the value of r 1 does not change within errors when increasing the sea quark mass from m s /20 to m s /5 [12].
Therefore, the same r 1 value can be used to set the scale for m l = m s /5 calculations.
We calculate pseudo-scalar meson correlators for different heavy quark masses using random color wall sources [25]. This reduces the statistical errors in our analysis by an order of magnitude compared the previous study [10]. Here we consider several quark masses in the region between the charm and bottom quark, namely m h = m c , 1.5m c , 2m c , 3m c , 4m c and m b . This helps us to study the running of α s at low energies and provides additional cross-checks on the error analysis.
The bare quark mass m c0 that corresponds to the physical charm quark mass has been determined in Ref. [10]. We use the same values of m c0 in this work. For the three largest values of β we determined the bare charm quark mass by requiring that the mass of η c meson obtained on the lattice in physical units agrees with the corresponding PDG value. We determine the b quark mass at each value of gauge coupling β by performing calculations at several values of the heavy quark mass near the b quark mass and linearly interpolating to find the quark mass at which the pseudo-scalar mass is equal to the physical mass of η b from PDG. In these calculations we use both random color wall sources and corner wall sources. In Appendix A we summarize the gauge ensembles used in our study, the number of measurements for each β as well as the values of the bare charm and bottom quark masses.

III. MOMENTS OF QUARKONIUM CORRELATORS AND QCD PARAMETERS
We consider the moments of the pseudo-scalar quarkonium correlator, which are defined as Here j 5 =ψγ 5 ψ is the pseudo-scalar current and m h0 is the lattice heavy quark mass. To take into account the periodicity of the lattice of temporal size N t the above definition of the moments can be generalized as follows: The moments G n are finite for n ≥ 4 (n even) since the correlation function diverges as t −4 for small t. Furthermore, the moments G n do not need renormalization because the explicit factors of the quark mass are included in their definition [38]. They can be calculated in perturbation theory in MS scheme Here µ is the MS renormalization scale. The scale at which the MS heavy quark mass is defined can be different from µ [39], though most studies assume µ m = µ. The coefficient g n (α s (µ), µ/m h ) is calculated up to 4-loop, i.e. up to order α 3 s [40][41][42]. Given the lattice data on G n one can extract α s (µ) and m c (µ) from the above equation. However, as it was pointed out in Ref. [38] it is more practical to consider the reduced moments where G (0) n is the moment calculated from the free correlation function. The lattice artifacts largely cancel out in these reduced moments. Our numerical results on R n and some of the relevant ratios, e.g. R 8 /R 10 , R 8 /m c0 are given in Appendix A. From the tables one can clearly see that the statistical errors are tiny and can be neglected. The dominant errors in our calculations are the errors due to finite size of the lattices and the errors induced by mistuning of the heavy quark mass. We estimated the finite size errors in the free theory and included them in the analysis. The free theory estimate provides an upper bound on the finite size errors. It is known that the finite size effects in the interacting theory are much smaller [25,27]. To estimate the errors due to the mistuning of the heavy quark mass we performed interpolation of R n in the heavy quark mass and examined the changes in R n when the value of the heavy quark mass is varied by one sigma. These systematic errors are summarized in Appendix A.
It is straightforward to write down the perturbative expansion for R n : From the above equations it is clear that R 4 as well as the ratios R 6 /R 8  of m h needs to be taken into account, which increases the uncertainty of the perturbative result [39].
There are also non-perturbative contribution to the moments proportional to gluon condensate [43].
We included this contribution at tree level using the value α s π G 2 = 0.006 ± 0.012 (10) from the analysis of τ decay [44].
In order to extract α s and the heavy quark masses from R n continuum extrapolation needs to be performed. Since tree-level lattice artifacts cancel out in the reduced moments R n we expect that discretization errors are proportional to α n s (am h0 ) j . Therefore, we fitted the lattice spacing dependence of R 4 , R n /m h0 , n ≥ 6 and of the ratios R 6 /R 8 and R 8 /R 10 with the form Here for α s we use the boosted lattice coupling defined as where g 2 0 = 10/β is a bare lattice gauge coupling and u 0 is an averaged link valuable defined by the plaquette u 4 0 = TrU /3. We performed fits using up to N = 3 and J = 5 in Eq. (11). We use data corresponding to am h0 < 1.1 to avoid uncontrolled cutoff effects. We note that the radius of convergence of the series expansion in am h0 for the free theory is π/2 [45], and thus, our upper limit on am h0 is well within the radius of convergence of the expansion. The coefficients of higher order terms in α s are not well constrained by the data and therefore in many cases we assumed that they are equal to a certain fraction of the lowest order coefficient. The fit range in (am h0 ) was also varied. We find that for larger values of m h and larger range of lattice spacings more terms in Eq. (11) had to be included, as expected. The lattice spacing (cutoff) dependence of R 4 turned out to be the most complicated. This is not completely surprising as R 4 being the lowest moment is most sensitive to short distance physics. Here we had to use up to five order polynomial in am h0 and at least two powers of α s to describe the lattice data on R 4 . Simpler fit forms only worked for the lowest mass and very small value of the lattice spacings. For the ratios R 6 /R 8 and R 8 /R 10 we also had to use high order polynomials in (am h0 ) 2 though the leading order in α s turned out to be sufficient. On the other hand the lattice spacing dependence of the ratios R n /m h0 are described by the leading order (α s a 2 ) or leading order plus next-to-leading order (α s (a 2 + a 4 )) forms even for large values of m h . To demonstrate these features we show sample continuum extrapolations for the moments in Fig. 2  increases with decreasing a 2 . Therefore, if there are no data points at small a the continuum limit may be underestimated. The leading order fit only works for the four smallest lattice spacings but agrees with the fit that uses fifth order polynomial in (am h0 ) 2 with N = 2, and extends to the whole range of the lattice spacings. Thus, the additional three lattice spacings included in this study are important for cross-checking the validity of continuum extrapolation, although the correct continuum result for R 4 can be obtained without these additional data points. For a-dependence of the ratio R 8 /R 10 we see the opposite trend, the slope decreases at small a. Not having lattice results at small lattice spacing may lead to an overestimated continuum result. This difference in the two quantities provides an important cross-check for the determination of α s . In the right panel of Fig. 2 we show the fits using fourth (solid line) and third order (dashed line) polynomials in am h0 and leading order in α s . We see that the two fits give very similar result and we find that higher terms in α s do no have a big impact here. Because high order polynomials are needed for extrapolations, continuum results for R 4 , R 6 /R 8 and R 8 /R 10 could only be obtained for m h ≤ 3m c .
For larger values of the quark masses we simply do not have enough data satisfying am h0 < 1.1 to perform the continuum extrapolations. The lattice spacing dependence of R n /m h0 , n ≥ 6 is well described by the next-to-leading order form for all quark masses as can be seen for example in the right panel of Fig. 3. For R 6 /m c0 even leading a 2 dependence is sufficient to describe the data, cf. left panel of Fig. 3. Including higher order terms in α s has no effect in this case. For  (8) 1.5m c 1.228(2) 1.0895(11) 1.0403 (10) 2.0m c 1.194(2) 1.0791 (7) 1.0353 (5) 3.0m c 1.158(6) 1.0693 (10) (13) 2.0m c 0.5584(35) 0.5156 (17) 0.4972(17) 3.0m c 0.3916(23) 0.3647 (19) 0.3527(20) 4.0m c 0.3055(23) 0.2859 (12)   results for R 4 , R 6 /R 8 and R 8 /R 10 are shown in Table II. In Table III we give our continuum results for R n /m h0 for n ≥ 6. These two tables present the main result of this study.
As will be discussed in the following section using the continuum results on the reduced moments R n and their ratios presented in Tables II and III one can obtain the strong coupling constant as well as the values of the heavy quark masses, and may perform many important cross-checks.
However, before discussing the determination of α s and the quark masses let us compare our continuum results for the moments and their ratios with other lattice determinations. In Fig. 4 Figure 5 we compare our continuum results on R 6 , R 8  in R n we assume that the coefficient of the unknown α 4 s term varies between −5r n3 and +5r n3 , i.e. r n4 = ±5 × r n3 . This error estimate should be sufficiently conservative. We also take into account the non-perturbative contribution to the reduced moments at tree level according to Ref. [43] and the value of the gluon condensate given in Eq. (10). In Table IV we give our results for α s (µ = m h ) for different quark masses. We see that both the perturbative uncertainty as well as the uncertainty due to the gluon condensate drastically decrease with increasing m h . The α s values determined from R 6 /R 8 and R 8 /R 10 have much larger perturbative uncertainties than the ones from R 4 . As one can see from the table the strong coupling constants obtained from R 4 , R 6 /R 8 and R 8 /R 10 do not always agree with each other. Only for m h = 2m c we find that all three determination agree within errors. We usually find that α s determined from R 8 /R 10 is systematically low. To obtain our final estimate of α s (µ) for µ = m c , 1.5m c , 2m c and 3m c we performed weighted average of the results obtained from R 4 , R 6 /R 8 and R 8 /R 10 . These are given in fifth column of Table IV. The uncertainty of this averaged α s values was determined such that it agrees with all individual α s extraction within the estimated errors.
Using the continuum results for R n /m h0 , n ≥ 6, given in Table III together with the corresponding perturbative expression for R n , n ≥ 6, and the averaged value of α s (m h ) given in the fifth column of Table IV we obtain the values of m h in MS scheme at µ = m h . These are presented in Table V. The heavy quark masses obtained from R 6 , R 8 and R 10 agree with each other very well.
Therefore, we calculated the corresponding weighted average to obtain our final estimates and the errors for the heavy quark masses. These are given in the last column of  MS from the value of the coupling at µ = m h we use the implicit scheme [47]. We also calculated the Λ n f =3 MS in explicit scheme and the small differences between the two schemes have been treated as systematic errors.  Table IV to a constant the χ 2 /d f is larger than two. One way to obtain a χ 2 /d f of about one for this fit is to assume that the error of α s (2m c ) has been underestimated by a factor two. However, it is not obvious how to justify such an assumption. Therefore, to obtain our final estimate for Λ n f =3 MS we take a weighted average of the data in the last column of Table IV and use the spread around this central value as our (systematic) error: The too low value of Λ n f =3 MS obtained for m h = 2m c is of some concern. One could imagine that the continuum extrapolation of R 4 at this quark mass is not reliable and the corresponding α s should not be considered in the analysis. If we determine α s (2m c ) using only the results for R 6 /R 8 and R 8 /R 10 we obtain a value for α s (2m c ), which is only one tenth sigma smaller than the value in Table IV. Finally if we take α s (2m c ) only from R 6 /R 8 we obtain Λ  Table III. RunDeC package [47] and first match α s at the charm threshold which we choose to be 1.5 GeV and then at the bottom quark threshold, which we choose to be 4.7 GeV and finally evolve α s to µ = M Z . We get α s (M Z , n f = 5) = 0.1161 (12) .
Here we also included the uncertainty in Λ

V. CONCLUSION
In this paper we calculated the moments of quarkonium correlators for several heavy quark masses in 2+1 flavor QCD using HISQ action. From the moments of quarkonium correlators we extracted the strong coupling constant and the heavy quark masses. Our main results are given in Tables II and III and by Eqs. (13), (14), (15) and (17). We improved and extended the previous 2+1 flavor HISQ analysis published in Ref. [10]. We drastically reduced the statistical errors on the moments by using random color wall sources, extended the calculations to smaller lattice spacings and considered several values of the heavy quark masses in the region between the charm and bottom quark mass. The calculations of the reduced moments at several heavy quark masses enabled us to map out the running of the coupling constant at low energy scales and better control the systematic errors due to the truncation of the perturbative series. Having many lattice spacings enabled us to study the cutoff dependence of the moments in detail and control the errors due continuum extrapolation. Furthermore, to have a handle on the uncertainties due to continuum extrapolations the strong coupling constant and the quark masses have been extracted using different moments and/or their ratios. The heavy quark masses in MS scheme obtained from various moments and different valence quark masses are remarkably consistent with each other. On the other hand there is some tension in the determination of the strong coupling constant performed using different quantities and different valence heavy quark masses. To deal with this tension the errors on the Λ-parameter have been chosen accordingly.
Evolving the low energy determination of α s to µ = M Z we obtain the value α s (M Z , n f = 5) = 0.1161 (12), which agrees with the previous result [10] but has larger error. Our result for α s is lower than many lattice QCD determinations as well as the PDG average. On the other hand it agrees with the determination of α s static QQ energy [9].
From the sixth, eighth and tenth moments we determined the charm and bottom quark masses.
Our results on the heavy quark masses agree well with the previous 2+1 flavor HISQ determination [10] but have smaller errors. We also found that our results agree well with other lattice determinations that are based on various approaches. Thus, our analysis suggests that lattice determination of the heavy quark masses is under control. [25] B. Chakraborty, C. T. H. Davies [32], while the lower part to lattices from Ref. [12].
Tables VII-XVIII the results for the same β but different light sea quark masses agree within the statistical errors. Thus, the use of heavier than the physical light sea quark masses has no effect on our results.

Appendix B: Details of continuum extrapolations
In this Appendix we discuss further details of the continuum extrapolations. As mentioned in the main text we performed a variety of continuum extrapolations using Eq. (11) and keeping different number of terms in the sum. In doing so we varied the fit interval such the χ 2 /d f is close to or below one. Fewer terms in Eq. (11) usually means more restricted interval in am h0 .
We first discuss the continuum extrapolation of the fourth reduced moment, R 4 , which is the most challenging. Here we find that including terms up to N = 2 in Eq. (11)  11) could be treated as free fit parameter as we have many data points for relatively small am h0 to have stable fits. We also consider fits where c 2 j /c 1 j was fixed to some value between −4 and −5.5 but the continuum extrapolated R 4 value did not change much. For m h = 3m c the value of the coefficient of α 2 s term was varied between −4 and −5.5, but setting the coefficient to zero did not change the result much. In Fig. 7 we show continuum estimates for R 4 obtained from different fits. The central values of these estimates show some scattering, although most of them agree within errors. We take the weighted average of these estimates to obtain our final continuum result for R 4 , which is shown as the solid horizontal line in Fig. 7. We then assign an error to this continuum result. The size of the error band is defined in such way that all individual fits agree with the final continuum estimate or it represents the typical error of the fits, whichever is the largest. The errors on the final continuum estimate are represented by dashed vertical bands in Fig. 7. We estimated the finite volume errors on the reduced moments using the free theory result, see above. As one can see from Fig. 2  In the 4th through 7th columns we list the moments with statistical error, finite size error and error due to the uncertainty of the quark mass.
the three finest lattices. To exclude the possibility that finite volume errors are underestimated we also performed fits omitting the data points corresponding to the three smallest lattice spacings for m h = m c , where finite volume effects are the largest. We find that doing so does not affect the final continuum estimate within errors.
A similar analysis has been performed also for the ratios of the reduced moments, R 6 /R 8 and Here including terms up to N = 2 in Eq. (11) turned out to be less important and good fits could be obtained already with N = 1 in the entire am h0 interval. The corresponding results are shown in Fig. 8 and Fig. 9.