Charmonium properties from lattice QCD + QED: hyperfine splitting, J/ψ leptonic width, charm quark mass and a c μ
D. Hatton, C. T. H. Davies, B. Galloway, J. Koponen, G. P. Lepage, A. T. Lytle
CCharmonium properties from lattice QCD + QED: hyperfine splitting,
J/ψ leptonicwidth, charm quark mass and a cµ D. Hatton, ∗ C. T. H. Davies, † B. Galloway, J. Koponen, G. P. Lepage, and A. T. Lytle (HPQCD collaboration) ‡ SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK High Energy Accelerator Research Organisation (KEK), Tsukuba 305-0801, Japan Laboratory for Elementary-Particle Physics, Cornell University, Ithaca, New York 14853, USA INFN, Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Roma RM, Italy
We have performed the first n f = 2 + 1 + 1 lattice QCD computations of the properties (massesand decay constants) of ground-state charmonium mesons. Our calculation uses the HISQ action togenerate quark-line connected two-point correlation functions on MILC gluon field configurationsthat include u/d quark masses going down to the physical point, tuning the c quark mass from M J/ψ and including the effect of the c quark’s electric charge through quenched QED. We obtain M J/ψ − M η c (connected) = 120.3(1.1) MeV and interpret the difference with experiment as theimpact on M η c of its decay to gluons, missing from the lattice calculation. This allows us todetermine ∆ M annihiln η c =+7.3(1.2) MeV, giving its value for the first time. Our result of f J/ψ =0.4104(17) GeV, gives Γ(
J/ψ → e + e − )=5.637(49) keV, in agreement with, but now more accuratethan experiment. At the same time we have improved the determination of the c quark mass,including the impact of quenched QED to give m c (3 GeV) = 0.9841(51) GeV. We have also usedthe time-moments of the vector charmonium current-current correlators to improve the lattice QCDresult for the c quark HVP contribution to the anomalous magnetic moment of the muon. We obtain a cµ = 14 . × − , which is 2.5 σ higher than the value derived using moments extracted fromsome sets of experimental data on R ( e + e − → hadrons). This value for a cµ includes our determinationof the effect of QED on this quantity, δa cµ = 0 . × − . I. INTRODUCTION
The precision of lattice QCD calculations has beensteadily improving for some time and is now approach-ing, or has surpassed, the 1% level for multiple quan-tities. Good examples are the masses and decay con-stants of ground-state pseudoscalar mesons [1]. The me-son masses can be used to tune, and therefore determine,quark masses. The decay constants can be combinedwith experimental annihilation rates to leptons to deter-mine elements of the Cabibbo-Kobayashi-Maskawa ma-trix. The accuracy of modern lattice QCD results meansthat sources of small systematic uncertainty that couldappear at the percent level need to be understood. Atthis level QED effects, i.e. the fact that quarks carry elec-tric as well as color charge, come into play. A naive ar-gument that such effects could be O ( α QED ) would implya possible 1% contribution. One key driver for the lat-tice QCD effort to include QED effects has been that ofcalculations of the hadronic vacuum polarisation (HVP)contribution to the anomalous magnetic moment of themuon, a µ . New results are expected soon from the Muon g − ∗ [email protected] † [email protected] ‡ the precision of a few percent for a µ and the systematicuncertainties, for example from neglecting QED effects,have become a major focus of attention (see, for exam-ple [4–6]).QED effects can have large finite-volume correctionswithin a lattice QCD calculation because the Coulombinteraction is long-range. For electrically neutral correla-tion functions, such as the ones needed for a calculationof the HVP, and that we study here, this is much lessof an issue (and we demonstrate this). Since α ns α QED isnot very different in size from α QED at hadronic scales,calculations must be fully nonperturbative in QCD. Aconsistent calculation must also allow for the retuning ofquark masses needed when QED effects are included.Here we examine the properties of ground-state char-monium mesons more accurately than has been possiblein previous lattice QCD calculations . Since it is possibleto obtain statistically very precise results for the charmo-nium system it is a good place to study small systematiceffects from QED and other sources. We are able to seesuch effects in our results and quantify them. We include u , d , s and c quarks in the sea for the first time and haveresults on gluon field configurations with physical u/d sea quarks. We also analyse the impact of including anelectric charge on the valence c quarks (only). This ap-proximation, known as ‘quenched QED’, should capture For a different kind of lattice QCD calculation that maps out thespectrum more completely, but paying less attention to ground-states, see [7]. a r X i v : . [ h e p - l a t ] M a y the largest QED effects and enable us to see how muchof a difference QED makes. For the vector ( J/ψ ) me-son properties we improve on earlier results by using anexact method to renormalise the lattice vector current(both in QCD and QCD+QED). To tune the c quarkmass we use the J/ψ meson mass, taking into accountthe retuning that is required when QED is switched on.We find the impact of this to be of similar size to themore direct QED effects.The quantities that we focus on here are the massesand decay constants of the ground-state η c and J/ψ mesons, the c quark mass and the contribution of the c vacuum polarisation to a µ .The correlation functions that we calculate in our lat-tice QCD and QCD+QED calculations are ‘connected’correlation functions i.e. they are constructed from com-bining charm quark propagators from the source to thesink. We do not include diagrams in which the c and c quarks annihilate to multiple gluons and hence hadrons.We can gain some insight into the effect this annihilationchannel has on the meson masses by looking at the mesonwidths, which are twice the imaginary part of the massand are dominated by hadronic channels. The J/ψ hasa tiny width of 93 keV [1] but the pseudoscalar η c canannihilate to two gluons, allowing it to mix with otherflavour-singlet pseudoscalars, and it has a width of 32MeV [1]. The annihilation channel might then be ex-pected to have a larger impact on the η c mass and leadlattice QCD calculations of the mass from connected cor-relators to disagree with experiment. The only way toachieve the O (1 MeV) accuracy required to see this is todetermine the mass difference between the J/ψ and η c (the hyperfine splitting). Any shift in the η c mass willhave a much larger (by a factor of 30) relative effect inthis splitting. Previous lattice QCD calculations of thehyperfine splitting from connected correlation functionshave not been accurate enough to see a significant differ-ence with experiment. Here, for the first time, we can seesuch a difference because we have very good control bothof discretisation effects and sea-quark mass effects andcan also determine the impact of QED on this quantity.A further place in which QED effects need to be quan-tified, given our accuracy, is that of the determinationof the c quark mass. We do this by tuning our resultsto the experimental J/ψ meson mass with and withoutelectric charge on the valence c quarks and also determinethe small change in the mass renormalisation factor, Z m ,needed to convert to the standard MS mass.Our analysis of the vector charmonium correlationfunctions with a completely nonperturbative renormal-isation of the vector current allows much improved ac-curacy in the determination of the leptonic decay widthof the J/ψ . Using the same correlators, we determinethe charm quark portion of the hadronic vacuum polar-isation contribution to the anomalous magnetic of themuon, a HVP,c µ , along with the impact of quenched QEDon this quantity. We can compare this to phenomenolog-ical results from R ( e + e − → hadrons). The paper is laid out as follows: • Section II describes the lattice QCD calculationand the inclusion of quenched QED; • Section III describes our determination of the hy-perfine splitting; • Section IV, the c quark mass; • Section V, the
J/ψ and η c decay constants; • Section VI the time-moments of c vector-vector cor-relators and the c quark hadronic vacuum polarisa-tion contribution to a µ ; • Section VII gives our conclusions.Each section of results includes a description of the pureQCD calculation followed by a determination of the im-pact of quenched QED on the result and then a discus-sion subsection including comparison to both experimentand previous lattice QCD calculations, where applicable.Finally Section VII collects up all of our results and sum-marises our conclusions.
II. LATTICE SETUP
We perform calculations on a total of 15 gluon fieldensembles that include the effects of light (up/down withthe same mass), strange and charm quarks in the sea( n f = 2 + 1 + 1) all using the HISQ action [8] and gen-erated by the MILC collaboration [9, 10]. These includelattices at 6 different β (the bare QCD coupling) valuescorresponding to 6 different sets of lattice spacing val-ues with the finest lattice reaching a spacing of ∼ . u/d masses at the physical point on all but the finesttwo lattice spacings. We also employ three ensembleswith shared parameters except for their spatial extent inunits of the lattice spacing L s . These ensembles allow usto investigate finite volume effects in our QED analysis.The gluon action on these ensembles is improved so thatdiscretisation errors through O ( α s a ) are removed [12].Parameters for the ensembles are given in Table I.On these gluon field configurations we calculate propa-gators for valence c quarks by solving the Dirac equationfor a source consisting of a set of Gaussian random num-bers across a timeslice (a random wall source). We usemultiple time sources per configuration to improve sta-tistical accuracy. The number of configurations used andthe number of time sources is given in Table I. The ta-ble also gives the valence c quark masses in lattice units,which may differ from those in the sea because they aretuned more accurately. This will be discussed further TABLE I. Details of the lattice gluon field ensembles and calculation parameters used. The lattice spacing is determined fromthe Wilson flow parameter, w [13], with w /a values given in column 2. The lattice spacing can be determined in fm by using w = 0 . f π ). L s and L t are the lattice spatial and temporal extents in lattice units. am val c is thevalence c quark mass with (cid:15) Naik the corresponding Naik parameter (see text). The column headed N cfg × N t shows both thenumber of configurations used in the pure QCD calculation and the number of time sources for propagators per configurations. N cfg , QED refers to the number of configurations (and time sources) used in QCD+QED calculations. Sets 1-3 will be referredto as very coarse ( a ≈ a ≈ a ≈ a ≈ a ≈ .
045 fm) and 15 as exafine ( a ≈ β w /a L s L t am sea l am sea s am sea c am val c (cid:15) Naik N cfg × N t N cfg , QED × N t × × ×
163 5.80 1.1367(5) 36 48 0.00235 0.0647 0.831 0.863 -0.3670 1000 × × ×
166 6.00 1.4029(9) 32 64 0.00507 0.0507 0.628 0.650 -0.2378 1000 × ×
167 6.00 1.4029(9) 40 64 0.00507 0.0507 0.628 0.650 -0.2378 - 220 ×
168 6.00 1.4149(6) 48 64 0.00184 0.0507 0.628 0.643 -0.2336 1000 × × × × × × × × × × below. The HISQ action [8] includes an improved dis-cretisation of the covariant derivative in the Dirac equa-tion. This removes tree-level a discretisation errors bythe addition of a 3-link ‘Naik’ term to the symmetric1-link difference. For heavy quarks the coefficient of theNaik term is adjusted from 1 to 1+ (cid:15) Naik to remove ( am ) errors at tree-level [8]. A closed-form expression for (cid:15) interms of the tree-level quark pole mass is given in [15]along with the formula for the tree-level quark pole massin terms of the bare mass. Table I gives the values of (cid:15) that we use.We combine charm quark and antiquark propagatorsto calculate two types of quark-line connected two-pointcorrelation functions: pseudoscalar and vector. Theground state of the pseudoscalar correlation function cor-responds to the η c meson and the vector correlation func-tion, to the J/ψ . When using staggered quarks, as here,the different spin structures are implemented using po-sition dependent phases in the operators at source andsink. The two-point ‘Goldstone’ pseudoscalar ( γ ⊗ γ inspin-taste notation) correlation functions are simply con-structed from quark propagators S ( x,
0) from the originto x as C ( t ) = 14 (cid:88) x (cid:104) Tr( S ( x, S † ( x, (cid:105) , (1)where the factor of 4 accounts for the taste multiplicitywith staggered quarks. For the vector correlation func-tions we use a local vector operator (spin-taste γ i ⊗ γ i ).The correlation functions then combine S ( x,
0) with apropagator made from patterning the source with a phase( − x i and inserting ( − x i at the sink timeslice as we tie the propagators together and sum over spatial sites.Our vector correlation functions average over all spatialpolarisations, i , for improved statistical precision. Notethat we do not calculate any quark-line disconnected cor-relation functions.The HISQ local vector current is not conserved and re-quires renormalisation. For this purpose we use the RI-SMOM momentum subtraction scheme implemented onthe lattice as discussed in [16]. In [16] it was shown that,because of the Ward-Takahashi identity, these renormal-isation factors do not suffer any contamination by non-perturbative artefacts (condensates) and can therefore besafely used in calculations such as those presented here.The quenched QED correction to the RI-SMOM vectorcurrent renormalisation was also given in [16] and shownto be tiny ( ∼ . Z V values only differ from 1 at the1% level and quenched QED provides a small correctionto this difference from 1). Here we use the Z V valuesfrom [16] at a scale µ of 2 GeV. We will also demonstrate(see Section V) that using µ = 3 GeV gives the same re-sult as it must for a Z V that correctly matches the latticeto continuum physics.Since we make use of an ensemble (set 14) with a finerlattice spacing than those studied in [16] we have directlycalculated the value of Z V on set 14 at µ = 2 GeV in ad-dition. We have, however, only used a small number ofconfigurations (6) in that calculation due to the compu-tational limitation of the stringent Landau gauge fixingrequired. We therefore double the statistical uncertaintyfor Z V on that ensemble. This has very little impact onour final results as the Z V uncertainty is small. See Ap-pendix A for a discussion of our Z V values, where we alsoderive a Z V value for set 15.In order to tune the mass of the valence charm quarkwe use bare charm mass values on each lattice that pro-duce a J/ψ mass equal to the experimental value (both inpure QCD and in QCD+QED). We choose the
J/ψ hererather than the η c because the relatively large width ofthe η c means that annihilation effects that we are notincluding could lead to small (order 0.1% ) uncertaintiesin the mass. This is mentioned in Section I and will bediscussed further in Section III. We measure our valence c mass mistunings as the difference between our lattice J/ψ mass and the experimental average value of 3.0969GeV (with negligible uncertainty) [1]. The two panelsof Fig. 1, where the horizontal line is the experimentalvalue, show that our mistunings are well below 0.5%.These mistunings are allowed for in our final fits.
A. Two-point correlator fits
We fit the two-point correlation functions describedabove as a function of the time separation, t , betweensource and sink. The aim is to extract the energies(masses) and amplitudes (giving decay constants) of theground-state mesons in each channel. However it is im-portant to allow for the systematic effect of excited statesthat are present in the correlation functions and can af-fect the ground-state values if they are not taken intoaccount. We do this by fitting the correlators to sumsof exponentials associated with each energy eigenvalueand using Bayesian priors to constrain the (ordered) ex-cited states in the standard way [17]. The pseudoscalarcorrelators are fit to C P ( t ) = (cid:88) i A Pi f ( E Pi , t ) , (2) f ( E, t ) = e − Et + e − E ( T − L t ) . The vector correlators require a more complicated formbecause of the presence of opposite parity states as aresult of the use of staggered quarks: 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) . (3)We cut out the correlator values at low values of t ,below some value t min (5 −
10) where excited state con-tamination is most pronounced. We also use a standardprocedure (see Appendix D of [18]) to avoid underesti-mating the low eigenvalues of the correlation matrix andhence the uncertainty.The fit parameters that we need from Eq. 2 and 3 arethe mass of the ground-state ( E P and E V ) and the am-plitude ( A P and A V ). From the amplitude we determinethe decay constant, see Section V. . . . . . am c . . . . M J / ψ [ G e V ] . . . . . am c . . . . M J / ψ [ G e V ] QCD+QEDQCD
FIG. 1. Top panel: Pure QCD
J/ψ masses on all 15 sets fromTable I using the valence masses in that Table. Two errorbars are shown. The error can be broken into parts that areuncorrelated between different sets and the contribution fromfixing the lattice spacing from the physical value of w whichis correlated. The outer error bar shows the full uncertaintyand the inner bar the uncertainty without the contributionfrom w . Bottom panel: The J/ψ masses on sets 2, 6, 10and 12 with and without the inclusion of quenched QED.On each set the same valence mass (and lattice spacing) wasused for both pure QCD and QCD+QED but the points havebeen separated on the x -axis for clarity. The error bars arethe same as for the top panel, but note that here there is acorrelation of the uncertainty from w /a for the QCD andQCD+QED results on each set. The QCD+QED results areabove the pure QCD results in every case. B. Including quenched QED
We perform calculations in both lattice QCD and inlattice QCD with quenched QED. By quenched QED wemean that we include effects from the valence quarks hav-ing electric charge but we neglect effects from the elec-tric charge of the sea quarks. To include quenched QEDeffects we generate a random momentum space photonfield A µ ( k ) in Feynman gauge for each QCD gluon fieldconfiguration. This choice of gauge simplifies the genera-tion of the photon field as the QED path integral weighttakes the form of a Gaussian with variance 1 / ˆ k whereˆ k = 2sin( k µ / L formulation [19]. A µ is then Fourier transformedinto position space. We have checked that these Feyn-man gauge A µ fields produce the plaquette and averagelink expected from O ( α QED ) perturbation theory (read-ily obtained from O( α s ) calculations in lattice QCD withWilson glue [20, 21]). These gauge fields are exponenti-ated as exp( ieQA µ ) to give a U(1) field which is thenmultiplied into the QCD gauge links before HISQ smear-ing. Q is the quark electric charge in units of the protoncharge e .Quenched QED affects only the valence quarks; the seaquarks remain uncharged. Hence there is no coupling ofQED effects to purely gluonic quantities. This meansthat the Wilson flow parameter w /a , measured on eachensemble, is unchanged. The physical value of w that isused to determine the lattice spacing on each ensemblewas determined in [14] by matching the decay constantof the π meson, f π , in lattice QCD to that obtained fromexperiment. The experimental value of f π is obtainedfrom measurement of the rate for π → (cid:96)ν [ γ ] decay where[ γ ] indicates that the rate is fully inclusive of additionalphotons. The rate obtained is then adjusted to removeelectromagnetic and electroweak corrections and to give a‘purely leptonic rate’ corresponding to weak annihilationat lowest order in the absence of QED [22]. Combiningthis with a determination of | V ud | from nuclear β decay [1]gives an experimental value of f π ≡ f expt π which is a ‘pureQCD’ value, albeit that for a physical π + meson. Thedominant uncertainty in f expt π is that from the remaininguncertainty in the electromagnetic corrections to the ex-perimental rate, mainly from the hadronic-structure de-pendent contributions to the emission of additional pho-tons. This is set at 0.1% in [22].Because f expt π is a pure QCD quantity it can be usedto set the lattice spacing in lattice QCD in a way thatshould be minimally different for lattice QCD+QED .Small differences might still be expected between the lat-tice QCD f π and the experimental value from the waythat the quark masses are tuned in a pure QCD scenario.The lattice QCD calculation in [14] used m u = m d andtuned the average mass, m l , to the experimental mass ofthe π , which is the mass that both neutral and charged π mesons have in the absence of QED, up to quadraticcorrections in the u − d mass difference. An uncertaintywas included in the π mass to allow for these corrections,taking an estimate from chiral perturbation theory [24].We expect the impact of such effects to be tiny, well be-low 0.1%. They are at the same level as potential effectsfrom the sea and would therefore be only possible to pindown with a calculation that included the impact of hav-ing electrically charged quarks in the sea. Indeed, f π cannot readily be calculated in lattice QCD+QEDbecause of infrared QED effects from an electrically charged π + .Calculations have been done that confirm the size of radiativecorrections to f π , however [23]. These expectations are backed up by recent latticeQCD+QED results [6] that used the Ω baryon mass todetermine w . The impact of QED for the sea quarkswas included to first order in α QED . No effects linearin m u − m d are expected in M Ω because, like f π and M π above, it is symmetric under u ↔ d interchange.Strong-isospin breaking effects were therefore ignored.The impact of QED in the sea on w M Ω was foundto be O (0 . f π analysis) was O (0 . w using M Ω from [6] agreeswell with the result using f π from [14], although the un-certainties in both cases are completely dominated bythose from the pure QCD, isospin-symmetric part of thecalculation.From this we conclude that the largest effect fromadding QED to lattice QCD will come from the inter-action between the electric charges of the valence quarksand to study this effect we can compare lattice QCD plusquenched QED with pure lattice QCD using the samevalue of the lattice spacing, determined from f π , in bothcalculations.An estimate of the size of corrections from quenchedQED (often simply referred to as QED in what follows)in charmonium systems can be obtained by studying theeffect on the J/ψ mass. The bottom panel of Fig. 1shows the
J/ψ mass for the same valence c mass for bothQCD+QED and pure QCD calculations on sets 2, 6, 10and 12. The QCD+QED and QCD results at the samelattice spacing are separated on the x -axis for clarity.All points share a correlated uncertainty (the outer errorbar) from w and this dominates the uncertainty. Theuncorrelated error is shown by the smaller inner errorbar. Note that the points at the same lattice spacing arealso correlated through their w /a value. The shift ofthe mass in QCD+QED compared to pure QCD is verysmall, at the level of 0.1%, and is upwards.When discussing QCD+QED and pure QCD calcu-lations of some quantity X we will use the notation X [QCD + QED] and X [QCD] respectively. We will of-ten consider the ratio of the two for which we will usethe shortened notation R (0)QED [ X ]. R will refer to the‘bare’ ratio defined using the same bare quark mass am c in both QCD+QED and pure QCD calculations. R QED will refer to the final QED-renormalised ratio which in-cludes the impact of retuning the c quark mass to givethe experimental J/ψ mass in both the QCD+QED andpure QCD cases. So R [ X ] ≡ X [QCD + QED] X [QCD] (cid:12)(cid:12)(cid:12)(cid:12) fixed am c (4) R QED [ X ] ≡ X [QCD + QED] X [QCD] (cid:12)(cid:12)(cid:12)(cid:12) fixed M J/ψ . As shown in Fig. 1 the bare c quark mass has to re-adjusted downwards for QCD+QED relative to pureQCD. C. Fitting strategy
We have results in pure QCD for all of the sets in Ta-ble I and QCD+QED results on a subset of ensembles.To be able to simultaneously account both for the ‘direct’effects of QED and for the effects of valence c mass mis-tuning, which may be similarly sized, we choose to fit allof this data in a single fit for each quantity we consider.The generic form of the fit we use for a quantity X is X ( a , Q ) = x (cid:20) (cid:88) i =1 c ( i ) a ( am c ) i + (5) c m, sea δ sea ,udsm { c a , sea (Λ a ) + c a , sea (Λ a ) } + c c, sea δ sea ,cm + c c, val δ val ,cm + α QED Q { c QED + (cid:88) p =1 c ( p ) aQ ( am c ) p + c val ,Q δ val ,cm } (cid:21) . Here Q is the valence quark electric charge (in unitsof e ) used in the calculation and is therefore 0 in pureQCD. The pure QCD value of this fit at the physicalpoint (the continuum limit with quark masses set to theirphysical values) is x . The value including quenched QEDcorrections is x [1 + α QED Q c QED ]. Note that the factorof α QED multiplying the QED part of the fit function isthere so that the fit parameters are order 1; this is not a perturbative expansion in α QED . All orders of α QED are included along with α QED α s terms. The only piecesmissing from our calculation and fit are α QED α s termsfrom interaction with electrically charged sea quarks.We now describe each of the terms in Eq. (5) in turn.The ( am c ) i terms on the first line account for discretisa-tion effects. Because we are dealing with heavy quarkshere, the scale of discretisation effects can be set by m c and will typically be larger than those for light-quarkquantities. Since m c > Λ QCD any discretisation effectsset by scale Λ
QCD will simply appear as m c -scale discreti-sation effects with a small coefficient.The terms on the second line allow for mistuning ofthe sea u/d and s masses and discretisation effects inthat mistuning (we shall see that those are important forthe hyperfine splitting). The total of the mistuning ofthe sea masses is defined as in [25]: δ sea ,udsm = 2 m sea l + m sea s − m phys l − m phys s m phys s . (6) m phys s is taken from [25] or, where not available on thefinest lattices, calculated from the tuned c quark massand the m c /m s ratio given in [25]. The value of Λ in thediscretisation effects multiplying the sea-quark mistuningis taken as 1 GeV ( ∼ m c ).The effect of mistuning the charm mass in the sea isincluded in the third line of Eq. (5) using δ sea ,cm = m sea c − m phys c m phys c . (7) The values of m phys c are taken from [25]. Although thisused a slightly different tuning method the differences arenegligible for this purpose. We have tested that includingdiscretisation effects for this term has no effect on the fit.Mistuning in the valence mass is accounted for through δ val ,cm on the third line. We define this as δ val ,cm = M J/ψ − M expt J/ψ M expt J/ψ . (8)The experimental value of the J/ψ mass is 3.0969 GeV [1]with negligible uncertainty.In order to make use of the correlations between ourQCD+QED and pure QCD results on the same gluonfield configurations we perform simultaneous fits to thecorrelators in each case. The fits then capture the corre-lations and we can propagate them to the fit of Eq. (5).At the same time it allows us to determine the ratio ofQCD+QED to pure QCD for the quantities that we willstudy. We will give results for these ratios in the sectionsthat follow.The fit form of Eq. (5) has been constructed such thatthe coefficients (apart from x ) are expected to be of order1. We therefore use priors of 0 ± x for which we take a prior width on its expectedvalue of 20% (the prior mean for x depends on the quan-tity being fitted). III. HYPERFINE SPLITTINGA. Pure QCD
The hyperfine splitting, ∆ M hyp , is calculated on eachensemble as the difference of the vector and pseudoscalarground state masses, in lattice units, divided by the lat-tice spacing. The results for aM η c and aM J/ψ and theirdifference are given in Table II for the pure QCD case.Although the pseudoscalar and vector correlators on eachconfiguration are correlated the fit outputs for the vectorcorrelator dominate the uncertainties and so the correla-tions have very little effect as a result.The pure QCD results are plotted in Fig. 2 along withthe fit form of Eq. (5) for the pure QCD case (i.e. Q = 0).Note the small range of the y -axis - this is possible for ourresults because we have a highly-improved quark actionwith small discretisation errors. Since all tree-level a errors have been removed, the shape of the curve reflectsthe fact that higher-order a and a errors are visible;discretisation errors of this kind are present in all for-malisms of course, but often they are hidden below muchlarger a effects and consequently overlooked. Note alsothe clear dependence on the light sea quark mass seenon the finest lattices. To pin down the value of the va-lence mass mistuning parameter, c val , we include resultsat deliberately mistuned c quark masses (see Table II).These are not shown in the Figure but are included inthe fit. The result for the hyperfine splitting in the pure TABLE II. The η c and J/ψ masses and their difference( a ∆ M hyp ) in pure QCD on each set in lattice units. Thevalues of am val c are those given in Table I except for two caseswith a deliberately mistuned c quark mass: set 6 denotedby a * where am c = 0 .
643 and set 14 denoted by a † where am c =0.188. The pseudoscalar and vector correlator fits havebeen performed separately and the correlations between aM η c and aM J/ψ have therefore been ignored because they have lit-tle impact.Set aM η c aM J/ψ a ∆ M hyp ∗ † . . . . . am c ) M J / ψ − M η c [ M e V ] FIG. 2. The charmonium hyperfine splitting as a functionof lattice spacing on pure QCD ensembles. Our results arefrom ensembles including u , d , s and c quarks in the sea withvarying u/d average quark mass. The raw lattice results, forwell-tuned c quark masses, are given by symbols with errorbars. The error bars include both statistical uncertainties andthose from determining the lattice spacing from w , whichare correlated between values. The results on each group ofensembles with approximately the same lattice spacing aregiven the same symbol. Within these groups, the results gofrom right to left as the u/d quark mass changes from m s / y -axis; thisis the results of discretisation effects being so small for theHISQ action. Results at mistuned c masses are not plottedbut are included in the fit. The fit line is the output of thefit from Eq. (5) at physical quark masses and with Q = 0.The orange cross gives our result in the continuum limit forphysical quark masses. The black cross gives the experimentalaverage result [1]. TABLE III. QCD+QED η c and J/ψ masses and hyperfinesplitting presented as the ratio of the QCD+QED result tothe pure QCD one on that set. Correlations between thecalculations in the QCD+QED and pure QCD cases are usedin the determination of the ratio and result in the very highstatistical accuracy obtained. Note that the ratio is calculatedfor the same am c value in the two cases i.e. the ratio givenhere does not include the impact of retuning the c quark mass.Set R [ M η c ] R [ M J/ψ ] R [∆ M hyp ]2 1.000450(26) 1.000750(27) 1.0086(10)6 1.0008335(59) 1.0010742(81) 1.00774(28)10 1.0011861(54) 1.0014044(76) 1.00739(26)12 1.0015755(48) 1.001787(11) 1.00750(33) . . . am c ) . . . . R Q E D [ ∆ M h y p ] FIG. 3. The fractional effect of quenched QED on the char-monium hyperfine splitting plotted against ( am c ) . The frac-tional effect is determined at the same am c value, i.e. it doesnot include c quark mass retuning effects. The dashed line ishorizontal and shows the weighted average value. The resultsshow the precision that can be obtained by capitalising onthe correlations between QCD+QED and QCD. This enablesa clear demonstration that the impact of quenched QED hereis not dependent on the lattice spacing. QCD case in the continuum limit and for physical quarkmasses is 118.6(1.1) MeV, which is higher than the ex-perimental average value, as is clear in Fig. 2. In orderto understand what this means, we need to quantify allpossible sources of small systematic effects in our calcu-lation, including those from QED.
B. Impact of Quenched QED
The fractional direct effect of quenched QED on the η c and J/ψ masses and the hyperfine splitting are givenin Table III. The correlation between the QCD+QEDand the pure QCD results enables very high statisticalaccuracy to be obtained in the ratio. The inclusion ofquenched QED shifts both the η c and J/ψ masses up by O (0 . am c value. Although these mass shifts are small, there is adifference between the shift for the J/ψ and that for η c .
025 0 .
030 0 .
035 0 . /L s . . . R Q E D [ M p ] p = J/ψp = η c FIG. 4. The fractional shift in the
J/ψ and η c masses fromthe inclusion of quenched QED plotted as a function of 1 /L s on sets 5, 6 and 7. The dashed lines are horizontal lines atthe weighted average values. and so the inclusion of quenched QED also changes thehyperfine splitting. The impact here is more substan-tial, 0.7%, because the hyperfine splitting is so muchsmaller. The size of the direct QED effect on the hyper-fine splitting can be simply estimated by replacing C F α s by Q α QED in a potential model estimate of the splitting.This gives a fractional effect of α QED / (3 α s ), consistentwith what we find.The values of R [∆ M hyp ] are plotted against ( am c ) in Fig. 3. This shows that the results are consistent acrossall lattice spacings and thus discretisation effects in thisratio are smaller than for the masses themselves.Finite-volume effects are an issue in general for QEDcorrections to meson masses but we expect them to besmall for the electrically neutral and spatially small char-monium mesons that we study here. In [26] it is shownthat the finite volume expansion for electrically neutralmesons starts at O (1 /L s ). In Fig. 4 we compare resultsfor the fractional effect of QED on the J/ψ and η c as afunction of 1 /L s . This calculation is done on sets 5, 6 and7 (see Table I) which differ only in their spatial extent.We see no finite-volume effects to well within 0.01%, andwe therefore ignore such effects.Our results including both pure QCD and QCD+QEDare shown in Fig. 5, plotted against ( am c ) . The fit curvefrom Eq. (5) at physical quark masses is also shown.The fit has a χ / dof of 0.59 and gives a hyperfine split-ting in the continuum limit at physical quark masses of119.6(1.1) MeV. Taking the (correlated) ratio betweenthe physical value of the full QCD+QED fit and the phys-ical value from the fit at Q = 0 (i.e. the pure QCD result)we obtain R QED [∆ M hyp ] = 1.00804(43). This ratio nowdoes include the effect of retuning the c quark mass to ob-tain the experimental J/ψ mass when quenched QED isincluded. This retuning requires a reduction of the bare c quark mass by O (0 . O (0 . . . . . . am c ) M J / ψ − M η c [ M e V ] FIG. 5. The charmonium hyperfine splitting as a functionof lattice spacing, including both QCD+QED and pure QCDpoints. The red open triangles are the same lattice resultsas in Fig. 2. The additional QCD+QED points are given asblue open hexagons. The green fit band is the output of thefit from Eq. (5), but now with the c quark electric charge, Q = 2 /
3. The orange cross and orange band gives our resultin the continuum limit for physical quark masses. The blackcross and black band gives the experimental average result [1].TABLE IV. Error budget for our final result for the charmo-nium hyperfine splitting including quenched QED corrections.The uncertainties shown are given as a percentage of the finalresult. The largest uncertainties are clearly from the deter-mination of the lattice spacing. ∆ M hyp a → w /a w There is an additional pure QED contribution to the
J/ψ mass that has not been included yet since it is quark-line disconnected. This comes from a diagram in whichthe cc annihilate into a photon which converts back into cc . The contribution of this diagram is8 πα QED Q M J/ψ | ψ (0) | , (9)where ψ (0) is the nonrelativistic J/ψ wavefunctionequal (in the normalisation being used here) to f J/ψ (cid:112) M J/ψ / √ M J/ψ − M η c (connected) = 120 . .
1) MeV . (10)The error budget for our hyperfine splitting result isgiven in Table IV. We follow Appendix A of [28] forthe meaning of the uncertainties contributing to the er-ror budget. The majority of the uncertainty is associ-ated with the lattice spacing determination, either fromthe correlated w uncertainty or the individual w /a un-certainties. This is not surprising because the hyper-fine splitting is sensitive to uncertainties in the deter-mination of the lattice spacing for the reasons discussedin [29]. We have separated out the uncertainty arisingfrom the pure QCD data and the R [∆ M hyp ] valuesfrom Table III which we label ‘Pure QCD Statistics’ and‘QCD+QED Statistics’ in Table IV. The sea mistuninguncertainty comes from the c m coefficients in Eq. 5 andthe valence mistuning uncertainty from the c val and c val ,Q coefficients. The a → c a and c aQ coefficients.Our final result is for the charmonium hyperfine split-ting determined from quark-line connected correlationfunctions in QCD and including the impact of QED ef-fects, through explicit calculation in quenched QED. Weexpect the effect of further QED effects in the sea tobe negligible compared to our 1% uncertainty. The onlysignificant Standard Model effect then missing is that ofquark-line disconnected diagrams in which the cc annihi-late to gluons. We expect this effect to be much larger forthe η c than for the J/ψ so a comparison of our result forthe hyperfine splitting to experiment can yield informa-tion on the size and sign of the annihilation contributionto the η c mass. This is discussed in the next subsection. C. Discussion: Hyperfine Splitting
The experimental average value of the hyperfine split-ting (113.0(5) MeV) from the Particle Data Group(PDG) [1] is calculated as the difference of the sepa-rate averages for the
J/ψ and η c masses. The differentexperimental results contributing to the PDG averageof the two masses are shown in Fig. 6. For the J/ψ mass the average is dominated by the most recent re-sult from KEDR [41]. There are only three experimen-tal results used in these analyses that can independentlyproduce values for the hyperfine splitting. These are theKEDR experiment [34, 41] and two LHCb analyses indifferent channels [31, 32]. The LHCb result in [32] usedthe η c (2S) → pp decay while the analysis of [31] used η c (1S) → pp . In the comparison plot of Fig. 7 [31] isreferred to as LHCb15 and [32] as LHCb17.Fig. 7 shows a comparison of lattice QCD results forthe charmonium hyperfine splitting along with the PDGaverage value and separate experimental values that mea-sured this splitting. Previous calculations on gluon fieldconfigurations that included n f = 2 + 1 flavours ofsea quarks by HPQCD [29] and by Fermilab/MILC [38]both obtained values above the experimental average, al-though only by just over one standard deviation.The result we present here is substantially more precisethan these earlier studies and for the first time displaysa significant, 6 σ , difference from the experimental aver- . . . . . M η c [GeV] BES3 BES3BABAR BABARKEDRLHCb15LHCb17 LHCb17BELL18 . . . . . . M J/ψ [GeV]
SPECE760OLYAKEDR
FIG. 6. Comparison of different experimental results forthe
J/ψ and η c masses along with the PDG average values.The η c results represent a recent subset of those used in thePDG average. The most recent result is from BELLE (de-noted BELL18) [30]. There are three different determinationsfrom LHCb [31–33], two of which also measured the hyperfinesplitting. We include a KEDR measurement [34], two fromdifferent BaBar analyses [35] and two from BESIII [36, 37]. age, clearly showing that the lattice result lies above theexperimental one. We interpret this as the effect of ig-noring annihilation to gluons in the calculation of the η c mass. From the comparison of our results to experimentwe conclude that these annihilation effects increase the η c mass by 7.3(1.2) MeV, where the uncertainty is dom-inated by that from the lattice calculation.Previous analyses of this issue have given mixed re-sults. In NRQCD perturbation theory [8] we can relatethe shift in the η c mass to its total (hadronic) widththrough the perturbative amplitude for cc → gg → cc atthreshold [42]. Then∆ M η c = Γ η c (cid:18) − π + O ( α s , v /c ) (cid:19) (11)= 31 . × (cid:0) − .
195 + O ( α s , v /c ) (cid:1) . The leading term here gives − η c and other flavour-singlet pseudoscalar mesons. Since0
105 110 115 120 125 M J/ψ − M η c [MeV] KEDRLHCb15LHCb17HPQCD12Brice˜no et al 12 χ QCD14Fermilab/MILC 19This work: pure QCDThis work: QCD+QED
FIG. 7. Comparison of different lattice results for the char-monium hyperfine splitting and separate experimental resultsthat measure this difference, as well as the PDG average (pinkband). The PDG average is obtained from taking the differ-ences of the PDG
J/ψ and η c masses (see Fig. 6) rather thanonly from experiments that directly measure the splitting.The squares give lattice QCD results from calculations thatinclude n f = 2+1 flavours of quarks in the sea. The hexagonsgive results that include n f = 2 + 1 + 1 flavours of sea quarks,including the results we present here at the top of the plot.All lattice QCD results have had uncertainties from neglect-ing η c annihilation removed so that we might expect somedifference between them and experiment (see text), but pre-vious lattice QCD results have not been accurate enough tosee this. The Fermilab/MILC result used Fermilab c quarkson gluon field configurations with asqtad sea quarks [38] andthe previous HPQCD result [29] used HISQ quarks on thesame asqtad-sea ensembles. Brice˜no et al [39] used a modifi-cation of the Fermilab approach known as relativistic heavyquarks on the n f = 2 + 1 + 1 HISQ sea quark ensembles thatwe use here. The χ QCD result [40] used overlap quarks ongluon field configurations including n f = 2 + 1 domain-wallsea quarks. these are lighter than the η c this mixing could give apositive correction to the η c mass. Direct lattice QCDdetermination of the effect, by calculating the appropri-ate quark-line disconnected correlation functions, has sofar not proved possible. This is because the lighter statesthat are introduced by the mixing make it very hard topin down a small effect on the mass of a particle, the η c ,which is so much further up the spectrum in this chan-nel. An estimate of the mass shift of +1–4 MeV wasobtained in the quenched approximation in which thismixing is not possible but where mixing with glueballscould happen instead [43].Our result for the hyperfine splitting, by its accuracy,provides for the first time a clear indication of the size ofthe impact of the η c annihilation to gluons on its mass:∆ M annihln η c = +7 . .
2) MeV . (12) IV. DETERMINATION OF m c A. Pure QCD
In [44] we showed that it is possible to determinethe strange and charm quark masses accurately in lat-tice QCD using an intermediate momentum-subtractionscheme. By intermediate we mean that the mass renor-malisation factor to convert the tuned bare lattice quarkmass to the momentum-subtraction scheme is calculatedon the lattice. The conversion from the momentum-subtraction to the final preferred MS scheme is carriedout using QCD perturbation theory in the continuum. Todo this accurately it is important to use a momentum-subtraction scheme that has only one momentum scale, µ . This means that the squared 4-momentum on eachleg of the vertex diagram, from which the mass renor-malisation factor is calculated, is µ . The RI-SMOMscheme [45] used in [44] is such a scheme. A further im-portant point is that the mass renormalisation factor willbe contaminated by nonperturbative (condensate) arte-facts through its nonperturbative calculation on the lat-tice. To identify and remove these artefacts (that appearas inverse powers of µ ) requires calculations at multiplevalues of µ and a fit to the results, as discussed in [44].Below we briefly summarise the procedure followedin [44]:1. Determine the tuned bare quark mass and latticespacing at physical sea quark masses for each setof gluon field configurations at a fixed β value. Wedo this following Appendix A of [25].2. Calculate the mass renormalisation factor, Z SMOM m ,that converts the lattice quark masses to the RI-SMOM scheme for each β value at multiple valuesof µ . We thereby obtain the quark mass in theRI-SMOM scheme at scale µ .3. Convert the mass to the MS scheme at scale µ us-ing a perturbative continuum matching calculation.We denote this conversion factor by Z MS / SMOM m ( µ ).4. Run all the MS quark masses at a range of scales µ to a reference scale of 3 GeV using the four loopQCD MS β function. We denote these running fac-tors r (3 GeV , µ ).5. Fit all of the results for the MS mass at 3 GeV to afunction that allows for discretisation effects andcondensate contamination, which begins at 1 /µ with the nongaugeinvariant (cid:104) A (cid:105) condensate.6. Obtain from the fit the physical value for the quarkmass in the MS scheme at 3 GeV with condensatecontamination removed.Here we provide three small updates to [44]. We firstlist them and then discuss them in more detail below.The three updates are:1 TABLE V. Lattice spacing values and tuned c quark massesat physical sea quark masses for each group of ensembles at afixed β value, denoted by their name in column 1 (see Table I).The lattice spacing value is given in units of w in column 2and the c quark mass, fixed from the J/ψ meson mass, isgiven in GeV units in column 3. The first uncertainty on themass is uncorrelated between lattice spacing values, and thesecond is correlated. w /a m tuned c [GeV]coarse 1.4055(33) 1.0524(10)(30)fine 1.9484(33) 0.9736(10)(30)superfine 3.0130(56) 0.8973(10)(30)ultrafine 3.972(19) 0.8592(20)(30) .
000 0 .
005 0 .
010 0 . a [fm ]5 GeV 4 GeV 3 GeV2 GeV 2.5 GeV a → , δ sea m → h JJ i [1408.4169]0 . . . . . m c ( G e V ) [ G e V ] FIG. 8. m c (3 GeV) extrapolated to the continuum with afit form that allows for condensate terms that behave likeinverse powers of the renormalisation scale µ . This plot is anupdated version of the upper section of Fig. 10 in [44], withadded data points on ultrafine lattices (set 14) and a retunedbare c quark mass fixed from the J/ψ , rather than η c , meson.The ultrafine points have slightly mistuned µ values comparedto the corresponding lines (see text).
1. we improve the uncertainty in the tuning of thebare lattice c quark mass by using the J/ψ massrather than the η c ;2. we include results from a finer ensemble of lattices(set 14) to provide even better control of the con-tinuum limit;3. we use the new 3-loop-accurate SMOM to MSmatching factors, Z MS / SMOM m , calculated in [46, 47]to reduce the perturbative matching uncertainty.The first update is to change how the tuning of the barecharm quark mass is done. In [44] bare charm masseswere used that had been tuned to the experimental η c mass, adjusted to allow for estimates of missing QED(from a Coulomb potential) and gluon annihilation ef-fects (from perturbation theory). A 100% uncertaintywas included on the adjustment [25]. Now that we areexplicitly including quenched QED it makes more sense to have a tuning process that uses an experimental me-son mass with no adjustments. We also want to use thesame tuning process for both the pure QCD case andthe QCD+QED case to allow for a clear comparison andone that reflects the procedures that would be followedin a complete QCD+QED calculation. This means thatwe should use the J/ψ meson mass, as we have done inSection III. The
J/ψ mass is more accurate experimen-tally than that of the η c and the J/ψ has a much smallerwidth, implying little effect on the mass from its 3-gluonannihilation mode. The impact of
J/ψ annihilation to asingle photon is a sub-1 MeV shift to the mass which isa 0.02% effect, so negligible.Using our
J/ψ meson masses and following the pro-cedure of [25] we obtain tuned bare c masses for each β value. These are given in Table V along with the w /a values corresponding to physical sea quark masses at that β value, which are also updated slightly from [44]. Theseslight changes in w /a lead to small adjustments in the µ values relative to those given in Table IV of [44]. Thisis accounted for when we run the Z m to the correct ref-erence scale in MS.The second update is to include results from the ultra-fine β = 7 . w /a value for physical sea quark masses isgiven in Table V. We have also calculated new Z SMOM m ( µ )values on set 14. These are given in Appendix A.The third update is to add the α s correction to theSMOM to MS conversion factor, Z MS / SMOM m , for the massrenormalisation. This correction was recently calculatedin continuum perturbative QCD [46, 47]. For n f = 4, ashere, the α s correction is a small effect ( 0.2%), continu-ing the picture seen at O ( α s ) and O ( α s ) and consistentwith the uncertainty taken from missing it in [44].Once we have determined results for m c (3 GeV) at avariety of lattice spacing values using the SMOM inter-mediate scheme at a variety of µ values, we need to fit theresults to determine m c (3 GeV) in the continuum limit.We do this following our previous calculation [44], wherethe fit function is given in Eq. (26). The fit includes dis-cretisation effects and condensate artefacts in the latticecalculation of Z SMOM m . In [44] we included a term in thefit, c α α s ( µ ) (with α s in the MS scheme) to allow for thethen-missing α s term in the SMOM to MS conversion.Here we replace that term with c α α s ( µ ) since the con-version is now calculated through α s and included in ourresults. We take a prior value on c α of 0 . ± .
5. Thisallows for the coefficient of the α s term in the conversionfactor to be 4 times as large as the α s coefficient.The updated fit to our results for m c in the MS schemeat 3 GeV is shown in Fig. 8. The fit has a χ / dof of0.71. The error budget for this calculation is shown inTable VI. Most of the entries are very similar to thosein [44]. The contribution due to the continuum extrap-olation has, unsurprisingly, dropped a little, as has theuncertainty from the missing higher order terms (here α s ) in the SMOM to MS conversion. The correlated tun-ing uncertainty comes largely from the uncertainty in the2 TABLE VI. Error budget (in %) for the calculation of thecharm quark mass in the MS scheme at a scale of 3 GeVusing RI-SMOM as an intermediate scheme. The listed con-tributions have the same meaning as those in [44] except thatwe use r here for the running factor rather than R and wehave an additional one labelled ‘QED effects’ which comesfrom the continuum extrapolation shown in Fig. 10. m c (3 GeV) a → α s term 0.10Condensate 0.21 m sea effects 0.00 Z MS / SMOM m and r Z SMOM m m tuned m tuned µ error from w physical value of w used to fix the lattice spacing. Wewill be able to reduce that uncertainty in future by im-proving the determination of the value of w .Our updated pure QCD result is m c (3 GeV) QCD = 0 . . (13)Running this down to a scale equal to the mass gives m c ( m c ) QCD = 1 . . (14)These results improve on and supersede the value in [44]. B. Impact of Quenched QED
To include quenched QED effects in the determinationof the c quark mass we must determine both the barequark mass and the mass renormalisation factor withquenched QED switched on.We include quenched QED effects for Z m in the RI-SMOM scheme in the same way as that described forthe vector current renormalisation in [16]. This involvesthe generation of U(1) fields in Landau gauge to remainconsistent with the Landau gauge QCD configurationsused in the pure QCD calculation. When convertingfrom the RI-SMOM scheme for Z m to the MS schemeit is also necessary to include QED effects in the per-turbative matching factor. We can evaluate the QEDcontribution to Z MS / SMOM m at order α QED by multiply-ing the known coefficient of α s [45] by (3/4) Q to give thecoefficient of α QED . The impact is very small ( < . O ( α QED ) term we do include forthe RI-SMOM to MS matching are given in Table VII.Note that the inclusion of QED in the lattice correlation
TABLE VII. Table giving factors needed for the determina-tion of the quark mass in a calculation including quenchedQED for different µ values and lattice spacings (denoted byset numbers). The fractional QED correction to Z SMOM m isgiven in the third column, the QED component of the RI-SMOM to MS matching for each µ in the fourth column andthe factor giving the QED mass running from µ to a referencescale of 3 GeV in the fifth column.Set µ [GeV] R QED [ Z SMOM m ] Z MS / SMOM m, QED r QED (3 GeV, µ )5 2 1.001200(83) 0.999872 0.99937210 2 1.001516(35) 0.999872 0.99937212 2 1.001853(83) 0.999872 0.9993725 2.5 1.000827(31) 0.999873 0.9997185 3 1.000540(15) 0.999873 -10 3 1.000851(11) 0.999873 -12 3 1.001308(18) 0.999873 -10 4 1.0005001(21) 0.999873 1.00044612 4 1.0009331(34) 0.999873 1.000446 .
000 0 .
005 0 .
010 0 .
015 0 . a [fm ]1 . . . . . C m FIG. 9. The factor C m that forms part of the QED effect inthe lattice to MS mass renormalisation, as defined in Eq. (15),plotted against the square of the lattice spacing. The factthat C m , as shown here, has almost no µ or a dependencedemonstrates that the µa dependence seen in R QED [ Z MS m ]from columns 3 and 4 of Table VII is simply that expectedfrom perturbation theory. functions is fully nonperturbative and it is only the con-tinuum matching and running that are simply done to O ( α QED ).To assess the QED impact on the tuned bare c mass weuse the QCD+QED J/ψ masses given in Table III andshown in Fig. 1. As we have corrected the pure QCD de-termination of m c to be tuned to the experimental J/ψ mass this is the tuning we will use for the QCD+QEDcase as well. The fractional shift in am c required to ob-tain the correct J/ψ mass after QED has been included(which we denote R QED [ am c ]) can be evaluated from thefractional change in the J/ψ mass. R QED [ am c ] is thefractional change in am c required to return the J/ψ massto the value it had in pure QCD (i.e. the experimentalvalue) once QED is switched on. Thus the increase in
J/ψ mass seen with QED must be compensated by a3 .
000 0 .
005 0 .
010 0 . a [fm ] 2 GeV3 GeV4 GeV0 . . . . . R Q E D [ m c ( G e V ) ] FIG. 10. QED correction to the charm quark mass in the MSscheme at a scale of 3 GeV. The different µ values, shown asdifferent colours and shapes, have all been run to 3 GeV andonly differ by discretisation and condensate effects. The redpoint on the left is the result for R QED [ m c (3 GeV)] returnedby the fit of Eq. (16). corresponding reduction in am c . From the deliberatelymistuned am c values in Table II we see that the fractionalchange in the c mass is 1.5 times larger than the changeseen in the meson mass. We can use this factor, alongwith R [ M J/ψ ] values from Table III and the requiredchange of sign in the shift, to determine the retuningof the quark masses in QCD+QED. We therefore take R QED [ am c ] on coarse, fine and superfine lattices (sets 6,10 and 12) to be 0.99840(8), 0.99790(4) and 0.99734(9)respectively. We increase the uncertainty on R QED [ am c ]compared to that for R [ M J/ψ ] by a factor of 5 toallow for the uncertainty in our conversion factor of 1.5above.We calculate the ratio of Z m in the SMOM schemein QCD+QED to pure QCD ( R QED [ Z SMOM m ]) followingthe methods of [44]. The calculations are carried out atmultiple values of the renormalisation scale µ and ex-trapolated to zero valence quark mass. Results are givenin Table VII. Notice that these values are larger than 1and so compensate to a large extent for the changes inthe tuning of the bare lattice am c induced by QED. Thisreflects the fact that most of the shift is an unphysicalultraviolet self-energy effect. We then convert from theSMOM scheme to MS by making use of the ratio of theperturbative conversion factor for QCD+QED to pureQCD, i.e. the O ( α QED ) piece of Z MS / SMOM m also given inTable VII. This is less than 1 but only by a very smallamount.Multiplying these two factors together gives the ra-tio of the lattice to MS mass renormalisation factors forQCD+QED to that for pure QCD, i.e. R QED [ Z MS m ]. FromTable VII, multiplying columns 3 and 4, we can see that R QED [ Z MS m ] varies with µ and with lattice spacing over arange of about 0.001. In perturbation theory we expect R QED [ Z MS m ] to consist of a power series in α QED and α s multiplied by constants and powers of logarithms of aµ . The leading logarithm at each order can be derived fromthe anomalous dimensions of the mass, allowing us towrite [48] R QED [ Z MS m ] = 1 + C m − α QED Q π log( µ a ) . (15)Here C m is a constant, up to discretisation effects andhigher order terms multiplying powers of log( aµ ). Fig-ure 9 plots our results for C m . These show very littlevariation with a and µ , confirming that the dependenceon a and µ of R QED [ Z MS m ] is almost entirely captured byEq. (15).Once the impact of QED on the c mass in the MSscheme at scale µ is obtained, as above, we then needto allow for QED effects in the running of the massesfrom µ to the reference scale of 3 GeV. This is done bymultiplication by a factor r QED calculated in O ( α QED )perturbation theory and given in Table VII. These num-bers are also very close to 1. The α s α QED term could inprinciple have some impact here but it is very small andwe neglect it [48].Multiplying R QED [ am c ] and R QED [ Z MS m ] together al-lows us to determine the ratio of the c quark mass in theMS scheme at 3 GeV from QCD+QED to that in pureQCD, i.e. R QED [ m c (3 GeV)]. The values that we havefor this ratio come from results at multiple values of µ and multiple values of the lattice spacing. To determinethe physical ratio in the continuum limit with condensatecontamination removed (in this case QED corrections toQCD condensates) we need to fit the results to a similarform to that used in [44]. We use R QED [ m c (3 GeV , µ, a )] = R QED [ m c (3 GeV)] × (16) (cid:34) α QED Q (cid:88) i =1 c ( i ) a ( a (1 GeV)) i (cid:35) × (cid:20) α QED Q (cid:18) (cid:88) j =1 c ( j ) µ a ( aµ ) j + (cid:88) n =1 α s ( µ ) c ( n )cond (1 GeV) n µ n (cid:19)(cid:21) . The first term on the second line of Eq. (16) accountsfor discretisation effects in R QED [ am c ]; a scale of 1 GeVis chosen in these effects as this is close to the c quarkmass. The term multiplying this (on the bottom twolines) models the a and µ dependence of the QED cor-rection to Z m . This includes discretisation effects of theform ( aµ ) i and terms to model condensate contribu-tions, starting at 1 /µ . The priors on all coefficientsare taken as 0 . ± .
0, except for the physical result, R QED [ m c (3 GeV)], for which we take prior 1.00(1).The lattice QCD results for R QED [ m c (3 GeV)] and thefit output are shown in Figure 10. The fit has a χ / dof of0.87 and returns a physical value of R QED [ m c (3 GeV)] of0.99823(17). We conclude that the impact of quenchedQED is to lower the c quark mass, m c (3 GeV) by a tinyamount: 0.18(2)%.4 .
20 1 .
25 1 .
30 1 .
35 1 . m c ( m c , n f = 4)(GeV) HPQCD HISQ RI-SMOMHPQCD HISQ RI-SMOMHPQCD HISQ JJcFNAL/MILC/TUM HISQ MRSETMC twisted mass RI-MOM
PDG u,d,s,c sea + qQEDu,d,s,c sea
FIG. 11. Comparison of lattice QCD results for m c that in-clude u , d , s and c quarks in the sea. The top two results arethe ones from this paper. Our QCD + quenched QED result isgiven in Eq. (17). Our pure QCD result, Eq. (13), supersedesour earlier result in [44]. The Fermilab/MILC/TUM resultis from [49] and uses a method based on charm-light mesonmasses. The ‘HPQCD HISQ JJc’ result is from [25] and usescurrent-current correlator techniques. These three resultsagree to better than 1%. The ETMC result is from [50] anduses the RI-MOM intermediate scheme. The grey band givesthe ± σ uncertainty band from the Particle Data Group [1]. We obtain our final answer for the c quark massin QCD+QED by multiplying R QED [ m c (3 GeV)] byour pure QCD result for m c (3 GeV). This gives theQCD+QED result of m c (3 GeV) QCD+QED = 0 . . (17)Running down to the scale of the mass with QCD+QEDgives: m c ( m c ) QCD+QED = 1 . , (18)very close to the pure QCD value at this scale. This is thefirst determination of the c quark mass to include QEDeffects explicitly, rather than estimate them phenomeno-logically as has been done in the past. The uncertaintyachieved here of 0.5% is smaller than the 0.6% from [44]because we have reduced several sources of uncertainty,mainly those from the extrapolation to the continuumlimit and from missing higher order terms in the SMOMto MS matching. C. Discussion: m c Figure 11 gives a comparison of lattice QCD results for m c . We plot the masses at the scale of the mass, as isconventional even though this is a rather low scale. Werestrict the comparison to results that were obtained on gluon field configurations including u , d , s and c quarksin the sea. The top result is our value from Eq. (17) thatexplicitly includes a calculation of the impact of quenchedQED on the determination of the quark mass.The top three results in the pure QCD section of thefigure all include an estimate of, and correction for, QEDeffects. These corrections are made, however, by allow-ing for ‘physical’ QED effects such as those arising fromthe Coulomb interaction between quark and antiquarkin a meson. They do not allow for the QED self-energycontribution which is substantial. Although a large partof this is cancelled by the impact of QED on the massrenormalisation, a consistent calculation has to includeboth effects, as we have done here.An important point about Figure 11 is that the topthree pure QCD results all have uncertainties of less than1% and agree to better than 1%, using completely differ-ent methods. This implies a smaller uncertainty on m c than the 1.5% allowed for by the Particle Data Group [1].This impressive agreement is not changed by our new re-sult including quenched QED because, as we have shown,the impact of this is at the 0.2% level. V. J/ψ
AND η c DECAY CONSTANTS
The decay constant of the
J/ψ , f J/ψ , is defined fromthe matrix element between the vacuum and a
J/ψ mesonat rest by (cid:104) | ψγ µ ψ | J/ψ (cid:105) = f J/ψ M J/ψ (cid:15) µ , (19)where (cid:15) µ is the component of the polarisation of the J/ψ in the direction of the vector current. In terms of theground state amplitude, A V , and mass, M V ( ≡ E V ),obtained from the fit of Eq. (3) to the charmonium vectorcorrelator it is (in lattice units) f J/ψ = Z V (cid:115) A V M V . (20) Z V is the renormalisation factor required to match thelattice vector current to that in continuum QCD if a non-conserved lattice vector current is used (as here). Wediscuss the renormalisation of vector currents using in-termediate momentum-subtraction schemes in [16] andwe will make use of the results based on the RI-SMOMscheme here (see Section II). Note that there is no ad-ditional renormalisation required to get from the RI-SMOM scheme to MS because the RI-SMOM schemesatisfies the Ward-Takahashi identity [16].The partial decay width of the J/ψ to an (cid:96) + (cid:96) − pair( (cid:96) = e, µ ) is directly related to the decay constant. Atleading order in α QED and ignoring ( m (cid:96) /M J/ψ ) correc-tion terms, the relation isΓ( J/ψ → (cid:96) + (cid:96) − ) = 4 π α , eff ( M J/ψ ) Q c f J/ψ M J/ψ , (21)5where Q c is the electric charge of the charm quark inunits of the charge of the proton. Note that the for-mula contains the effective coupling, α QED , eff evaluatedat the scale of M J/ψ but without including the effectof the
J/ψ resonance in the running of α QED to avoiddouble-counting [51].Experimental values of Γ(
J/ψ → e + e − ) are obtainedby mapping out the cross-section for e + e − → e + e − and e + e − → hadrons through the resonance region [52] or byusing initial-state radiation to map out this region via e + e − → µ + µ − γ [53]. In either case initial-state radi-ation and non-resonant background must be taken careof [54, 55]. A cross-section fully inclusive of final-state ra-diation is obtained; interference between initial and final-state radiation is heavily suppressed [56]. The resonanceparameter determined by the experiment is then the ‘full’partial width [55, 57],Γ (cid:96)(cid:96) = Γ (0) (cid:96)(cid:96) | − Π | (22)where Γ (0) is the partial width to lowest order in QEDand Π is the photon vacuum polarisation. The effectof the vacuum polarisation is simply to replace α QED in the lowest-order QED formula for the width with α QED , eff ( M ), as we have done in Eq. (21).The experimental determination of Γ (cid:96)(cid:96) is accurate to2% for the J/ψ [1]. This allows us to infer a decayconstant value from experiment, accurate to 1%, usingEq. (21). f expt J/ψ = (cid:18) M J/ψ πQ c (cid:19) / Γ / e + e − α QED , eff (23)= 40 . / Γ / e + e − α QED , eff . Using the experimental average of Γ e + e − = 5.53(10)keV [1], and α QED , eff ( M J/ψ ) = 1 / . f expt J/ψ = 406 . . .
5) MeV . (24)The first uncertainty comes from the experimental uncer-tainty in Γ and the second is an O ( α QED /π ) uncertaintyfor higher-order in QED terms, for example from final-state radiation, in the connection between Γ and f inEq. (21). Note that using α QED of 1/137 would increasethis number by 2.3% (9 MeV).This experimental value can then be compared to ourlattice QCD results for a precision test of QCD. Here weimprove on HPQCD’s earlier calculation [29] by workingon gluon field configurations that cover a wider range oflattice spacing values and with sea u/d quark masses nowdown to their physical values. In addition we now include c quarks in the sea and have a more accurate determi-nation of the vector renormalisation factor Z V [16]. Wewill also test the impact on f J/ψ of the c quark’s electriccharge. TABLE VIII. Columns 2 and 3 give the results in lattice unitsfor the
J/ψ and η c decay constants respectively in pure QCDon the ensembles listed in Table I. The values of am val c arethose given in Table I except for two cases with a deliber-ately mistuned c quark mass: set 6 denoted by a * where am c = 0 .
643 and set 14 denoted by a † where am c =0.188.The results for f J/ψ do not include the multiplication by Z V needed to normalise them (Eq. (20)). Columns 4 and 5 givethe electromagnetic corrections for the J/ψ and η c decay con-stants, as the ratio of the QCD+QED result to that in pureQCD. Again, the electromagnetic corrections for f J/ψ do notinclude the corrections to Z V . Z V values are given in Ta-ble IX.Set af J/ψ /Z V af η c R [ f J/ψ /Z V ] R [ f η c ]1 0.43370(55) 0.37659(18) - -2 0.42346(48) 0.370332(91) 1.00410(64) 1.00294(50)3 0.4163(11) 0.366127(57) - -4 0.29411(21) 0.268331(61) - -6 0.28835(15) 0.263727(60) 1.00341(37) 1.00326(13)6 ∗ † The decay constant of the pseudoscalar η c meson isdetermined from our pseudoscalar correlators (of spin-taste γ ⊗ γ ) using the ground-state mass and amplitudeparameters from the correlator fit, Eq. (2): f η c = 2 m c (cid:115) A P ( M P ) . (25)Note that this is absolutely normalised and no Z factor isrequired. Because the η c does not annihilate to a singleparticle there is no experimental process from which wecan directly determine f η c . Nevertheless it is a usefulquantity to calculate for comparison to f J/ψ and to fillout the picture of these hadronic parameters from latticeQCD [59]. Again we will improve on HPQCD’s earliercalculation [60] as discussed above for the
J/ψ . A. Pure QCD
The second column of Table VIII gives our results forthe (unnormalised) values of af J/ψ in pure QCD on 12 ofthe sets from Table I. We multiply af J/ψ /Z V by the valueof Z V and convert to physical units using the inverselattice spacing. Z V values are taken from [16] except fora new value calculated here for β = 7 . Z V results. The Z V values are very6 TABLE IX. Vector current renormalistion constants, Z V ( µ ),using the RI-SMOM scheme at µ = 2 GeV (column 2) and µ =3 GeV (column 3) in pure QCD for each β value correspondingto a group of ensembles in Table I. Column 4 gives the QEDcorrection to Z V at 2 GeV in the form of the ratio of theQCD+QED value to that of pure QCD. Most of these valuesare taken from [16] although the Z V value at β = 7 and theQED correction at β = 5 . β Z V (2 GeV) Z V (3 GeV) R QED [ Z V (2 GeV)]5.80 0.95932(18) - 0.999544(14)6.00 0.97255(22) 0.964328(75) 0.999631(24)6.30 0.98445(11) 0.977214(35) 0.999756(32)6.72 0.99090(36) 0.98702(11) 0.999831(43)7.00 0.99203(108) 0.99023(56) - . . . . . am c ) . . . . . f J / ψ [ G e V ] FIG. 12. The
J/ψ decay constant calculated on the ensem-bles of Table I in pure QCD and plotted against the squareof the bare c quark mass in lattice units. The different redshapes correspond to different groups of ensembles with sim-ilar lattice spacing and the error bars shown include the fulluncertainty on the points. Points at mistuned m c are notplotted but are included in the fit. The green curve marksour extrapolation to the physical point, where the black crossshows the result determined from the experimental averagefor Γ e + e − from Eq. (24). precise and so have little impact on the uncertainty in f J/ψ .Our results for f J/ψ for the pure QCD case are shownin Fig. 12 plotted against the square of the lattice spacing(in units of the bare c quark mass). Clear dependence onthe lattice spacing is seen. This dependence comes fromthe amplitudes of the two-point correlators; the latticespacing dependence of Z V contributes very little to it.We also plot in Fig. 12 the results of the fit using Eq. (5).The priors for the fit are as given in Section II C with theprior on the physical value of f J/ψ (i.e. x in Eq. (5)) of0.4(1). The χ / dof of the fit is 0.43. The agreement withthe result derived from experiment can clearly be seen.We obtain an f J/ψ value in pure QCD of f J/ψ,
QCD = 409 . .
6) MeV . (26)We will discuss this result further in Section V C. . . . . . am c ) . . . . . f J / ψ [ G e V ] FIG. 13. The continuum extrapolation of f J/ψ using Z V (2GeV) (results and fit curve in red, same values as in Fig. 12)and Z V (3 GeV) (results and fit curve in blue). The continumextrapolated results agree in the two cases as they should.A black cross shows the experimental average result fromEq. (24). We have used vector current renormalisation factors, Z V , in the RI-SMOM scheme at a scale of 2 GeV. The µ dependence of Z V should just be the result of discretisa-tion effects and results for the physical quantity, f J/ψ , us-ing different renormalisation scales µ should agree in thecontinuum limit. Here we verify that this is the case us-ing Z V (2 GeV) and Z V (3 GeV) results from Table IX [16].There is no 3 GeV result on the very coarse lattices since µa would be too large. The comparison for f J/ψ using µ =2 and 3 GeV is shown in Fig. 13 for the pure QCDcase. The difference between the two values of µ is barelyvisible. The values at µ = 3 GeV give a continuum limitresult of f J/ψ = 408.7(1.8) MeV, in good agreement withthat at µ = 2 GeV in Eq. (26) but slightly less accurate.The χ / dof of the fit was 0.45.Our results for af η c , the η c decay constant in latticeunits, are given in the third column of Table VIII forthe pure QCD case. After conversion to physical units,they are plotted in Figure 14. The curve is similar tothat for f J/ψ but with somewhat smaller discretisationeffects. We also plot the results of performing the samefit as for f J/ψ using Eq. (5). The χ / dof of the fit is 0.88giving a result for the decay constant in pure QCD of f η c , QCD = 397 . .
0) MeV . (27)This agrees well with the earlier HPQCD value on n f =2 + 1 gluon field configurations [60] of 0.3947(24) GeVbut has half the uncertainty. In [60] the effects fromneglecting the charm quark in the sea are estimated tobe O (0 . f J/ψ to f η c in pure QCD. A lot of the discretisation effects cancelin the ratio, as is evident in comparing this figure toFigures 12 and 14. Systematic uncertainties, for examplefrom the determination of the lattice spacing, are also7 . . . . . am c ) . . . . . f η c [ G e V ] HPQCD [1008.4018]
FIG. 14. The η c decay constant calculated on the ensemblesof Table I in pure QCD and plotted against the square of thebare c quark mass in lattice units. The different red shapescorrespond to different groups of gluon field ensembles withsimilar lattice spacing. The error bars on each point are thefull uncertainty, including correlated uncertainties from, forexample, the determination of the lattice spacing. The greencurve shows our fit and extrapolation to the physical point.The black cross gives the earlier HPQCD result on n f = 2 + 1gluon field configurations from [60]. . . . . . am c ) . . . . . . f J / ψ / f η c FIG. 15. The ratio of
J/ψ to η c decay constant determinedin pure QCD, plotted against the square of the bare c quarkmass in lattice units. The different red shapes correspond todifferent groups of gluon field ensembles with similar latticespacing. The error bars on each point are the full uncertainty,including correlated uncertainties from, for example, the de-termination of the lattice spacing. The green curve shows ourfit and extrapolation to the physical point. reduced. The shape of the curve again, as in the hyperfinesplitting case, reflects the fact that we have successfullyreduced sources of a error to the point where a and a are visible.We fit the ratio to the same fit as before (Eq. (5)) witha prior on the physical value of 1.0(1). The fit has a χ / dof of 0.62 and returns a physical value for the decayconstant ratio in pure QCD of f J/ψ,
QCD f η c , QCD = 1 . . (28) TABLE X. Error budget for the
J/ψ and η c decay constantsas a percentage of the final answer. f J/ψ f η c a → Z V w /a w . . . am c ) . . . . R Q E D [ f p ] p = J/ψp = η c FIG. 16. The fractional QED correction to the
J/ψ and η c de-cay constants as a function of lattice spacing. The horizontaldashed lines mark the weighted average of the points. Thus we see that the
J/ψ decay constant is nearly 3%larger than that of the η c with an uncertainty of 0.2%.Table X gives the error budget for our final valuesof f J/ψ and f η c , both for the pure QCD case and theQCD+QED case to be discussed in Section V B. Thecontributions from different sources are very similar be-tween the two decay constants. It is clear from this thatthe dominant sources of error are related to the determi-nation of the lattice spacing, as for the hyperfine split-ting.The error budget presented here for the J/ψ decay con-stant is markedly different from that of [29]. There thedominant contribution to the error was from the vectorrenormalisation constant, Z V , obtained using a matchingbetween lattice time moments and high order perturba-tive QCD [61]. Here that error is substantially reducedby using the Z V values obtained in lattice QCD fullynonperturbatively in the RI-SMOM scheme [16]. B. Impact of Quenched QED
Including quenched QED effects into our calculationsallows us to determine the effect on the
J/ψ and η c decayconstants of the electric charge of the valence c quarks.8 .
025 0 .
030 0 .
035 0 . /L s . . . . R Q E D [ f p ] p = J/ψp = η c FIG. 17. The volume dependence of the fractional QED effecton the
J/ψ and η c decay constants measured on sets 5-7. Thedashed lines are horizontal and indicate the weighted averageof the points. There is no observable finite volume effect atthe level of our statistical uncertainties. . . . . . am c ) . . . . . f J / ψ [ G e V ] FIG. 18. The
J/ψ decay constant calculated on the ensem-bles of Table I in pure QCD (red points) and including alsoquenched QED (blue points) plotted against the square ofthe bare c quark mass in lattice units. The green curve marksour extrapolation to the physical point, where the black crossshows the experimental average result from Eq. (24). Because the
J/ψ and η c are electrically neutral particles,there is no long-distance infrared component to causeproblems (as there is for f π + ) and we can simply proceedto determine the decay constants after the addition ofthe QED field as we did in the pure QCD case.The fractional QED effect on the J/ψ and η c decayconstants at fixed bare c quark mass in lattice units isgiven in Table VIII. We see a 0.3% increase, offset slightlyby the change in Z V in the J/ψ case. The fractionalQED effect on Z V is given in Table IX. The fractionalQED effect at fixed bare mass is plotted in Figure 16.We see that the effect is similar for the J/ψ and η c in thecontinuum limit and shows very little dependence on thelattice spacing.The volume dependence of the fractional QED effect isshown in Figure 17 on sets 5–7. We find that the effectis negligible well below the 0.1% level. . . . . . am c ) . . . . . f η c [ G e V ] HPQCD [1008.4018]
FIG. 19. The η c decay constant calculated on the ensem-bles of Table I in pure QCD (red points) and including alsoquenched QED (blue points) plotted against the square of thebare c quark mass in lattice units. The green curve marks ourextrapolation to the physical point. The black cross gives theearlier HPQCD result on n f = 2+1 gluon field configurationsfrom [60]. We now combine our QCD+QED results with our pureQCD results and the full fit of Eq. (5), which takes intoaccount the retuning of the quark mass needed whenquenched QED is included. Figure 18 shows our pureQCD results, QCD+QED results and fit curve for the
J/ψ decay constant. We obtain f J/ψ,
QCD+QED = 410 . .
7) MeV . (29)This is a 0.2% increase over the value in pure QCD(Eq. (26)) because retuning reduces the quark massand offsets some of the impact of quenched QEDseen in Table VIII. A more accurate statement isthat the final fractional effect from quenched QED is R QED [ f J/ψ ] =1.00188(36).A very similar picture is seen for f η c in Figure 19. Weobtain f η c , QCD+QED = 398 . .
0) MeV . (30)The final fractional effect from quenched QED is then R QED [ f η c ] =1.00166(25).Finally, in Figure 20 we plot results for the ratio of J/ψ to η c decay constants and show the fit curve extrapolatedto the continuum limit. This gives f J/ψ,
QCD+QED f η c , QCD+QED = 1 . . (31)This is almost the same as the pure QCD result. C. Discussion : f J/ψ and f η c Figure 21 compares our new pure QCD andQCD+QED results to previous results including n f =2 + 1 flavours of sea quarks for f η c [60] and f J/ψ [29].9 . . . . . am c ) . . . . . . f J / ψ / f η c FIG. 20. The ratio of
J/ψ to η c decay constants calculatedon the ensembles of Table I in pure QCD (red points) andincluding also quenched QED (blue points) plotted againstthe square of the bare c quark mass in lattice units. Thegreen curve marks our extrapolation to the physical point. .
38 0 .
39 0 .
40 0 .
41 0 . f (GeV) HPQCD HISQHPQCD HISQ 2010/12HPQCD HISQ
J/ψη c u,d,s,c sea + qQEDu,d,s,c seau,d,s sea FIG. 21. A comparison of our new results for f η c and f J/ψ with earlier lattice QCD results, also by HPQCD, on gluonfield configurations that include n f = 2 + 1 flavours of quarksin the sea. The results labelled ‘HPQCD HISQ’ are from thispaper and the results labelled ‘HPQCD HISQ 2010/12’ arefrom [29, 60]. The grey bands are the ± σ bands from ournew QCD+QED results. There is good agreement. These earlier calculations alsoused HISQ quarks but our new results are more accurate,particularly for f J/ψ because of the use of more accuratevalues of Z V [16].There have also been calculations that use n f = 2flavours of sea quarks. It is harder to make a comparisonto these results because it is not clear what the system-atic error is from not including at least the s quarks inthe sea, and no uncertainty is included for this. The cal-culation of [62] uses twisted-mass quarks on n f = 2 gluon . . . . . . Γ e + e − (keV) BES III 16HPQCD HISQKEDR 18PDGu,d,s,c sea + qQEDexperiment
FIG. 22. A comparison of the width for
J/ψ decay to e + e − implied by our new results for f J/ψ to that obtained fromtwo recent experiments. The point labelled ‘KEDR 18’ isfrom [52] and the point labelled ‘BES III 16’ is from [53].The grey bands are the ± σ bands from the Particle DataGroup average [1]. field configurations and obtains f η c =387(7)(2) MeV and f J/ψ =418(8)(5) MeV. The calculation of [63] uses cloverquarks on n f = 2 CLS gluon-field configurations to give f η c =387(3)(3) MeV and f J/ψ =399(4)(2) MeV. The n f = 2 results for f η c agree with each other and havea central value about 2 σ below ours. The σ here is thatfrom the n f = 2 results since our uncertainty is muchsmaller. The n f = 2 results for f J/ψ are compatible witheach other and with our result, again at 2 σ .As discussed in Section V, the J/ψ decay constant isthe hadronic quantity that is needed to determine therate of
J/ψ annihilation to leptons via Eq. (21). Ourresult for Γ using Eq. (21) with f J/ψ from Eq. (29) isΓ(
J/ψ → e + e − ) = 5 . . (32)The first uncertainty is from our lattice QCD+QED re-sult for f J/ψ and the second uncertainty allows for a rela-tive O ( α QED /π ) correction to Eq. (21) from higher-ordereffects.Figure 22 compares this width Γ( J/ψ → e + e − ) to re-sults from experiment. Recent experimental results fromKEDR [52] and BES III [53] are shown along with theParticle Data Group average [1] (as a grey band). Fig-ure 22 shows good agreement between our result and theexperimental values shown, as well as the experimentalaverage.The lattice QCD result is now more accurate than theexperimental values.0 VI. VECTOR CORRELATOR MOMENTS AND a cµ With new results expected from the Fermilab g − R ( e + e − → hadrons), it is also important tocompare lattice QCD results to these, disaggregated byflavour where possible.The first calculation of the quark-line connected c -quark contribution to the HVP, a cµ , was given in [64]using results for the time-moments of vector charmoniumcurrent-current correlators calculated in [29]. The timemoments are defined by G n = Z V (cid:88) t t n C V ( t ) , (33)where C V ( t ) is the vector current-current correlator and Z V is the vector current renormalisation factor, dis-cussed in Section V. Note that t ∈ {− L/ , − L/ , . . . , L/ } .The even-in- n time-moments for n ≥ q = 0 of the renormalised vacuumpolarisation function [65], ˆΠ( q ) ≡ Π( q ) − Π(0), by G n = ( − n/ ∂ n ∂q n q ˆΠ( q ) (cid:12)(cid:12)(cid:12)(cid:12) q =0 . (34)This means that ˆΠ( q ) can be reconstructed, using Pad´eapproximants, from the G n [64] and fed into the integralover q that yields the quark-line connected HVP con-tribution to a µ [66]. Only time-moments of low momentnumber are needed to give an accurate result for a cµ be-cause the integrand is dominated by small q . We willgive results for the four lowest moments, n = 4 , , , c quarks in the sea and calculate, rather than estimate,the impact of the leading QED effects.ˆΠ( q ) and hence a cµ can also be determined from ex-perimental results for R c ≡ R ( e + e − → cc → hadrons) asa function of squared centre-of-mass energy, s . This canbe done using inverse- s moments M k ≡ (cid:90) dss k +1 R c ( s ) . (35) R c is obtained from the full e + e − rate from just belowthe c threshold upwards by subtracting the backgroundcontribution from u , d , and s quarks perturbatively, seee.g. [61]. . . R [ G / ]0 . . R [ G / ]0 . . R [ G / ]0 . . R [ G / ]0 . . . . . . . . /L s . . R [ a cµ ] FIG. 23. A study of the volume dependence of the electromag-netic correction to the first four time moments of the vectorcurrent-current correlator and to a cµ on sets 5-7. There is noobservable dependence on the lattice spatial extent, L s , ascan be judged by comparison to the dashed horizontal linesat the weighted averages of the points. The relationship between G k +2 and M k is then G k +2 = (2 k + 2)! M k π Q . (36)A comparison of our correlator time-moments calculatedon the lattice and extrapolated to the continuum limitto the inverse- s moments determined from experiment isequivalent to a test of the agreement of the results for a cµ in the two cases. A. Vector correlator moments: Pure QCD andQCD+QED results
Table XI gives our raw results for the time-momentsof the same vector current-current correlators from whichwe have determined the mass and decay constant of the
J/ψ meson in Sections III and V. Notice that the sta-tistical uncertainties are tiny. The correlators make useof a local vector current that must be renormalised asdiscussed in Section V. The results in Table XI are cal-culated before renormalisation and are given in latticeunits. The quantity that is tabulated is (cid:18) G n Z V (cid:19) / ( n − . (37)We take the ( n − Z V values at µ = 2 GeV given in Table IXthat were used for f J/ψ in Section V.Table XI also gives the result of including quenchedQED as the ratios R for each rooted moment (atfixed am c ). These values are all very slightly less than 1,1 TABLE XI. Time-moments of charmonium vector current-current correlators calculated on the ensembles in Table I. The resultstabulated are values of ( G n /Z V ) / ( n − in lattice units along with their statistical uncertainties, given for n = 4, 6 , 8 and 10in columns 2, 3, 4 and 5. Uncertainties are statistical only. Note that the results for different moments are correlated becausethey are determined from the same correlation functions. The results marked with a ∗ and † are for deliberately mistuned am c values as detailed in the caption to Table II. Also included are the ratios at fixed am c , denoted R , of QCD+QED to QCDresults for each rooted moment on a subset of ensembles.Set n = 4 n = 6 n = 8 n = 10 R [ n = 4] R [ n = 6] R [ n = 8] R [ n = 10]1 0.389670(40) 0.949791(62) 1.410524(75) 1.815497(88) - - - -2 0.396283(22) 0.961260(35) 1.425498(42) 1.833868(49) 0.999954(26) 0.999910(17) 0.999858(15) 0.999810(15)3 0.400779(15) 0.969045(24) 1.435671(28) 1.846369(33) - - - -4 0.511194(12) 1.164351(19) 1.701040(26) 2.184698(34) - - - -6 0.5206344(85) 1.181180(14) 1.724311(19) 2.214708(24) 0.9998455(15) 0.9997169(11) 0.9995987(10) 0.9994908(11)6 ∗ † by up to 0.1% for n =4, and 0.2% for n =10. We can alsotest the finite-volume dependence of the quenched QEDeffect using sets 5, 7 and 8 and the results are shownin Fig. 23. There is no visible volume dependence in theQED effect on time-moments at the level of our statisticaluncertainties. This is to be expected, as seen for the J/ψ mass and decay constant in Sections III and V, since thevector current being used here is electrically neutral.To fit the results as a function of lattice spacing it isconvenient to work with the dimensionless combination: M J/ψ ( G n ) / ( n − , (38)using our M J/ψ masses from Table II. This reduces theuncertainty coming from the value of the lattice spacing.At the same time it also reduces the sensitivity to mis-tuning of the valence c quark mass. We fit the quantitydefined in Eq. (38) as a function of lattice spacing, allow-ing for quark-mass mistuning effects of both valence andsea quarks to derive results in the physical continuumlimit. We do this as before using the fit function given inEq. (5). Values in the physical continuum limit are thendivided by M J/ψ from experiment to obtain our final re-sults for G / ( n − n in GeV − . In doing this we allow anuncertainty of 0.7 MeV in M J/ψ from annihilation to aphoton (see Section III) since this effect is not includedin our results.Fig. 24 plots the results for the rooted moments mul-tiplied by M J/ψ as a function of ( am c ) . Also shown isthe fit result from Eq. (5). Only the pure QCD latticeresults are shown for clarity; those including the effectof quenched QED are very close to them. The fit resultplotted is that for the QCD+QED case.Table XII gives our QCD+QED results for the 4th, TABLE XII. QCD+QED results in the physical continuumlimit for the first four time-moments (column 2) comparedwith the results extracted from experiment in column 3 [67].Agreement within 2 σ is seen for all except the 4th moment,but the lattice QCD results are much more accurate. Column4 gives the effect of quenched QED as a ratio of the physicalresults in QCD+QED to those in pure QCD. n G / ( n − n ( G exp .n ) / ( n − R QED [ G / ( n − n ]GeV − GeV − − χ / dof of 0.62 for an svdcut of 1 × − .Column 3 of Table XII gives the results derived fromexperimental data for R ( e + e − → hadrons) and Γ (cid:96)(cid:96) andmass for the J/ψ and ψ (cid:48) by [67] for comparison (seeSection VI). We have converted these results into thequantity that we calculate according to Eq. (36). Wesee agreement within 2 σ for n ≥ σ tension at n = 4. Theresults given from [67] correspond to their standard se-lection of experimental datasets. Results shift by ± σ for other selections. Note that the lattice QCD resultsare now much more precise than those determined fromexperiment. The comparison between lattice QCD andexperiment will be discussed further in Section VI B.Column 4 of Table XII gives the final impact of2 TABLE XIII. Error budget for G / ( n − n as a percentage ofthe final answer. G / G / G / G / a → Z V M J/ψ w w /a quenched QED, including the quark mass retuning, onthe moments. We see that the ratios, R QED [ G / ( n − n ],are greater than 1, in contrast to the R of Table XIthat are less than 1. As discussed in Section II B, theinclusion of QED means that the c quark mass must beretuned downwards. The rooted time-moments are ap-proximately inversely proportional to the quark mass andso they increase under this retuning, more than offsettingthe direct effect of QED seen in R . . . . . . am c ) . . . . . . . . . . . ( n t h m o m e n t) / ( n − ) × M J / ψ FIG. 24. The four lowest time moments and their extrap-olation to a = 0. The symbols give results for the rootedmoment multiplied by M J/ψ with different symbols denotingdifferent groups of ensembles with similar lattice spacing. Thediffferent colours pick out the different moments, from n = 4at the bottom to n = 10 at the top of the plot. Uncertain-ties are too small to be visible. The extrapolation for all themoments are performed simultaneously including correlationsbetween moments determined from the same current-currentcorrelators. Only the results from the pure QCD calcula-tion at well-tuned c quark masses are shown for clarity. Thedashed lines give the QCD+QED fit curve using Eq. (5) andthe coloured crosses mark the result in the continuum limitat physical quark masses. Table XIII gives the error budget for the G n . The totaluncertainty in all cases is below 0.2%. .
305 0 .
310 0 .
315 0 .
320 0 . ( G ) / (GeV − ) HPQCD HISQ 2012K¨uhn et al 2007Dehnadi et al 2011HPQCD HISQu,d,s,c sea + qQEDu,d,s seapheno. .
660 0 .
665 0 .
670 0 .
675 0 . ( G ) / (GeV − ) JLQCD DW 2016K¨uhn et al 2007HPQCD HISQ 2012Dehnadi et al 2011HPQCD HISQu,d,s,c sea + qQEDu,d,s seapheno.
FIG. 25. Comparison of our new result to those of pre-vious lattice QCD calculations for the 4th and 6th time-moments (appropriately rooted) of the charmonium vectorcurrent-current correlator. Our new result obtained here on n f = 2 + 1 + 1 gluon field configurations and including theeffect of quenched QED is given at the top (blue cross).HPQCD’s 2012 result on n f = 2 + 1 gluon field configura-tions with HISQ valence c quarks is marked ‘HPQCD HISQ2012’ [29]. JLQCD’s 2016 result using n f = 2 + 1 domain-wall quarks is marked ‘JLQCD DW 2016’ for the 6th mo-ment only. We also compare to results (open red squares)denoted ‘pheno.’ that are derived from experimental data forthe cross-section for e + e − to hadrons as a function of centre-of-mass energy, by determining the cc component. The pointsplotted come from [61] and [67]. Open orange circles show al-ternative selections of datasets from [67]; the upper value isfor the ‘maximal’ set and the lower value for the ‘minimal’set. B. Discussion: Vector correlator moments
Our new results for the time-moments of vectorcurrent-current correlators improve significantly on ear-lier lattice QCD calculations and are now more accuratethan results derived from experiment.Figure 25 compares our new results for the 4th and 6thmoments to earlier lattice QCD results. Comparison for3the 8th and 10th moments gives a very similar pictureand so is not shown. The first lattice QCD calculation ofthe time-moments of vector charmonium current-currentcorrelators was given by HPQCD in [29] using HISQ va-lence quarks on n f = 2 + 1 asqtad gluon field configu-rations. Our new results have an uncertainty almost tentimes smaller than these. The error budget in the earlierresults was dominated by the uncertainty in Z V from theuse of continuum perturbation theory in the matchingfactors and there was also a sizeable uncertainty fromthe lattice spacing. These uncertainties have been enor-mously reduced here and in addition we no longer needan uncertainty from missing QED effects. We also show acomparison with the subsequent results from the JLQCDcollaboration [68] using domain-wall quarks on n f = 2+1gluon field configurations. JLQCD do not give a valuefor the 4th moment because of discretisation effects intheir formalism (tree-level O ( a )). The dominant uncer-tainties in their results are from statistics and from thevalue of t / used to fix the lattice spacing, on which theyhave a 2% uncertainty. Good agreement is seen for all ofthe lattice QCD results.Two of the most recent results from phenomenologi-cal determinations of the moments [61, 67] are also com-pared in Fig. 25. The results from [67] include exper-imental datasets for the inclusive cross-section that areboth older and newer than those used in [61]. Resultsfrom [67]’s ‘standard’ selection of datasets were given inTable XII and are shown in Fig. 25 in red. We alsoshow, in orange, the results from the ‘maximal’ set (allexperimental information available at that point) and the‘minimal’ set (datasets that are needed to cover the full √ s range from 2 GeV to 10.5 GeV without gaps, keepingthe most accurate results). Note that the resonance pa-rameters are the same for all selections. We see that thevariation with dataset selection covers almost 1 σ for the4th moment, but much less for the 6th moment. This isalso reflected in the differences between [67] and [61].These phenomenological analyses must subtract the‘non-charm’ background from experimental results for R ( e + e − → hadrons) to leave R c for Eq. (35). R c isdefined to be the result from diagrams with a charmquark loop connected to a photon at both ends [61]i.e. the quark-line connected vector current-current cor-relator that we study on the lattice. The subtractedbackground includes QED effects for the non-charm andsinglet (quark-line disconnected) contributions. The re-mainder, R c then includes the QED effects associatedwith the cc loop. The dominant source of uncertaintyin R c comes from the charmonium resonance ( J/ψ and ψ (cid:48) ) region and is set by the uncertainty in Γ ee for thesestates. The fractional uncertainty is approximately thesame for all moments [61, 67]. When the ( n − n . Good agreement is seen between the phenomenologicalresults and our new lattice results for n = 6, 8 and 10, al-though the lattice results are systematically at the upper TABLE XIV. Values of a cµ on the ensembles of Table I andthe direct quenched QED correction on a subset of those en-sembles. Those marked with a ∗ and † are at deliberatelymistuned c masses (see caption to Table II). The uncertain-ties quoted are correlated through the value of M J/ψ (for allensembles, see text) and Z V (for ensembles at a given β ).Set a cµ × R (cid:2) a cµ (cid:3) ∗ † end of the phenomenological range. The largest discrep-ancy is a 2.8 σ tension for the 4th moment between us andthe results of [67] for their minimal selection of datasets.The tension is 2.4 σ for the standard selection, and below2 σ for the maximal selection and for the results of [61].The σ here is that for the phenomenological results sincethe lattice uncertainty is much smaller. Because the 4thmoment dominates the determination of a cµ , this tensionbetween lattice QCD+QED and some of the phenomeno-logical results carries over to a cµ , to be discussed in thenext section.The time-moments can also be used to determine avalue for m c by comparing to O ( α s ) continuum QCDperturbation theory and this was the focus of [61, 67].We do not do this here because the scale of α s is ratherlow in these determinations meaning that uncertaintiesfrom missing higher-order corrections can be substantial.We prefer instead the method of [25], which enables ahigher scale to be used in the perturbation theory. Wehave checked, however that the m c value that wouldbe obtained from the time-moments is consistent withboth [25] and the value given in Section IV. C. a cµ : Pure QCD and QCD+QED results To calculate the quark-line connected HVP contribu-tion to a µ from c quarks, a cµ , we can either use the physi-cal results for the vector current-current correlator time-moments discussed in the previous subsection or we cancalculate a cµ on each lattice ensemble and perform a fitas a function of lattice spacing to extrapolate directly tothe continuum limit. We will do the latter here.The values of a cµ on each lattice are given in Table XIV.4 . . . . . am c ) . . . . a c µ × − moments FIG. 26. Extrapolation to the continuum physical point of theconnected charm HVP contribution to the anomalous mag-netic moment of the muon. Different symbols denote resultson groups of ensembles with similar lattice spacing. Results atdeliberately mistuned c quark masses are not plotted but areincluded in the fit. The red points correspond to pure QCD,the light blue points to QCD+QED and the dashed green fitcurve plotted is that for QCD+QED. The continuum result(red cross) is compared to the result (open black square) ob-tained by calculating a cµ from the individually extrapolatedtime-moments in Section VI A. These are determined from the time-moments rescaled bythe
J/ψ mass on each lattice : M J/ψ G / ( n − n M expt J/ψ . (39)As discussed in Section VI A rescaling by M J/ψ reducesthe lattice spacing uncertainty and the impact of mis-tuning the c quark mass. This was used for lighter quarkmasses in [69] (see also [70]). The reduced effect of mis-tuning is clear from comparing the mistuned results inTable XIV to those in Table XI.Table XIV also gives the direct effect of quenched QEDthrough the ratio R [ a cµ ]. Because the rescaled mo-ments have less sensitivity to the c quark mass these num-bers are larger than 1 (unlike the results in Table XI) andreflect more closely the final impact of QED on a cµ . Weobserve no finite-volume dependence for the quenchedQED corrections to a cµ , as for the correlator time mo-ments. This is shown in Fig. 23.The results from Table XIV are shown in Fig. 26 alongwith our standard fit of the form given in Eq. (5). The fithas a χ / dof of 0.44. This fit obtains the physical values a cµ = 14 . × − , QCD (40) a cµ = 14 . × − , QCD + QEDalong with R QED [ a cµ ] = 1 . σ ) with those obtained using the extrapo-lated values for the moments from Table XII and calcu-lating a cµ in the continuum, as is seen in Fig. 26.The error budget for our final value of a cµ is given inTable XV. The largest uncertainties come from the deter- TABLE XV. Error budget for a cµ from our fit to a cµ values oneach ensemble. a cµ a → Z V w /a w M exp .J/ψ . . . . . a cµ × HPQCD HISQETMC twisted mass 2017HPQCD HISQ 2014BMWc stout stagg. 2017u,d,s,c seau,d,s sea
FIG. 27. Comparison of lattice QCD results (not includingQED) for the connected c quark HVP contribution to a µ , a cµ . Results are divided according to the number of sea quarkflavours included in the gluon field configurations on whichthe calculation was done. The first result, labelled ‘HPQCDHISQ 2014’ is from [25] using values of time-moments deter-mined in [29] using HISQ valence quarks on gluon field con-figurations including n f = 2 + 1 flavours of asqtad sea quarks.The other results all include n f = 2 + 1 + 1 flavours of seaquarks. The result labelled ‘ETMC twisted mass 2017’ usesthe twisted mass formalism [71] and that labelled ‘BMW stoutstagg. 2017’ a smeared staggered quark action [72]. Our newresult (from Eq. (40)) labelled ‘HPQCD HISQ’ agrees with,but is much more accurate than, these earlier results. mination of the lattice spacing, although this uncertaintyis much reduced by our rescaling with M J/ψ . D. Discussion: a cµ Our new result for the pure QCD case (Eq. (40)) iscompared to earlier lattice QCD results with realistic seaquark content in Fig. 27. Our new result (the top pointon the plot) is much more accurate than the earlier re-sults. With respect to the first calculation of a cµ that alsoused HISQ quarks [25] we have greatly reduced the pre-vious dominant sources of uncertainty from Z V and thedetermination of the lattice spacing.5 .
000 0 .
005 0 .
010 0 .
015 0 . a [fm ]0 . . . . . a c µ × − FIG. 28. Comparison of the discretisation effects in a cµ in ourresults from Table XIV (red open symbols) and those fromthe BMW collaboration in [72] (filled blue circles) that use aless highly-improved staggered quark action. The filled bluecircles only give estimates of the position of the BMW datapoints and do not indicate the statistical uncertainties whichare much smaller than the size of the points. With respect to results using other formalisms, we giveone figure to demonstrate the control of discretisation ef-fects that is possible with the HISQ formalism. Figure 28compares the approach to the continuum limit of our re-sults (from Table XIV) with results from BMWc [72],for which the continuum extrapolated result is shown inFig. 27. The points plotted from [72] are estimates of thepositions read from Figure S4, and do not include anyindication of statistical uncertainties. The BMWc stoutstaggered quark action has O ( a ) discretisation errors attree-level since it uses an unimproved derivative in itsversion of the Dirac equation. Figure 28 shows that theprice to be paid for not improving the discretisation is avery large discretisation effect. This is particularly evi-dent when working with heavier quarks such as charm,since it means that the dominant (Λ a ) effects have Λ ≈ ∼ ∼ a ∼ a cµ is +0.214(19)% (Eq. 40). This is a shift in a cµ inthe presence of the dominant QED effect of δa cµ = +0 . × − . (41)We can compare this to the result obtained by ETMCin [5] following work on the QED effect on the renormal-isation of quark bilinears in [23]. The ETMC value is+0 . × − , 2.8 σ smaller than ours.The two calculations of the quenched QED effect aredifferent. Ours is a direct calculation of the quenchedQED effect, retuning the valence quark mass throughdetermination of a meson mass in the usual way. The . . . . . a cµ × Dehnadi et al 2011HPQCD HISQK¨uhn et al 2007u,d,s,c sea + QEDpheno.
FIG. 29. Comparison of our lattice QCD+QED result (bluecross) for the connected c quark HVP contribution to a µ , a cµ , with results determined from experimental information.The red open squares, denoted ‘pheno.’, use the charmoniummoments from R e + e − given in [61, 67] to determine a cµ . Theorange open circles give the alternative maximal (upper) andminimal (lower) inclusive cross-section datasets from [67]. ETMC calculation is perturbative in quenched QED andfixes the valence quark masses so that they agree in theMS scheme at 2 GeV in QCD+QED and pure QCD. Ourresults in Section IV show that the quark mass in the MSscheme is lower in QCD+QED than in QCD (at scalesabove m c ) when the quark mass is tuned so that M J/ψ agrees with experiment in the two cases. This leads us toexpect a larger result for the impact of quenched QED on a cµ with our tuning. Once full lattice QCD+QED calcu-lations are underway tuning of the quark masses will bedone through matching of meson masses to experiment.Figure 29 includes a comparison of our result for a cµ ,now in QCD+QED, with values obtained using the mo-ments determined from J/ψ and ψ (cid:48) properties and theinclusive cross-section for e + e − → hadrons, removingcontributions from quarks other than c . The results forthese moments from [61, 67] were discussed and com-pared to our results in Section VI B. In Fig. 29 we haveconverted the moments into a result for a cµ for compar-ison. As with the moments, we see a 2 . σ tension withthe [67] result using the standard selection of datasets( a cµ = 14 . × − ), and a larger tension with theminimal selection of datasets. There is agreement within2 σ for the maximal selection of datasets from [67] andwith [61]. Our result may then provide a pointer to theselection of R e + e − datasets for the phenomenological de-termination and/or indicate an issue with the perturba-tion theory used to subtract the u/d/s contribution toobtain R c ( s ).A more recent determination of the complete HVPcontribution to a µ using results from R e + e − is givenin [51]. Although the c component is not separatedout, the contribution from the J/ψ resonance is given6as 6 . × − . We can also readily determine thecontribution to a cµ from the J/ψ resonance alone since inthat case [73] G n = n ! f J/ψ M nJ/ψ . (42)Using our value for f J/ψ from Eq. (29) and the experi-mental value for M J/ψ [1] gives a cµ ( J/ψ ) = 6 . × − . (43)This is in good agreement with the result determinedfrom the experimental J/ψ parameters, but more accu-rate, reflecting the fact that our result for Γ e + e − in Fig. 22agrees with experiment but has smaller uncertainty.The J/ψ contribution provides almost half of a cµ - weconclude that it is the rest of the contribution, from theinclusive cross-section above the resonance region, thatcauses the tension between our results for a cµ and thatfrom R e + e − for some selections of experimental datasets.The tension then amounts to 7(3)% of this non-resonantcross-section. VII. CONCLUSIONS
We have performed the first n f = 2+1+1 lattice QCDcomputations of the properties of ground-state charmo-nium mesons. These have been done using the HISQaction to calculate quark-line connected two-point corre-lation functions on gluon field configurations that include u/d quark masses going down to the physical point. Thesmall discretisation effects in the HISQ action and highstatistics achievable have given us good control over boththe continuum and chiral extrapolations and has enabledus to obtain smaller uncertainties than previous latticeQCD computations of these properties (including previ-ous calculations by HPQCD). At the same time we haveimproved the tuning of the bare c quark mass to updatethe value of m c . From the same correlators that we usefor the masses and decay constants we have also derivedan improved result for the c quark HVP contribution tothe anomalous magnetic moment of the muon.The precision possible for c quark correlators with theHISQ action makes it possible to determine the impactof the c quark’s electric charge. We do this directly andnonperturbatively in quenched QED by multiplying anappropriate U(1) field into our gluon field configurations.We tune the bare c quark mass so that the J/ψ massagrees with experiment in both QCD+QED and QCDand we calculate mass and vector renormalisation con-stants in the RI-SMOM scheme in both cases, perform-ing a full analysis as a function of µ to determine the c quark mass in the MS scheme.Here we collect our final QCD+QED results (fromEqs. (10), (17), (29), (30), (40)) before discussing each in turn: M J/ψ − M η c = 0 . m c (3 GeV) = 0 . f J/ψ = 0 . f η c = 0 . a cµ = 14 . × − . (44)Error budgets are given in Tables IV, VI, X and XVrespectively.The precision of our result for the charmonium hy-perfine splitting allows us to resolve, for the first time,the sign and magnitude of the anticipated difference be-tween the lattice and experimental results arising fromthe fact that we do not include quark-line disconnectedcorrelation functions. We take this difference to be theeffect of the η c decay to two gluons which is prohibited inthe lattice calculation, and conclude that ∆ M annihln η c =+7 . .
2) MeV.The effect of QED on the hyperfine splitting is fairlysubstantial (1.4%) and the largest effect that we observehere. This has 3 components that all act in the samedirection: a direct effect of 0.7%, 0.1% from retuning the c quark mass and 0.6% from J/ψ annihilation to a photonthat we add by hand.Our updated value of the charm quark mass (in the MSscheme at a scale of 3 GeV) includes the effect of QEDon the bare c quark mass, tuned so that M J/ψ matchesexperiment, and on the mass renormalisation constant Z m determined on the lattice using the intermediate RI-SMOM scheme. The addition of results at a finer lat-tice spacing has improved the uncertainty slightly overHPQCD’s earlier result [44].The impact of QED on m c is small since QED effectstend to cancel between the retuning needed and changesin Z m . At a scale of 3 GeV we see a -0.18(2)% effect,falling towards zero at a scale of m c .Our result for the the J/ψ decay constant is the mostprecise to date and acts as a subpercent test of QCD.The gain in precision from the 2012 HPQCD calculationof [29] is a result of the use of a more accurate renormali-sation of the vector current [16] as well as gluon field con-figurations with a wider range of sea u/d quark massesand lattice spacing values. We can use our result for f J/ψ to determine a value for the width for
J/ψ decay toa lepton-antilepton pair (repeating Eq. (32)):Γ(
J/ψ → e + e − ) = 5 . . (45)This is more accurate than the current average of exper-imental results [1].The impact of quenched QED on f J/ψ is +0.2%, sincethe retuning of the c quark mass offsets some of the di-rect effect. The effect on f η c is almost the same so thatthe ratio of f J/ψ to f η c remains the same. This ratio isdetermined here to be 1.0289(19), so that it is definitelygreater than 1; this was not completely clear from earliercalculations.7Our results for the time-moments of the charmoniumvector current-current correlators also provide a new levelof accuracy for these quantities, improving by a factor of10 over the first such calculations in [29]. Our resultsare given in Table XII. We see some tension for the low-est (4th) moment with phenomenological results derivedfrom R ( e + e − → hadrons) in [67] when particular selec-tions of experimental datasets are made.We use the time-moments to derive the connected c quark HVP contribution to the anomalous magnetic mo-ment of the muon, a cµ . Although this is not a large partof the total HVP contribution and so improving its un-certainty has little impact on the full HVP contribution,nevertheless it is a piece that can be calculated very ac-curately in lattice QCD and provides a test case for com-parison of lattice calculations and a comparison with phe-nomenology.Our result for a cµ improves the accuracy by a factorof 3 over earlier lattice QCD values [25, 71, 72]. Com-parison of our result for a cµ to that determined fromphenomenology can be divided into contributions fromnarrow resonances and from the continuum e + e − → cc .The contribution from the J/ψ agrees well between ourresult and phenomenology, with our result being moreaccurate (see Section VI D). This reflects the situationdescribed above for f J/ψ and Γ (cid:96)(cid:96) . The total a cµ derivedfrom the time-moments determined from experimentaldata on R ( e + e − → hadrons) when the component froma cc loop is separated out shows some tension with our re-sults, depending on which experimental datasets are usedabove the resonance region. Our central value is higher,tending to reduce by a small amount the discrepancy in a µ between existing experiment and the Standard Model.We should stress, however, that more complete determi-nations of a µ , for example in [51], do not separate out a cµ and so we cannot make a direct comparison to them.We also determine the impact of quenched QED on a cµ in a scheme in which the c quark mass is tuned from M J/ψ in both QCD+QED and pure QCD. We find (repeatingEq. (41)) δa cµ = +0 . × − . (46)This is a 0.2% effect and dominated by the impact ofretuning the c quark mass.Our result for a cµ has an accuracy of 0.3%. Sub-0.5%uncertainty is the aim for lattice QCD calculations of thefull HVP contribution to a µ . We have shown that this ispossible, for a small piece of the HVP, at least. Acknowledgements
We are grateful to the MILC collaboration for the useof their configurations. We are also grateful for the useof MILC’s QCD code. We have modified it to gener-ate quenched U(1) gauge fields and incorporate thoseinto the quark propagator calculation as described here.We thank P. Knecht for contributions at the start of the .
000 0 .
005 0 .
010 0 .
015 0 . a [fm ]0 . . . Z l o c V ( S M O M ) FIG. 30. Results for the local vector current renormalisationfactor in the RI-SMOM scheme compared to the fit form givenin Eq. A1. All of the data points shown are included in thefit but only the 2 and 3 GeV fit lines are drawn. The datapoints not on the fit lines correspond to µ = 2.5 GeV (lightorange) and 4 GeV (light green). The dashed line shows thevalue of a on the exafine lattices that allows us to determinea Z V value there. project, C. McNeile for the calculation of w /a and A.Keshavarzi, D. Miller, D. Nomura, D. Smaranda and T.Teubner for useful discussions. Computing was done onthe Darwin supercomputer at the University of Cam-bridge High Performance Computing Service as part ofthe DiRAC facility, jointly funded by the Science andTechnology Facilities Council, the Large Facilities Capi-tal Fund of BIS and the Universities of Cambridge andGlasgow. 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: Z V and Z m determination on the finestlattices We can use the results of [16] to check the consistencyof the Z V value on set 14 that we present here and toobtain a value on set 15. We do this by fitting the RI-SMOM results for the local vector current in Table III of[16] to the expected functional form. This form is basedon a power series in α s evaluated in the MS scheme ata scale of 1 /a , which is the perturbative lattice-to-MSrenormalisation factor, remembering that the SMOM toMS renormalisation factor is exactly 1 in this case. Inaddition we must allow for possible discretisation effectsthat depend on aµ . The form we use is: Z loc V (SMOM)( a, µ ) = 1 + i =4 ,j =3 (cid:88) i,j =1 (cid:20) c i + d ij (cid:16) aµπ (cid:17) j (cid:21) α is . (A1)This is very similar to the approach adopted in AppendixB of [74] using results for Z V from the determination ofform factors between two identical mesons at rest. As8 TABLE XVI. Values of Z SMOM m obtained in the RI-SMOMscheme on set 14 in pure QCD. Results are given at two valuesof µ , listed in lattice units in column 1. They correspondapproximately to 2 GeV and 3 GeV. The correlation betweenthe two numbers is 0.122. aµ Z SMOM m there, we fix the α s coefficient to its known perturbativevalue of -0.1164(3).Fig. 30 shows the Z V data from [16] as hexagons,coloured according to their µ value. The fit lines for2 and 3 GeV are also shown. The fit has a χ / dof of0.93. We note that the results from set 14 are includedin this fit. The results from set 14, however, also agreewith the fit result when they are not included in the fit.As discussed in Section IV µ was slightly mistuned onset 14; the true values are 2.04 and 2.98 GeV rather than2 and 3 GeV. At this small lattice spacing the variationin Z V with µ (which is a discretisation effect) is small enough that this small mistuning can be neglected.Note that the underlying perturbation theory ofEq. (A1) should agree with that obtained from the fitin [74] and it does. The discretisation effects are differ-ent in the two cases, of course. For RI-SMOM there is amomentum scale µ which we take as 2 GeV (although itcan be taken as much smaller for Z V [16]) and this setsthe size of discretisation effects as is clear from Figure 30.Using the fit of Eq. (A1) we may also extract Z V val-ues at finer lattice spacings where it may not be practicalfor us to perform direct calculations due to the compu-tational cost of Landau gauge fixing. This includes theexafine lattices of set 15 in Table I with a lattice spacingadjusted for sea-mass mistuning of 0.032 fm. The valueof Z V from the fit for these lattices is 0.99296(21) at µ = 2 GeV and 0.99186(18) at 3 GeV.Results for the mass renormalisation factor from thelattice to the RI-SMOM scheme, Z SMOM m ( µ ), determinedon ultrafine set 14 are given in Table XVI at µ = 2 GeVand 3 GeV. These are calculated in the same way as thatdiscussed in [44]. [1] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D98 , 030001 (2018).[2] J. Grange et al. (Muon g-2), (2015), arXiv:1501.06858[physics.ins-det].[3] G. W. Bennett et al. (Muon g-2), Phys. Rev. Lett. ,161802 (2004), arXiv:hep-ex/0401008 [hep-ex].[4] P. Boyle, V. G¨ulpers, J. Harrison, A. J¨uttner, C. Lehner,A. Portelli, and C. T. Sachrajda, JHEP , 153 (2017),arXiv:1706.05293 [hep-lat].[5] D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo,and S. Simula, Phys. Rev. D99 , 114502 (2019),arXiv:1901.10462 [hep-lat].[6] S. Borsanyi et al. , (2020), arXiv:2002.12347 [hep-lat].[7] G. K. Cheung, C. O’Hara, G. Moir, M. Peardon, S. M.Ryan, C. E. Thomas, and D. Tims (Hadron Spectrum),JHEP , 089 (2016), arXiv:1610.01073 [hep-lat].[8] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P.Lepage, J. Shigemitsu, H. Trottier, and K. Wong(HPQCD, UKQCD), Phys. Rev. D75 , 054502 (2007),arXiv:hep-lat/0610092 [hep-lat].[9] A. Bazavov et al. (MILC), Phys. Rev.
D87 , 054505(2013), arXiv:1212.4768 [hep-lat].[10] A. Bazavov et al. , Phys. Rev.
D98 , 074512 (2018),arXiv:1712.09262 [hep-lat].[11] C. Bernard and D. Toussaint (MILC), Phys. Rev. D ,074502 (2018), arXiv:1707.05430 [hep-lat].[12] A. Hart, G. M. von Hippel, and R. R. Horgan (HPQCD),Phys. Rev. D79 , 074008 (2009), arXiv:0812.0503 [hep-lat].[13] S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, et al. , JHEP , 010 (2012), arXiv:1203.4469 [hep-lat].[14] R. Dowdall, C. Davies, G. Lepage, and C. Mc-Neile (HPQCD), Phys.Rev.
D88 , 074504 (2013),arXiv:1303.1670 [hep-lat]. [15] C. Monahan, J. Shigemitsu, and R. Horgan, Phys. Rev.
D87 , 034017 (2013), arXiv:1211.6966 [hep-lat].[16] D. Hatton, C. T. H. Davies, G. P. Lepage, andA. T. Lytle (HPQCD), Phys. Rev.
D100 , 114513 (2019),arXiv:1909.00756 [hep-lat].[17] G. P. Lepage, B. Clark, C. T. H. Davies, K. Hornbostel,P. B. Mackenzie, C. Morningstar, and H. Trottier,
Lat-tice field theory. Proceedings, 19th International Sym-posium, Lattice 2001, Berlin, Germany, August 19-24,2001 , Nucl. Phys. Proc. Suppl. , 12 (2002), arXiv:hep-lat/0110175 [hep-lat].[18] R. J. Dowdall, C. T. H. Davies, R. R. Horgan, G. P.Lepage, C. J. Monahan, J. Shigemitsu, and M. Wingate,Phys. Rev.
D100 , 094508 (2019), arXiv:1907.01025 [hep-lat].[19] M. Hayakawa and S. Uno, Prog. Theor. Phys. , 413(2008), arXiv:0804.2044 [hep-ph].[20] P. Weisz, Nucl. Phys. B , 1 (1983).[21] A. Hart, R. Horgan, and L. Storoni, Phys. Rev. D ,034501 (2004), arXiv:hep-lat/0402033.[22] J. Rosner, S. Stone, and R. Van de Water (Particle DataGroup review), in Phys. Rev. D98 , 030001 (2018).[23] M. Di Carlo, D. Giusti, V. Lubicz, G. Martinelli, C. T.Sachrajda, F. Sanfilippo, S. Simula, and N. Tantalo,Phys. Rev.
D100 , 034514 (2019), arXiv:1904.08731 [hep-lat].[24] G. Amoros, J. Bijnens, and P. Talavera, Nucl. Phys.
B602 , 87 (2001), arXiv:hep-ph/0101127 [hep-ph].[25] B. Chakraborty, C. T. H. Davies, B. Galloway, P. Knecht,J. Koponen, G. C. Donald, R. J. Dowdall, G. P. Lep-age, and C. McNeile, Phys. Rev.
D91 , 054508 (2015),arXiv:1408.4169 [hep-lat].[26] Z. Davoudi and M. J. Savage, Phys. Rev.
D90 , 054503(2014), arXiv:1402.6741 [hep-lat]. [27] A. Gray, I. Allison, C. T. H. Davies, E. Dalgic, G. P. Lep-age, J. Shigemitsu, and M. Wingate, Phys. Rev. D72 ,094507 (2005), arXiv:hep-lat/0507013 [hep-lat].[28] C. M. Bouchard, G. P. Lepage, C. Monahan, H. Na,and J. Shigemitsu, Phys. Rev.
D90 , 054506 (2014),arXiv:1406.2279 [hep-lat].[29] G. C. Donald, C. T. H. Davies, R. J. Dowdall, E. Follana,K. Hornbostel, J. Koponen, G. P. Lepage, and C. Mc-Neile, Phys. Rev.
D86 , 094501 (2012), arXiv:1208.2855[hep-lat].[30] Q. N. Xu et al. (Belle), Phys. Rev.
D98 , 072001 (2018),arXiv:1805.03044 [hep-ex].[31] R. Aaij et al. (LHCb), Eur. Phys. J.
C75 , 311 (2015),arXiv:1409.3612 [hep-ex].[32] R. Aaij et al. (LHCb), Phys. Lett.
B769 , 305 (2017),arXiv:1607.06446 [hep-ex].[33] R. Aaij et al. (LHCb), Eur. Phys. J.
C77 , 609 (2017),arXiv:1706.07013 [hep-ex].[34] V. V. Anashin et al. , Phys. Lett.
B738 , 391 (2014),arXiv:1406.7644 [hep-ex].[35] J. P. Lees et al. (BaBar), Phys. Rev.
D89 , 112004 (2014),arXiv:1403.7051 [hep-ex].[36] M. Ablikim et al. (BESIII), Phys. Rev. Lett. , 222002(2012), arXiv:1111.0398 [hep-ex].[37] M. Ablikim et al. (BESIII), Phys. Rev.
D86 , 092009(2012), arXiv:1209.4963 [hep-ex].[38] C. DeTar, A. S. Kronfeld, S.-h. Lee, D. Mohler, andJ. N. Simone (Fermilab Lattice, MILC), Phys. Rev.
D99 ,034509 (2019), arXiv:1810.09983 [hep-lat].[39] R. A. Briceno, H.-W. Lin, and D. R. Bolton, Phys. Rev.
D86 , 094504 (2012), arXiv:1207.3536 [hep-lat].[40] Y.-B. Yang et al. , Phys. Rev.
D92 , 034517 (2015),arXiv:1410.3343 [hep-lat].[41] V. V. Anashin et al. , Phys. Lett.
B749 , 50 (2015).[42] G. T. Bodwin, E. Braaten, and G. P. Lepage,Phys. Rev.
D51 , 1125 (1995), [Erratum: Phys.Rev.D55,5853(1997)], arXiv:hep-ph/9407339 [hep-ph].[43] L. Levkova and C. DeTar, Phys. Rev.
D83 , 074504(2011), arXiv:1012.1837 [hep-lat].[44] A. T. Lytle, C. T. H. Davies, D. Hatton, G. P. Lep-age, and C. Sturm (HPQCD), Phys. Rev.
D98 , 014513(2018), arXiv:1805.06225 [hep-lat].[45] C. Sturm, Y. Aoki, N. H. Christ, T. Izubuchi, C. T. C.Sachrajda, and A. Soni, Phys. Rev.
D80 , 014501 (2009),arXiv:0901.2599 [hep-ph].[46] B. A. Kniehl and O. L. Veretin, (2020), arXiv:2002.10894[hep-ph].[47] A. Bednyakov and A. Pikelner, (2020), arXiv:2002.12758[hep-ph].[48] A. Bednyakov, B. Kniehl, A. Pikelner, and O. Veretin,Nucl.Phys.B , 463 (2017), arXiv:1612.00660 [hep-ph].[49] A. Bazavov et al. (Fermilab Lattice, MILC, TUMQCD),Phys. Rev.
D98 , 054517 (2018), arXiv:1802.04248 [hep-lat].[50] N. Carrasco et al. (ETM), Nucl. Phys.
B887 , 19 (2014),arXiv:1403.4504 [hep-lat].[51] A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev.D , 114025 (2018), arXiv:1802.02995 [hep-ph]. [52] V. V. Anashin et al. (KEDR), JHEP , 119 (2018),arXiv:1801.01958 [hep-ex].[53] M. Ablikim et al. (BESIII), Phys. Lett. B761 , 98 (2016),arXiv:1604.01924 [hep-ex].[54] V. V. Anashin et al. , Phys. Lett.
B711 , 280 (2012),arXiv:1109.4215 [hep-ex].[55] J. P. Alexander, G. Bonvicini, P. S. Drell, R. Frey, andV. Luth, Nucl. Phys.
B320 , 45 (1989).[56] V. S. Fadin, V. A. Khoze, and A. D. Martin, Phys. Lett.
B320 , 141 (1994), arXiv:hep-ph/9309234 [hep-ph].[57] V. V. Anashin et al. , Proceedings, 6th Interna-tional Workshop on e+e- Collisions from Phi to Psi(PHIPSI09): Beijing, China, October 13-16, 2009 , Chin.Phys.
C34 , 836 (2010), arXiv:1110.0328 [hep-ex].[58] A. Keshavarzi, D. Nomura, and T. Teubner, private com-munication (2020).[59] B. Colquhoun, C. T. H. Davies, R. J. Dowdall, J. Kettle,J. Koponen, G. P. Lepage, and A. T. Lytle (HPQCD),Phys. Rev.
D91 , 114509 (2015), arXiv:1503.05762 [hep-lat].[60] C. T. H. Davies, C. McNeile, E. Follana, G. P. Lep-age, H. Na, and J. Shigemitsu, Phys. Rev.
D82 , 114504(2010), arXiv:1008.4018 [hep-lat].[61] J. H. Kuhn, M. Steinhauser, and C. Sturm, Nucl. Phys.
B778 , 192 (2007), arXiv:hep-ph/0702103 [HEP-PH].[62] D. Beˇcirevi´c, G. Duplanˇci´c, B. Klajn, B. Meli´c,and F. Sanfilippo, Nucl. Phys.
B883 , 306 (2014),arXiv:1312.2858 [hep-ph].[63] G. Bailas, B. Blossier, and V. Mor´enas, Eur. Phys. J.
C78 , 1018 (2018), arXiv:1803.09673 [hep-lat].[64] B. Chakraborty, C. T. H. Davies, G. C. Donald,R. J. Dowdall, J. Koponen, G. P. Lepage, andT. Teubner (HPQCD), Phys. Rev.
D89 , 114501 (2014),arXiv:1403.1778 [hep-lat].[65] I. Allison et al. (HPQCD), Phys. Rev.
D78 , 054513(2008), arXiv:0805.2999 [hep-lat].[66] T. Blum, Phys. Rev. Lett. , 052001 (2003), arXiv:hep-lat/0212018 [hep-lat].[67] B. Dehnadi, A. H. Hoang, V. Mateu, and S. M. Zebarjad,JHEP , 103 (2013), arXiv:1102.2264 [hep-ph].[68] K. Nakayama, B. Fahy, and S. Hashimoto, Phys. Rev. D94 , 054507 (2016), arXiv:1606.01002 [hep-lat].[69] B. Chakraborty, C. T. H. Davies, P. G. de Oliviera, J. Ko-ponen, G. P. Lepage, and R. S. Van de Water, Phys. Rev.