Up-, down-, strange-, charm-, and bottom-quark masses from four-flavor lattice QCD

We calculate the up-, down-, strange-, charm-, and bottom-quark masses using the MILC highly improved staggered-quark ensembles with four flavors of dynamical quarks. We use ensembles at six lattice spacings ranging from $a\approx0.15$~fm to $0.03$~fm and with both physical and unphysical values of the two light and the strange sea-quark masses. We use a new method based on heavy-quark effective theory (HQET) to extract quark masses from heavy-light pseudoscalar meson masses. Combining our analysis with our separate determination of ratios of light-quark masses we present masses of the up, down, strange, charm, and bottom quarks. Our results for the $\overline{\text{MS}}$-renormalized masses are $m_u(2~\text{GeV}) = 2.130(41)$~MeV, $m_d(2~\text{GeV}) = 4.675(56)$~MeV, $m_s(2~\text{GeV}) = 92.47(69)$~MeV, $m_c(3~\text{GeV}) = 983.7(5.6)$~MeV, and $m_c(m_c) = 1273(10)$~MeV, with four active flavors; and $m_b(m_b) = 4195(14)$~MeV with five active flavors. We also obtain ratios of quark masses $m_c/m_s = 11.783(25)$, $m_b/m_s = 53.94(12)$, and $m_b/m_c = 4.578(8)$. The result for $m_c$ matches the precision of the most precise calculation to date, and the other masses and all quoted ratios are the most precise to date. Moreover, these results are the first with a perturbative accuracy of $\alpha_s^4$. As byproducts of our method, we obtain the matrix elements of HQET operators with dimension 4 and 5: $\overline{\Lambda}_\text{MRS}=555(31)$~MeV in the minimal renormalon-subtracted (MRS) scheme, $\mu_\pi^2 = 0.05(22)~\text{GeV}^2$, and $\mu_G^2(m_b)=0.38(2)~\text{GeV}^2$. The MRS scheme [Phys. Rev. D97, 034503 (2018), arXiv:1712.04983 [hep-ph]] is the key new aspect of our method.


I. INTRODUCTION
Quark masses are fundamental parameters of QCD. They must be known accurately for precise theoretical calculations within the Standard Model, especially for testing whether quarks receive mass via Yukawa couplings to the Higgs field. Because of confinement, the quark masses can be defined only as renormalized parameters of the QCD Lagrangian. Thus, they must be determined by comparing theoretical calculations of an appropriate set of observables to experimental measurements of those observables. Lattice QCD makes it possible to calculate in a nonperturbative way simple observables, such as hadron masses.
To determine the quark masses in lattice QCD, one needs to tune the bare lattice quark masses such that a suitable set of hadron masses coincide with their experimental values.
The resulting bare masses must be renormalized, preferably to a regularization independent scheme, such as the recently introduced minimal renormalon subtracted (MRS) mass [1]. One approach is to use lattice perturbation theory, but multiloop calculations are difficult so, in practice, nothing more than two-loop matching [2][3][4] is available in the literature. Another is to use nonperturbative renormalization to, for example, momentumsubtraction [5,6] or finite-volume [7,8] schemes. Finally, one can use lattice gauge theory to obtain quantities in continuum QCD and apply multiloop continuum perturbative QCD to extract the quark masses. An example of the latter is the analysis of quarkonium correlators [9]. In practice, no regularization-independent scheme is in such common use as the modified minimal subtraction (MS) scheme [10] of dimensional regularization, so we shall use MS to quote results.
Our method studies how a heavy-light meson mass depends on the mass of its heavy valence (anti)quark [11][12][13]. Like the quarkonium correlators, our approach requires only continuum perturbation theory. On the other hand, the binding energy of a heavy-light meson is of order Λ QCD , so it is necessary to use heavy quark effective theory (HQET) to separate long-and short-distance scales. In this way, we can obtain the masses of the charm and bottom quarks and, at the same time, HQET matrix elements [13]. Because this analysis uses as inputs the bare masses of the up, down, and strange quarks-tuned to reproduce the pion and kaon masses [14], it also yields the renormalized masses of these quarks.
Following Ref. [13], our analysis is based on the HQET formula for the heavy-light meson mass [15] M H ( * ) = m h + Λ + where M H ( * ) is the pseudoscalar (vector) meson mass; m h is the heavy-quark mass; and Λ, µ 2 π , and µ 2 G (m h ) are matrix elements of HQET operators with dimension 4 and 5. The last three correspond to the energy of the light quarks and gluons, the heavy quark's kinetic energy, and the spin-dependent chromomagnetic energy, with coefficient d H = 1 for pseudoscalar mesons and d H * = − 1 3 for vector mesons. The chromomagnetic operator has an anomalous dimension, known to three loops [16], so µ 2 G (m h ) depends logarithmically on the mass m h . The strategy is to use lattice QCD to compute M H ( * ) as a function of m h and fit Eq. (1.1) to distinguish the terms on the right-hand side including, in principle, higher orders in 1/m h [13].
The utility of Eq. (1.1) rests on the definition of the quark mass m h . In HQET, the natural definition is the pole mass (also known as the on-shell mass). Although the pole mass is infrared finite [17] and gauge independent [17,18] at every order in perturbation theory, its value is ambiguous when all orders are considered [19,20]. At large orders, the coefficients of the self energy grow factorially, and a possible interpretation via Borel summation is obstructed by a series of renormalon singularities [19,20]. This behavior is a manifestation of the strongly-coupled long-range gluon field that, remarkably, appears in perturbation theory. Note that because M H is unambiguous, the ambiguity in m h must be canceled by those in Λ, µ 2 π , and higher-dimension terms. 1 To address this problem, some of us introduced the minimal renormalon-subtracted (MRS) mass in a companion paper [1]. It is defined by Eq. (2.24) of Ref. [1], where m = m MS (m MS ), the r n are the coefficients relating the MS mass to the pole mass, R n denote their asymptotic behavior, and J MRS (m), which is defined in Eqs. (2.25) and (2.26) of Ref. [1], is the unambiguous part of the Borel sum of R n α n+1 s . To compute m MRS one uses the known behavior of the R n [21][22][23], including their overall normalization [23]. In deriving Eq. (1.2), Ref. [1] puts the leading renormalon ambiguity into a specific quantity of order Λ MS , denoted δm, and transfers it from m h to Λ. Below we write m h,MRS and Λ MRS to denote the unambiguous definitions of m h and Λ in the MRS scheme.
A second feature of our technique may seem almost trivial. In Eq. where am r is the bare mass (in lattice units) of staggered fermions, and the subscript r labels a reference mass; see Sec. III. Owing to the remnant chiral symmetry of staggered fermions, the first factor in Eq. (1.3a) is 1+O(a 2 ). In Eq. (1.3b), the factors are, respectively, a convenient fit parameter, the factor to run from scale µ to m h , the quantity in the big parentheses in Eq. (1.2), and the ratio of the freely chosen heavy-quark lattice mass to the reference mass. Equation (1.3) plays a key role: with r n for MS in Eq. (1.2), the first factor in Eq. (1.3b) is in the MS scheme; with J MRS removing the leading renormalon ambiguity, the product on the right-hand side of Eq. (1.3b) is indeed the MRS mass. By taking m r = 0.4m s (the so-called p4s approach), the analysis yields m s as well as the heavy-quark masses m c and m b . The third important feature of our work is a data set with a wide range of lattice spacing, heavy-quark mass, and light valence and sea masses. These data, which were generated in a companion project to compute the B-and D-meson decay constants [14], are very precise, with statistical errors of 0.005-0.12%. It is very challenging to take advantage of the statistical power and parameter range of the data set. In this paper, we use heavymeson rooted all-staggered chiral perturbation theory [24] (HMrASχPT) to describe the dependence of the heavy-light pseudoscalar meson masses on the light mesons. To make possible a fit to lattice data, Ref. [1] combined the next-to-leading-order HMrASχPT with the MRS mass to write heavy-light meson masses as a function of lattice spacing and heavy-1 By forming the spin average, 1 4 (M H + 3M H * ), and spin difference, M H * − M H , it is easy to see that spin-independent and spin-dependent ambiguities are distinct. and light-quark masses. The fit function, by construction, has the correct nonanalytic form in the chiral and HQET limits. Here, it is extended with enough analytic terms to mimic higher-order corrections and obtain a good fit.
A preliminary report of this analysis can be found in Ref. [51]. Instead of the MRS mass, at that time we used the renormalon-subtracted (RS) mass [22], which also subtracts the leading renormalon ambiguity but at the same time introduces a factorization scale ν f . In principle, the MS masses emerging from Eqs. (1.3) and (1.1) should not depend on ν f , but we found more dependence than one would like. Moreover, it turns out to be necessary to introduce three scales in all, ν f < µ < m h , with µ being used for α s [51]. For that reason, we prefer the MRS over the RS mass.
This paper is organized as follows. Section II contains a description of the lattice-QCD simulations, focusing on the way we eliminate the lattice scale in favor of physical units. In Sec. III, we present our function of quark masses and lattice spacing that describes masses of heavy-light pseudoscalar mesons. In Sec. IV, we perform a combined-correlated fit to the meson masses; the fit is then extrapolated to the continuum and interpolated to physical values of the light quark masses. In Sec. V, we present our final results for the masses of the strange, charm and bottom quarks as well as quark mass ratios m c /m s , m b /m s , and m b /m c . Combining our results with our separate determination of the quark-mass ratios m u /m d and m s /m l , where m l ≡ 1 2 (m u + m d ), we also report the up-and down-quark masses. In addition, we present our lattice-QCD determinations of Λ MRS , µ 2 π , and µ 2 G (m b ) as well as flavor splittings and low-energy constants of heavy-meson chiral perturbation theory. Section VI compares our main results with work in the literature and offers some remarks on further work. An appendix gives the correlation matrices of the MRS masses of the charm and bottom quarks with HQET matrix elements, and of the charm-quark mass and quark-mass ratios.

II. SIMULATIONS SUMMARIZED
The lattice data used in this work come from the same correlation functions used to determine leptonic decay constants of charmed and b-flavored mesons in a companion paper [14]. For a full description of the simulation, the reader should consult Ref. [14]. Here we provide a brief summary.
We employ a data set that includes ensembles with five values of lattice spacings ranging from approximately 0.12 fm to 0.03 fm, enabling good control over the continuum extrapolation. Ensembles at a sixth lattice spacing, approximately 0.15 fm, are used only to estimate the continuum extrapolation error. The data set includes ensembles with the light (up-down), strange, and charm sea masses close to their physical values ("physical-mass ensembles") at all but the smallest lattice spacing, 0.03 fm. The data set also includes ensembles where either the mass of light sea quarks is heavier than in nature, or the mass of the strange sea quark is lighter than in nature, or both. As in Ref. [14], we set the scale of the lattice spacing a with a two-step procedure that uses the value of f π from the Particle Data Group (PDG), f π,PDG = 130.50 (13) MeV [52], combined with the so-called p4s method.
The first step in the scale-setting procedure takes f π,PDG to set the overall scale on each physical-mass ensemble. On these ensembles, we tune the valence light, strange, and charmed quark masses to reproduce the pion, kaon, and D s -meson masses. Then we calculate M p4s and f p4s , which are the mass and decay constant of a pseudoscalar meson with both valence-quark masses set equal to m p4s ≡ 0.4m s . We then form the ratio R p4s ≡ f p4s /M p4s and take the continuum limit of f p4s and R p4s . These values and those of the quark mass ratios are then used as inputs to the second step of the procedure, which we call the p4s method. In the p4s method, the values of am p4s and af p4s are calculated on a given physical-mass ensemble, with a = 0, by adjusting the valence-quark mass until af p4s /aM p4s equals the physical-mass continuum limit of R p4s : In the p4s method, all ensembles at the same bare gauge coupling, β = 10/g 2 0 , as a given physical-mass ensemble are chosen to have the same lattice spacing a and the same am p4s . This choice is known as a mass-independent scale-setting scheme.
At a ≈ 0.03 fm, we have only a 0.2m s ensemble, so this procedure cannot be carried out. In this case, we rely on the derivatives with respect to a, which are given in Ref. [50].

III. CONSTRUCTION OF THE FIT FUNCTION
In this section, we discuss in detail how to construct a function of quark masses and lattice spacing that describes masses of heavy-light pseudoscalar mesons. To this end, we use three effective field theories (EFTs), HQET, and HMrASχPT, as mentioned already, and the Symanzik effective theory of cutoff effects [34,53,54]. We start with the merger of HQET and HMrASχPT [1] and incorporate generic lattice-spacing dependence, as well as higher-order terms in HQET and HMrASχPT. Putting everything together, we obtain an EFT fit function for masses of heavy-light pseudoscalar mesons.
A. Leading-order χPT Let us start with fixing our notation for quark masses associated with lattice ensembles with 2+1+1 flavors of quarks. We use m l , m s , and m c to denote the simulation masses of the light (up-down), strange, and charm quarks, respectively; without the primes, we use m l = 1 2 (m u + m d ), m s , and m c to denote the correctly tuned masses of the corresponding quarks; last, we use m q to denote a generic light quark mass. Further, we use H where B 0 is the low energy constant (LEC) in the relation m 2 π = B 0 (m u + m d ) between the pion mass and the quark mass; d H ( * ) = 1 (− 1 3 ) for pseudoscalar (vector) mesons; λ 1 and λ 1 are LECs that appear in (continuum) heavy-meson chiral perturbation theory (HMχPT) [55]; and δM H ( * ) x is the one-loop corrections to the mass of the H ( * ) x meson in HMrASχPT [1]. The arguments of M H ( * ) x and δM H ( * ) x in Eq. (3.1) correspond to the light valence-quark mass; the set of three light sea-quark masses, which are not necessarily tuned to their physical values; and the lattice spacing a. As usual for a one-loop χPT result, δM Hx contains a term nonanalytic as m 2 π → 0 (a "chiral log"). For the pseudoscalar mesons with (2 + 1) light flavors in the sea, we have.
where the indices S and Ξ run over light sea-quark flavors and meson tastes, respectively; M Sx,Ξ is the mass of the pseudoscalar meson with taste Ξ and flavors S and x; ∆ * is the lowest-order hyperfine splitting; δ Sx is the flavor splitting between a heavy-light meson with light quark of flavor S and one of flavor x; g π is the H-H * -π coupling; δ A and δ V are the taste-breaking hairpin parameters; a 2∆ is the mean-squared pion taste splitting; and λ a 2 and λ a 2 are parameters in SχPT related to taste breaking in meson masses. Definitions of the residue functions R [n,k] j , the sets of masses in the residues, and the chiral-log function K 1 at infinite and finite volumes are given in Ref. [1] and references therein. The expression for δM H * x is also given in Ref. [1], but because we have lattice data only for pseudoscalar mesons, it is not needed here.
In Eq. (3.1), we set so that in the continuum limit the usual expression is recovered for physical values of sea-quark masses and m x = m q . With this choice for C, the values that we obtain for Λ MRS , µ 2 π and µ 2 G (m h ) are readily applicable for calculations in HQET. 2 In this work, we set m q = 1 2 (m u + m d ), and we report Λ MRS , µ 2 π and µ 2 G (m h ) for this choice.
At this stage, the fit parameters are m r,MS (µ = 2 GeV) via Eq. (1.3), Λ MRS , the kinetic energy µ 2 π , the chromomagnetic energy µ 2 G (m b ) from which we obtain µ 2 G (m h ) as in Eq. (3.6) below, and the LECs λ 1 , λ 1 , λ a 2 , and λ a 2 . Ideally, one would have data for both pseudoscalarand vector-meson masses, and then one could set up separate fits for spin-independent and spin-dependent terms. In this work, however, only the pseudoscalar masses are available. The experimental masses of the B * and B mesons can be used to estimate which neglects contributions to the hyperfine splitting suppressed by a power of 1/m b . The chromomagnetic operator has an anomalous dimension, however, so we obtain using the three-loop relation [16] for the Wilson coefficient C cm (m h ). For four active flavors, As discussed in Sec. I and Ref. [1], the matrix elements of HQET suffer in general from ambiguities related to renormalon singularities, although the ambiguities cancel in observables such as the meson mass. For instance, the ambiguity in Λ cancels the leading-renormalon ambiguity in the pole mass. By construction, only the leading renormalon is removed to define the MRS mass. In principle, renormalon ambiguities in µ 2 π and µ 2 G (m h ) remain. In practice, numerical investigation indicates that the subleading infrared renormalon of the pole mass is small [1], which implies that the corresponding renormalon ambiguity in µ 2 π is not large. Moreover, the leading spin-dependent renormalon in µ 2 G is suppressed by a further power of 1/m h .

B. Higher-order terms in χPT
Because we have very precise data with statistical errors of 0.005-0.12%, we can anticipate that NLO χPT is not enough to fully describe the quark-mass dependence, especially for data with m x near m s . We therefore extend the function given in Eq. (3.1) by adding higherorder analytic corrections in powers of light quark masses and in inverse powers of the heavy quark mass. For the expansion in inverse powers of the heavy-quark mass, we introduce the dimensionless variable with Λ HQET = 600 MeV. Then the natural size of coefficients of the 1/m h corrections is of order 1. For expansion in light quark masses, following Refs. [14,50], we define dimensionless quark masses, which are natural expansion parameters in χPT: where q can be either the valence or sea light quarks. For simplicity, we drop the primes on the simulation x q s. The quark masses in the formula for δM Hx can also be expressed in We include all mass-dependent analytic terms at order x 2 q by adding to the expression for M Hx in Eq. (3.1). With f π to set the overall scale of these higher-order terms, the coefficients q i become of order 1 or less. We also include all mass-dependent analytic terms at order x 3 q , namely In practice, one can expect the terms without x x to be less important, but we keep all of them for consistent power counting.
To improve the expansion in inverse powers of the heavy quark, we add with three fit parameters ρ i to the right-hand side of Eq. (3.1). We also add w h and w 2 h corrections to the LECs λ 1 , λ 1 and g π ; and w h corrections to the fit parameters q i in Eq. (3.10).
The heavy quark mass also affects the hyperfine splitting ∆ * and the flavor splitting δ Sx in Eq. (3.2). Although we could express these quantities in terms of µ 2 G (m h ) and λ 1 , we exploit the experimental values for the hyperfine splittings and flavor splittings in the D and B systems to calculate ∆ * and δ Sx for different quark masses. See our companion paper on decay constants [14] for details.
We now discuss the effects of mistuning in the sea charm-quark mass m c . The effects can be divided into two parts: the effects on the pole mass (and, hence, the MRS mass) and the effects on the effective theory after the charm quark is integrated out. The former effects are taken into account in calculating the MRS mass from the MS mass; cf. Eq. (3.24). We treat the latter effects as in Ref. [14]. We use Λ where the extra fit parameter k 1 describes higher-order corrections to Eq. (3.13). We must also include generic lattice artifacts in our analysis. Taste-breaking discretization errors from staggered fermions are already included in Eq. (3.2). In addition to these effects, various discretization errors, from gluons for example, must be taken into account. We include the leading lattice artifacts for Λ MRS by replacing where Λ is the scale of generic discretization effects, set to 600 MeV in this analysis. The factor of α s in the second-order term arises because the HISQ action is tree-level improved to order a 2 . Note that Λ MRS is not affected by heavy-quark discretization errors. As discussed in the Appendices of Ref. [14], at leading order (LO) in HQET, heavy-quark discretization errors only affect the normalization of the heavy-quark state. Thus, Λ MRS and also λ 1 , λ 1 and g π at leading order in 1/m h are free of heavy-quark discretization errors. For λ 1 we replace 16) where the c 3 term is added to incorporate effects of heavy-quark discretization errors. We incorporate similar corrections for λ 1 and g π . Finally, we add α s (aΛ) 2 and α s (am h ) 2 corrections to µ 2 π and µ 2 G (m b ), and α s (aΛ) 2 corrections to the parameters q i in Eq. (3.10).

C. Heavy-quark mass
Although the MRS mass is the key to our interpretation of the HQET mass formula, as indicated in Eq. (1.3) we arrange the fit to yield the MS mass. For am 1, the relation between the MS and bare masses is where a in the denominator is set from the scale setting quantity (here f p4s , as described in Sec. II). With staggered fermions, there is no additive mass renormalization, and to eliminate tree-level discretization errors from Eq. (3.17), we take the mass am to be the tree-level pole mass. 3 Taking the ratio between two masses 4 where the dots stand for higher-order terms in a 2 and α s . In fact, each higher order in α s is also multiplied by a quantity of order a 2 , as stated in the Introduction. In this analysis, we set the reference-quark mass m r to m p4s ≡ 0.4m s and the scale of the MS scheme to µ = 2 GeV. Thus, m p4s,MS (2 GeV) is a free parameter left to be determined in the fit to lattice data; cf. Eq. (1.3b).
To incorporate further heavy-quark discretization effects into Eq. (3.18), we multiply the right-hand side of Eq. (1.3b) by where the dimensionless coefficients k n are free fit parameters, and We multiply am h by a factor of 2/π so that the parameters k n become of order 1, based on the radius of convergence of various tree-level formulas for the HISQ action; see Appendix A of Ref. [14]. Because (am p4s /am c ) 2 ≈ 0.001, the effects of a nonzero value of am p4s are negligible compared with the heavy-quark discretization effects. To incorporate generic lattice-spacing dependence into our analysis, we additionally multiply the right-hand side of Eq. (1.3b) by 1 +c 1 α s (aΛ) 2 +c 2 (aΛ) 4 +c 3 (aΛ) 6 . (3.21) To complete our approach to introducing m p4s,MS (2 GeV) via m h,MRS , we must describe the calculation of the second and third factors in Eq. (1.3b). The second factor simply uses the anomalous dimension to run from µ = 2 GeV to the self-consistent scale where with four active flavors [57] C(πu) = u 12/25 1 + 1.01413u + 1.
The coefficient of u 4 is obtained from the five-loop results for the quark-mass anomalous dimension [57] and beta function [58]. Finally, the third factor in Eq. (1.3b) is simply the relation derived in Ref. [1], which at the four-loop level reads where the r n are known through order α 4 s [59,60]; the R n depend only on the coefficients of the beta function [21,23] up to an overall normalization, which is given in Ref. [23]; the function mJ MRS (m) = J MRS (m) appears in the definition of the MRS mass [1]; and ∆m (c) contains the contribution from the charm sea quark. Because the nonzero mass of the charmed sea quark cuts off the infrared region that is the origin of factorial growth in the r n [61], we subtract the renormalon with three massless active quarks and lump the charmed loops' contributions into ∆m (c) [62].
The detailed formulas for J MRS (m) and ∆m (c) can be found in Ref. [1]. The crucial aspects of Eq. (3.24) for the fits of the next section is that the renormalon-subtracted perturbative coefficients are small: r n − R n = (−0.1106, −0.0340, 0.0966, 0.0162) for n = (0, 1, 2, 3) and three active flavors. The Borel resummed renormalon is computed from a function with a convergent expansion in 1/α s . (In fact, our implementation of one of the factors in J MRS uses the convergent expansion until it saturates to numerical precision.)

D. Summary formulas
In summary, we fit our data for aM (m h , m x ; {m l , m l , m s }; a) to where F is the fit function and f p4s is in the continuum limit. From the preceding subsections [with free fit parameters in blue (arXiv)]: where y = (aΛ) 2 and w h = Λ HQET /m h,MRS . The HMrASχPT self energy δM Hx depends on f , λ a 2 , λ a 2 ,g π , δ V , and δ A , as well as ∆ * and taste-independent δ Sx . The breved quantities areΛ Ref. [1] am 0h am 0,p4s sim × (3.28) k n x n h × 1 +c 1 α s y +c 2 y 2 +c 3 y 3 .
Thus, there are 61 free parameters, 4 parameters [f , g π , µ 2 G (m b ), and, in m h,MRS /m h , R 0 ] with external priors, and 2 hairpin parameters (δ V and δ A ) from light-meson χPT. ∆ * and δ Sx introduce 2 parameters each that are, however, frozen to reproduce PDG hyperfine and flavor splittings. The total number of fit parameters is 67 (compared with 60 for the decay-constant fit [14]).
In Eq.

IV. EFT FIT TO DETERMINE THE QUARK MASSES
In Sec. III, we have constructed a function with 67 fit parameters that is motivated by EFTs. Here, we use this function to perform a correlated fit to partially-quenched data at five lattice spacings, from a ≈ 0.12 fm to ≈ 0.03 fm, and at several values of the light sea-quark masses. A sixth lattice spacing, a ≈ 0.15 fm, is used to check discretization errors but is not included in the base fit. At the coarsest lattice spacings, we only have data with two different values for valence heavy-quark mass: m h = m c and m h = 0.9m c , where m c is the simulation value of sea charm-quark mass in each ensemble. It is close to but not precisely equal to the physical charm mass m c because of tuning errors. We include data with 0.9m c ≤ m h ≤ 5m c subject to the condition am h < 0.9, which is chosen to avoid large lattice artifacts. For every valence heavy quark, we use several light valence quarks with masses m l m x m s ; on ensembles with the mass of the strange sea quark close to its physical value, m x /m s takes values in a subset of {0.036, 0.1, 0.2, 0.4, 0.6, 1.0} (in several cases the whole set). In the base fit, we obtain the meson masses from fits to two-point correlators with three pseudoscalar states and two opposite-parity states, which we denote "3+2" below. To investigate the error arising from excited state contamination, we also use meson-mass data from (2+1)-state fits.
The values of the bare masses corresponding to the light and strange quarks are taken from combinations of the physical pion and kaon masses, as discussed in Refs. [14] and [50]. Similarly, the physical charmed and bottom quarks are defined so that the D s -and B s -meson masses take their physical values. Because the gauge-field ensembles omit electromagnetism, we need to subtract electromagnetic effects from the experimentally measured masses, which means introducing a specific scheme to do so. We identify M QCD π 0 = M expt π 0 and adjust m l accordingly. Values of the charged and neutral kaon masses in pure QCD, are obtained by subtracting the electromagnetic effects via the quantities (M 2 K 0 ) γ and . Here (M 2 K 0 ) γ is the electromagnetic contribution to the squared mass of the neutral kaon. The quantity parametrizes higher-order corrections to Dashen's theorem: In this paper, we use the most recent values from the MILC Collaboration [63]: Our scheme is the one introduced for u and d quarks in Refs. [64,65]. It defines the isospin limit in the presence of electromagnetism to be the point at which the masses of both the uū and dd pseudoscalar mesons (neglected quark-line-disconnected contributions) are equal to M QCD π 0 . The scheme is extended naturally to the s quark using the fact that mass renormalization for staggered quarks is multiplicative [63]. Numerically, the scheme dependence predominantly affects M 2 K 0 γ and has relatively little influence on the value of .
Using Eqs. (4.2) and (4.3), m s is tuned to obtain As in Ref. [14], we tune m c and m b with the phenomenological formula [66][67][68] M expt amounts to a specific QED renormalization scheme for the heavy-quark mass. Another choice, for example, would be to subtract the leading QED contribution to the self-energy of the heavy quark, which is proportional to e 2 h . Finally, we set (aM Hs /af p4s ) sim = M QCD Hs /f p4s to find the physical am c and am b on each ensemble. In the one-loop χPT result, Eq. (3.2), finite-volume effects enter through the function K 1 . Because the numerical evaluation of those effects is time-consuming, our base fit, as well as various alternative fits that we employ to estimate or check statistical and systematic errors, use the infinite-volume version of K 1 . The finite-volume correction is determined only in a single fit at the end of the analysis. Cross terms between finite-volume and other systematic errors are missed with this approach, but they are negligible.
We use a constrained fitting procedure [69] with priors set as follows. For the main objectives of the analysis, we choose extremely wide priors: 0 ± 6 GeV for both m p4s,MS (2 GeV) and Λ MRS , and (0±1)Λ 2 HQET for µ 2 π . Several other parameters are set from external considerations. As discussed in Sec. III A, the value of µ 2 G (m b ) should be close to the B * -B hyperfine splitting; following Ref. [70], we set the prior distribution of µ 2 G (m b ) to (0.35 ± 0.07) GeV 2 . For the LECs that appear at LO in HMrASχPT and are common for both decay constants and meson masses, we use the same prior constraints as in our work on decay constants [14]: where a 2∆ is related to the differences in squared pion masses, as discussed in Ref. [14]. For the LECs λ 1 and λ 1 , we use wide priors of (0 ± 2) GeV −1 , which are 10 times larger than what can be extracted from the flavor splittings of B or D mesons, namely λ 1 ≈ 0.2 GeV −1 (see, for example, Ref. [71]). Similarly, for the dimensionless LECs λ a 2 and λ a 2 , we use priors of 0 ± 10, which are much wider than the expected size of order 1.
For the overall normalization of R n , which is denoted by N m in Ref. [22] and N in Refs. [21,23], we use R 0 = 0.535 ± 0.010 (4.7) for a theory with three massless active quarks [23]. In the fits reported in this section, we use Eq. (4.7) to provide a prior for R 0 .
Finally, the remaining parameters, which are dimensionless, are given the prior 0 ± 1. The calculation of the MRS mass relies on having a precise estimate for the strong coupling. In this paper, we use α MS (5 GeV; n f = 4) = 0.2128 (25), (4.8) which has been obtained by HPQCD Collaboration [72] for four active flavors. This value corresponds to α MS (m Z ; n f = 5) = 0.11822 (74), whereas the PDG quotes α MS (m Z ; n f = 5) = 0.1181(11) with a somewhat larger uncertainty. The advantage of Eq. (4.8) is that it has been determined on a subset of the same ensembles used here, so it is consistent to use it with our lattice-QCD data. We use the mean value in our base fit, and we introduce an uncertainty associated with α MS by varying its value by 1σ. To run the coupling constant to the scale µ, we use the QCD beta function at five-loop order accuracy [58] and integrate the differential equation numerically. In Sec. V, we comment on how the results would change using the PDG's estimate of the uncertainty in α MS . In general, our data for the meson masses are more precise than the data for scale-setting quantities. We incorporate the latter uncertainties as follows. Let us use af p4s and am p4s to denote the p4s quantities computed from light mesons at each lattice spacing and Σ p4s to denote their covariance matrix. We introduce two fit parameters at each lattice spacing, af p4s,opt and am p4s,opt , to represent optimized values for the p4s quantities under the influence of the heavy-light data. We then employ the so-called penalty trick [73] to take into account the uncertainties in af p4s and am p4s . Thus, we add δχ 2 = af p4s − af p4s,opt am p4s − am p4s,opt (Σ p4s ) −1 af p4s − af p4s,opt am p4s − am p4s,opt (4.9) to our χ 2 function, where the sum is over all lattice spacings. Because data at 5 different lattice spacings enter the base fit, 10 additional parameters are required. The optimized values for the scale setting quantities are then obtained simultaneously in the EFT fit. Given the size of errors in our data, the bias discussed in Ref. [73] is negligible. Altogether we have 384 lattice data points and 77 parameters in our base fit: 67 parameters in the EFT fit function and 10 parameters for optimized values of scale-setting quantities. The fit returns a correlated χ 2 data /dof = 320/307, giving a p value of p = 0.3. Figures 1 and 2 illustrate the base fit at the four (five) lattice spacings for the physical mass (0.2m s ) ensembles and in the continuum limit. The valence light mass m x is tuned to m s : the graphs illustrate a snapshot for heavy-strange meson masses. We plot the heavy-strange meson mass or the difference of the meson mass and theh-antiquark MRS mass versus the continuum limit of the h-quark MRS mass (in Fig. 1) or its reciprocal (in Fig. 2). Data points with open symbols to the right (left) of the dashed vertical line of the corresponding color in Fig. 1 (Fig. 2) are omitted from the fit because they have am h > 0.9. In the continuum extrapolation the masses of sea quarks are set to the physical (correctly tuned) quark masses m l , m s and m c , while at nonzero lattice spacing the masses of the sea quarks take their simulation values.
The width of the fit lines in Figs. 1 and 2 show the statistical error coming from the fit, which is only part of the total statistical error, since it does not include the statistical errors in the inputs of the light quark masses and the lattice scale. Furthermore, the statistical error reported by the fit is sensitive to numerical errors in computing the fit parameters' covariance matrix. For a robust determination of the total statistical error of each output determined directly from the fit to the lattice data. Moreover, the fit function evaluated at zero lattice spacing and physical sea-quark masses yields the meson masses as a function of the valence heavy and light quark masses; see, e.g., Figs. 1 and 2. Figure 3 shows the stability of our final results for MS quark mass ratios m b /m c , m c /m s and m b /m c ; for masses of strange, charm and bottom quarks; and for the HQET matrix element Λ MRS . We test the systematic error in the continuum extrapolation by repeating the fit after either adding in the coarsest (a ≈ 0.15 fm) ensembles or omitting the finest (a ≈ 0.03 fm) ensemble. These changes are shown in Fig. 3 and seen to have no significant effect, so we consider these tests to be cross-checks. The meson-mass data in our base fit are obtained from the (3+2)-state fits to two-point correlators. To investigate the error arising from excited state contamination, we repeat the EFT fit with meson-mass data from the (2+1)-state fits to two-point correlators. As seen in Fig. 3, the effects from this change are small too. Because we have no other handle on systematic errors due to excited states, we take the difference between the results from the two types of correlator fits as an estimate of this uncertainty.
We now turn to effects from truncating perturbative QCD in the relation between quarkmass definitions and the beta function. As explained with Eq. (1.3), the MRS mass connects the MS mass of the h quark to the heavy-light-meson mass H x . By design, the fit yields m p4s,MS (2 GeV), and we use the continuum limit of am h /am p4s to convert to m h,MS (2 GeV). We then use Eqs. (3.22) and (3.23) to calculate m h and Eq. (3.24) to calculate m h,MRS . The beta function and quark-mass anomalous dimension are known at five loops [57,58], and the pole mass at four loops [59,60].
To monitor the errors from truncating perturbative QCD, we rerun the analysis with fewer orders in Eqs. (3.23) and (3.24) and in the beta function without, however, changing C cm , the Wilson coefficient for µ 2 G , Eq. (3.7). Figure 4 shows the stability of our results for MS quark mass ratios m b /m c , m c /m s and m b /m c ; for masses of strange, charm and bottom quarks; and for the HQET matrix element Λ MRS , as the order of perturbation theory is increased. In Fig. 4, we denote by O(α n s ) a fit that includes n orders beyond the leading terms in Eqs. sensitive to the truncations in the perturbative-QCD relations; essentially, these ratios are the continuum limit of the corresponding bare masses. For the quark masses and the HQET matrix elements, one finds good convergence, within the statistical errors, as the order of α s in the perturbative expressions is increased. Based on this observation, we do not introduce any additional systematic error associated with truncation in perturbative-QCD results. Note that the truncation effects are negligible because the renormalon-subtracted perturbative coefficients in the MRS mass are all very small. If one employs the RS mass [22], for instance, the truncation error for the bottom quark mass would be about 10 to 20 MeV, depending on the details.
Our data prefer an overall coefficient of the one-loop HMrASχPT contribution, δM Hx , namely g 2 π /f 2 , well below the prior width of the product of g 2 π and 1/f 2 in Eqs. (4.6a) and (4.6b). In our base fit, the posterior for this product is g 2 π /f 2 = 3.8 GeV −2 (for D systems) while the prior value is (14 ± 7) GeV −2 . To investigate the effects of treating g 2 π and 1/f 2 as free parameters, we consider two alternative fits. First, we fix g π = 0.45, one sigma below its nominal value, and 1/f 2 = 1/f 2 K . Second, we fix g π = 0, which is equivalent to fitting to a polynomial in the quark masses. In Fig. 5, we label the first of these "g π = 0.45" and the second "g π = 0". As one can see, the quark-mass ratios, quark masses themselves, and the HQET matrix element Λ MRS do not change significantly under these variations. Consequently, we do not introduce any additional systematic error associated with our treatment of g π and 1/f 2 .
Our full error budget for quark-mass ratios, quark masses, and the HQET matrix element Λ MRS is given in Table I. The row labeled "statistics and EFT fit" lists the uncertainty reported by the Bayesian fit, which incorporates associated systematic effects of extrapolation. There are further systematic effects not captured in the EFT fit. The excited-state contamination in two-point correlator fits is explained above. Our method of estimating the systematic error associated with the tuned quark masses and scale setting quantities is similar to that in Ref. [14]. We correct for (exponentially small) finite-volume effects using the finite-volume version of the NLO χPT for the heavy-light mesons, and using NLO or NNLO χPT for the light-quark and scale-setting inputs following Ref. [50]. Residual finite-volume effects from higher orders in χPT are estimated, as in Ref. [14], as 0.3 times the calculated finite-volume correction. The nonequilibration of topological charge in our finest ensembles causes small finite-volume effects that are not exponentially suppressed [74]. Although this error is negligible for masses of heavy-light mesons, even with our high statistics, we include the shift expected from Ref. [74] as a systematic error. Note that, despite the fact that we take the full topological shift as the associated error, these errors all round to zero at the precision shown in Table I. The uncertainties stemming from the omission of electromagnetism are discussed in detail below. Finally, our results have uncertainties from the parametric inputs α s , given in Eq. (4.8), and f π,PDG = 130.50(13) MeV [52]. Table II shows a breakdown of the uncertainties from matching a pure-QCD calculation such as this to QCD+QED. Briefly, "K + -K 0 splitting" is the uncertainty in connecting the K + -K 0 splitting to that of π + -π 0 , stemming from , and "K 0 mass" refers to the uncertainty in the electromagnetic contribution to the neutral kaon mass, M 2 K 0 γ . These two effects are negligible compared with the other sources of uncertainty. In this work, we choose a specific scheme [63,64] for the electromagnetic contribution to the neutral kaon masses; other works, for example the FLAG report [75], choose other schemes. Changing 17%. When using these ratios and the strange-quark mass in a setting that ignores the subtleties of the QED scheme, one may wish to incorporate an additional uncertainty of ±0.17%. Another uncertainty comes from the estimates of the electromagnetic correction to the heavy-light meson mass, described above with Eq. (4.5). It is denoted "H x mass" in Table II. For the associated error, we take the difference between results obtained with and without the electromagnetic shift. Results for heavy-quark masses depend on the chosen QED quarkmass scheme. As discussed above, we do not subtract any part of the QED self-energy.

V. RESULTS
In this section, we collect the results that stem from the EFT fits described in the previous two sections. These fall into four categories: quark masses themselves and their ratios, HQET matrix elements, flavor splittings in the D and B systems, and LECs of HMχPT. We emphasize again that our final results for quark masses depend on our prescription for calculating QCD-only meson masses; cf. the discussions around Eq. (4.5) and about Table II. where m l is again the average of the up-and down-quark masses. To obtain these results, we take the large side of the asymmetric uncertainties reported in Ref. [14], namely, m s /m l = 27.178 (47)  The quark masses given above are for four active flavors. The mass of the bottom quark with five active flavors can be calculated from [76] where n l = n f − 1 and µ = m b (n f ) . Setting n f = 5, we obtain (5.15) or, adding all errors in quadrature, m b = 4195 (14) MeV. The five-flavor mass can be run from m to higher scales using the five-loop anomalous dimension [57] and beta function [58] with n f = 5. For completeness, we run the b mass to 10 GeV, finding m b,MS (10 GeV; n f = 5) = 3665(11) stat (1) syst (1) αs (1) f π,PDG MeV, (5.16) in which the α s uncertainty has become very small. This value for µ 2 G (m b ) cannot be considered as a pure lattice-QCD determination because, as discussed in Sec. III A, the prior for µ 2 G (m b ) ∼ 0.35(7) GeV 2 comes from the B-meson hyperfine splitting. The definition of µ 2 π used here still has a renormalon ambiguity of order Λ 2 QCD , although it is expected to be small [1,77]. In any case, the result in Eq. (5.22) cannot be directly compared with results in the "kinetic" scheme [78,79], where µ 2 π ≈ µ 2 G is expected [80] and roughly holds [12,81]. We checked whether our χ 2 function could be consistent with such an outcome by starting the fit at µ 2 π = 0.35 MeV 2 , but found the same minimum as in Eq. (5.22). We also have tried changing the prior for µ 2 π from (0 ± 0.36) GeV 2 to (0.35 ± 0.36) GeV 2 , in which case χ 2 is minimized for µ 2 π = 0.09(16) GeV 2 and µ 2 G (m b ) = 0.39(1) GeV 2 , where the errors are statistical only here. To compare Eq. (5.21) with the RS scheme at a given factorization scale ν f , one can use [1] Λ RS (ν f ) = Λ MRS + J MRS (ν f ), (5.24) with the function J MRS given in Eq. (2.37) of Ref. [1]. Setting ν f = 1 GeV, we find The uncertainty associated with α s is larger here than for Λ MRS , because J MRS (ν f ) in Eq. (5.24) depends on α s (ν f ). Our result for Λ RS (1 GeV) agrees with Λ RS (1 GeV) = 659 MeV (no error quoted) [22] and 623 MeV (after rough conversion of a result in the "RS " scheme) [82], which are obtained from the B-meson mass and the RS mass for the bottom quark. For future phenomenological studies, Table III in the Appendix provides the correlation matrix of the MRS masses of the charm and bottom quarks with the HQET matrix elements Λ MRS , µ 2 π and µ 2 G (m b ).

C. Flavor splittings
We use the D s -and B s -meson masses as experimental input to set the c-and b-quark masses. Comparing the output of the fit at m x = m d with m x = m s , we obtain the flavor splittings In these combinations of meson masses, the leading-order electromagnetic contributions cancel. The last uncertainty here stems from the significant changes found in the alternate fits with g π fixed to 0.45 or to 0 (the polynomial fit).
In a similar vein, we can set the quark masses to m x = m l = m s = 0 to obtain the SU(3) chiral limit of charmed and b-flavored mesons, or set m x = m l = 0 and leave m s = m s to obtain the SU (2)

D. Low-energy constants in HMχPT
Reference [74] uses the LECs λ 1 and λ 1 obtained in this work. In particular, the values used are those for D mesons, which come from the simple, polynomial analysis without chiral expressions, i.e., g π = 0:λ where the errors are statistical only, which suffices for Ref. [74]. Here, the breve is a reminder that finite-mass corrections to the LECs in the HMχPT Lagrangian are included. From the experimental data for the flavor splittings of D mesons, one findsλ 1 ≈ 0.2 GeV −1 [71].

VI. SUMMARY, COMPARISONS, AND OUTLOOK
The results presented in Sec. V show that the new HQET-based method, developed here and in Ref. [1], is both qualitatively and quantitatively successful. The qualitative success relies on the clean separation of scales provided by HQET with the MRS definition of the heavy-quark mass, while the quantitative success relies on the high statistics of the MILC Collaboration's HISQ ensembles [25][26][27], all 24 of which have been employed here. Also relevant to the success of the method is the availability of the order-α 5 s perturbation theory for the running of the quark mass [57] and strong coupling [58], and the order-α 4 s coefficient linking the MS mass to the pole mass and, hence, the MRS mass [59,60]. These features are not (yet) shared by other determinations of quark masses using lattice QCD. Although the HQET method separates the heavy-quark scale from the QCD scale, mass ratios determined in the course of this work and Ref. [14] yield results for all quarks except the top quark.
Our results for heavy-quark masses m c and m b are compared with other results in the literature in Fig. 6. Both panels show the most recent lattice-QCD calculation with a complete error budget from each combination of method and collaboration. For nonlattice calculations, we also show the most recent result from each method and-or collaboration, but include only those with perturbative-QCD accuracy of order-α 3 s matching and, if needed, order-α 4 s running. As noted in Sec. V, the parametric uncertainty in α s is one of our largest uncertainties, but, thanks to the MRS mass, higher-order perturbative corrections are likely to be negligible compared with this and our statistical uncertainty; cf. Fig. 4 and Eq. (5.19).
For m c , the overall agreement is very good, and our result's uncertainty is about the same as those from charmonium correlators and (continuum) perturbative QCD, using either lattice [72] or experimental [90,94] data as input. (References [90] and [94] differ in the moments used.) The difference between our result for m c,MS (3 GeV) and the recent update from Chetyrkin et al. [90] is 0.9σ. For m b , the overall agreement is good. The difference between our result and those of Narison [93], Bodenstein et al. [102], Chetyrkin et al. [103], and Penin and Zerf [101] is 1.3σ, 1.6σ, 1.6σ, and 1.7σ, respectively. Such discrepancies among 19 independent results, especially given the importance of systematic uncertainties in all determinations, should not be seen as alarming.
It is noteworthy that for m c = m c,MS (m c,MS ) the result of Ref. [90] is more precise than ours, while for m c,MS (3 GeV) ours is more precise. In both cases, the error bar runs as dictated by the quark-mass anomalous dimension and beta function. In addition, the orderα s coefficient is proportional to [ln(3 GeV/µ) + c]. For the relation between m and m MRS , c > 0, so the first-order α s error vanishes for some µ > 3 GeV. On the other hand, for the  [72]; ETM 14 [84]; Maezawa and Petreczky 16 [85]; RBC/UKQCD 14 [105]; BMW 10 [106]; HPQCD 10 [88]; and MILC 09 [107].
relation between m and moments of the charmonium correlator, c < 0, so the first-order α s error vanishes for some µ < 3 GeV. We therefore also provide light, strange, and charm masses at 3 GeV: In contexts beyond the Standard Model, one needs the masses-that is the Yukawa coupling to the Higgs field-at scales of 100 GeV or higher. Table IV in the Appendix provides the correlation matrix for our charm-quark mass at 3 GeV and quark-mass ratios.
Our results for light-quark masses are compared with other results from lattice QCD in Fig. 7. As above, both panels show the most recent lattice-QCD calculations with a complete error budget from each combination of method and collaboration. As can be seen from the plots, and similar comparisons of m u and m d , ours are the most precise results to date. Here the precision stems from very precise quark-mass ratios from the pseudoscalar meson spectrum, together with the overall scale of quark masses from the EFT fit. Consequently, the results inherit an uncertainty due to α s , which is largest except in the cases of m d,MS (2 GeV) and m u,MS (2 GeV), which have larger statistical and electromagnetic systematic uncertainties from m u /m d .
As compelling as these results are, they could be improved in several ways. First, because the EFT fit controls systematics, the statistical error (after propagation through the fit) is often the second-largest source of uncertainty, so, as usual, having more data would reduce the error. The additional data need not be more precise per se: the right panel of Fig. 2 suggests that finer lattice spacings will be needed. Second, because the other dominant uncertainty is the parametric error of α s , it would be interesting to carry out a simultaneous determination of α s and the quark masses, for example in a combined analysis of heavylight meson masses and quarkonium correlators. Such an analysis would output m c , m b , and α s with their correlations, which would be very convenient for determining Higgs-boson branching ratio in the Standard Model and extensions thereof. Third, QCD+QED simulations would eliminate the scheme dependence arising from the matching of QCD+QED to pure QCD. Finally, the ideal determination of the matrix elements µ 2 π and µ 2 G , and analogous quantities that enter at order 1/m 2 Q and higher, would require computing heavy-light vector mesons on the lattice, in addition to the pseudoscalar mesons studied here. In particular, this would make possible a pure lattice result for µ 2 G , without making use of the experimental information on the B-meson hyperfine splitting.