{'fraction_non_alphanumeric': 0.0442401855510695, 'fraction_numerical': 0.01803968731208659, 'mean_word_length': 4.07244008714597, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-Precision c and b Masses, and QCD Coupling from Current-Current Correlators in Lattice and Continuum QCD

We extend our earlier lattice-QCD analysis of heavy-quark correlators to smaller lattice spacings and larger masses to obtain new values for the c mass and QCD coupling, and, for the first time, values for the b mass: m_c(3GeV,n_f=4)=0.986(6)GeV, alpha_msb(M_Z,n_f=5)=0.1183(7), and m_b(10GeV,n_f=5)=3.617(25)GeV. These are among the most accurate determinations by any method. We check our results using a nonperturbative determination of the mass ratio m_b(mu,n_f)/m_c(mu,n_f); the two methods agree to within our 1% errors and taken together imply m_b/m_c=4.51(4). We also update our previous analysis of alpha_msb from Wilson loops to account for revised values for r_1 and r_1/a, finding a new value alpha_\msb(M_Z,n_f=5)=0.1184(6); and we update our recent values for light-quark masses from the ratio m_c/m_s. Finally, in the Appendix, we derive a procedure for simplifying and accelerating complicated least-squares fits.


I. INTRODUCTION
Precise values for the QCD coupling α MS and the quark masses are important for high-precision tests of the Standard Model of particle physics. In a recent paper we showed how to use realistic lattice QCD simulations to extract both the coupling and the charm quark's mass m c from zero-momentum moments of correlators built from the c quark's (UV cutoff-independent) pseudoscalar density operator m c ψ c γ 5 ψ c [1]. In this paper we refine our previous analysis and extend it to include other quark masses, up to and including the b-quark mass. As a result our coupling constant and mass determinations from these correlators are among the most accurate by any method.
Low moments of heavy-quark correlators are perturbative and several are now known through O(α 3 s ) in perturbation theory (that is, four-loop order) [2][3][4][5][6]. Moments of correlators built from the electromagnetic currents can be estimated nonperturbatively, using dispersion relations, from experimental data for the electron-positron annihilation cross section, σ(e + e − → γ * → X). Accurate values for both the c and b masses can be obtained by comparing these perturbative and nonperturbative determinations of the moments (for a recent discussion see [7]).
In our earlier paper we showed that heavy-quark correlator moments are easily and accurately computed nonperturbatively using lattice QCD simulations, in place of experimental data, provided: 1) the electromagnetic current is replaced by the pseudoscalar density multiplied by the bare quark mass; 2) the discretization of the quark action has a partially conserved axial vector current (PCAC); and 3) the discretization remains accurate when applied to heavy quarks. In our simulations we use the HISQ discretization of the quark action, which is a highly corrected version of the standard staggeredquark action [8]. It has a chiral symmetry (PCAC) and has been used in a wide variety of accurate simulations involving c quarks [8][9][10][11][12].
Here we show that the HISQ action can be pushed to still higher masses -indeed, very close to the b masson new lattices, from the MILC collaboration [13], with the smallest lattice spacings available today (0.044 fm). Currently most high-precision lattice work on b physics relies upon nonrelativistic effective field theories, like NRQCD [10, 12,14,15]. In this paper we show how to obtain accurate b physics using the fully relativistic HISQ action on these new lattices.
In what follows, we first review how the QCD coupling and quark masses are extracted from heavy-quark correlators, in Section II. Then in Section III we describe our lattice QCD simulations and discuss in detail the chief systematic errors in our simulation results. In Section IV we describe our fitting procedure and the results of our analysis of the heavy-quark correlators. We check our calculation using a different, nonperturbative method to determine m b /m c in Section V. We then, in Section VI, update our previous calculation of the QCD coupling from Wilson loops to compare with our new result from the correlators. We summarize our findings in Section VII and compare our results with work by others. There we also update our recent calculations of the light-quark masses from the c mass. In the Appendix we present a powerful simplification for complicated least-squares fits that can greatly reduce the computing required for fits. We use this technique in dealing with finite-a errors in our analysis.

II. HEAVY-QUARK CORRELATOR MOMENTS
Following our earlier paper [1], we focus on correlators formed from the pseudoscalar density of a heavy quark, j 5 = ψ h γ 5 ψ h : x (am 0h ) 2 0|j 5 (x, t)j 5 (0, 0)|0 (1) where m 0h is the heavy quark's bare mass (from the lattice QCD lagrangian), t is euclidean and periodic with period T , and the sum over spatial positions x sets the total three momentum to zero. In our earlier paper we examined only c quarks; here we will consider a range of masses between the c and b masses. While we have written this formula for use with the lattice regulator, it is important to note that the correlator is UV-finite because we include the factors of am 0h . Consequently lattice and continuum G(t)s are equal when t = 0 up to O((am h ) m ) corrections, which vanish in the continuum limit.
The moments of G(t) are particularly simple to analyze: where, on our periodic lattice, Low moments emphasize small ts and so are perturbative; and moments with n ≥ 4 are UV-cutoff independent. Therefore G n = g n (α MS (µ), µ/m h ) (am h (µ)) n−4 for small n ≥ 4, where m h (µ) is the heavy quark's MS mass at scale µ, and the dimensionless factor g n can be computed using continuum perturbation theory. Again following our previous paper, we introduce reduced moments to suppress both lattice artifacts and tuning errors in the heavy quark's mass [16]: where G (0) n is the moment in lowest-order, weak-coupling perturbation theory, using the lattice regulator, and m η h is the (nonperturbative) mass of the pseudo-Goldstone hh boson. The reduced moments can again be written in terms of continuum quantities: and r n is obtained from g n (Eq. (4)) and its value, g n , in lowest-order continuum perturbation theory: for n ≥ 6.
Our strategy for extracting quark masses and the QCD coupling relies upon lattice simulations to determine nonperturbative values for the R n , using simulation results for am η h /am 0h . We then compare this simulation "data" to the continuum perturbation theory formulas (Eq. (6)). That is, we find values for α MS (µ) and z(µ/m h , m η h ) that make lattice and continuum results agree for small n ≥ 4. The function z(µ/m h , m η h ) can then be combined with experimental results for m ηc and m η b to obtain masses for the c and b quarks: Parameter µ sets the scale for α MS in the perturbative expansions of the r n . An obvious choice for this parameter is µ = m h since the quark mass, together with n, sets the momentum scale in our correlators. As noted in our previous paper, however, perturbation theory is somewhat more convergent if we use larger µs in the cquark case. Consequently here we take µ/m h = 3, which is approximately what we did in our previous paper.
The mass and coupling determinations were done separately in our previous paper. Here we extract them simultaneously, to guarantee consistency between results. Also in our previous paper we considered only heavyquark masses near the c mass. Here we explore a variety of masses ranging from just below the c mass to just below the b mass. This allows us to obtain a value for b-quark's mass.

A. Simulation Results
The gluon-configuration sets we use were created by the MILC collaboration. The relevant simulation parameters are listed in Table I. Given a lattice spacing, the QCD action is specified completely by the values of the bare coupling constant and the bare quark masses. In our analyses we set the u and d quark masses equal; this approximation results in negligible errors (≪ 1%) for the quantities studied in this paper. It is too costly to simulate QCD at the correct value for the u/d mass; we typically use masses that are 2-5 times too large and extrapolate to values that give I: Parameter sets used to generate the gluon configurations analyzed in this paper. The lattice spacing is specified in terms of the static-quark potential parameter r1 = 0.3133 (23) fm; values for r1/a are from [13]. The bare quark masses are for the ASQTAD formalism and u0 is the fourth root of the plaquette. The spatial (L) and temporal (T ) lengths of the lattices are also listed, as are the number of gluon configurations (N cf ) and the number of time sources (Nts) per configuration used in each case. Sets with similar lattice spacings are grouped. the correct mass for the π 0 -meson. We tune the strange quark mass to give the correct mass for the (fictitious) η s meson [10]. The c and b masses are tuned to give correct masses for the η c and η b mesons, respectively. It is convenient in QCD simulations to specify a value for the bare coupling constant and then extract the value of the lattice spacing from the simulation. We set the lattice spacing using MILC results for r 1 /a, computed from the heavy-quark potential, and [10] The MILC configurations include vacuum polarization contributions from only the lightest three quark flavors, using the ASQTAD discretization. Vacuum polarization effects from the heavier c and b quarks are easily incorporated into our final results for quark masses and the QCD coupling using perturbation theory.
We computed heavy-quark correlators (Eq. (1)) using the HISQ discretization [8] for a variety of bare heavyquark masses am 0h on the MILC gluon configurations. Our results for the reduced moments R n with n = 4-18 are given in Table II.
In Table II, we also give masses am η h from the simulations for the pseudo-Goldstone meson made from two heavy quarks. These were computed using singleexponential fits to G(t) for the middle 30% of ts on the lattice for all configurations except the two smallest lattice spacings where we used only 8% of the ts. We have less statistics for the two finest lattice spacings and consequently the fits did not work as well for these. We increased the statistical errors on our results for am η h by factors of 1.4 and 2 for the next-to-finest and finest lattice spacings (sets 8 and 9), respectively, to account for this. The statistical errors here are very small and have only a small impact on our final results. We also verified our results with multi-exponential fits in every case.

B. Systematic Errors
As discussed above, our goal is to find values for α MS (µ) and z(µ/m h , m η h ) (Eq. (7)) that make the theoretical results from perturbative QCD agree, to within statistical and systematic errors, with Monte Carlo simulation "data" for the reduced moments. We simultaneously analyze results for all of our lattice spacings and most of our masses, and for moments with 4 ≤ n ≤ 10. We focus on these particular moments for our final results since their perturbation theory is known to third order.
Systematic errors are larger here than statistical errors, which contribute less than 0.3%. We discuss the most important sources of systematic error in this section.

m h Extrapolations
We need the m η h dependence of the mass-ratio function z(µ/m h = 3, m η h ) in order to extract c and b masses from our simulation (using Eq. (9)). We parameterize this dependence as follows: where the z j s are determined in our fit. This is an expansion in the QCD scale, which we take to be divided by m η h /2, which we use as a proxy for the quark mass. The expansion is adequate for the range of quark masses used in our analysis, where (2Λ/m η h ) 2 ranges approximately between 1/m 2 η b = 0.01 and (1/m ηc ) 2 = 0.1; the singular point m h = 0 is infinitely far away in this parameterization. In our fits we keep terms only through order N z = 4, but, as we discuss later, our results are unchanged by additional terms. On dimensional grounds, we assume a priori that the coefficients are

Finite-Lattice Spacing Errors
Discretization errors are of order (am h ) 2i α s for i ≥ 1. We model these by where: fit function R n (µ, m η h , a, N am ) has the double expansion  (8) 1.5234 (9) 1.3739 (7) 1.3069 (6) 1.2657 (6) 1.2366 (6) 1.2145 (6) 1.1969 (5)  the c (n) ij s are determined in our fit, R cont n is given by Eq. (6), and again we use m η h /2 in place of the quark mass. This expansion allows for finite-a corrections involving (am η h /2) 2 , (aΛ) 2 , and cross terms, with m η h -dependent coefficients. We assume a priori that c (n) which implies smaller a dependence for larger ns. This is expected (and obvious in our simulation data) since the reduced moments become more infrared as n increases.
The exact functional form of the n dependence has little effect on our results, as we show later. In our fits we take N z = 4. While low orders suffice for the 2Λ/m η h expansion, expansion parameter am η h /2 ranges between 0.3 and 1.1, and higher orders are necessary, especially given our tiny statistical errors. We find that our fit results don't converge well unless N am is larger than 10-20. Also we have difficulty getting good fits if we include data with am η h > 1.95 from Table II. The am η h /2 expansion may not converge for these last cases and therefore we exclude such data from our final analysis.
The fit function has many more fit parameters c (n) ij than we have simulation data points when N am is so large. This does not cause problems in (Bayesian) constrained fits since the parameters' priors (Eq. (17)) are included in the fit as extra data [17]. Each parameter has a prior and therefore we always have more data than parameters.
It is, however, very time consuming to fit a function with so many fit parameters. Although it is not essential for our analysis, there is a trick that greatly accelerates this kind of fit. The idea is to fit a modified moment R latt andN am ≪ N am . The modified moment is fit with the much simpler formula (simpler sinceN am ≪ N am ) where R n (. . .) is again given by Eq. (15). To evalu-ateR latt n from Eq. (18), we treat the coefficients c (n) ij with i >N am as new data with means and standard deviations specified by the prior, Eq. (17). Uncertainties coming from the c (n) ij s are combined in quadrature with the statistical error in R latt n to obtain a new error estimate forR latt n (but leaving the central value unchanged). In effect we are increasing the error in the reduced moment to account for high-order (am η h /2) 2i terms omitted from the fit formula Eq. (19). By choosingN am ≪ N am , most of the am η h /2 terms are incorporated intoR latt n (Eq. (18)), where they are inexpensive, and relatively few end up in the fit functionR n (. . .) (Eq. (15)), where they add parameters to the fit and increase its cost. Note that the new errors introduce correlations betweenR latt n s computed with different lattice spacings or quark masses, since the same c (n) ij s are used for all as and m η h s. These correlations are important and need to be preserved in the fit.
Our procedure, whereby terms are moved out of the fitting function and incorporated into new (correlated) errors in the Monte Carlo fit data, is generally useful. Somewhat remarkably, final fit results are completely (or almost completely) independent of the number of terms that are transferred when fits are linear (or almost linear) in the associated parameters. (The general theorem from which this result follows is proven in the Appendix.) Consequently, in our analysis here, we can take N am very large -say, N am = 80 -and still have very fast fits by keepingN am very small. With N am = 80 we find, for example, that settingN am = 0 inR n (. . .) (no terms) gives essentially identical results for our quark masses and coupling as settingN am = 30 (140 terms), even though the latter fit requires 22 times more computing. We used this procedure, withN am = 0, for most of our testing and development in this project.

Truncated Perturbation Theory
The perturbative part, of the reduced moments is known at best through third order. We present coefficients r nj through j = 3 in Table III [2][3][4][5][6]; the values for n = 4-10 are exact, while r n3 is estimated for the others. In our fits we include higherorder terms by treating the coefficients of these terms as fit parameters with prior for any coefficient that hasn't been computed in perturbation theory. We set N pth = 6 since then contributions from still higher orders should be less than 0.1% (and setting N pth = 8 doesn't change our results). The perturbative coefficients for µ/m h = 1 (Table III) are small and relatively uncorrelated from order-to-order. This is less true for µ/m h = 3, which is where we wish to work (see Section II), because of log(µ/m h ) m terms. In order to capture these effects, we use renormalization group equations to express the r nj (3) coefficients (for all j ≤ N pth ) in terms of the r nj (1) coefficients and log(µ/m h ), and substitute the results from Table III for j ≤ 3 and from the prior (Eq. (21)) for j > 3. This procedure generates (correlated) priors for the unknown coefficients at µ/m h = 3 that properly account for renormalization-group logarithms.

α MS Evolution
As discussed above, we fix the ratio of µ/m h (µ) in our analysis. This means that the renormalization scale µ varies over a wide range of values for the different m h s we use. The coupling constant α MS (µ) used in the perturbative expansions for the r n s is specified at µ = 5 GeV by fit parameter α 0 , with prior where The prior corresponds to α MS (M Z ) = 0.118(3) -a very broad range, which means that the prior has little impact on our final fit results. The coupling value at any scale µ = 5 GeV is obtained by integrating (numerically) the QCD evolution equation for α MS (µ) starting with value α 0 at scale 5 GeV. We use the MS beta function through sixth order in α MS , where β 0 . . . β 3 are known from perturbation theory and β 4 is taken as a fit parameter with prior where σ β is the root-mean-square average of β 0 . . . β 3 [18,19]. We include this last term to estimate the uncertainties in our final results caused by unknown terms in the beta function.
Our simulations include vacuum polarization effects from only the three lightest quarks. We use perturbation theory, together with the c and b masses that come out of our analysis, to incorporate vacuum polarization effects from the heavier quarks into our final results for the masses and QCD coupling (using formulas from [20,21] to add the c and b quarks at scales µ = m c and m b , respectively).

Nonperturbative Condensates
As discussed in our previous paper, nonperturbative effects dominate the reduced moments when n is large. The dominant nonperturbative contribution, which is from the gluon condensate, is quite small, however, for the range of ns and quark masses we use here. We correct for it by replacing where d n is computed to leading order in perturbation theory [22] with m h = m h (m h ), which we approximate by m η h /2.27. We take which covers the range of most current estimates [23]. The correction factor in Eq. (26) adds (slightly) to the error in R latt n (and introduces new correlations between different moments, since the same α s G 2 /π is assumed for every moment, lattice spacing and quark mass).

Finite Volume Errors
We expect small errors due to the fact that our simulation lattices are only about 2.5 fm across. We allow for the possibility of finite-volume errors by replacing where ∆R pth n is the finite volume error in leading-order perturbation theory and The true finite-volume errors are expected to be smaller, because of quark confinement, than the perturbative errors that we use to model them here. We verified this by running two sets of simulations that were identical except for the spatial volume (gluon configuration sets 4 and 5 in Table I). The differences between the two simulations are smaller than our statistical errors, but the statistical errors are much smaller than our estimate above. Our error estimate here is very conservative, but has negligible impact on our final results.

Sea-Quark Masses
The sea-quark masses used in our simulations are not exactly correct. To correct for this we replace where δm l and δm s are the errors in the u/d and s masses (see [10] for more details), respectively, and This correction introduces (correlated) errors into the R latt n s that are of order 0.5-1%. Direct comparison of results from configuration sets 6 and 7 (or 1-2 and 3-4) in Table I suggests that sea-quark mass effects are no larger than 0.1%, so our error estimate is conservative.
We have only included the leading dependence on the sea-quark mass, which comes from nonperturbative (chiral) effects. Quadratic terms from perturbation theory and other nonperturbative sources are negligible.

IV. ANALYSIS AND RESULTS
We have computed reduced moments for 30 different sets of lattice spacing, lattice volume and quark masses (Table II). To extract quark masses and the QCD coupling, we fit moments with 4 ≤ n ≤ 10 from 22 of these parameter sets (the ones with am η h ≤ 1.95) -88 pieces of simulation data in all. In this section we first describe the fitting method used to extract the masses and coupling, and then we review our results.

A. Constrained Fits
We analyze all four R n s for all 22 parameter sets simultaneously using a constrained fitting procedure based upon Bayesian ideas [17]. In this procedure we minimize an augmented χ 2 function of the form where: the R latt n come from Table II with  • lattice spacings a i for each gluon configuration set, with priors specified by simulation results for r 1 /a (Table I) and the current value for r 1 (Eq. (10)); • values for am η h i , with priors specified by our simulation results (Table II).
The renormalization scales µ i are obtained from the ratio µ/m h = 3, simulation results for m η h , and Eq. (7). We take N am = 30 for our final results.

B. Results
We fit our simulation data for the reduced moments R latt n (Table II) using fit function R n (. . .) (Eq. (15)) with N am = 30, as discussed in the previous section. The best-fit values for parameters z j give us the mass-ratio function z(µ/m h = 3, m η h ) (Eq. (7)), which we plot in Figure 1. We also show our simulation results there for R latt n /r n , together with the best-fit lines for each lattice spacing. Results are shown for the three moments that depend upon z, 5 different lattice spacings, and quark masses ranging from below the c mass almost to the b mass. The simulation data were all fit simultaneously, using the same functions z(3, m η h ) and α MS (µ) (with µ = 3m η h /(2z)) for all moments. The fits are excellent, with χ 2 /88 = 0.19 for the 88 pieces of simulation data we fit. The solid line, plus gray error envelope, shows the a = 0 extrapolation obtained from our fit. This is compared with simulation results for Rn/rn for n = 6, 8, 10 from our 5 different lattice spacings, together with the best fits (dashed lines) corresponding to those lattice spacings. Dashed lines for smaller lattice spacings extend further to the right. The points marked by an "x" are for the largest mass we tried (last line in Table II); these are not included in the fit because amη h is too large. Finite-a errors become very small for the larger-n moments, causing points from different lattice spacings to overlap.
Evaluated at m ηc = 2.985(3) GeV [24], the massratio function is z(3, m ηc ) = 1.507 (7). Combining this with Eq. (9) and perturbation theory, we can obtain the following results for the MS c-quark mass at different scales: Similarly at m η b = 9.395(5) GeV [25], the mass-ratio function is z(3, m η b ) = 1.296 (8), and we obtain the following results for the MS b-quark mass at different scales: Note that the ratio m b (µ, n f )/m c (µ, n f ) is independent of µ and n f . We obtain the following result for this The other important output from our fit is a value for the parameter α 0 ≡ α MS (5 GeV, n f = 3) = 0.2034 (21).
To compare with other determinations of the coupling, we add vacuum polarization corrections from the c and b quarks, using the masses above, and evolve to the Zmeson mass [18][19][20][21]: α MS (M Z , n f = 5) = 0.1183 (7).  Figure 2 shows how consistent our simulation results are with the theoretical curve for α MS (µ, n f = 3) corresponding to our value for α 0 . For this figure we extracted values for α MS from each R n separately by dividing out the a 2 dependence and z(3, m η h ) using our best-fit parameters, and then solving for α MS by matching with perturbation theory for r n . (In our fit, of course, we fit all R n s simultaneously to obtain a single α MS for all of them.) The dominant sources of error for our results are listed in Table IV. The largest uncertainties come from: extrapolations to a = 0, especially for quantities involving b quarks; unknown higher-order terms in perturbation theory, especially for quantities involving c quarks; statistical fluctuations; extrapolations in the heavy quark mass, especially for quantities involving b quarks; and uncertainties in static-quark parameters r 1 /a and r 1 .The pattern of errors is as expected in each case. The nonperturbative contribution from the gluon condensate is negligible except for m c , again as expected; and errors due to mistuned sea-quark masses, finite volume errors, and uncertainties in MS coupling and mass evolution are negligible (< 0.05%).
The a 2 extrapolations of our data are not large. This is illustrated for m h ≈ m c in Figure 3, which shows the a 2 dependence of the reduced moments. The smallest two lattice spacings are sufficiently close to a = 0 that the extrapolation is almost linear from those points. The a = 0 extrapolated values we obtain here for the R n agree to within (smaller) errors with those in our previous paper: here we get 1.282(4), 1.527(4), 1.373(3), 1.304(2) with n = 4, 6, 8, 10, respectively, for the masses used in the figure.
We tested the stability of our analysis in several ways: • Vary perturbation theory: We chose µ = 3m h in order to keep scales large and α MS (µ) small. Our results are quite insensitive to µ, however. Choosing µ = m h , for example, shifts none of our results by more than 0.2σ, and leaves all errors unchanged except for m c (3), where the error increases by a third. Taking µ = 9m h shifts results by less than 0.4σ, and reduces the m c error by a third, leaving others only slightly reduced. Adding more terms to the perturbative expansions (N pth = 6 → 8) also has essentially no effect on the results. The prior for the unknown perturbative coefficients (Eq. (21)) is twice as wide as suggested by our simulation results (using the empirical Bayes criterion [17]); we choose the larger width to be conservative.
• Include more/fewer finite-a corrections: We set N am = 30 for our results above. Using N am = 15 gives results that differ by less than 0.5σ for m b and much less for the other quantities. Much larger N am s can be tested easily using the trick described in Section III B 2. For example, replacing R latt n bȳ R latt n (Eq. (18)) with N am = 80 andN am = 30 gives results that are essentially identical to those above. As discussed above, takingN am = 0 with the same N am also gives the same results and is 22 times faster (see the appendix for further discussion).
• Change n dependence of finite-a corrections: Replacing the n-dependent prior for the expansion coefficients (Eq. (17)) by the n-independent prior 0 ± 0.5 causes changes that are less than 0.3σ. The width of the original prior is optimal according to the empirical Bayes criterion -that is, it is the width suggested by the size of finite-a deviations observed in our simulation data.
• Add more/fewer Λ/mη h terms in z: Increasing the number of terms in the expansion for z from N z = 4 to 6 changes nothing by more than 0.1σ. Decreasing to N z = 3 also has no effect. Again the width of the prior is optimal according to the empirical Bayes criterion.
• Include more/fewer moments: Keeping all mo-ments 4 ≤ n ≤ 18 changes nothing by more than 0.5σ and reduces errors slightly for everything other than m b , where the errors are cut almost in half: m b (10) = 3.623(15) GeV or m b (m b ) = 4.170(13) GeV, both for n f = 5. We continue to restrict ourselves to moments with n ≤ 10 because these are the only moments for which we have exact third-order perturbation theory. Keeping just n = 4, 6 gives almost identical results for m c and α MS , with almost the same errors, but doubles the error on m b .
• Omit simulation data: The coarsest two lattice spacings (configuration sets 1-5) affect our results only weakly. Leaving these out shifts no result by more than 0.5σ and leaves errors almost unchanged. Leaving out the smallest lattice spacing, however, increases errors significantly (almost double for α MS ), while still shifting central values by less than 0.5σ.
• Add large masses: Including cases with am η h > 1.95 from Table II leads to poor fits. The excluded data, however, do not deviate far from the best-fit lines. For example, the points marked with an "x" in Figure 1 are for the largest mass we studied, corresponding to m η h = 9.15 GeV (last line in Table II). Although am η h is too large for this case to be included in our fit, the values of R n /r n are only slightly below the fit results.
It is possible to extract the ratio of quark masses m b /m c directly, without using the moments and without using perturbation theory. This provides an excellent nonperturbative check on our results from the moments.
Ratios of quark masses are UV-cutoff independent and therefore the ratio of MS masses for any µ and n f , where m 0b and m 0c are the bare quark masses in the lattice quark action that give correct masses for the η c and η b , respectively. We obtain accurate mass ratios from this relationship by extrapolating to a = 0. We used such a method recently to determine m c /m s [11].
Here we have to modify our earlier method slightly because we cannot reach the b-quark mass directly, but rather must simultaneously extrapolate to the b mass and the continuum limit. This is most simply done by determining the functional dependence of the ratio on the η h mass and the lattice spacing. The ratio of MS masses is then given by the experimental masses of the η c and η b and the equation: It might seem simpler to fit m 0h directly, rather than the ratio w; but using w significantly reduces the m η h dependence (and therefore our extrapolation errors), and also makes our results quite insensitive to uncertainties in our values for the lattice spacing. We parameterize function w with an expansion modeled after the one we used to fit the moments: where, as for the moments, Coefficients c ij and w n are determined by fitting function w(m η h , a) to the values of 2am 0h /(am η h ) from Table II. The fit also determines the parameters Z m (a), one for each lattice spacing, which account for the running of the bare quark masses between different lattice spacings. The finite-a dependence is smaller here than for the moments, because the η h is nonrelativistic (finite-a errors are suppressed by additional powers of v/c [8]), and the variation with m η h stronger (twice that of z(3, m η h )). So here we use priors c ij = 0 ± 0.05 (44) w n = 0 ± 4 Z m (a) = 1 ± 0.5 with N w = 8. We again take N am = 30, although identical results are obtained with N am = 15.
Our fit results are illustrated by Figure 4 which plots the ratio m 0h /m η h divided by m 0c /m ηc for a range of η h masses. Our data for different lattice spacings is compared with our fit, and with the a = 0 limit of our fit (solid line). The fit is excellent, with χ 2 /22 = 0.42 for the 22 pieces of data we fit (we again exclude cases with am η h > 1.95). Using the η c and η b masses from Section IV B, and Eq. (41) with the best-fit values for the parameters, we obtain finally m 0b m 0c → 4.49(4) as a → 0 (45) which agrees well with our result from the moments (Eq. (36)). The solid line shows the a = 0 extrapolation obtained from our fit. This is compared with simulation results for our 4 smallest lattice spacings, together with the best fits (dashed lines) corresponding to those lattice spacings. The point marked by an "x" is for the largest mass we tried (last line in Table II); this was not included in the fit because amη h is too large.

VI. α MS FROM WILSON LOOPS
In a recent paper [26], we presented a very accurate determination of the QCD coupling from simulation results for Wilson loops. Here we want to compare those results to the value we obtain from heavy-quark correlators. First, however, we must update our earlier analysis to take account of the new value for r 1 [10] given in Eq. (10) and improved values for r 1 /a [13] given in Table I. (The Wilson-loop paper uses some additional configuration sets: from Table II in that paper, sets 1, 6, 9, and 11 whose new r 1 /as are 1.813 (8) This agrees very well with the result in the earlier paper, α MS (M Z ) = 0.1183(8), but has a slightly smaller error, as expected given the smaller error in r 1 . This new value also agrees well with our very different determination from heavy-quark correlators (Eq. (38)). A breakdown of the error into its different sources can be found in Table IV of [26] (reduce the r 1 and r 1 /a errors in that table by half to account for the improved values used here).

VII. CONCLUSIONS
In this paper, we improve significantly on our previous determinations of the QCD coupling and c-quark mass from heavy-quark correlators. This is principally due to the inclusion of a new, smaller lattice spacing in our analysis. We also generated results for a variety of quark masses near m c , allowing us to interpolate more accurately to the physical value of m c . New third-order perturbation theory makes R 10 as useful now as R 4 , R 6 , and R 8 were in the earlier paper. Finally, in this paper, we fit multiple moments simultaneously, determining consistent values simultaneously for both the QCD coupling and the quark masses for all moments. Previously we examined each moment or ratio of moments independently, extracting m c s or α MS s independently of each other. Our  (7), agree well with our older results of 0.986(10) GeV and 0.1174 (12), respectively [1].
The much heavier b quark is usually analyzed using effective field theories like NRQCD or the static-quark approximation. By using very small lattice spacings and the very highly improved HISQ discretization for the heavy quarks, we are able to extend our analysis almost to the b-quark mass, using the same relativistic discretization that we use for c and lighter quarks. A 1.5% extrapolation of z(3, m h ), from the largest m η h used in our fits to m η b , gives us a new, accurate determination of the b-quark mass, m b (10 GeV, n f = 5) = 3.617(25) GeV.
This calculation demonstrates the utility of the HISQ formalism for studying b quarks on lattices that are computationally accessible today. This represents a breakthrough for b physics on the lattice since far greater precision becomes possible when all quarks are treated using the same formalism, and that formalism is relativistic and has a chiral symmetry. Even better would be to work right at the b mass, as opposed to extrapolating from nearby; this would require a lattice spacing of order 0.03 fm.
Both of our new c and b masses agree well with non-lattice determinations from vector-current correlators and experimental e + e − collisions. A recent analysis of the continuum data gives [7]  Function z(µ/m h , m η h ) is a by-product of our analysis. It relates the MS quark mass m h (µ) to the η h mass (Eq. (7)). We show our result again in Figure 6 for µ = 3m h , as well as for µ = m h and µ = m h /2, which we obtain by evolving perturbatively from µ = 3m h . The latter two curves are relatively flat, and the last surprisingly close to 1 for most masses.
Questions have been raised about the way perturbation theory is used in analyzing the perturbative parts of the moments [27]. Like [7] we favor using larger scales than m c for c-quark correlators, but, as we have shown, our results are quite insensitive to µ over a broad range. Furthermore, the fact that our results, from pseudoscalar-density correlators, agree so well with the continuum results, from vector-current correlators, is also compelling evidence that perturbation theory is being handled correctly. We also find consistent results from several different moments, which is only possible if perturbation theory is working well. Compare, for example, Figure 7 for the moments, as a function of m η h , with the plots of R n /r n in Figure 1. Figure 7 shows very different m η h behavior, at the 10-20% level, for different moments R n ; Figure 1, where the perturbative part r n is divided out, shows behavior that is almost momentindependent.
An additional check on our use of perturbation theory comes from the close agreement between our perturbative result for the ratio m b /m c of MS masses (Eq. (36)) and our nonperturbative result for the ratio of HISQ masses (Eq. (45)). These should be and are equal to within our 1% errors. Taken together they suggest a composite result of: The validity of our perturbative analyses is further supported by the close agreement between the QCD coupling we get from the heavy-quark correlators, α MS (M Z ) = 0.1183 (7), and that obtained from Wilson loops, 0.1184 (6). These are radically different methods for determining the coupling. The first relies upon a continuum quantity, extrapolated to a = 0, and continuum perturbation theory. The second relies upon quantities that are highly sensitive to the UV cutoff (π/a) but are analyzed to all orders in the cutoff using lattice perturbation theory. Systematic errors are almost completely different in the two cases. The fact that they agree to within our 0.6% uncertainties is highly nontrivial evidence that perturbative and other potential errors are understood.
Our coupling values also agree well with determinations from non-lattice methods. Figure 8 summarises recent results that were included in a world average by Bethke [28]. The world average result, 0.1184 (7), was dominated by our previous determination from the Wilson loop analysis. The average excluding our result was 0.1186 (11), which also agrees well. Including our new results into a new error-weighted world average gives α MS (M Z ) = 0.1184(4).
Our new c mass is the most accurate currently available. With it we can improve slightly on our recent determination of light quark masses using an accurate value for m c /m s , 11.85 (16), derived completely nonperturbatively from lattice calculations [11]. Our new c mass, which becomes 1.093(6) GeV when converted to n f = 3 at 2 GeV, implies:  Figure 9. Agreement is excellent, but our uncertainties are much smaller in every case, and by an order of magnitude for the strange and light quarks.

Particle Data Group's 2009 values in
Finally we note that the consistency between quark masses from lattice and non-lattice analyses, and between couplings from heavy-quark correlators and Wilson loops provides further evidence that taste-changing interactions in the HISQ and ASQTAD quark formalisms are understood and vanish as a → 0. While early concerns about the validity of these formalisms have been largely addressed both by formal arguments [13,[30][31][32][33][34] and by extensive empirical studies [8-11, 26, 35-39], it remains important to test the simulation technology of lattice QCD with increasing precision, given the growing importance of lattice results for phenomenology.
In Section III B 2 we used a trick to simplify our fits by, in effect, transferring fit terms from the fit function into the errors of the fit data. This trick can greatly speed up complicated fits. Here we present a formal derivation of this procedure for three increasingly complicated situations.

A. Linear Least Squares -Exact Data
Assuming we know D values y i for a quantity y which can be expressed as a power series in x, we wish to obtain a best fit for the first F unknown coefficients c n . The c n are then our random variables. If we are able to make reasonable estimates for their means and standard deviations σ n , in the absence of additional information, maximizing entropy suggests a Gaussian prior of For simplicity, we assume throughout that the c n are uncorrelated and have a prior mean of zero; extending to more general cases is straightforward.
If we knew all coefficient values, then the data y i would be completely determined, with Bayes' theorem P (c|y) ∝ P (y|c)P (c) (55) allows us to convert this into a distribution for c given the data y.
If we are only interested in fitting a subset of coefficients c n< with n < F , we integrate over the remaining c n> , giving We replace the delta function by its Fourier representation, integrate over first the c n> , then the Fourier variables, to obtain (det σ 2 ∆ ) −1/2 e −∆y·(2σ 2 ∆ ) −1 ·∆y . Here is the discrepancy between the measured y i and the portion of the series to be kept in the fit, the dot product sums over the D data points, and The correlation matrix σ 2 ∆ is independent of c < (so the determinant is constant), and is the same as one would compute directly by using c m> c n> c> = σ 2 n δ mn .
Finally, we fit c < by minimizing χ 2 , which includes these correlations and is augmented by the remaining c < priors. Because the distribution is Gaussian, the c < at their minima are equal to their average values. The correlation matrix σ 2 ∆ properly accounts for correlations in the discrepancy, due to the neglected terms, between y and the portion of the series retained. If F terms are kept in the series, σ ∆ is O(x F ), enforcing agreement between y and the finite series to this order, as appropriate. It also suggests an alternative but equivalent approach. We may define new random (rather than exact) versions of y, whose correlation matrix is σ 2 ∆ , by moving the c > terms to the left side of Eq. (52). Using the truncated series as a model for these random data, straightforward application of Bayes' theorem [17] again implies the distribution in Eq. (57).
One useful consequence is that, as long as we include the correlations for the c > , we may arbitrarily reduce the number of coefficients c < retained, even to as few as one, and still obtain the same minimization values. To see this, note that to compute a particular c n< , we could start with the full distribution and integrate over all cs. The integral over c > produces P (c < |y), which we then use in the integral over c < ; the result will be the same regardless of where the dividing line is set, as long as it does not include c n< . (We could even include in σ 2 ∆ terms of order less than n.) Because averaging and minimization give the same result, the minimization value for c n< will also remain unchanged. This is also true of the c n< error. While the result is the same, reducing the number of terms in the series to fit can significantly improve the fitting time.

B. Fits to Nonlinear Functions -Exact Data
We now consider fitting to the data y i a general function g i (c n ) not necessarily linear in the parameters c, and where we assume y i = g i (c n ) exactly for properly chosen c n . Now P (y|c) ∝ D−1 i=0 δ(y i − n g i (c n )) . (62) Combining with the prior P (c) and integrating over the c > gives P (c < |y). If our estimate of prior means is good, expanding g around c > = 0 should give a reasonable approximation; an expansion to first order gives a Gaussian. More specifically, defining and and integrating over c > in this Gaussian approximation gives as before P (c < |y) ∝ e − n < c 2 n < /2σ 2 n × (65) (det σ 2 ∆ (c < )) −1/2 e −∆y·(2σ 2 ∆ (c<)) −1 ·∆y , but with This is again the correlation one would compute directly for ∆y i ∆y j c> after expanding g to first order in c > . We have not expanded in c < , so σ 2 ∆ depends on c < , the determinant in front is not constant, and the dependence of ∆y on c < is not in general linear. In practice, however, we will often further approximate the distribution by setting the c < to their prior means in σ 2 ∆ (c < ) before minimization.
Because g(c < ) is nonlinear, c < from minimization can differ slightly from c < , and due to approximations made, can vary somewhat with the number of terms retained.

C. Fits to Data with Intrinsic Statistical Errors
Finally we consider the most general case, in which the data y contribute intrinsic statistical uncertainties in addition to those associated with the truncated series. If we measure a range of values for y with an average y and correlation matrix σ 2 y , then for sufficiently large samples we expect a Gaussian distribution P ( y |c) ∝ e −( y −g(c))·(2σ 2 y ) −1 ·( y −g(c)) (67) rather than the delta function above. Combining with the prior P (c) gives P (c| y ).
The resulting correlation matrix is a combination of true statistical and neglected series contributions, with as one would obtain by including both sources of uncertainty in computing ∆y i ∆y j directly. With no statistical fluctuations in y, σ 2 y = 0, and it reduces to the previous result. When σ 2 y is nonzero but small, σ 2 ∆ still makes an important contribution.

D. Application to this Paper
We used the technique described here in much of our testing and tuning (but not for our final results) to speed up the (am η h /2) 2 fit. As described in Section III B 2, we kept corrections through order N am = 80 but moved all butN am ≪ N am out of the fit function and into the errors for the reduced moments. If we setN am = 3, for example, our fit to the R n simulation data changes from Figure 7 to Figure 10. The smallN am means that each point in Figure 10 has much larger error bars, coming from (am η h /2) 2 terms moved into the R n s. The final fit results, however, are almost identical in both cases (to within less than 0.1σ), with the same errors. Note that the R n errors in Figure 10 are highly correlated, which is why the fit curve passes through the central value for each point. As discussed above these correlations are essential if results are to be independent of the value ofN am .