# Determination of \overline{m}_b/\overline{m}_c and \overline{m}_b from n_f=4 lattice QCD+QED

D. Hatton, C. T. H. Davies, J. Koponen, G. P. Lepage, A. T. Lytle

DDetermination of m b /m c and m b from n f = 4 lattice QCD + QED D. Hatton and C. T. H. Davies ∗ SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK

J. Koponen

Helmholtz Institute Mainz, Johannes Gutenberg-Universit¨at Mainz, 55099 Mainz, Germany

G. P. Lepage † Laboratory for Elementary-Particle Physics, Cornell University, Ithaca, New York 14853, USA

A. T. Lytle

Department of Physics, University of Illinois, Urbana, IL 61801, USA (HPQCD collaboration) ‡ (Dated: February 22, 2021)We extend HPQCD’s earlier n f = 4 lattice-QCD analysis of the ratio of MS masses of the b and c quark toinclude results from ﬁner lattices (down to 0.03 fm) and a new calculation of QED contributions to the mass ratio.We ﬁnd that m b ( µ ) /m c ( µ ) = 4 . at renormalization scale µ = 3 GeV. This result is nonperturbative.Combining it with HPQCD’s recent lattice QCD + QED determination of m c (3 GeV) gives a new value for the b -quark mass: m b (3 GeV) = 4 . GeV. The b -mass corresponds to m b ( m b , n f = 5) = 4 . GeV.These results are the ﬁrst based on simulations that include QED.

PACS numbers: 11.15.Ha,12.38.Aw,12.38.Gc

I. INTRODUCTION

Accurate masses for heavy quarks are important for QCD-phenomenology generally, but they will be particularly im-portant for high-precision searches for new physics in Higgsdecays [1]. In this paper we present a new result for the ra-tio m b /m c of the MS masses of the b and c quarks. Ouranalysis of the mass ratio is completely nonperturbative. Thisis in contrast to lattice-QCD determinations of the separatequark masses, which need QCD perturbation theory to relate MS masses to lattice quantities. Thus the (nonperturbative)mass ratio provides a nontrivial check on (perturbative) de-terminations of the separate masses. The ratio can also becombined with recent accurate determinations of the c -quarkmass to obtain new results for the b -quark mass.Lattice simulations of b quarks are complicated by thequark’s large mass, which leads to large lattice-spacing errorswhen the b quarks are described by the Dirac equation (as op-posed to, say, NRQCD [2]). We address this problem by us-ing a Highly Improved Staggered-Quark discretization of theDirac equation (HISQ [3]) that is also highly efﬁcient, mak-ing simulations at very small lattice spacings feasible. Ourprevious analysis of the mass ratio [4] used lattices with spac-ings down to 0.06 fm but still required an extrapolation in thequark mass to reach m b . Here we reduce the lattice spacingto 0.03 fm, where am b ≈ . , which allows us to simulateat the b mass. Lattice spacing errors at m b are less than 1%on our ﬁnest lattice, and we are able to remove most of that ∗ [email protected] † [email protected] ‡ error by extrapolating from results covering a range of latticespacings and heavy-quark masses m h .Our new result is accurate to about 0.25%, so it becomesimportant to include QED effects. We recently analyzed theQED contributions to the c quark’s mass [5]. Here we adaptthe methods from our earlier paper to provide the ﬁrst resultsfor QED contributions to m b /m c and m b . Here and in ourearlier paper we use the quenched QED approximation, whichomits contributions from photons coupling to sea quarks. Thequenched approximation should capture the bulk of the QEDcorrection in mesons whose valence quarks are both heavy;contributions from sea quarks are expected to be an order ofmagnitude smaller [5].In Section II we describe our general strategy and the latticeQCD simulations we employed. In Section III, we extract avalue for m b /m c using results from simulations without QED.We then add QED effects in Section IV. We summarize ourresults for the mass ratio in Section V and combine them withHPQCD’s recent c -quark mass to obtain a new result for the b -quark mass. II. LATTICE QCD SIMULATIONS

We use gluon conﬁguration sets generated on a variety oflattices (by the MILC collaboration [8]), with n f = 4 ﬂavorsof HISQ sea quark and lattice spacings ranging from 0.09 fmto 0.03 fm. These sets are described in Table I. The u and d quark masses are set equal to m (cid:96) ≡ ( m u + m d ) / ; cor-rections to this approximation are quadratic in the light-quarkmasses for our analysis, and so are negligible. We include re-sults where the light-quark masses in the sea are tuned closeto their physical value, but we also include results with muchlarger light-quark masses in the sea. Results from these last a r X i v : . [ h e p - l a t ] F e b TABLE I. Gluon conﬁguration sets used in this paper. Sets are grouped by approximate lattice spacing, with lattice spacings of 0.09 fm (Sets 1and 2), 0.06 fm (Sets 3 and 4), 0.045 fm (Set 5), and 0.03 fm (Set 6). Lattice spacings are determined from the values shown for the Wilsonﬂow parameter w /a [6] where w = 0 . fm [7]. The sea quark masses are given in lattice units for u/d quarks ( am (cid:96) ), s quarks ( am s ),and c quarks ( am c ). Tuned values for the lattice c masses are also given (in GeV). These masses are adjusted to give correct masses for eitherthe η c or J/ψ mesons (Eq. (3)). The c masses are tuned using slope dm c /dm cc , which is the same (within errors) for the η c and J/ψ . Thespatial and temporal sizes of the lattices, L and T , are also listed.Set w /a am sea (cid:96) am sea s am sea c ˜ m tuned c ( η c ) ˜ m tuned c ( J/ψ ) d ˜ m c /dm cc L/a T /a simulations are unphysical but are easily corrected [4]. In-creasing the light-quark mass in the sea signiﬁcantly reducesthe cost of our analysis at small lattice spacings (becausesmaller lattice volumes are used).Ignoring QED for the moment, the ratio of the b and c MS masses equals the ratio of the corresponding bare quarkmasses used in the lattice Lagrangian, up to corrections thatvanish in the continuum limit [9]: m b ( µ ) m c ( µ ) = m tuned b m tuned c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) latt + O ( α s ( π/a ) a ) (1)where the bare masses are tuned so that the QCD simulationsreproduce the experimental results for meson masses. This re-lationship between the MS and lattice quark masses is nonper-turbative and independent of the MS renormalization scale µ .Pseudoscalar and vector meson masses from our simula-tions are listed in Table II for a variety of (valence) heavy-quark masses, ranging on the ﬁnest lattices (Sets 5 and 6)approximately from the c mass to the b mass. The analysismethods for extracting these masses (and most of the results)come from [5, 10].The quark masses am h in Table II are what is used in theHISQ Lagrangian. The a ˜ m h masses are corrected to removetree-level ( am h ) n errors (in the pole mass) through order n =10 [3, 11]: a ˜ m h ≡ am h (cid:32) − (cid:16) am h (cid:17) + 14722240 (cid:16) am h (cid:17) + 456448537600 (cid:16) am h (cid:17) − (cid:16) am h (cid:17) (cid:33) . (2)We write the expansion as powers of am h / becausethis makes the leading coefﬁcients roughly the same size(about / ). The correction is − % at am h = 0 . , whichis the largest value we use.We give results in Table I for the tuned bare c mass for eachof the conﬁgurations. In each case we adjust the c mass so asto reproduce the continuum value for either the η c mass or the J/ψ mass: m cont η c = 2 . m cont J/ψ = 3 . . (3)Here we have subtracted 7.3(1.2) MeV from the experimentalvalue for m η c [12] to account for the fact that we are not in-cluding contributions from cc annihilation in our simulations;this correction is determined in [5]. The analogous correc-tion to the J/ψ mass is negligible, but we have subtracted0.7(2) MeV from the mass to account for cc annihilation to aphoton; this correction is estimated perturbatively in [5]. Weextrapolate the c masses to their correct values using a ˜ m tuned c = a ˜ m c − (cid:0) am lat cc − am cont cc (cid:1) d ˜ m c dm cc (4)where m cc is either the η c or the J/ψ mass and the slopes (seeTable I) are estimated from splines ﬁt to the entries in Table II.

III. m b /m c WITHOUT QED

We could tune the lattice b mass the same way we tuned m c ,but we have only a few simulation results near the b and thesehave signiﬁcant ( am h ) n errors. Instead we will use the datain Table II to deﬁne functions that relate the ratio of quarkmasses to the pseudoscalar ( P ) or vector ( V ) masses: m Phh = f Phh ( m h /m c ) m Vhh = f Vhh ( m h /m c ) , (5)where f Phh (1) ≡ m cont η c f Vhh (1) ≡ m cont J/ψ . (6)Given these functions, we then obtain two estimates for m b /m c by solving each of the equations f Phh ( m b /m c ) = m cont η b f Vhh ( m b /m c ) = m contΥ (7) TABLE II. Lattice QCD results for the ground-state pseudoscalar andvector hh mesons in lattice units: am Phh and am Vhh , respectively.Results are given for each conﬁguration sets (Table I) and a variety ofbare quark masses am h and corrected masses a ˜ m h (in lattice units).The uncertainties in the meson masses are negligible compared withother errors and have no impact on the ﬁnal results.Set am h a ˜ m h am Phh am Vhh for m b /m c , where m cont η b = 9 . m contΥ = 9 . . (8)The two mass ratios should agree. Here we account for the bb annihilation contribution to the η b mass by adding an ex-tra error of ± MeV to the experimental result [12]; this esti-mate is based on NRQCD perturbation theory and the meson’swidth [3]. The analogous contribution to the Υ mass is negli-gible, as is Υ annihilation via a photon.In what follows, we ﬁrst describe our lattice-QCD analysisof f Phh and f Vhh , and then discuss the results.

A. Analysis

We determine the f hh functions (Eq. (5)) by ﬁtting the me-son masses am hh from Table II to functions of the followingform, am hh × (1 ± σ u ) = af hh ( r )+ am hh (cid:0) δ a + δ sea uds + δ sea c (cid:1) , (9) where r ≡ a ˜ m h a ˜ m tuned c ξ m ( ˜ m tuned c , δm sea uds ) ξ m ( ˜ m h , δm sea uds ) (10)is the effective ratio of quark masses m h /m c . Pseudoscalarand vector mesons are ﬁt separately, to determine each of f Phh and f Vhh . We describe each element of the ﬁt function in turn:• We increase the fractional error on each value of am hh from Table II to ± σ u . The am hh errors listed in thetable are very small. It is impossible to ﬁt the almostsix signiﬁcant digits in these data with a model as sim-ple as we use here. So we increase the fractional erroron each value to σ u , which is then a measure of the partof the variation in the data that is unexplained by ourmodel. The σ u errors are uncorrelated from one am hh to another. We use the same value for σ u for every datapoint and adjust its size to maximize the Bayes Factorfrom the ﬁt [13]. For the parameters and model usedhere, we ﬁnd that σ u = 0 . , (11)which means that our model explains the individual datapoints to within ± . %. A simpler model wouldhave a larger σ u : for example, σ u more than doubles ifthe ξ m factors in ratio r are dropped. Note that the sta-tistical errors listed in Table II can be neglected when σ u is included.• We parameterize the f hh functions as splines [14] with6 knots evenly spaced from r = 1 to . ( ≈ m b /m c ),inclusive. The ﬁt parameters are the function values atthe knots. These functions are linear up to correctionsof order v /c ∼ . – 0.3, where v is the typical ve-locity of the heavy quarks in the meson. Therefore weuse the following priors for the values at the knots with r > : f hh ( r knot ) = 1 . × (cid:16) m cc + r knot − . m bb − m cc ) (cid:17) , (12)where m cc and m bb are the continuum masses ofthe pseudoscalar/vector mesons composed of c and b quarks, respectively (Eqs. (3) and (8)). At r = 1 ,we require f hh ( r = 1) = m cc . (13)We choose 6 knots to maximize the Bayes Factor fromthe ﬁt. Results obtained using 5 or 7 knots agree wellwith those from 6 knots. Doubling the width of the pri-ors has no effect on our results.• The ξ m factors in the mass ratio r rescale the quarkmasses to correct for detuned values of the light seaquarks. From [4], ξ m ( m h , δm sea uds ) = 1 + g m ( m h /m c ) ζ δm sea uds m s (14)where δm sea uds ≡ (cid:88) q = u,d,s (cid:0) m sea q − m tuned q (cid:1) (15)is the difference between the masses used in the sim-ulation and their tuned values. The tuned masses aredetermined from the tuned c -quark mass using resultsfor m c /m s and m s /m (cid:96) from [15]. The priors for the ﬁtparameters are g m = 0 . ζ = 0 . , (16)which come from ﬁts described in [4].• The largest simulation errors are from the discretiza-tion. These are suppressed by α s ( π/a ) in order a be-cause we are using the HISQ formalism [3]. Beyondthis order they are suppressed either by α s ( π/a ) or by v /c , since we have removed the tree-level a n errorsin the quark masses using Eq. (2). The ﬁt can’t dis-tinguish easily between α s and v /c since both arearound 0.2 for our data, so we include only an α s cor-rection, modeled after Eq. (2): δ a ≡ α s ( π/a ) (cid:88) n =1 f na ( r ) (cid:32) a ˜ m h (cid:33) n . (17)Here functions f na ( r ) are 6-knot splines with priors atthe knots (same locations as above) of f na ( r knot ) = 0 . . (18)Terms beyond n = 3 have no effect on the ﬁt results.The splines allow for m h dependence in the a n correc-tions.• We include a corrections to ξ m since δm sea uds is largefor some of our conﬁguration sets: δ sea uds = α s ( π/a ) f sea uds ( r ) δm sea uds m s (cid:32) a ˜ m h (cid:33) . (19)where function f sea uds ( r ) is again a 6-knot spline, nowwith priors at the knots of f sea uds ( r knot ) = 0 . . (20)where the width is chosen to be somewhat larger thansuggested by g m above.• We also include a correction to ξ m from detuned c -quark masses in the sea. This correction should be smallbecause of heavy-quark decoupling [4] — the momen-tum transfers in the heavy-quark mesons are too smallto produce cc pairs efﬁciently. We include the correc-tion δ sea c = f sea c ( r ) δm sea c m c (21) m hh (GeV)12345 m h / m c m η b m η c m Υ m J/ψ

FIG. 1. Lattice QCD (without QED) results for m h /m c plottedversus the hh meson masses. Values for m h /m c are corrected asin Eq. (10). The lines, which vary in thickness, show results fromthe best-ﬁt values for the functions f P/Vhh ; the line thickness showsthe σ uncertainty in these functions. These functions can be re-constructed from information in the Appendix. Separate results areshown using the pseudoscalar masses m Phh (top line, squares) andthe vector masses m Phh (bottom line, circles). Differ colors indicatedifferent conﬁguration sets, with Sets 6 (brown) and 5 (purple) hav-ing the largest masses, followed by Sets 4 (red), 3 (green), 2 (orange)and 1 (blue), in that order. where δm c ≡ m c − m tuned c , and f sea c ( r ) is a 6-knotspline with f sea c ( r knot ) = 0 . . (22)We chose the width to maximize the Bayes Factor fromthe ﬁt.The ﬁt parameters are the values of the coefﬁcient functions(splines) at the knots, together with g m and ζ from the ξ m fac-tors (Eq. (14)). We use the lsqfit Python module to do theﬁts [16, 17].

B. Results

The functions f P/Vhh obtained from the ﬁts described inthe previous section (see Appendix) are plotted in Fig. 1, to-gether with the data from Table II. Fig. 2 shows that the model(Eq. (9)), with best-ﬁt values for the ﬁt parameters in the cor-rections (on the right-hand side), reproduces the data withinerrors. The difference between the lattice results with andwithout corrections is − . % for the highest quark masson the ﬁnest lattice (Set 6).We can use functions f P/Vhh to extract values for the ratio of χ is less useful as a measure of goodness-of-ﬁt here because we adjust σ u to give a good ﬁt. χ per degree of freedom was 0.9 for the pseudoscalardata (26 points) and 0.8 for the vector data (25 points). . . . . a ˜ m h − . . . ( d a t a − m o d e l ) / d a t a FIG. 2. Relative difference between the data for am hh from Table IIand the model in Eq. (9) with best-ﬁt values for the ﬁt parameters.Results are shown for both pseudoscalar (squares, offset right) andvector (circles, offset left) mesons, where data points for the twomesons are offset slightly in opposite directions to improve visibility.TABLE III. Contributions to the total ( σ ) error in m b /m c fromQCD simulations (without QED), as a percentage of the mean value.Results are given for determinations using pseudoscalar mesons( m Phh ) and vector mesons ( m Vhh ), and for the weighted average ofthese results. The dominant errors come from the extrapolation tozero lattice spacing and from uncertainties in the lattice spacing. Ad-ditional errors are from residual uncertainties taken in the ﬁt data( σ u ), uncertainties in ξ m used to correct for unphysical sea-quarkmasses, uncertainties in the η c and η b masses, and tuning uncertain-ties in the sea-quark masses. The error budgets are the same whenQED is included aside from an additional uncertainty of 0.03% asso-ciated with the QED corrections. m b /m c [ m Phh ] m b /m c [ m Vhh ] m b /m c [avg]( am h ) → w , w /a σ u g m , ζ m cc m bb ( am h ) δm sea uds → δm sea c → d ˜ m c /dm cc MS masses by solving Eqs. (7). We obtain m b /m c = (cid:40) . from the am Phh . from the am Vhh , (23)independent of renormalization scale. The two estimatesagree to within . %. The weighted average, taking ac-count of correlations, is m b /m c = 4 . . (24)We tabulate the leading uncertainties in our two results inTable III. The error budgets are similar for the two mesons,and are dominated by uncertainties associated with discretiza-tion errors and the lattice spacings. TABLE IV. Ratio R of m hh with QED corrections to m hh withoutQED corrections, evaluated at the same quark mass m h . Resultsare shown for ground-state pseudoscalar and vector mesons analyzedon two conﬁguration sets. The quark’s QED charge is Q times theproton’s charge; results for Q = 2 / can be converted to Q = 1 / by replacing R with R − / .Set a ˜ m h Q R ( m Phh , Q ) R ( m Vhh , Q ) / / / / / / / / / IV. ADDING QED

Adding QED complicates the analysis of m b /m c becausethe quarks have different QED charges and therefore differentmass anomalous dimensions. Thus the nonperturbative rela-tion in Eq. (1) is only true up to O ( α QED ) corrections. Wedeal with this complication by introducing QED through tworatios R : m b ( µ ) m c ( µ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCDQED = R (cid:0) m b /m c , Q c,b = 0 → (cid:1) R (cid:0) m c ( µ ) , Q c = → (cid:1) × m b m c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCD . (25)Here R (cid:0) m c ( µ ) , Q c = → (cid:1) ≡ m c ( µ ) with Q c = m c ( µ ) with Q c = (26)is the ratio of the c mass in a theory with c -quark charge tothe mass in a theory with c -quark charge . Similarly, R (cid:0) m b /m c , Q c,b = 0 → (cid:1) ≡ m b /m c with Q c,b = m b /m c with Q c,b = 0 , (27)where the c and b charges are equal ( Q c = Q b ) in each case(and so the ratio is µ independent). In every case, the quarkmasses are tuned to reproduce the continuum meson massesin Eqs. (3) and (8). Either the pseudoscalar or vector mesonscan be used; they give the same results to within the precisionneeded here. Only the ﬁrst of the R factors (Eq. (26)) dependson the MS renormalization scale µ ; we take µ = 3 GeV, fol-lowing [5]. We approximate full QED by quenched

QED,where only the valence quarks carry electric charge. This isexpected to be the dominant contribution in O ( α QED ) and ismuch less costly to analyze. The techniques we use for intro-ducing QED into simulations are standard and are describedin [5, 10].The R factor for the c mass (Eq. (26)) is the most importantand can be inferred from our earlier result [5]: R (cid:0) m c (3 GeV) , Q c = 0 → (cid:1) = 0 . . (28) m h /m c R ( m h / m c , Q = / ) FIG. 3. Ratio R ( m h /m c , Q = ) is plotted versus m h /m c . It isthe ratio of m h /m c computed with QED charge Q = to the resultwithout QED ( Q = 0 ), where the quark masses are tuned to give thesame results for m Phh (bottom, red) or m Vhh (top, blue). Results areshown from conﬁguration Sets 1 (squares) and 3 (circles). Errors aresmaller than the plot symbols. The blue and red shaded areas showthe ± σ ﬁts to the data (Eq. (32)). R is quadratic in Q c to better than 0.01% so the QED correc-tion ( R − ) required to go from charge to charge is threequarters that required to go from 0 to : R (cid:0) m c (3 GeV) , Q c = → (cid:1) = 0 . . (29)The other R factor, for m b /m c , is expected to be muchcloser to one for two reasons: the QED corrections for the b and c masses are similar and tend to cancel in the ratio;and the charges Q c,b = are smaller (and the QED effectis quadratic in the charge). To estimate the effect, we cal-culated the ratio R of meson masses m hh with and without Q c,b = QED, holding the quark masses constant, for two ofour conﬁguration sets; our results are in Table IV. This quan-tity can be related to the R factor for m b /m c by re-expressingthe R -factor in terms of lattice masses, using Eq. (1) (since Q c = Q b ), and writing it as R ( m b /m c , Q ) = 1 + δ ˜ m Qb ˜ m b − δ ˜ m Qc ˜ m c + O (cid:0) δ ˜ m (cid:1) , (30)where δ ˜ m Qc,b are the quark mass shifts needed to hold the me-son masses constant when QED is added to the simulation.The mass shifts can be calculated for different heavy-quarkmasses m h from the R factors in Table IV: δ ˜ m Qh = (cid:0) − R ( m hh , Q ) (cid:1) m hh d ˜ m h dm hh . (31)Here the derivative is estimated for each conﬁguration by ﬁt-ting a cubic spline to the a ˜ m h values in Table II as a functionof the corresponding am hh values.Values for R ( m h /m c , Q = ) are plotted versus m h /m c in Fig. 3 for both pseudoscalar (below) and vector (above)mesons from the two conﬁguration sets. We ﬁt these data to asimple function suggested by QED perturbation theory: R = 1 + (cid:88) i =1 c i log i ( ˜ m h / ˜ m c ) + (cid:88) j =1 d j (cid:0) a ˜ m h / (cid:1) j (32) . . m b ( m b , n f = ) (GeV) HPQCD ’21 (HISQ)Fermilab/MILC/TUMQCD ’18Gambino et al ’17ETM ’16HPQCD ’14 (NRQCD b )HPQCD ’14 (HISQ) QCDQCD + QED

FIG. 4. Values for the MS mass of the b quark from lattice QCD sim-ulations with n f = 2 + 1 + 1 ﬂavors of sea quark. Results are shownfrom: HPQCD ’21 (this paper), Fermilab/MILC/TUMQCD [15],Gambino et al [18], ETM [19], HPQCD ’14 (NRQCD) [20], andHPQCD ’14 (HISQ) [4]. The gray band corresponds to the top re-sult (HPQCD ’21), the only one from simulations that include QED. with priors c i = 0 . and d j = 0 . . Extrapolating tothe b mass gives: R ( m b /m c , Q = ) = (cid:40) . from m Phh . from m Vhh . (33)The two results agree with each other, but the corrections aretoo small to affect our ﬁnal results signiﬁcantly. We use thelarger error in the error budgets for our ﬁnal result.Including both R factors, we arrive at new results for thequark mass ratio at µ = 3 GeV that include (quenched) QED: m b (3 GeV) m c (3 GeV) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCDQED = (cid:40) . from m Phh . from m Vhh . (34)These again agree with each other. The weighted average,which is our ﬁnal result, is: m b (3 GeV) m c (3 GeV) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCDQED = 4 . . (35)The error budgets for these ratios are the same as those in Ta-ble III, but with an additional uncertainty of 0.03% associatedwith the QED correction. Mass ratios for other values of therenormalization scale are readily calculated using QED per-turbation theory: m b ( µ ) m c ( µ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCDQED = (cid:16) µ (cid:17) α QED / π m b (3 GeV) m c (3 GeV) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCDQED . (36) R ( m b /m c , Q = ) = 1 . to leading order in QED perturbationtheory. Our results are close to this value but also include nonperturbativecorrections from QCD. The QED uncertainty is obtained by adding (in quadrature) the 0.013%uncertainty in Eq. (29), the 0.019% uncertainty in Eq. (33), and 0.017% forpossible corrections due to quenching QED (10% of the QED correction).

Here and elsewhere, we ignore the running of α QED since ithas a negligible effect compared with our errors.

V. CONCLUSIONS

In this paper we described a new calculation of the ratio ofthe MS masses of the b and c quarks: m b (3 GeV , n f = 4) m c (3 GeV , n f = 4) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) QCDQED = 4 . , (37)where n f is the number of quark ﬂavors in the sea. This is theﬁrst calculation of the mass ratio based on simulations thatinclude QED, and it makes no use of weak-coupling perturba-tion theory. Earlier analyses used phenomenological modelsto estimate QED corrections to the mass ratio, but the preci-sion of the most recent results requires a more accurate treat-ment, like the one described in this paper. QED increasedthe mass ratio by 0.17(3)% relative to our ratio without QED(Eq. (24)), which is almost equal to the standard deviation ofour ﬁnal result. Our new result is consistent at the σ level with earlier re-sults from n f = 4 simulations that did not include QED (andso are µ independent): m b m c = (cid:40) . Fermilab/MILC/TUMQCD [15] . HPQCD [4] . (38)In both cases the listed uncertainties include estimates of theQED effects. Our new result and HPQCD’s previous resultsare nonperturbative; the Fermilab/MILC/TUMQCD result re-lies upon perturbation theory (and Heavy Quark EffectiveTheory), although sensitivity to the perturbative contributionsmostly cancels in the ratio. The Fermilab/MILC/TUMQCDresult comes from simulations of heavy-light mesons ( D s and B s ) rather than the heavy-heavy mesons used here.In a recent paper, HPQCD presented a new value for the c mass that includes (quenched) QED effects as we do here: m c (3 GeV , n f = 4) (cid:12)(cid:12) QCDQED = 0 . . (39)Combining this result with our mass ratio gives a new resultfor the b -quark’s MS mass, m b (3 GeV , n f = 4) (cid:12)(cid:12) QCDQED = 4 . , (40) Note that the “QED correction” to a QCD-only analysis depends in detailon how parameters are set in the QCD-only simulation. Since QCD withoutQED is not the real world, it makes a difference, for example, which hadronmass is used to tune a quark mass; and the QED correction will differ fordifferent choices. Our QED correction is relative to the speciﬁc QCD-onlytheories deﬁned in Section III. The second error in the Fermilab/MILC/TUMQCD result is their estimateof residual QED uncertainties not included in their main analysis (andtherefore not included in the errors stated in their abstract) [15]. which is the ﬁrst based on simulations that include QED. Us-ing perturbation theory to run to the b mass gives m b ( m b ) (cid:12)(cid:12) QCDQED = (cid:40) . n f = 44 . n f = 5 , (41)where we now must include a factor from QED: Z QED m ( µ ) = (cid:0) µ/ (cid:1) − α QED / π . (42)with µ = m b . This new result for the b quark is compared withearlier results in Fig. 4. All of these results agree to withinerrors. APPENDIX

The f Phh ( r ) function plotted in Fig. 1 can be recreated fromthe its values at the knot locations r = m h /m c , f Phh ( r ) = . . . . . . at r = . . . . . . , (43)together with the correlation matrix for these values: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The analogous results from the vector mesons are f Vhh ( r ) = . . . . . . at r = . . . . . . (44)with correlation matrix: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . We use α MS (5 GeV , n f = 4) = 0 . from [4], together with5-loop results for the beta function and mass anomalous dimension, and4-loop results for adding a ﬂavor [21–28]. Acknowledgements

We are grateful to the MILC collaboration for the use oftheir conﬁgurations. We are also grateful for the use ofMILC’s QCD code. We have modiﬁed it to generate quenched U (1) [1] G. P. Lepage, P. B. Mackenzie, and M. E. Peskin, (2014),arXiv:1404.0319 [hep-ph].[2] A. Lee, C. Monahan, R. Horgan, C. Davies, R. Dowdall,and J. Koponen (HPQCD), Phys. Rev. D , 074018 (2013),arXiv:1302.3739 [hep-lat].[3] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage,J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD),Phys. Rev. D , 054502 (2007), arXiv:hep-lat/0610092.[4] B. Chakraborty, C. Davies, B. Galloway, P. Knecht, J. Koponen,G. C. Donald, R. J. Dowdall, G. P. Lepage, and C. McNeile,Phys. Rev. D , 054508 (2015), arXiv:1408.4169 [hep-lat].[5] D. Hatton, C. Davies, B. Galloway, J. Koponen, G. P. Lep-age, and A. Lytle (HPQCD), Phys. Rev. D , 054511 (2020),arXiv:2005.01845 [hep-lat].[6] S. Borsanyi et al. , JHEP , 010 (2012), arXiv:1203.4469 [hep-lat].[7] R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. McNeile,Phys. Rev. D , 074504 (2013), arXiv:1303.1670 [hep-lat].[8] A. Bazavov et al. (MILC), Phys. Rev. D , 054505 (2013),arXiv:1212.4768 [hep-lat].[9] C. Davies, C. McNeile, K. Wong, E. Follana, R. Horgan,K. Hornbostel, G. P. Lepage, J. Shigemitsu, and H. Trottier,Phys. Rev. Lett. , 132003 (2010), arXiv:0910.3102 [hep-ph].[10] D. Hatton, C. T. H. Davies, J. Koponen, G. P. Lepage, and A. T.Lytle, (2021), arXiv:2101.08103 [hep-lat].[11] E. McLean, C. T. H. Davies, A. T. Lytle, and J. Koponen, Phys.Rev. D , 114512 (2019), arXiv:1904.02046 [hep-lat].[12] P. Zyla et al. (Particle Data Group), PTEP , 083C01(2020).[13] See Section 5.2 on the Empirical Bayes criterion in G. P.Lepage, B. Clark, C. Davies, K. Hornbostel, P. Mackenzie,C. Morningstar, and H. Trottier, Nucl. Phys. B Proc. Suppl. , 12 (2002), arXiv:hep-lat/0110175.[14] We use the monotonic spline described in M. Steffen, Astron.Astrophys. , 443 (1990).[15] A. Bazavov et al. (Fermilab Lattice, MILC, TUMQCD), Phys.Rev. D , 054517 (2018), arXiv:1802.04248 [hep-lat].[16] G. P. Lepage, lsqfit v. 11.7, https://github.com/gplepage/lsqfit , doi:10.5281/zenodo.4037174 (2020).[17] The splines are implemented using the gvar Python mod-ule: G. P. Lepage, gvar v. 11.9.1, https://github.com/gplepage/gvar , doi:10.5281/zenodo.4290884 (2020).[18] P. Gambino, A. Melis, and S. Simula, Phys. Rev. D , 014511(2017), arXiv:1704.06105 [hep-lat].[19] A. Bussone et al. (ETM), Phys. Rev. D , 114505 (2016),arXiv:1603.04306 [hep-lat].[20] B. Colquhoun, R. J. Dowdall, C. T. H. Davies, K. Horn-bostel, and G. P. Lepage, Phys. Rev. D , 074514 (2015),arXiv:1408.5768 [hep-lat].[21] P. A. Baikov, K. G. Chetyrkin, and J. H. K¨uhn, Phys. Rev. Lett. , 082002 (2017), arXiv:1606.08659 [hep-ph].[22] F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren, and A. Vogt,JHEP , 090 (2017), arXiv:1701.01404 [hep-ph].[23] T. Luthe, A. Maier, P. Marquard, and Y. Schroder, JHEP ,166 (2017), arXiv:1709.07718 [hep-ph].[24] P. A. Baikov, K. G. Chetyrkin, and J. H. K¨uhn, JHEP , 076(2014), arXiv:1402.6611 [hep-ph].[25] T. Luthe, A. Maier, P. Marquard, and Y. Schr¨oder, JHEP ,081 (2017), arXiv:1612.05512 [hep-ph].[26] Y. Schroder and M. Steinhauser, JHEP , 051 (2006),arXiv:hep-ph/0512058.[27] K. G. Chetyrkin, J. H. Kuhn, and C. Sturm, Nucl. Phys. B ,121 (2006), arXiv:hep-ph/0512060.[28] T. Liu and M. Steinhauser, Phys. Lett. B746