Bottomonium precision tests from full lattice QCD: hyperfine splitting, Υ leptonic width and b quark contribution to e^+e^- \rightarrow hadrons
D. Hatton, C. T. H. Davies, J. Koponen, G. P. Lepage, A. T. Lytle
BBottomonium precision tests from full lattice QCD: hyperfine splitting, Υ leptonicwidth and b quark contribution to e + e − → hadrons. 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: March 1, 2021)We calculate the mass difference between the Υ and η b and the Υ leptonic width from latticeQCD using the Highly Improved Staggered Quark formalism for the b quark and including u , d , s and c quarks in the sea. We have results for lattices with lattice spacing as low as 0.03 fmand multiple heavy quark masses, enabling us to map out the heavy quark mass dependence anddetermine values at the b quark mass. Our results are: M Υ − M η b = 57 . . .
0) MeV (where thesecond uncertainty comes from neglect of quark-line disconnected correlation functions) and decayconstants, f η b = 724(12) MeV and f Υ = 677 . .
7) MeV, giving Γ(Υ → e + e − ) = 1 . b quark contribution to σ ( e + e − → hadrons) determined from experiment. Moments 4–10provide a 2% test of QCD and yield a b quark contribution to the anomalous magnetic moment ofthe muon of 0.300(15) × − . Our results, covering a range of heavy quark masses, may also beuseful to constrain QCD-like composite theories for beyond the Standard Model physics. I. INTRODUCTION
Weak decay matrix elements calculated in lattice QCDare critical to the flavour physics programme of over-determining the Cabibbo-Kobayashi-Maskawa (CKM)matrix to find signs of new physics. For that programmethe weak decays of b quarks are particularly importantsince they give access to the least well known CKM el-ements, V ub and V cb . These CKM matrix elements canbe determined either using exclusive decay channels andlattice QCD form factors or inclusive decay channels andmeasured spectral shape functions. There continues to besome tension between the exclusive and inclusive deter-minations [1] that needs further improvements to bothapproaches to resolve. On the lattice QCD side thismeans developing improved approaches to B meson weakdecay matrix elements, such as [2, 3], but also providingmore stringent tests of lattice QCD results in b physicsto make sure that sources of systematic error are un-der full control. Here we provide three such tests usingbottomonium correlation functions; the ground-state hy-perfine splitting (the mass difference between the Υ and ∗ [email protected] † [email protected] ‡ η b ), the Υ leptonic width and the b quark contributionto R ( e + e − → hadrons). The latter two, being electro-magnetic processes, can be compared with experimentfree from CKM uncertainties. We obtain the most accu-rate results to date for these quantities and are able toinclude the effect of the b quark’s electric charge in thecalculation for the first time.We use the Highly Improved Staggered Quark (HISQ)discretisation of the quark action for these calculations.The HISQ action was developed in [4] to have small dis-cretisation errors with the leading errors, quadratic in thelattice spacing, removed. This makes the action partic-ularly good for heavier quarks when discretisation errorsappear as powers of the quark mass in lattice units, whichcan be relatively large. This action enabled the first ac-curate lattice calculations in charm physics [5–8]. Morerecently it has been used to achieve sub-1% accuracy inthe charmonium hyperfine splitting, J/ψ leptonic width, m c and c quark vacuum polarisation contribution to theanomalous magnetic moment of the muon [9]. The cal-culation used a range of lattice spacing values from 0.15fm to 0.03 fm with u , d , s and c quarks in the sea andincluded the effect of the c quark’s electric charge [9].Here we will extend this latter calculation to bottomo-nium. Because the b quark mass is much larger than thatof c , we need fine lattices to reach the b with a quarkmass in lattice units, am b < a r X i v : . [ h e p - l a t ] F e b sation errors. Our strategy, known as the heavy-HISQapproach [10], is to perform calculations for a range ofmasses between c and b on lattices with a range of finelattice spacings. We can then map out the dependenceon the heavy quark mass both of the quantity being cal-culated and its discretisation errors. This enables us todetermine a physical result at the b quark mass.This approach has been very successful for decay con-stants and spectroscopy for heavy-light ( B , B s and B c )mesons [10–12] and is now being used for the form factorsfor B meson weak decays [2, 3]. Here we will apply thisapproach to the Υ for the first time.There are alternative nonrelativistic approaches thatcan be used at the b quark mass on coarser lattices;see [13] for the determination of the Υ and Υ (cid:48) leptonicwidths using lattice NonRelativistic QCD. The nonrel-ativistic expansion of the Hamiltonian and the currentsthat appear in matrix elements gives systematic uncer-tainties from missing higher-order relativistic correctionsand from renormalisation of the lattice current to matchthe continuum current. These uncertainties hinder testsat high precision.In contrast the HISQ action is relativistic and theHISQ vector current can be matched accurately and fullynonperturbatively to that in continuum QCD [14]. Aswe will show below, this enables us to improve the lat-tice QCD accuracy on the bottomonium hyperfine split-ting to better than 5% and to achieve percent-level pre-cision on the Υ and η b decay constants and on mo-ments that parameterise the b quark contribution to R ( e + e − → hadrons).Our results cover the range of heavy quark masses from c to b and we will give results for decay constant to massratios over this range. These could be useful both fortuning of phenomenological models of QCD but also asconstraints on QCD-like composite theories for beyondthe Standard Model physics.The paper is organised as follows. In the next Sectionwe give details of the lattice calculation we perform. Thisincludes a general description of the fits that we use todetermine the heavy mass dependence of the quantitiescalculated. We then present results for the bottomoniumhyperfine splitting in Section III, the decay constants forboth the η b and Υ in Section IV and the time-moments ofthe vector current-current correlators in Section V. Eachsection includes a description of the calculation and thena discussion subsection with comparison to experimentand previous lattice QCD results. Section IV on decayconstants also includes plots of decay constant to massratios and the ratio of vector to pseudoscalar decay con-stants over the quark mass range from c to b . We thengive our conclusions in Section VI. II. LATTICE CALCULATION
We use ensembles of lattice gluon field configurationsprovided by the MILC collaboration [15] at values of the lattice spacing, a ≈ α s a -improved discretisation of the gluon action [16] and in-clude the effect of u , d , s and c quarks in the sea with theHISQ formalism [4]. The u and d masses are taken to bethe same and we denote this mass m l . For most of theensembles we have unphysically heavy u/d quarks with m l /m s ≈ m l and lattice spacing values 0.09 fm and 0.06fm. We expect sea quark mass effects to be small for theΥ because it has no valence light quarks. However, ananalysis of such effects is needed for accurate results.Table I lists the parameters of the ensembles that weuse. The lattice spacing is determined in terms of theWilson flow parameter w [17]. On these ensembles wecalculate quark propagators from random wall sourcesusing the HISQ action and with a variety of masses, m h ,from that of the c quark upwards. The valence heavyquark masses that we use on each ensemble are listed inTable II. The value of (cid:15) Naik used in the coefficient of theNaik term in the HISQ action [4] is taken as the tree-levelfunction of the quark mass given in [6].We combine the quark propagators into (connected)meson correlation functions for both pseudoscalar ( η h )and vector ( φ h ) mesons, using the local γ and γ i operators converted to appropriate form for staggeredquarks [8, 9]. Note that we do not include quark-line dis-connected corrrelation functions that take account of theheavy quark/antiquark annihilation to gluons. We ex-pect the effect of the disconnected correlation functionsto be very small in the heavyonium system. In [9] ourresult for the mass difference between J/ψ and η c mesonswas accurate enough, for the first time, to see a differ-ence with experiment of 7.3(1.2) MeV. We concluded thatthis was the effect of the missing disconnected correlationfunction on the η c mass. Here we will test for a similareffect on the η b .On the coarsest two ensembles we use 16 time-sourceson each gluon field configuration for high statistics; wetake 8 time-sources on the other ensembles. We use atleast 100 configurations on each ensemble for a good sta-tistical sample. In generating the very fine lattice (set 6in Table I) a slow evolution in Monte Carlo time of thetopological charge was observed [15], so that the ensem-ble does not explore many topological sectors. However,it has been shown that the impact of this on calculationsfor heavy mesons is negligible [18].We fit the correlation functions from each ensemble us-ing a multi-exponential constrained fit [21] and followingthe method in [9]. The fit form used for the pseudoscalarcorrelators as a function of t , the time separation betweensource and sink, is C P ( t ) = (cid:88) i A Pi f ( E Pi , t ) , (1)and the vector fit form is C V ( t ) = (cid:88) i (cid:16) A Vi f ( E Vi , t ) − ( − t A V,oi f ( E V,oi , t ) (cid:17) . (2) TABLE I. Sets of MILC configurations [15] used here withHISQ sea quark masses, m sea l ( l = u/d ), m sea s and m sea c givenin lattice units. The lattice spacing is given in units of w [17];the physical value of w was determined to be 0.1715(9) fmfrom f π [19]. Sets 1 and 2 are ‘fine’ ( a ≈ .
09 fm), sets 3 and4 are ‘superfine’ ( a ≈ .
06 fm), set 5 ‘ultrafine’ ( a ≈ . a ≈ .
03 fm). The final two columnsgive the extent of the lattice in each spatial direction ( L s ) andtime ( L t ).Set label w /a am sea l am sea s am sea c L s L t Here f ( E, t ) = e − Et + e − E ( L t − t ) . (3)The term that oscillates in time in the vector case re-sults from the use of staggered quarks. E is the massof the lowest lying state (either pseudoscalar or vector)and A is related to the meson decay constants. Theground-state pseudoscalar meson we will denote as η h and the vector as φ h . We fit the correlation functionsfor all masses on a given ensemble simultaneously (withtwo exceptions, see Table II). This means that the cor-relations between results for different masses are carriedthrough the rest of the calculation. The correlations be-tween the φ h and η h correlators are safely neglected asthe uncertainty in the φ h results dominates that for the η h . Results for the ground-state mesons are listed in Ta-ble II.We also have a limited amount of data which includesthe effects of quenched QED (electrically charged valencequarks, but not sea quarks). This allows us to assess theimpact of QED and appropriately account for it in ourerror budgets. As in [9] we use photon fields in Feynmangauge in the QED L formalism [22]. Our quenched QEDcalculations [9, 20] used a valence quark electric chargeof 2/3 e (i.e. the charge on a c quark), where e is themagnitude of the charge on an electron. We can use theseresults to determine the electromagnetic correction forthe b electric charge of − (1 / e . Given the smallness of α QED we take QED corrections to be linear in the quarkcharge squared, Q , and simply rescale the effect of QEDby a factor of 1/4 from that for Q = (2 / e . Results aregiven in Table III in the form of the ratio, R , of resultsin QCD+QED to those in pure QCD at a fixed valueof the valence quark mass in lattice units. We see therethat the impact of QED is tiny but visible. We showedin [9] that QED finite-volume effects were negligible forthe electrically neutral charmonium mesons. This willcontinue to be true for the heavier mesons that we studyhere and so we ignore such effects.As we showed in [9], fixing the lattice spacing from w and f π , as we have done, means that QED corrections tothe lattice spacing should be at the sub-0.1% level (com-ing from QED effects in the quark sea). We can thereforecompare QCD plus quenched QED to pure QCD usingthe same value of the lattice spacing (i.e. that from Ta-ble I). QED affects the tuning of the lattice quark masses,however. We use the simple and natural scheme of tun-ing the b quark mass in both the QCD+QED and pureQCD cases so that the Υ mass has the physical value. Wecan then use our results to determine the impact of QEDon the quantities we study, taking the renormalisation ofthe quark mass into account. We will give that informa-tion, after fitting, as the renormalised ratio R which isthe ratio of the QCD+QED result to that in pure QCDwhen both theories have a b quark mass separately tunedso that the Υ mass takes the experimental value in bothcases.For each quantity that we study we must fit our re-sults, in physical units, as a function of heavy quark massand lattice spacing to determine results in the continuumlimit at the physical b quark mass. We will use the φ h mass as a physical proxy for the heavy quark mass andthen the physical point is defined by the φ h mass becom-ing equal to that of the Υ.In previous studies of heavy meson masses and decayconstants using the heavy-HISQ method the HPQCD col-laboration have used fit forms that capture the heavymass dependence as a polynomial in the inverse H s or η h mass [10, 11]. In the case of heavy-light mesons this formis justified by the heavy quark effective theory (HQET)expansion. In the case of heavyonium the HQET expan-sion is not valid but the same form may still be expectedto work as a Taylor expansion over a finite region in m h .Here, however, we choose to use a fit form that is moreagnostic with regards to the dependence on the heavyquark mass. We achieve this by using cubic splines be-tween specified knot positions . We do not expect toneed many knots because the quantities we study hereshould be smooth monotonic functions of M φ h in thecontinuum limit at physical sea quark masses. The fitfunction for our lattice results also needs to include de-pendence on the lattice spacing and the mistuning of seaquark masses. Both of these effects can also depend onthe heavy quark mass ( M φ h ) through smooth monotonicfunctions and so we also include cubic splines in theirdescription.We use fits of the following form for the pure QCD We use splines that are monotonic between knots (Steffensplines [23]).
TABLE II. Results in lattice units for the masses of the ground-state pseudoscalar meson, η h , and ground-state vector meson, φ h , for valence heavy quark masses in lattice units listed in column 2, for the ensembles listed in column 1. Results come fromsimultaneous fits to all heavy quark masses on a given ensemble, except for the cases marked with an asterisk [20]. These useddifferent random numbers for the sources and so are not correlated with the other results for that ensemble. Column 5 gives themass difference in lattice units between the φ h and η h , column 6 the η h decay constant and column 7 the raw (unrenormalised) φ h decay constant. The required Z V factors are taken from [14].Set am h aM η h aM φ h a ∆ M hyp af η h af φ h /Z V ∗ ∗ e/
3, to a subset of the results of Table II presented as theratio, R , of the value in QCD+QED to that in pure QCD at fixed valence quark mass in lattice units.Set am h R [ aM η h ] R [ aM φ h ] R [ a ∆ M hyp ] R [ af η h ] R [ af φ h /Z V ]2 0.866 1.0002170(45) 1.0002637(14) 1.00307(25) 1.001255(68) 1.001343(62)3 0.274 1.0003937(42) 1.0004468(28) 1.001875(83) 1.00078(12) 1.000708(83)0.548 1.0003306(12) 1.00036876(24) 1.002882(80) 1.001127(35) 1.0009805(75) part of our fit: F ( a, M φ h )[QCD] = A (cid:2) F ( M φ h ) + G (1 /M φ h ) (4)+ (cid:88) i =1 G ( i )1 ( M φ h )( am h ) i + (cid:88) j =1 G ( j )2 ( M φ h )( a Λ) j + G ( M φ h )( a Λ) ( am h ) + (cid:88) k =0 G ( k )4 ( M φ h )( am h ) k δ l + G ( M φ h ) δ c (cid:3) F ( a, M φ h ) are the lattice QCD results in physical unitsof GeV for the hyperfine splitting and decay constants and GeV − for the time-moments. A is then a dimen-sionful number of a size commensurate with the size ofthe quantity being fitted (this will be given in each sec-tion) so that the rest of the fit in square brackets is di-mensionless. F is a dimensionless function of M φ h thatdiffers between the different quantities we examine (andwill also be given in each section). It is a very simplefunction of M φ h that captures the general trend in themass dependence, on top of which the corrections mod-elled by splines are relatively minor. G n denotes a cubicspline. All of the splines have different parameters butkeep the same positions for the knots. The splines arefunctions of M φ h , where M φ h is in GeV, except for thefirst spline, G , which provides the physical correctionsto F . We found that a spline in 1 /M φ h rather than M φ h gave a better χ for quantities like the hyperfine splittingwhich fall as M φ h grows, approximately as the inverse.We allow for two kinds of discretisation effects, those thatare set by the heavy quark mass am h and those that areindependent of the heavy quark mass and instead are setby a fixed scale that we call Λ. We take Λ to be 0.5 GeV.The G ( i )1 splines allow for heavy mass dependence in the am h discretisation effects and the G ( j )2 splines allow forheavy mass dependence in the a Λ discretisation effects.We also include a mixed term in ( am h ) ( a Λ) with spline G , although this has little impact on the fits.The last two lines of the fit form in Eq. (4) allow for seaquark mass mistunings. The light quark mass mistuningparameter δ l is defined as: δ l = 2 m sea l + m sea s − m phys l − m phys s m phys s . (5) m phys s and m phys l are taken to be the same as those usedin [9]. The charm quark mass mistuning parameter isdefined similarly: δ c = m sea c − m phys c m phys c . (6)We allow for discretisation effects and heavy mass depen-dence set by splines G ( k )4 within the light sea quark massmistuning term. We allow for heavy mass dependence inthe sea charm mass mistuning term through the spline G .For data that includes quenched QED effects we addextra terms to the fit function in Eq. (4). The full fittakes the form: F ( a, M φ h )[QCD + QED] = F ( a, M φ h )[QCD] + (7) Aα QED Q (cid:104) ˆ G (1 /M φ h ) + c QED ,am h ˆ G ( M φ h )( am h ) (cid:105) . ˆ G and ˆ G are additional spline functions.For all the fits considered here we use knots placed at { . , . , } GeV, taking knots at the beginning and endof the fit range and one in the middle. This is the opti-mal number of knots according to the Bayes Factor. Wehave checked that varying the sums over discretisationeffects so that the number of terms included changes by ± c F , in F andfor the values of the spline functions at each knot. We usepriors of 0 ± lsqfit python M φ h [GeV]0 . . . . M φ h − M η h [ G e V ] f-5f-physsf-5sf-phys uf-5ef-5PDG FIG. 1. The heavyonium hyperfine splitting as a functionof the vector heavyonium mass, M φ h . The points show ourlattice results from Tables II and III, with different symbolsdenoting different ensembles as in the legend. The errorsare dominated by uncertainties from a that are correlatedbetween the points. QCD+QED points are shown in cyan.They are not distinguishable from their pure QCD counter-parts, but are visible by being shown on top of these values.Cubic splines are used to fit the heavy mass dependence, asdescribed in the text. The fit evaluated at the physical pointand zero lattice spacing is given by purple dashed line witherror band. The experimental average value for the hyperfinesplitting [1] is plotted as the black cross at the physical Υmass. module [25] to do the fits, implementing the splines withthe gvar module [26].To obtain our final results, the fit function is evaluatedat lattice spacing equal to zero, sea quark masses tunedto their physical values and M φ h equal to M Υ . This istaken from experiment as 9.4603 GeV [1] with negligibleuncertainty. III. HYPERFINE SPLITTING
The hyperfine splitting, ∆ M hyp , is the difference inmass between the ground-state φ h and η h mesons. Thevalues for the hyperfine splitting on each ensemble for avariety of heavy quark masses are given in lattice unitsin column 5 of Table II. The separate η h and φ h massesare given in lattice units in columns 3 and 4. The impactof quenched QED at fixed valence quark mass is givenin Table III. The effect of QED is similar to that forcharmonium [9] but reduced because of the smaller elec-tric charge of the b quark. The direct effect of QED onthe hyperfine splitting is to increase it, through a QEDhyperfine effect which has the same sign as the QCDhyperfine effect. QED also increases the meson masses,however, and this requires a retuning of the bare quarkmasses downwards to match the same meson mass. Thisthen has an indirect effect, increasing the hyperfine split-ting by a very small amount.Our lattice results are plotted in Fig 1 as a functionof the vector heavyonium mass, M φ h . The points in-clude both pure QCD and QCD+QED values but theQCD+QED values are indistinguishable from pure QCDon this scale.To fit our results for the hyperfine splitting usingEq. (4) we take A = 0.1 GeV and the simple form for F , F = c (0) F + c (1) F (3 GeV) /M φ h . The coefficients c (0) F and c (1) F have prior values 0(1). Note that we multiplythe QED correction term in the fit (Eq. (7)) by a factorof 2 because of the size of the QED corrections that wesee in the results (Table III and [20]). That these priorwidths are very conservative can be judged from the val-ues and variation across Fig. 1. Evaluating the fit resultat zero lattice spacing, tuned quark masses and with M φ h equal to the Υ mass, we obtain the physical result for thebottomonium hyperfine splitting using connected corre-lation functions of: M Υ − M η b (connected) = 57 . .
3) MeV . (8)This is the QCD+QED value. For the ratio of theQCD+QED value to the pure QCD result, we obtain: R QED [∆ M hyp ] = 1 . . (9)Note that this is the ‘renormalised’ ratio with the b quarkmass tuned from the Υ in both QCD+QED and QCD.We see no significant impact of quenched QED at the0.2% level. The fit has a χ / dof of 0.73 using an SVDcut of 5 × − .Fig. 1 shows our fit curve as a function of M φ h inthe continuum limit for tuned sea quark masses. Thisgives useful physical insight into how the hyperfine split-ting falls as the quark mass increases. At the high massend of the plot we mark with a black cross the exper-imental average value [1] for the bottomonium system.We will discuss the comparison to experiment further inSection III A. We note that at the lower mass end ofthe curve we have results for charmonium. Our fit heredoes not include all of the charmonium results that wentinto [9] but gives a value for the charmonium hyperfinesplitting that is consistent (within 1 σ ) with [9] for thepure QCD case. The QED+QCD result here is too smallat the charmonium end of the fit curve because the QEDis being included with quark charge 1/3 e rather than thecorrect charm quark charge of 2/3 e .We will discuss in Section III A what the impact ofquark-line disconnected (but gluon-connected) correla-tion functions could be on the bottomonium hyperfinesplitting. For our charmonium calculation of [9] we in-cluded an estimate of the QED quark-line disconnectedcontribution to the hyperfine splitting coming from cc an-nihilation to a single photon, which then converts back to cc . The contribution of this to the charmonium hyperfinesplitting is 0.7 MeV, which was a little more than halfthe uncertainty in our result. The equivalent contribu-tion for the Υ here is much smaller, at 0.17 MeV, becauseof the smaller electric charge of the b quark. At a size ofone tenth of the uncertainty in our result in Eq. (8), thiswould then have negligible impact and we do not includeit. TABLE IV. Error budget for the hyperfine splitting and decayconstants as a percentage of the final answer. M Υ − M η b f Υ f η b statistics 2.40 0.77 0.38SVD cut 1.48 0.44 0.67 w w /a Z V - 0.29 - F G G G G G G G G A complete error budget for the bottomonium hyper-fine splitting is given in Table IV. Statistical uncertain-ties are divided between those arising from our 2-pointfits and those coming from the lattice spacing determina-tion, both correlated between ensembles ( w ) and uncor-related ( w /a ). The uncertainty from the 2-point fits isfurther divided in two. As already mentioned, the use ofan SVD cut is conservative and increases the uncertaintyin the fit output. We can calculate the contribution to anerror budget of both the data with and without the SVDcut applied to its correlation matrix. In the error bud-gets of Table IV we give the contribution from the datawith the original correlation matrix under the heading“statistics”. The additional contribution from the SVDcut is then defined as the square root of the difference ofthe squared contributions from the data with and with-out an SVD cut applied. The contributions from variousparts of the heavy mass dependence in Eqs. (4) and (7)are shown individually, labelled by the set of spline func-tions for that contribution.The fit parameters required to reproduce the physi-cal curve of the hyperfine splitting as a function of M φ h plotted in Fig. 1 are given in Table IX of Appendix A. A. Discussion: Hyperfine splitting
Our bottomonium hyperfine splitting result of Eq. 8is compared to earlier lattice QCD results in Fig. 2, go-ing back to the first lattice QCD calculation to includesea quarks [27]. Clearly the use of the heavy-HISQ ap-proach has allowed us to reduce the uncertainty signifi-cantly (by a factor of 3) relative to these earlier results.The earlier results all use non-relativistic actions, or ac-tions with non-relativistic input such as the Fermilab for-malism [36], for the b quarks. This leads to uncertaintiesfrom the normalisation of relativistic corrections to theHamiltonian, such as the σ · B term that is responsible
40 50 60 70 M Υ − M η b [MeV] HPQCD/UKQCD05FNAL/MILC09Meinel10RBC/UKQCD12HPQCD13/15This work FIG. 2. Comparison of lattice QCD determinations of thebottomonium hyperfine splitting. Our result from Eq. (8)is given by the top purple hexagon. Previous results (greensquares) come from: HPQCD/UKQCD using O ( v ) NRQCD b quarks and 2+1 flavours of asqtad sea quarks [27]; the Fermi-lab Lattice/MILC collaborations using the Fermilab formal-ism for the b quark and 2+1 flavours of asqtad sea quarks [28];S.Meinel using NRQCD b quarks with O ( v ) spin-dependentterms and 2+1 flavours of domain-wall sea quarks [29]; theRBC/UKQCD collaboration using the RHQ formalism for the b quark and 2+1 flavours of domain-wall sea quarks [30] andHPQCD using radiatively-improved NRQCD b quarks with O ( v ) spin-dependent terms and 2+1+1 flavours of HISQ seaquarks [31]. All of these results come from calculation ofconnected correlation functions and do not include an uncer-tainty from missing quark-line disconnected diagrams, exceptfor [31]. [31] includes the effect of these disconnected diagramsthrough the inclusion of 4-quark operators with coefficients,calculated in perturbation theory through O ( α s ). See thetext for discussion of the impact on the hyperfine splittingthrough η b annihilation to gluons. The red band is the PDGexperimental average [1]. The result for the hyperfine split-ting calculated here shows a clear improvement on previouslattice QCD results, as well as being the first to include QEDeffects. This improvement is in large part due to the elimina-tion of systematic uncertainties from the use of nonrelativisticactions which arise in previous calculations. for the hyperfine splitting. We avoid this uncertaintywith the HISQ action at the cost of having to calculateat multiple heavy quark masses rather than directly atthe b quark mass.As discussed in Sec. II, we have only computed con-nected correlators. This is also true for the earlier resultsexcept for that in [31]. This means that we are neglect-ing the contribution to the η b mass from its annihilationto gluons. This contribution can be related to the η b hadronic width using NRQCD perturbation theory [4]:∆ M η b = Γ η b (cid:18) − π + O ( α s , v /c ) (cid:19) . (10)Using the total width of the η b of 10(5) MeV [1] givesa shift to the η b mass from the leading order term of-1.0(5) MeV. This would result in an upward shift inthe hyperfine splitting of approximately 1 MeV, whichamounts to 0.5 σ for our result (Eq. (8)).
55 60 65 70 75 M Υ − M η b [MeV] BABAR02BABAR01CLEOBELLEThis work(lattice)This work PDG FIG. 3. Comparison of different experimental results for thebottomonium hyperfine splitting. The red band shows thePDG average of these experimental results [1]. The filledblue hexagon is our result (Eq. (11)) and is carried down-wards as the blue band. Note that our result here includesan uncertainty from the effect of η b annihilation missing fromour lattice calculation. There is some tension between the dif-ferent experimental results with our value favouring the mostrecent result from BELLE [32]. The result labelled CLEO isfrom [33], BABAR01 from [34] and BABAR02 from [35]. We recently showed, for the first time, that thisleading-order analysis fails in the case of the charmo-nium hyperfine splitting [9] where, with the improvedaccuracy we were able to achieve, it becomes clear thatthe lattice QCD+QED result is significantly higher thanthe experimental average. Assuming that this differenceis the result of the effect of η c annihilation missing fromthe lattice calculation, it seems that the leading-orderperturbative analysis is misleading in this case. Presum-ably missing higher-order terms in the perturbative anal-ysis or nonperturbative effects from mixing between the η c and other flavour-singlet pseudoscalar mesons [37], orboth combined, have a larger effect than the leading-order term and opposite sign. In the bottomonium casethe η b is considerably further from these lighter statesand so we may expect a much smaller effect from this.We also expect perturbation theory to be more reliable atthe higher energy associated with bottomonium states.We therefore allow an additional 1 MeV uncertainty forthe impact of η b annihilation on the hyperfine splittingand give a final result of M Υ − M η b = 57 . . .
0) MeV . (11)The first uncertainty is from the lattice calculation andthe second from missing quark-line disconnected contri-butions.The experimental average value for the bottomoniumhyperfine splitting (62.3 ± σ .Our result is also in agreement with the PDG averageto within 1.5 σ . We see no disagreement with the experi-mental result that would signal a larger contribution from η b annihilation than the 1 MeV that we have allowedabove. Indeed a shift upwards of our hyperfine split-ting result by 1 MeV, as suggested by leading-order per-turbation theory, would improve the agreement betweenlattice QCD and experiment, although the shift wouldnot be significant. In contrast, a shift downwards ofthe bottomonium hyperfine splitting by several MeV, aswe found for the charmonium hyperfine splitting, wouldcause tension with the experimental results.Finally we note that the high precision we are ableto achieve for the bottomonium hyperfine splitting is theresult of concentrating on the ground-state mesons with ahighly-improved relativistic action. For a more completepicture of the bottomonium spectrum, obtained on ananisotropic lattice with the Fermilab heavy quark actionand focussing on highly excited states see [38]. IV. Υ AND η b DECAY CONSTANTS
We define the vector heavyonium meson ( φ h ) decayconstant from the annihilation matrix element as (cid:104) | ψγ i ψ | φ h (cid:105) = f φ h M φ h (cid:15) i . (12)This means that we can determine the decay constantfrom our fits to the vector meson correlation functionsusing: af φ h Z V = (cid:115) A V E V , (13)where A V is the ground-state amplitude from a correla-tor fit of the form given in Eq. (2). Z V is the renormal-isation constant required to match the local vector cur-rent in lattice QCD to that of continuum QCD at eachvalue of the lattice spacing. We use Z V values calculatedin a nonperturbative implementation of the RI-SMOMscheme [14, 39, 40]. The pure QCD results for Z V forthe HISQ action are given in [9, 14]; we use values atscale µ = 2 GeV. Note that no additional matching fac-tor is required to reach MS from the RI-SMOM schemeand, because Z V has no anomalous dimensions, any µ dependence is purely a discretisation effect [14].The vector meson decay constant is the amplitude forannihilation of the valence quark/antiquark pair, into aphoton, for example. It is related to the experimentallymeasurable leptonic width by:Γ( φ h → e + e − ) = 4 π α e h f φ h M φ h (14)where e h is the quark electric charge (1/3 for b ). The α QED here is evaluated at the mass of the heavy quarkand is equal to 1/132.15 [41] at the b . We also compute the decay constant of the pseu-doscalar heavyonium meson, f η h . In terms of the pa-rameters of our correlator fit, Eq. (1) this is defined as: f η h = 2 m h (cid:115) A P ( E P ) . (15)Because the partially conserved axial current (PCAC)relation holds for HISQ quarks the pseudoscalar decayconstant is absolutely normalised and no Z factor is re-quired to match to the continuum regularisation of QCD.Since the pseudoscalar meson does not annihilate to asingle particle, there is no experimental decay processthat gives direct acces to the decay constant. Its value isnevertheless of interest for comparison to that of the cor-responding vector meson and other pseudoscalar mesons.The values of the decay constants, in lattice units, oneach ensemble and for each heavy mass are given in thesixth and seventh columns of Table II. The decay con-stants converted to GeV units, and renormalised in thecase of the vector decay constant, are plotted as a func-tion of the φ h mass in Fig. 4. The decay constants in-crease with increasing φ h mass. Discretisation effectsare clearly visible that cause the lattice results to peelaway from the physical curve upwards. The same effectwas seen previously for both heavy-light and heavyoniummesons [10, 11].We also show results in Fig. 4 that include the effect ofquenched QED. Those results are given in Table III as theratio of values in QCD+QED to those in pure QCD. Forthe decay constant of the φ h these ratios do not includethe impact of QED on the vector current renormalisationfactor, Z V . This was calculated in [14] for the case of aquark with electric charge 2 e/
3, again as a ratio of resultsin QCD+QED to those in pure QCD. These results aregiven in Table IV of [14], with further results in TableX of [9]. The ratio is within 0.05% of 1, as expected foran O ( α QED ) correction to a Z factor that is already veryclose to 1 for the HISQ action in pure QCD. Here we needresults for an electric charge of e/ µ = 2 GeV) and dividing thedifference from 1 by a factor of four.To fit our decay constant results as a function of latticespacing and heavy quark mass we again use the fit formof Eq. (4) but we use A = 0.7 GeV, as appropriate for thedecay constant values, and a different form for F to thatused for the hyperfine splitting case. The dependence ofthe decay constants on the heavy mass is approximatelylinear and so we choose F = c (0) F + c (1) F M φ h / (3 GeV),where c (0) F and c (1) F are fit parameters with prior valuesof 0 ±
1. We fit the φ h and η h decay constants simul-taneously, including the correlations between them andtake the same F for both since they are so close in value.The spline functions that map out the differences from F in physical heavy quark mass dependence and the de-pendence on the lattice spacing and sea quark massestake independent values in the two cases. The fit has a χ / dof value of 0.44 using an SVD cut of 1 × − . Weagain evaluate our fits at zero lattice spacing, physicalsea quark masses and with M φ h = M Υ to obtain thephysical bottomonium results.We obtain, for the Υ, f Υ = 677 . .
7) MeV (16)with R QED [ f Υ ] = 1 . . (17)For the η b , f η b = 724(12) MeV , (18)with R QED [ f η b ] = 1 . . (19)Again, QED effects are not discernible within our 0.1%uncertainties. At the charmonium end of our range ourresults agree within uncertainties with the values we ob-tained in [9], remembering that the calculation done hereis for an electric charge that does not match that of the c quark. The error budget for both decay constants isgiven in Table IV. The fit curves evaluated at zero lat-ting spacing and physical sea quark masses are plotted asa function of heavy quark mass (given by M φ h ) in Fig. 4.The fit parameters required to reproduce these physicalcurves of the decay constants as a function of M φ h aregiven in Table X of Appendix A.Given that the heavy mass dependence and discreti-sation effects in the vector and pseudoscalar decay con-stants are similar we can study the ratio of the two as afunction of the heavy mass to high precision. Our resultsfor the ratio are shown as a function of M φ h in Fig. 5. Aslow downward drift of the ratio is seen with increasing M φ h from a value slightly above 1 for c quarks to a valueslightly below 1 for b quarks.To obtain a physical result for the ratio we again usethe fit form of Eq. (4), now taking F to be a constant, c F , since the ratio is relatively flat, so that the splinefunctions handle all of the mass dependence. We take theprior value of c F to be 1(1), i.e. with a very conservativewidth. Since we expect a lot of systematic effects tocancel in this ratio (and Fig. 5 shows that they do) wehalve the prior widths on all of the correction terms inEq. 4 i.e. we take prior values on the function values atthe knots of 0.0(5). The fit has a χ / dof of 0.22 and noSVD cut is required. Evaluating the fit function at thephysical point gives f Υ f η b = 0 . , (20)and R QED (cid:20) f Υ f η b (cid:21) = 0 . . (21) M φ h [GeV]0 . . . . f φ h [ G e V ] f-5f-physsf-5 sf-physuf-5ef-5 M φ h [GeV]0 . . . . f η h [ G e V ] f-5f-physsf-5 sf-physuf-5ef-5 FIG. 4. Upper panel: The φ h decay constant plotted againstthe φ h mass. The symbols correspond to different gluon fieldensembles, as given in the legend (see Table I for a list). Pointsincluding quenched QED are shown in cyan, indistinguishablefrom pure QCD points underneath. The dashed line and er-ror band show the fit described in the text evaluated at zerolattice spacing and physical sea quark masses. Lower panel:The η h decay constant plotted against the φ h mass, symbolsand fit line as above. The total uncertainty in the ratio for the b is 1%, witha value clearly below 1. The fit curve evaluated at zerolattice spacing and physical sea quark masses is plotted asa function of M φ h in Fig. 5. The fit parameters requiredto reproduce this physical curve are given in Table XI ofAppendix A. A. Discussion : Decay constants
Figure 6 compares our result for the Υ decay constant, f Υ , to that of an earlier lattice QCD calculation on a sub-set of the same gluon field configurations used here butusing an improved NRQCD action for the b quarks [13].Clearly, we achieve a considerably improved uncertaintyover that of [13]. A large amount of the NRQCD uncer-tainty arises from the normalisation and O ( v /c ) im-provement of the NRQCD vector current, where v isthe nonrelativistic quark velocity. Here, since we usethe HISQ action which is relativistic and we have per-formed the vector current renormalisation to very highprecision previously [14], these sources of uncertainty areeffectively eliminated.0 M φ h [GeV]0 . . . . f φ h / f η h f-5f-physsf-5 sf-physuf-5ef-5 FIG. 5. The ratio of the vector to pseudoscalar heavyoniumdecay constants as a function of vector heavyonium mass. Atthe charmonium point the ratio is above 1. By the bottomo-nium point the ratio has shifted to be below 1. The symbolscorrespond to results on the different gluon field configura-tions listed in the legend with cyan points corresponding toQCD+QED. The line is the fit curve evaluated at the physicalpoint as a function of M φ h described in the text.
620 640 660 680 f Υ [MeV] PDGHPQCD (latticeNRQCD)This work(lattice) FIG. 6. A comparison of our result (filled blue hexagon) forthe decay constant of the Υ with HPQCD’s earlier latticeQCD result using NRQCD b quarks [13] (open green square).We also include the value inferred from the experimental lep-tonic decay width in Eq. (22) (pink open circle). Figure 6 also compares our result for f Υ to that ob-tained from the experimental average for the Υ leptonicwidth using Eq. (14). Using Γ(Υ → e + e − ) = 1.340(18)keV [1] gives f exptΥ = 689 . . .
8) MeV (22)where the first uncertainty comes from the experimentaluncertainty in Γ and the second allows for an O ( α QED /π )uncertainty from higher-order in QED terms in Eq. (14)coming, for example, from final-state radiation. Notethat using α QED of 1/137 here instead of 1/132.15 wouldincrease the experimental value of f Υ by 3.7% or 25 MeV.This is several times larger than either the experimentaluncertainty or our lattice QCD uncertainty.Figure 6 shows good agreement, within 1 σ , be-tween our lattice QCD result and that from experiment(eq. (22)). The experimental uncertainty is about half that from our lattice QCD result.Our result for f Υ can be converted into a determinationof the width for Υ decay to light leptons in the StandardModel using Eq. (14). This givesΓ(Υ → e + e − ) = 1 . O ( α QED /π )correction to Eq. (14) from higher-order QED effects.Our result for f η b can be compared to an earlierHPQCD lattice QCD result using HISQ quarks andthe heavy-HISQ approach on gluon field configura-tions including the effect of 2+1 flavours of asqtad seaquarks [11]. That work obtained a value f η b = 667(6)MeV, which is significantly lower (by 4 σ ) than our re-sult here. The discrepancy is most likely to result froma bias in the earlier results from not having values onlattices with spacings as fine as we do here. Anotherpossible source of the discrepancy is the fact that theearlier calculation did not include c quarks in the sea.Having more flavours of quarks in the sea results in aslower running of the strong coupling constant. Hence,using the language of potential models, we expect theCoulomb-like term in the heavy quark potential (of theform − α s ( r ) / (3 r )) to have a larger value for α s at theshort-distance scales to which the η b meson decay con-stant is sensitive. This corresponds to a deeper po-tential at short distances and a correspondingly larger‘wavefunction-at-the-origin’, which is the quantity in apotential model that translates approximately into thedecay constant. This effect could explain some of thediscrepancy but is unlikely to be large enough to explainit all. The calculations in [11] also used a different formto fit the lattice results as a function of heavy quark mass(in that case using as proxy M η h ). This consisted of mul-tiple powers of the inverse heavy quark mass multipliedby a leading function of the form ( M/M ) b where b wasallowed to float. We have checked that using that fitform here gives us results for f η b very consistent withour spline fits, so the discrepancy with [11] is not relatedto the form of the fit used.The ratio of vector to pseudoscalar decay constantsas a function of heavyonium mass provides a test of ourunderstanding of these mesons. In the language of po-tential models the heavyonium vector and pseudoscalarmesons differ only through spin-dependent relativisticcorrections to the central potential [42]. The size of rela-tivistic corrections fall as the heavy quark mass increasesand the mean squared-velocity of the heavy quarks fall.In the infinite quark mass limit pseudoscalar and vectorheavyonium mesons have the same mass and the samewavefunction-at-the-origin. The decay constants differ,however, by the matching factors that are needed torenormalise temporal axial and spatial vector currentsfrom this nonrelativistic framework to full continuumQCD. The ratio of the vector to pseudoscalar heavyo-nium decay constants would then be expected to becomethe ratio of the vector to temporal axial vector matching1factors in the heavy quark limit. The matching factorscome from high-momentum regions of phase-space andso can be calculated in QCD perturbation theory. An O ( α s ) matching calculation was done in [43] for spin-independent nonrelativistic QCD and gave the result Z V Z A = 1 − g π = 1 − α s π . (24)From this we conclude that we would expect the ratio ofvector to pseudoscalar decay constants to be below 1 forlarge heavy quark mass. Eq. (24) expects the differencefrom 1 to be O (5%), taking α s ≈ .
25, but this formulawill have corrections from higher orders in α s . A valuefor the ratio of 5% below 1 is very consistent with ourresults in Fig. 5, however.Very similar behaviour is seen for the ratio of vectorto pseudoscalar decay constants for heavy-light mesonsfrom lattice QCD calculations. The decay constant ofthe D ∗ s meson is found to be several percent larger thanthat of the D s [44–46] whereas that of the B ∗ s is a fewpercent below that of the B s [47] [45]. This behaviourcan be understood on the same basis as the argumentsfor heavyonium above. In the heavy-light case an α s calculation of the matching factors is available in the in-finite heavy quark mass limit [48]. The corrections tothe O ( α s ) formula for the ratio (which is the same as forheavyonium in Eq. (24)) are sizeable but have the same(negative) sign and so do not change the qualitative be-haviour of the difference of the ratio from 1.Having performed the fits of the previous subsectionswe now have physical values for the decay constants notonly at the b quark mass but at the full range of massesbetween the c and b quark masses. The physical curves asa function of meson mass in Figs. 4 and 5 could be used totune phenomenological QCD potential models, that oftendiffer markedly on features of heavyonium physics suchas details of the wavefunction even when reproducing thespectrum (see, for example [50–52]).They may also be useful beyond QCD. In [49] lat-tice QCD results across a range of masses were collectedwith the intention of providing useful information forphenomenologists studying strongly coupled beyond theStandard Model (BSM) theories. These theories are of-ten QCD-like but typically with heavier (relative to theconfinement scale) fundamental fermions than the lightquarks of QCD. Ref. [49] makes the point that informa-tion from lattice QCD calculations about how (for exam-ple) meson masses and decay constants depend on quarkmasses can be useful to constrain such BSM theories.This then requires lattice QCD results for quark massesnot at their physical values, as we have here. The latticeQCD results need to be presented in an appropriate waywith dimensionless combinations of decay constants andmasses on both axes. A convenient x -axis is the ratio ofpseudoscalar to vector meson mass. In [49] the squareof this quantity was used since the lattice QCD resultswere concentrated at light quark masses. Here, since wehave heavy quarks and the ratio of pseudoscalar to vector .
97 0 .
98 0 . M η h /M φ h M η h / f η h .
96 0 .
97 0 .
98 0 . M η h /M φ h . . . f φ h / M φ h FIG. 7. Upper plot: The physical ratio (evaluated at zero lat-tice spacing and with sea quark mass mistunings set to zero)of the mass to decay constant for the pseudoscalar heavyo-nium meson, η h as a function of the ratio of pseudoscalar tovector heavyonium masses. Lower plot: The physical ratio ofdecay constant to mass for the vector heavyonium meson, φ h ,(the quantity denoted f V in [49]) plotted against the sameratio of masses. meson masses is close to 1, we simply use the ratio.Dimensionless ratios are readily obtained for our rawlattice results using the values in Table II. Correlationscan be ignored because statistical uncertainties are sosmall. In the following we construct appropriate ratiosfrom our fit functions in the limit of zero lattice spacingand physical sea quark masses and do not include theraw lattice results in the figures, for clarity.One useful quantity [53] is the ratio of the pseudoscalarmeson mass and decay constant for a meson made ofquarks of degenerate mass (i.e. the ‘pions’ of the BSMmodel). Using the physical heavy mass dependence of f η h extracted from our fit we display the ratio of M η h and f η h as a function of the ratio of pseudoscalar to vec-tor meson masses in Fig. 7. Our results show values of M η h /f η h around 10, and continuing to rise, as the ratioof pseudoscalar to vector meson masses heads towards 1.Note that our definition of the pseudoscalar decay con-stant in Eq. (15) corresponds to the normalisation f π ≈
130 MeV.As discussed in [49] composite models of a dark sec-tor in which a ‘dark ρ ’ meson couples to ordinary mat-ter through a dark photon (e.g. [54]) need informationon the vector meson decay constant for an appropriate2range of fermion masses. The ratio of vector meson de-cay constant to vector meson mass is denoted f V in [49].In our convention for the vector meson decay constant(Eq. (12)) it is f φ h /M φ h . We plot f φ h /M φ h against thepseudoscalar to vector meson mass ratio in Fig. 7. Wesee that this ratio becomes small as the ratio of pseu-doscalar to vector meson masses heads towards 1. It alsohas relatively strong dependence on the mass ratio, sousing an approximately constant value (based, for exam-ple, on naive dimensional analysis) would not agree wellwith our results.We also note that accurate lattice QCD results areavailable at the ratio of pseudoscalar to vector mesonmasses of 0.673 which corresponds to ss mesons whenonly connected correlation functions are calculated. Thismeans that the pseudoscalar meson is not allowed to an-nihilate and mix with flavour-singlet mesons made fromlighter quarks, and likewise the vector decay to two pseu-doscalar mesons incorporating lighter quarks is not in-cluded. This is then the scenario that would matchthat required in a composite BSM scenario. For thiscase HPQCD calculates M η s /f η s = 3.801(16) [19] and f φ s /M φ s = 0.233(3) [55]. These results must connectsmoothly to the ones shown in Fig. 7 as the quark massis increased. V. VECTOR CURRENT-CURRENTCORRELATOR TIME-MOMENTS AND a bµ The ground-state vector heavyonium decay constant isdetermined by the amplitude of the state that dominatesthe correlator at large times and this can be connectedto experiment via the leptonic width, as we have seen.We can also calculate the time moments of the correla-tor. These depend on the behaviour of the correlator atshorter time distances and can also be connected to ex-perimental results [8, 56]. The moments of the vectorheavyonium correlator are defined by: G Vn = Z V (cid:88) ˜ t ˜ t n C φ h (˜ t ) (25)where ˜ t is lattice time symmetrised around the centre ofthe lattice, C φ h is the vector two-point correlation func-tion and Z V is the renormalisation factor for the heavy-onium vector current operator used.Results for ( G Vn /Z V ) / ( n − in lattice units on each ofour ensembles are given in Table V for n = 4 to 10. Thepower 1 / ( n −
2) is taken to reduce all the moments to thesame dimension. We take the Z V factor for the vectorcurrent to be the same one used for the leptonic widthabove [14]. Figure 8 then shows the physical results forthese moments as a function of M φ h .Results that include quenched QED corrections for asubset of ensembles are given in Table VI. These are givenas the ratio of the result in QCD+QED to that in QCDat fixed valence quark mass in lattice units. The valuesof R are very slightly below 1, as for charmonium [9]. M φ h [GeV]0 . . . ( t h m o m e n t) / [ G e V − ] f-5f-physsf-5 sf-physuf-5ef-5 M φ h [GeV]0 . . . ( t h m o m e n t) / [ G e V − ] f-5f-physsf-5 sf-physuf-5ef-5 M φ h [GeV]0 . . . . ( t h m o m e n t) / [ G e V − ] f-5f-physsf-5 sf-physuf-5ef-5 M φ h [GeV]0 . . . . ( t h m o m e n t) / [ G e V − ] f-5f-physsf-5 sf-physuf-5ef-5 FIG. 8. Results for the 4th, 6th, 8th and 10th time momentsof the heavyonium vector correlator plotted as a function of M φ h . The symbols correspond to different gluon field en-sembles, as given in the legend (see Table I for a list). Theerrors on the points are dominated by uncertainties from thedetermination of Z V that are correlated between the points.Points including quenched QED are shown in cyan, indistin-guishable from pure QCD points underneath. The dashed linewith purple error band displays our continuum/chiral fit, asdiscussed in the text. Values determined from experimentalresults for R e + e − (eq. (26)) are plotted as the black crossesat M φ h = M Υ [57]. The difference from 1 is even smaller here because ofthe smaller quark electric charge. Note that the vector3
TABLE V. Results in lattice units for time moments of the vector heavyonium correlator as defined in Eq. (25). We giveraw results here in which the vector current has not been renormalised and we also take the ( n − G Vn /Z V ) / ( n − for n =4 to 10.Set am h n = 4 n = 6 n = 8 n = 101 0.6 0.562768(11) 1.263877(18) 1.849015(25) 2.386587(33)0.8 0.4342507(59) 1.0316257(96) 1.525578(12) 1.967383(15)2 0.6 0.5628101(56) 1.2640282(83) 1.849341(10) 2.387158(12)0.8 0.4342571(32) 1.0316599(48) 1.5256620(57) 1.9675445(65)0.866 0.395358(44) 0.966054(54) 1.439945(54) 1.861205(52)3 0.274 1.070712(58) 2.27651(10) 3.35545(14) 4.37419(17)0.4 0.797940(92) 1.72744(16) 2.54252(21) 3.31432(25)0.5 0.665034(57) 1.466318(96) 2.15508(13) 2.80305(16)0.548 0.60993(22) 1.36481(25) 2.00869(24) 2.61113(24)0.6 0.569093(37) 1.282926(63) 1.886184(84) 2.44653(10)0.7 0.495371(26) 1.145651(43) 1.689190(58) 2.186270(72)0.8 0.436182(18) 1.037685(31) 1.538234(41) 1.989299(51)4 0.260 1.114660(44) 2.366266(78) 3.48827(11) 4.54699(14)0.4 0.798236(18) 1.728246(32) 2.544009(44) 3.316608(55)0.6 0.5691755(75) 1.283151(13) 1.886612(17) 2.447216(22)0.8 0.4362111(38) 1.0377647(63) 1.5383891(83) 1.989559(11)5 0.194 1.431378(91) 3.03675(16) 4.49434(22) 5.86769(29)0.4 0.808461(20) 1.757499(36) 2.597493(51) 3.396810(65)0.6 0.5722526(77) 1.292481(13) 1.904940(19) 2.476827(24)0.8 0.4371741(39) 1.0407975(65) 1.5447616(86) 2.000643(11)0.9 0.3876830(29) 0.9510722(49) 1.4215906(63) 1.8418824(80)6 0.138 1.91475(23) 4.06357(42) 6.02429(55) 7.86806(66)0.45 0.739093(19) 1.623965(33) 2.403967(45) 3.147423(57)0.55 0.621096(13) 1.389882(21) 2.052994(29) 2.679049(36)0.65 0.5342576(87) 1.222777(15) 1.806579(20) 2.349544(24)TABLE VI. Quenched QED corrections, for quark electric charge e/
3, to the time-moments given for a subset of the results inTable V for the ( n − G Vn /Z V . The results are given as the ratio, R , of the value in QCD+QEDto that in pure QCD at fixed valence quark mass in lattice units.Set am h n = 4 n = 6 n = 8 n = 102 0.866 0.999669(32) 0.999731(16) 0.999697(11) 0.999641(8)3 0.274 0.999774(30) 0.999692(25) 0.999646(24) 0.999622(24)0.548 0.999475(43) 0.999353(22) 0.999194(14) 0.999060(11) current is not renormalised in these raw results and QEDeffects in Z V must also be taken into account, as for thedecay constant [14]. These results are also plotted inFig. 8 as the cyan points. The impact of QED is notvisible.To fit the time-moment results as a function of latticespacing and heavy quark mass we again use the fit ofEq. (4), supplemented with QED effects in Eq. (7). Forthe time-moments we use F = ( c (0) F + c (1) F (3 GeV) /M φ h ),and the dimensionful parameter A is taken as 0.5 GeV − for every moment. The prior values on c (0) F and c (1) F aretaken to be 0 ± × − in all cases.The χ / dof of the fits, in order of increasing n , are 0.9,0.19, 0.26 and 0.4. The curves in Fig. 8 show the fitresults evaluated at zero lattice spacing and with tunedsea quark masses.Table VII gives our results for the time-moments eval-uated at the b quark mass in the continuum limit, with their total uncertainties. The corresponding error bud-get is given in Table VIII. In the next section we comparethese results to earlier lattice analyses and values deter-mined from experimental data for R ( e + e − → hadrons).We will also use the results to improve the determinationof the b quark contribution to the hadronic vacuum po-larisation term in the Standard Model determination ofthe anomalous magnetic moment of the muon. Column3 of Table VII gives the ratio of the QCD+QED resultto that in pure QCD for each moment. Again we arenot able to distinguish any impact of QED on the resultsat the level of our uncertainties (which range from 0.4%down to 0.1%). A. Discussion: time moments and a bµ In Fig. 9 we compare our results for the time-momentsto those of an earlier HPQCD calculation that used4
TABLE VII. Results for the time moments of the bottomo-nium vector current-current correlator obtained from evalu-ating our fit functions in the continuum limit at the b quarkmass. These are given in the second column for moment num-bers listed in the first column. The results extracted fromexperimental data in [57] are given in the third column forcomparison. The fourth column gives the quenched QED cor-rection to these moments, as a ratio of the value in QCD plusQED to that in pure QCD with a tuned b quark mass (toreproduce the Υ mass from experiment) in both cases. All ofthe ratios are consistent with 1.0. n G / ( n − n ( G exp .n ) / ( n − R QED (cid:104) G / ( n − n (cid:105) [GeV − ] [GeV − ]4 0.0905(23) 0.09151(31) 0.9996(38)6 0.1920(39) 0.19910(49) 0.9999(19)8 0.2934(55) 0.29964(55) 0.9999(13)10 0.3918(66) 0.39548(59) 0.9999(10)TABLE VIII. Error budget for the n th time-moment, G / ( n − n , as a percentage of the final answer. n w w /a Z V F G G G G G G G G . . . . n th moment) / ( n − [GeV − ] ExperimentalHPQCD NRQCD[1408.5768]This work n = 4 n = 6 n = 8 n = 10 FIG. 9. Comparison of different determinations of the fourlowest time moments of the bottomonium vector current-current correlator. The three determinations are, from thetop, this work, previous calculation by the HPQCD collabo-ration using NRQCD b quarks [13] and the values obtainedfrom experimental data on R ( e + e − → hadrons) in [57]. NRQCD b quarks [13]. For the NRQCD results the keysources of error were from the vector current normal-isation (using a method based on matching the time-moments to continuum perturbation theory) and fromthe lattice spacing dependence effects in the NRQCD ac-tion. Our uncertainties here are a considerable improve-ment (by over a factor of two) on the NRQCD results,because we have a very accurate vector current normali-sation and have results over a large range of lattice spac-ing values to control the lattice spacing dependence.Figure 9 shows that our results agree within 2 σ withthe values extracted for the q -derivative moments, M k ( n = 2 k + 2), of the b quark vacuum polarizationusing experimental values for R e + e − = σ ( e + e − → hadrons) /σ pt [57]. The appropriate normalisation ofthese results for the comparison to ours, is: G exp n = (cid:18) M exp k n !12 π e b (cid:19) / ( n − (26)Our results from lattice QCD have considerably largeruncertainties than those of the experimental values buttogether these results provide a further test of QCD atthe level of 2%.We may also use these time-moments to extract the b quark connected contribution to the leading orderhadronic vacuum polarisation contribution to the anoma-lous magnetic moment of the muon. This was done in [13]and, given that we have improved on the time-momentsof that work, we provide an update here. We obtain a bµ = 0 . × − (27)This agrees with the value in [13] with an improvementin uncertainty of a factor of 2.5. Since the b quark is soheavy, this is not a significant contribution to the anoma-lous magnetic moment of the muon [58]. VI. CONCLUSIONS
We have used the fully relativistic HISQ action to cal-culate the masses and decay constants of ground-statebottomonium mesons in lattice QCD including the ef-fects of u , d , s and c quarks in the sea. We have usedvery fine lattices and a range of heavy quark masses ateach lattice spacing to control the discretisation effects asa function of heavy quark mass along with the physicaldependence on the heavy quark mass of the quantitiesbeing studied. We have used a fit function with com-pletely generic dependence on the heavy quark mass ineach of its component pieces, capturing this dependencethrough cubic spline functions. Values for bottomoniumare obtained by evaluating the fit function at zero latticespacing with tuned sea quark masses and a valence quarkmass tuned to that of the b , defined to be the point at thewhich the Υ mass agrees with experiment. We have alsoincluded an analysis of the impact of the electric chargeof the valence b quarks on the quantities being studied.5The results given are from the QCD+QED fit but, in allcases, we find the impact of QED to be negligible at thelevel of our uncertainties.Our results yield the most precise, to date, lattice cal-culation of the bottomonium hyperfine splitting. We ob-tain the value (repeating Eq. (11)): M Υ − M η b = 57 . . .
0) MeV . (28)The first uncertainty is from our fit results (see errorbudget in Table IV) and the second uncertainty is froman estimate of missing quark-line disconnected contribu-tions that would affect the mass of the η b meson. Ourresult is in agreement with, but on the low side of, theexperimental average value [1]. It tends to favour themost recent experimental result obtained by the BELLEcollaboration [32], although uncertainties (both ours andfrom the experiment) are still too large to draw strongconclusions from this.We also provide the most precise lattice QCD deter-mination of the Υ decay constant, which can be used todetermine the Υ leptonic width. Our uncertainty of 1.5%is three times better than the previous lattice QCD cal-culation of [13]. The big advantage of using a relativisticformalism, as we do here, is that the vector current can benormalised very accurately and nonperturbatively [14].Our result (repeating Eq. (16)) is f Υ = 677 . .
7) MeV , (29)with error budget in Table IV. Using this result to obtainthe Υ leptonic width gives (repeating Eq. (23)):Γ(Υ → e + e − ) = 1 . . (30)The first uncertainty is from our result for f Υ and the sec-ond from possible O ( α QED /π ) corrections to the formulaconnecting decay constant and leptonic width (Eq. (14)).This is to be compared with the current experimental av-erage of 1.340(18) keV [1]. We see that our result is ingood agreement with experiment and our uncertainty isjust twice as large.The decay constant of the η b can also be accuratelycalculated with our approach. There is no experimen-tal decay rate that can be directly compared to this de-termination, but the value of f η b is important for ourphenomenological understanding of the relationships be-tween decay constants for different mesons. We obtain(repeating Eq. (18)) f η b = 724(12) MeV . (31)In particular, repeating Eq. (20), we find that f Υ f η b = 0 . , (32)i.e. less than 1. This is in contrast to the charmoniumcase where f J/ψ /f η c is larger than 1 [9]. Fig. 5 showshow the ratio of the decay constants for vector and pseu-doscalar heavyonium mesons varies with heavy quark mass. This is qualitatively similar to the behaviour seenfor the decay constants of heavy-light mesons [47]. Fi-nally, in Fig. (7) we plot the ratios of mass to decay con-stant for pseudoscalar and vector mesons as a function ofthe ratio of pseudoscalar to vector meson masses. Thesemay provide useful information to constrain these ratiosin QCD-like beyond the Standard Model scenarios.The low time moments of the bottomonium vectorcurrent-current correlator provide a further opportunityto compare lattice QCD results to experiment, where thematching inverse- s moments of the b -quark contributionto R ( e + e − → hadrons) can be determined. Our resultsfor the 4 th , 6 th , 8 th and 10 th time moments are givenin Table VII where they can be compared to the resultsobtained from experiment. Our uncertainties are 2% soprovide the most accurate test to date for these quanti-ties. The time moments can be used to determine the b quark contribution to the anomalous magnetic momentof the muon. We find (repeating Eq. (27)) a bµ = 0 . × − . (33)Together these results demonstrate how the propertiesof low-lying bottomonium states can be determined in afully relativistic calculation in lattice QCD and the gainsin precision that such an approach makes possible. Theresults given here also allow us to improve the fully non-perturbative determination of the ratio of quark masses, m b to m c . We will present this analysis separately. Acknowledgements
We are grateful to the MILC collaboration for theuse of their gluon field configurations and for the use ofMILC’s QCD code. We have modified the code to gen-erate quenched U(1) gauge fields and incorporate thoseinto the quark propagator calculation as described here.We are grateful to B. Galloway for contributions to thisproject at a very early stage, and to R. Horgan, C. Mc-Neile and J. Rosner for useful discussions. Computingwas done on the Darwin supercomputer at the Universityof Cambridge High Performance Computing Service aspart of the DiRAC facility, jointly funded by the Scienceand Technology Facilities Council, the Large FacilitiesCapital Fund of BIS and the Universities of Cambridgeand Glasgow. We are grateful to the Darwin support stafffor assistance. Funding for this work came from the Sci-ence and Technology Facilities Council and the NationalScience Foundation.
Appendix A: Reconstructing the heavy quark massdependence
We give here the fit parameters that enable our fit re-sults for the dependence on heavy quark mass of the hy-perfine splitting, decay constants and ratio of decay con-stants to be reconstructed. The pieces of Eq. (4) that give6
TABLE IX. Fit parameters for F and G for the fit of Eq. (4)to the hyperfine splitting as a function of the vector heavyo-nium mass, M φ h . The dimensionful constant, A is 0.1 GeVin this case and F = c (0) F + c (1) F × /M φ h . The meanand standard deviation for c (0) F and c (1) F and the values atthe 3 knot positions for G are: c (0) F = 0 . c (1) F =0 . G k = 0 . G k = 0 . G k = − . c (0) F c (1) F G k G k G k F and G for the fit of Eq. (4)to the vector decay constant of the vector heavyonium φ h meson (upper set) and the decay constant of the pseudoscalarheavyonium η h meson (lower set), both as a function of themass M φ h . The dimensionful constant, A is 0.7 GeV in thesecases and F = c (0) F + c (1) F × M φ h / c (0) F and c (1) F and the values at the 3 knot positions for G . The correlationmatrix for these 5 parameters is given underneath. Theseresults enable the red fit curves of both plots of Figure 4 tobe reconstructed within the errors given. f φ h c (0) F c (1) F G k G k G k f η h c (0) F c (1) F G k G k G k the physical curves in the continuum limit are F ( M φ h )and G (1 /M φ h ), multiplied by dimensionful constant, A (absent for the case of the ratio of decay constants). Weignore here the QED pieces of the fit; these have negligi-ble effect in all cases. F is a simple function with at most two parameters, c (0) F and c (1) F . G is a Steffen spline function [23] with 3knots at 2.5, 4.9 and 10.0 GeV in M φ h , so that in 1 /M φ h they are at 1/2.5, 1/4.9 and 1/10.0. Tables IX, X and XIgive the mean and standard deviation of c (0) F , c (1) F , and thevalues at the 3 knots of G : G k , G k and G k . This isfollowed underneath by the correlation matrix betweenthese parameters. The parameters are strongly corre-lated and this is why we give the values to 4 significantfigures. The splines can easily be implemented using the gvar Python module [26].
TABLE XI. Fit parameters for F and G for the fit of Eq. (4)to the ratio of vector to pseudoscalar heavyonium decay con-stants as a function of the mass of the vector heavonium me-son, M φ h . F is simply a constant, c (0) F in this case. The toprow gives the mean and standard deviation for c (0) F and thevalues at the 3 knot positions for G . The correlation matrixfor these 4 parameters is given underneath. These results en-able the red fit curve of Figure 5 to be reconstructed withinthe errors given. c (0) F G k G k G k [1] P. Zyla et al. (Particle Data Group), PTEP ,083C01 (2020).[2] E. McLean, C. Davies, J. Koponen, and A. Lytle, Phys.Rev. D , 074513 (2020), arXiv:1906.00701 [hep-lat].[3] J. Harrison, C. T. H. Davies, and A. Lytle (HPQCD),Phys. Rev. D , 094518 (2020), arXiv:2007.06957 [hep-lat].[4] E. Follana et al. (HPQCD Collaboration), Phys.Rev. D75 , 054502 (2007), arXiv:hep-lat/0610092 [hep-lat].[5] E. Follana, C. Davies, G. Lepage, and J. Shigemitsu(HPQCD Collaboration), Phys.Rev.Lett. , 062002(2008), arXiv:0706.1726 [hep-lat].[6] C. Davies, C. McNeile, E. Follana, G. Lepage, H. Na, et al. (HPQCD Collaboration), Phys.Rev.
D82 , 114504(2010), arXiv:1008.4018 [hep-lat].[7] H. Na, C. T. Davies, E. Follana, G. P. Lepage, andJ. Shigemitsu (HPQCD Collaboration), Phys.Rev.
D82 ,114506 (2010), arXiv:1008.4562 [hep-lat].[8] G. Donald, C. Davies, R. Dowdall, E. Follana, K. Horn-bostel, et al. (HPQCD collaboration), Phys.Rev.
D86 ,094501 (2012), arXiv:1208.2855 [hep-lat].[9] D. Hatton, C. Davies, B. Galloway, J. Koponen, G. Lep-age, and A. Lytle (HPQCD), Phys. Rev. D , 054511(2020), arXiv:2005.01845 [hep-lat].[10] C. McNeile, C. Davies, E. Follana, K. Hornbostel, andG. Lepage (HPQCD Collaboration), Phys.Rev.
D85 ,031503 (2012), arXiv:1110.4510 [hep-lat].[11] C. McNeile, C. Davies, E. Follana, K. Hornbostel,and G. Lepage, Phys. Rev. D , 074503 (2012),arXiv:1207.0994 [hep-lat].[12] A. Bazavov et al. , Phys. Rev. D , 074512 (2018),arXiv:1712.09262 [hep-lat].[13] B. Colquhoun, R. Dowdall, C. Davies, K. Hornbostel,and G. Lepage, Phys. Rev. D , 074514 (2015),arXiv:1408.5768 [hep-lat].[14] D. Hatton, C. Davies, G. Lepage, and A. Ly-tle (HPQCD), Phys. Rev. D , 114513 (2019),arXiv:1909.00756 [hep-lat].[15] A. Bazavov et al. (MILC), Phys. Rev. D87 , 054505(2013), arXiv:1212.4768 [hep-lat].[16] A. Hart, G. M. von Hippel, and R. R. Horgan (HPQCD),Phys. Rev.
D79 , 074008 (2009), arXiv:0812.0503 [hep-lat].[17] S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, et al. , JHEP , 010 (2012), arXiv:1203.4469 [hep-lat].[18] C. Bernard and D. Toussaint (MILC), Phys. Rev. D ,074502 (2018), arXiv:1707.05430 [hep-lat].[19] R. Dowdall, C. Davies, G. Lepage, and C. Mc-Neile (HPQCD), Phys.Rev. D88 , 074504 (2013),arXiv:1303.1670 [hep-lat].[20] D. Hatton, C. Davies, and G. Lepage, Phys. Rev. D ,094514 (2020), arXiv:2009.07667 [hep-lat].[21] G. P. Lepage et al. , Nucl. Phys. Proc. Suppl. , 12(2002), arXiv:hep-lat/0110175.[22] M. Hayakawa and S. Uno, Prog. Theor. Phys. , 413(2008), arXiv:0804.2044 [hep-ph].[23] M. Steffen, Astron. Astrophys. , 443 (1990).[24] R. Dowdall, C. Davies, R. Horgan, G. Lepage, C. Mona-han, J. Shigemitsu, and M. Wingate, Phys. Rev. D ,094508 (2019), arXiv:1907.01025 [hep-lat]. [25] G. P. Lepage, lsqfit v. 11.7, https://github.com/gplepage/lsqfit , doi:10.5281/zenodo.4037174 (2020).[26] G. P. Lepage, gvar v. 11.9.1, https://github.com/gplepage/gvar , doi:10.5281/zenodo.4290884 (2020).[27] A. Gray, I. Allison, C. Davies, E. Dalgic, G. Lepage, et al. (HPQCD Collaboration), Phys.Rev.
D72 , 094507 (2005),arXiv:hep-lat/0507013 [hep-lat].[28] T. Burch, C. DeTar, M. Di Pierro, A. X. El-Khadra, E. D.Freeland, S. Gottlieb, A. S. Kronfeld, L. Levkova, P. B.Mackenzie, and J. N. Simone, Phys. Rev.
D81 , 034508(2010), arXiv:0912.2701 [hep-lat].[29] S. Meinel, Phys. Rev.
D82 , 114502 (2010),arXiv:1007.3966 [hep-lat].[30] Y. Aoki, N. H. Christ, J. M. Flynn, T. Izubuchi,C. Lehner, M. Li, H. Peng, A. Soni, R. S. Van de Water,and O. Witzel (RBC, UKQCD), Phys. Rev.
D86 , 116003(2012), arXiv:1206.2554 [hep-lat].[31] R. J. Dowdall, C. T. H. Davies, T. Hammant, R. R.Horgan, and C. Hughes (HPQCD), Phys. Rev.
D89 ,031502 (2014), [Erratum: Phys. Rev.D92,039904(2015)],arXiv:1309.5797 [hep-lat].[32] R. Mizuk et al. (Belle), Phys. Rev. Lett. , 232002(2012), arXiv:1205.6351 [hep-ex].[33] G. Bonvicini et al. (CLEO), Phys. Rev. D , 031104(2010), arXiv:0909.5474 [hep-ex].[34] B. Aubert et al. (BaBar), Phys. Rev. Lett. , 161801(2009), arXiv:0903.1124 [hep-ex].[35] B. Aubert et al. (BaBar), Phys. Rev. Lett. , 071801(2008), [Erratum: Phys.Rev.Lett. 102, 029901 (2009)],arXiv:0807.1086 [hep-ex].[36] A. X. El-Khadra, A. S. Kronfeld, and P. B. Macken-zie, Phys.Rev. D55 , 3933 (1997), arXiv:hep-lat/9604004[hep-lat].[37] L. Levkova and C. DeTar, Phys.Rev.
D83 , 074504(2011), arXiv:1012.1837 [hep-lat].[38] S. M. Ryan and D. J. Wilson (Hadron Spectrum), (2020),arXiv:2008.02656 [hep-lat].[39] Y. Aoki et al. , Phys. Rev. D , 054510 (2008),arXiv:0712.1061 [hep-lat].[40] C. Sturm, Y. Aoki, N. Christ, T. Izubuchi, C. Sachra-jda, and A. Soni, Phys. Rev. D , 014501 (2009),arXiv:0901.2599 [hep-ph].[41] A. Pivovarov, Phys. Atom. Nucl. , 1319 (2002),arXiv:hep-ph/0011135.[42] C. Davies, Lect. Notes Phys. , 1 (1998), arXiv:hep-ph/9710394.[43] B. Jones and R. Woloshyn, Phys. Rev. D , 014502(1999), arXiv:hep-lat/9812008.[44] G. Donald, C. Davies, J. Koponen, and G. Lepage, Phys.Rev. Lett. , 212002 (2014), arXiv:1312.5264 [hep-lat].[45] V. Lubicz, A. Melis, and S. Simula (ETM), Phys. Rev.D , 034524 (2017), arXiv:1707.04529 [hep-lat].[46] Y. Chen, W.-F. Chiu, M. Gong, Z. Liu, and Y. Ma(chiQCD), (2020), arXiv:2008.05208 [hep-lat].[47] B. Colquhoun, C. Davies, R. Dowdall, J. Kettle, J. Ko-ponen, G. Lepage, and A. Lytle (HPQCD), Phys. Rev.D , 114509 (2015), arXiv:1503.05762 [hep-lat].[48] S. Bekavac, A. Grozin, P. Marquard, J. Piclum, D. Sei-del, and M. Steinhauser, Nucl. Phys. B , 46 (2010),arXiv:0911.3356 [hep-ph]. [49] T. DeGrand and E. T. Neil, Phys. Rev. D , 034504(2020), arXiv:1910.08561 [hep-ph].[50] W. Kwong, J. L. Rosner, and C. Quigg,Ann.Rev.Nucl.Part.Sci. , 325 (1987).[51] E. J. Eichten and C. Quigg, Phys.Rev. D52 , 1726 (1995),arXiv:hep-ph/9503356 [hep-ph].[52] E. Eichten, S. Godfrey, H. Mahlke, and J. L. Rosner,Rev.Mod.Phys. , 1161 (2008), arXiv:hep-ph/0701208[hep-ph].[53] Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky,and J. G. Wacker, Phys. Rev. Lett. , 021301 (2015),arXiv:1411.3727 [hep-ph]. [54] K. Harigaya and Y. Nomura, Phys. Rev. D , 035013(2016), arXiv:1603.03430 [hep-ph].[55] B. Chakraborty, C. Davies, G. Donald, J. Koponen, andG. Lepage (HPQCD), Phys. Rev. D , 074502 (2017),arXiv:1703.05552 [hep-lat].[56] I. Allison et al. , Phys.Rev. D78 , 054513 (2008),arXiv:0805.2999 [hep-lat].[57] J. H. Kuhn, M. Steinhauser, and C. Sturm, Nucl.Phys.
B778 , 192 (2007), arXiv:hep-ph/0702103 [HEP-PH].[58] T. Aoyama et al. , Phys. Rept.887