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

We extend HPQCD's earlier $n_f=4$ lattice-QCD analysis of the ratio of $\overline{\mathrm{MSB}}$ masses of the $b$ and $c$ quark to include results from finer lattices (down to 0.03fm) and a new calculation of QED contributions to the mass ratio. We find that $\overline{m}_b(\mu)/\overline{m}_c(\mu)=4.586(12)$ at renormalization scale $\mu=3$\,GeV. This result is nonperturbative. Combining it with HPQCD's recent lattice QCD$+$QED determination of $\overline{m}_c(3\mathrm{GeV})$ gives a new value for the $b$-quark mass: $\overline{m}_b(3\mathrm{GeV}) = 4.513(26)$GeV. The $b$-mass corresponds to $\overline{m}_b(\overline{m}_b, n_f=5) = 4.202(21)$GeV. These results are the first based on simulations that include QED.


I. INTRODUCTION
Accurate masses for heavy quarks are important for QCDphenomenology generally, but they will be particularly important for high-precision searches for new physics in Higgs decays [1].In this paper we present a new result for the ratio m b /m c of the MS masses of the b and c quarks.Our analysis of the mass ratio is completely nonperturbative.This is in contrast to lattice-QCD determinations of the separate quark masses, which need QCD perturbation theory to relate MS masses to lattice quantities.Thus the (nonperturbative) mass ratio provides a nontrivial check on (perturbative) determinations of the separate masses.The ratio can also be combined with recent accurate determinations of the c-quark mass to obtain new results for the b-quark mass.
Lattice simulations of b quarks are complicated by the quark's large mass, which leads to large lattice-spacing errors when the b quarks are described by the Dirac equation (as opposed to, say, NRQCD [2]).We address this problem by using a Highly Improved Staggered-Quark discretization of the Dirac equation (HISQ [3]) that is also highly efficient, making simulations at very small lattice spacings feasible.Our previous analysis of the mass ratio [4] used lattices with spacings down to 0.06 fm but still required an extrapolation in the quark mass to reach m b .Here we reduce the lattice spacing to 0.03 fm, where am b ≈ 0.6, which allows us to simulate at the b mass.Lattice spacing errors at m b are less than 1% on our finest lattice, and we are able to remove most of that error by extrapolating from results covering a range of lattice spacings and heavy-quark masses m h .
Our new result is accurate to about 0.25%, so it becomes important to include QED effects.We recently analyzed the QED contributions to the c quark's mass [5].Here we adapt the methods from our earlier paper to provide the first results for QED contributions to m b /m c and m b .Here and in our earlier paper we use the quenched QED approximation, which omits contributions from photons coupling to sea quarks.The quenched approximation should capture the bulk of the QED correction in mesons whose valence quarks are both heavy; contributions from sea quarks are expected to be an order of magnitude smaller [5].
In Section II we describe our general strategy and the lattice QCD simulations we employed.In Section III, we extract a value for m b /m c using results from simulations without QED.We then add QED effects in Section IV.We summarize our results for the mass ratio in Section V and combine them with HPQCD's recent c-quark mass to obtain a new result for the b-quark mass.

II. LATTICE QCD SIMULATIONS
We use gluon configuration sets generated on a variety of lattices (by the MILC collaboration [8]), with n f = 4 flavors of HISQ sea quark and lattice spacings ranging from 0.09 fm to 0.03 fm.These sets are described in Table I.The u and d quark masses are set equal to m ≡ (m u + m d )/2; corrections to this approximation are quadratic in the light-quark masses for our analysis, and so are negligible.We include results where the light-quark masses in the sea are tuned close TABLE I. Gluon configuration sets used in this paper.Sets are grouped by approximate lattice spacing, with lattice spacings of 0.09 fm (Sets 1 and 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 flow parameter w0/a [6] where w0 = 0.1715(9) fm [7].The sea quark masses are given in lattice units for u/d quarks (am ), s quarks (ams), and c quarks (amc).Tuned values for the lattice c masses are also given (in GeV).These masses are adjusted to give correct masses for either the ηc or J/ψ mesons (Eq.( 3)).The c masses are tuned using slope d mc/dmcc, which is the same (within errors) for the ηc and J/ψ.The spatial and temporal sizes of the lattices, L and T , are listed, as are the number of configurations used in our analysis (the two numbers for Set 1 are for quark masses am h below and above 0.5; the three numbers for Set 2 are for the pseudoscalar correlators, the vectors with mass with mass below 0.5 and the vectors with mass above 0.5).The three polarizations were averaged for vectors.Eight time sources were used on each configuration except for Set 6 where four were used.to their physical value, but we also include results with much larger light-quark masses in the sea.Results from these last simulations are unphysical but are easily corrected [4].Increasing the light-quark mass in the sea significantly reduces the cost of our analysis at small lattice spacings (because smaller 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 quark masses used in the lattice Lagrangian, up to corrections that vanish in the continuum limit [9]: where the bare masses are tuned so that the QCD simulations reproduce the experimental results for meson masses.This relationship between the MS and lattice quark masses is nonperturbative and independent of the MS renormalization scale µ.Pseudoscalar and vector meson masses from our simulations are listed in Table II for a variety of (valence) heavyquark masses, ranging on the finest lattices (Sets 5 and 6) approximately from the c mass to the b mass.The analysis methods for extracting these masses (and most of the results) come from [5,10].We use multi-exponential fits to calculate the masses.Fig. 1 compares the result from our fit with the effective mass values at various times for the correlator closest to the η b mass on our finest lattice.
The quark masses am h in Table II are what is used in the HISQ Lagrangian.The a mh masses are corrected to remove tree-level (am h ) 2n errors (in the pole mass) through order 2n = 10 [3,11]: We write the expansion as powers of am h /2 because this makes the leading coefficients roughly the same size (about 1/2).The correction is −2% at am h = 0.9, which is the largest value we use.
We give results in Table I for the tuned bare c mass for each of the configurations.In each case we adjust the c mass so as to reproduce the continuum value for either the η c mass or the J/ψ mass: Here we have subtracted 7.3(1.2) MeV from the experimental value for m ηc [12] to account for the fact that we are not including contributions from cc annihilation in our simulations; this correction is determined in [5].The analogous correction to the J/ψ mass is negligible, but we have subtracted 0.7(2) MeV from the mass to account for cc annihilation to a photon; this correction is estimated perturbatively in [5].We extrapolate the c masses to their correct values using where m cc is either the η c or the J/ψ mass and the slopes (see Table I) are estimated from splines fit to the entries in Table II.
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 these have significant (am h ) 2n errors.Instead we will use the data in Table II to define functions that relate the ratio of quark masses to the pseudoscalar (P ) or vector (V ) masses: where TABLE II.Lattice QCD results for the ground-state pseudoscalar and vector hh mesons in lattice units: am P hh and am V hh , respectively.Results are given for each configuration sets (Table I) and a variety of bare quark masses am h and corrected masses a mh (in lattice units).The uncertainties in the meson masses are negligible compared with other errors in our analysis and so have no impact on the final results.Most of these results are from [5,10].  1.The effective mass plotted versus time for the am h = 0.65 pseudoscalar correlator from configuration Set 6.For clarity, the plot includes only every third point.The orange band and dotted line show the corresponding mass (Table II) obtained from a multiexponential fit [5,10].The error in the fit result (orange band) is almost entirely statistical in origin.In particular, possible biases due to excited states are completely negligible (50× smaller), as is typical in fits for heavy-quark ground-state masses.

Set
Given these functions, we then obtain two estimates for m b /m c by solving each of the equations for m b /m c , where The two mass ratios should agree.Here we account for the bb annihilation contribution to the η b mass by adding an extra error of ±1 MeV to the experimental result [12]; this estimate is based on NRQCD perturbation theory and the meson's width [3].The analogous contribution to the Υ mass is negligible, as is Υ annihilation via a photon.
In what follows, we first describe our lattice-QCD analysis of f P hh and f V hh , and then discuss the results.

A. Analysis
We determine the f hh functions (Eq.( 5)) by fitting the meson masses am hh from Table II to functions of the following form, where is the effective ratio of quark masses m h /m c .Pseudoscalar and vector mesons are fit separately, to determine each of f P hh and f V hh .We describe each element of the fit 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 the table are very small.It is impossible to fit the almost six significant digits in these data with a model as simple as we use here.So we increase the fractional error on each value to σ u , which is then a measure of the part of the variation in the data that is unexplained by our model.The σ u errors are uncorrelated from one am hh to another.We use the same value for σ u for every data point and adjust its size to maximize the Bayes Factor from the fit [13].For the parameters and model used here, we find that which means that our model explains the individual data points to within ±0.025%.A simpler model would have a larger σ u : for example, σ u more than doubles if the ξ m factors in ratio r are dropped (but gives consistent results within the larger errors).Note that the statistical errors listed in Table II can be neglected when σ u is included.
• We parameterize the f hh functions as splines [14] with 6 knots evenly spaced from r = 1 to 4.6 (≈ m b /m c ), inclusive.The fit parameters are the function values at the knots.These functions are linear up to corrections of order v2 /c 2 ∼ 0.1-0.3,where v is the typical velocity of the heavy quarks in the meson.Therefore we use the following priors for the values at the knots with r > 1: where m cc and m bb are the continuum masses of the pseudoscalar/vector mesons composed of c and b quarks, respectively (Eqs.( 3) and ( 8)).At r = 1, we require We choose 6 knots to maximize the Bayes Factor from the fit.Results obtained using 5 or 7 knots agree well with those from 6 knots, with similar or smaller errors.Doubling the width of the priors has no effect on our results.
• The ξ m factors in the mass ratio r rescale the quark masses to correct for detuned values of the light sea quarks.From [4], where is the difference between the masses used in the simulation and their tuned values.The tuned masses are determined from the tuned c-quark mass using results for m c /m s and m s /m from [15].The priors for the fit parameters are which come from fits described in [4].
• largest simulation errors are from the discretization.These are suppressed by α s (π/a) in order a 2 because we are using the HISQ formalism [3].Beyond this order they are suppressed either by α s (π/a) or by v 2 /c 2 , since we have removed the tree-level a 2n errors in the quark masses using Eq. ( 2).The fit can't distinguish easily between α s and v 2 /c 2 since both are around 0.2 for our data, so we include only an α s correction, modeled after Eq. ( 2): Here functions f n a 2 (r) are 6-knot splines with priors at the knots (same locations as above) of Terms beyond n = 3 have no effect on the fit results; keeping just the n = 1 term gives the same final results but with errors that are 25% smaller.The splines allow for m h dependence in the a 2n corrections.
• We include a 2 corrections to ξ m since δm sea uds is large for some of our configuration sets: where function f sea uds (r) is again a 6-knot spline, now with priors at the knots of where the width is chosen to be somewhat larger than suggested by g m above.Omitting this correction has negligible effect on our final results.
• We also include a correction to ξ m from detuned cquark masses in the sea.This correction should be small because of heavy-quark decoupling [4] -the momentum transfers in the heavy-quark mesons are too small to produce cc pairs efficiently.We include the correction where δm c ≡ m c − m tuned c , and f sea c (r) is a 6-knot spline with f sea c (r knot ) = 0.00(1).
We choose the width to maximize the Bayes Factor from the fit.Omitting this correction has negligible effect on our final results.
The fit parameters are the values of the coefficient functions (splines) at the knots, together with g m and ζ from the ξ m factors (Eq.( 14)).We use the lsqfit Python module to do the fits [16,17].

B. Results
The functions f P/V hh obtained from the fits described in the previous section (and detailed in the Appendix) are plotted in Fig. 2, together with the data from Table II.Fig. 3 shows that the model (Eq.( 9)), with best-fit values for the fit parameters in the corrections (on the right-hand side), reproduces the data within errors. 1 The difference between the lattice results with  (10).The lines, which vary in thickness, show results from the best-fit values for the functions f P/V hh ; the line thickness shows the 1σ uncertainty in these functions.These functions can be reconstructed from information in the Appendix.Separate results are shown using the pseudoscalar masses m P hh (top line, squares) and the vector masses m P hh (bottom line, circles).Different colors indicate different configuration sets, with Sets 6 (brown) and 5 (purple) having the largest masses, followed by Sets 4 (red), 3 (green), 2 (orange) and 1 (blue), in that order.Error bars are smaller than the plot symbols.II and the model in Eq. ( 9) with best-fit values for the fit parameters.Results are shown for both pseudoscalar (squares, offset right) and vector (circles, offset left) mesons, where data points for the two mesons are offset slightly in opposite directions to improve visibility.and without corrections is −0.69(23)% for the highest quark mass on the finest lattice (Set 6).
We can use functions f P/V hh to extract values for the ratio of MS masses by solving Eqs.(7).We obtain m b /m c = 4.578(12) from the am P hh 4.578 (15) from the am V hh , independent of renormalization scale.The two estimates agree to within 0.01(23)%.The weighted average, taking ac- Results are given for determinations using pseudoscalar mesons (m P hh ) and vector mesons (m V hh ), and for the weighted average of these results.The dominant errors come from the extrapolation to zero lattice spacing and from uncertainties in the lattice spacing.Additional errors are from residual uncertainties taken in the fit data (σu), uncertainties in ξm used to correct for unphysical sea-quark masses, uncertainties in the ηc and η b masses, and tuning uncertainties in the sea-quark masses.The error budgets are the same when QED is included aside from an additional uncertainty of 0.03% associated with the QED corrections.
We tabulate the leading uncertainties in our two results in Table III.The error budgets are similar for the two mesons, and are dominated by uncertainties associated with discretization errors and the lattice spacings.Doubling the widths of any of the priors associated with these uncertainties has negligible effect on the central values from our fits (< σ/3), and only doubling the discretization priors (Eq.( 18)) has an appreciable impact on the final uncertainties, as expected from Table III.Omitting results from the coarsest lattices (Sets 1 and 2) has negligible effect on our results (< σ/10).Omitting results from the finest lattice (Set 6) increases the final uncertainties significantly (by factors of 5-6) because there is then insufficient data at large masses to constrain the 6-knot splines used in the fit function; reducing the number of knots decreases the errors by a third.In either case the results agree with our final results within errors.
Finally, as discussed in [5], we expect errors from the finite lattice volume and strong-isospin breaking (m u = m d ) in the sea to be less than 0.01% and so negligible here.We have verified this for the meson masses (using configuration Sets 3A-3B and 5-7 from [5]).The Wilson flow parameter w 0 should be similarly insensitive, and we have verified this to the level of our statistical errors (0.03%) for w 0 /a.See [5] for further details.

IV. ADDING QED
Adding QED complicates the analysis of m b /m c because the quarks have different QED charges and therefore different TABLE IV.Ratio R0 of m hh with QED corrections to m hh without QED corrections, evaluated at the same quark mass m h .Results are shown for ground-state pseudoscalar and vector mesons analyzed on two configuration sets.The quark's QED charge is Q times the proton's charge; results for Q = 2/3 can be converted to Q = 1/3 by replacing R0 with 1 + (R0 − 1)/4.mass anomalous dimensions.Thus the nonperturbative relation in Eq. ( 1) is only true up to O(α QED ) corrections.We deal with this complication by introducing QED through two ratios R:

Set
Here is the ratio of the c mass in a theory with c-quark charge 2 3 to the mass in a theory with c-quark charge 1  3 .Similarly, 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 quark masses are tuned to reproduce the continuum meson masses in Eqs. ( 3) and (8).Either the pseudoscalar or vector mesons can be used; they give the same results to within the precision needed here.Only the first of the R factors (Eq.( 26)) depends on the MS renormalization scale µ; we take µ = 3 GeV, following [5].We approximate full QED by quenched QED, where only the valence quarks carry electric charge.This is expected to be the dominant contribution in O(α QED ) and is much less costly to analyze.The techniques we use for introducing QED into simulations are standard and are described in [5,10].
The R factor for the c mass (Eq.( 26)) is the most important and can be inferred from our earlier result [5]: R is quadratic in Q c to better than 0.01% so the QED correction (R − 1) required to go from charge 1 3 to charge 2 3 is three 3 ) is plotted versus m h /mc.It is the ratio of m h /mc computed with QED charge Q = 1  3 to the result without QED (Q = 0), where the quark masses are tuned to give the same results for m P hh (bottom, red) or m V hh (top, blue).Results are shown from configuration Sets 1 (squares) and 3 (circles).Errors are smaller than the plot symbols.The blue and red shaded areas show the ±1σ fits to the data (Eq.( 32)).
quarters that required to go from 0 to 2  3 : The other R factor, for m b /m c , is expected to be much closer 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 = 1 3 are smaller (and the QED effect is quadratic in the charge).To estimate the effect, we calculated the ratio R 0 of meson masses m hh with and without Q c,b = 1  3 QED, holding the quark masses constant, for two of our configuration sets; our results are in Table IV.This quantity can be related to the R factor for m b /m c by re-expressing the R-factor in terms of lattice masses, using Eq.(1) (since Q c = Q b ), and writing it as where δ mQ c,b are the quark mass shifts needed to hold the meson masses constant when QED is added to the simulation.The mass shifts can be calculated for different heavy-quark masses m h from the R 0 factors in Table IV: Here the derivative is estimated for each configuration by fitting a cubic spline to the a mh values in Table II as a function of the corresponding am hh values.
Values for R(m h /m c , Q = 1 3 ) are plotted versus m h /m c in Fig. 4 for both pseudoscalar (below) and vector (above) mesons from the two configuration sets.We fit these data to a simple function suggested by QED perturbation theory: R = 1 + (33) The two results agree with each other, but the corrections are too small to affect our final results significantly. 2 Doubling the fit priors leaves the results unchanged.We use the larger error in the error budgets for our final result.
Including both R factors, we arrive at new results for the quark mass ratio at µ = 3 GeV that include (quenched) QED: These again agree with each other.The weighted average, which is our final result, is: The error budgets for these ratios are the same as those in Table III, but with an additional uncertainty of 0.03% associated with the QED correction. 3Mass ratios for other values of the renormalization scale are readily calculated using QED perturbation theory, where the additional QED correction is negligible compared to our errors for typical values of µ.Here and elsewhere we ignore the running of α QED and O(α QED α s ) corrections since they are also negligible compared with our errors.

V. CONCLUSIONS
In this paper we described a new calculation of the ratio of the MS masses of the b and c quarks: where n f is the number of quark flavors in the sea.This is the first calculation of the mass ratio based on simulations that 2 R(m b /mc, Q = 1 3 ) = 1.00059 to leading order in QED perturbation theory.Our results are close to this value but also include nonperturbative corrections from QCD. 3 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% for possible corrections due to quenching QED (10% of the QED correction).Gambino et al [18], ETM [19], HPQCD '14 (NRQCD) [20], and HPQCD '14 (HISQ) [4].The gray band corresponds to the top result (HPQCD '21), the only one from simulations that include QED.
include QED, and it makes no use of weak-coupling perturbation theory.Earlier analyses used phenomenological models to estimate QED corrections to the mass ratio, but the precision of the most recent results requires a more accurate treatment, like the one described in this paper.QED increased the mass ratio by 0.17(3)% relative to our ratio without QED (Eq.( 24)), which is almost equal to the standard deviation of our final result. 4ur new result is consistent at the 1σ level with earlier results from n f = 4 simulations that did not include QED (and so are µ independent): .578(8)(10) Fermilab/MILC/TUMQCD [15] 4.528(54) HPQCD [4].
(38) In both cases the listed uncertainties include estimates of the QED effects. 5Our new result and HPQCD's previous results are nonperturbative; the Fermilab/MILC/TUMQCD result relies upon perturbation theory (and Heavy Quark Effective Theory), although sensitivity to the perturbative contributions mostly cancels in the ratio.The Fermilab/MILC/TUMQCD result 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: where we now include an evolution factor from QED: with µ = m b (which shifts the result by less than 0.02% and so is negligible).This new result for the b quark is compared with earlier results in Fig. 5.All of these results agree to within errors.

APPENDIX
The f P hh (r) function plotted in Fig. 2 can be recreated from the its values at the knot locations r = m h /m c ,

FIG. 2 .
FIG.2.Lattice QCD (without QED) results for m h /mc plotted versus the hh meson masses.Values for m h /mc are corrected as in Eq.(10).The lines, which vary in thickness, show results from the best-fit values for the functions f

FIG. 3 .
FIG.3.Relative difference between the data for am hh from TableIIand the model in Eq. (9) with best-fit values for the fit parameters.Results are shown for both pseudoscalar (squares, offset right) and vector (circles, offset left) mesons, where data points for the two mesons are offset slightly in opposite directions to improve visibility.

15 )
from m V hh .

TABLE III .
Contributions to the total (1σ) error in m b /mc from QCD simulations (without QED), as a percentage of the mean value.