Charmonium properties from lattice QCD + QED: hyperfine splitting, $J/\psi$ leptonic width, charm quark mass and $a_{\mu}^c$

We have performed the first $n_f = 2+1+1$ lattice QCD computations of the properties (masses and decay constants) of ground-state charmonium mesons. Our calculation uses the HISQ action to generate quark-line connected two-point correlation functions on MILC gluon field configurations that include $u/d$ quark masses going down to the physical point, tuning the $c$ quark mass from $M_{J/\psi}$ and including the effect of the $c$ quark's electric charge through quenched QED. We obtain $M_{J/\psi}-M_{\eta_c}$ (connected) = 120.3(1.1) MeV and interpret the difference with experiment as the impact on $M_{\eta_c}$ of its decay to gluons, missing from the lattice calculation. This allows us to determine $\Delta M_{\eta_c}^{\mathrm{annihiln}}$ =+7.3(1.2) MeV, giving its value for the first time. Our result of $f_{J/\psi}=$ 0.4104(17) GeV, gives $\Gamma(J/\psi \rightarrow e^+e^-)$=5.637(49) keV, in agreement with, but now more accurate than experiment. At the same time we have improved the determination of the $c$ quark mass, including the impact of quenched QED to give $\overline{m}_c(3\,\mathrm{GeV})$ = 0.9841(51) GeV. We have also used the time-moments of the vector charmonium current-current correlators to improve the lattice QCD result for the $c$ quark HVP contribution to the anomalous magnetic moment of the muon. We obtain $a_{\mu}^c = 14.638(47) \times 10^{-10}$, which is 2.5$\sigma$ higher than the value derived using moments extracted from some sets of experimental data on $R(e^+e^- \rightarrow \mathrm{hadrons})$. This value for $a_{\mu}^c$ includes our determination of the effect of QED on this quantity, $\delta a_{\mu}^c = 0.0313(28) \times 10^{-10}$.


I. INTRODUCTION
The precision of lattice QCD calculations has been steadily improving for some time and is now approaching, or has surpassed, the 1% level for multiple quantities. Good examples are the masses and decay constants of ground-state pseudoscalar mesons [1]. The meson masses can be used to tune, and therefore determine, quark masses. The decay constants can be combined with experimental annihilation rates to leptons to determine elements of the Cabibbo-Kobayashi-Maskawa matrix. The accuracy of modern lattice QCD results means that sources of small systematic uncertainty that could appear at the percent level need to be understood. At this level QED effects, i.e. the fact that quarks carry electric as well as color charge, come into play. A naive argument that such effects could be O(α QED ) would imply a possible 1% contribution. One key driver for the lattice QCD effort to include QED effects has been that of calculations of the hadronic vacuum polarisation (HVP) contribution to the anomalous magnetic moment of the muon, a µ . New results are expected soon from the Muon g − 2 experiment at Fermilab [2] which aims to clarify the observed tension between experiment and Standard Model theory seen by the Brookhaven E821 experiment [3]. Current lattice QCD calculations have reached the precision of a few percent for a µ and the systematic uncertainties, for example from neglecting QED effects, have become a major focus of attention (see, for example [4][5][6]).
QED effects can have large finite-volume corrections within a lattice QCD calculation because the Coulomb interaction is long-range. For electrically neutral correlation functions, such as the ones needed for a calculation of the HVP, and that we study here, this is much less of an issue (and we demonstrate this). Since α n s α QED is not very different in size from α QED at hadronic scales, calculations must be fully nonperturbative in QCD. A consistent calculation must also allow for the retuning of quark masses needed when QED effects are included.
Here we examine the properties of ground-state charmonium mesons more accurately than has been possible in previous lattice QCD calculations 1 . Since it is possible to obtain statistically very precise results for the charmonium system it is a good place to study small systematic effects from QED and other sources. We are able to see such effects in our results and quantify them. We include u, d, s and c quarks in the sea for the first time and have results on gluon field configurations with physical u/d sea quarks. We also analyse the impact of including an electric charge on the valence c quarks (only). This approximation, known as 'quenched QED', should capture the largest QED effects and enable us to see how much of a difference QED makes. For the vector (J/ψ) meson properties we improve on earlier results by using an exact method to renormalise the lattice vector current (both in QCD and QCD+QED). To tune the c quark mass we use the J/ψ meson mass, taking into account the retuning that is required when QED is switched on. We find the impact of this to be of similar size to the more direct QED effects.
The quantities that we focus on here are the masses and 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 lattice QCD and QCD+QED calculations are 'connected' correlation functions i.e. they are constructed from combining charm quark propagators from the source to the sink. 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 annihilation channel has on the meson masses by looking at the meson widths, which are twice the imaginary part of the mass and are dominated by hadronic channels. The J/ψ has a tiny width of 93 keV [1] but the pseudoscalar η c can annihilate to two gluons, allowing it to mix with other flavour-singlet pseudoscalars, and it has a width of 32 MeV [1]. The annihilation channel might then be expected to have a larger impact on the η c mass and lead lattice QCD calculations of the mass from connected correlators to disagree with experiment. The only way to achieve the O(1 MeV) accuracy required to see this is to determine the mass difference between the J/ψ and η c (the hyperfine splitting). Any shift in the η c mass will have a much larger (by a factor of 30) relative effect in this splitting. Previous lattice QCD calculations of the hyperfine splitting from connected correlation functions have not been accurate enough to see a significant difference with experiment. Here, for the first time, we can see such a difference because we have very good control both of discretisation effects and sea-quark mass effects and can also determine the impact of QED on this quantity.
A further place in which QED effects need to be quantified, given our accuracy, is that of the determination of the c quark mass. We do this by tuning our results to the experimental J/ψ meson mass with and without electric charge on the valence c quarks and also determine the small change in the mass renormalisation factor, Z m , needed to convert to the standard MS mass.
Our analysis of the vector charmonium correlation functions with a completely nonperturbative renormalisation of the vector current allows much improved accuracy in the determination of the leptonic decay width of the J/ψ. Using the same correlators, we determine the charm quark portion of the hadronic vacuum polarisation contribution to the anomalous magnetic of the muon, a HVP,c µ , along with the impact of quenched QED on this quantity. We can compare this to phenomenological results from R(e + e − → hadrons).
The paper is laid out as follows: • Section II describes the lattice QCD calculation and the inclusion of quenched QED; • Section III describes our determination of the hyperfine splitting; • Section IV, the c quark mass; • Section V, the J/ψ and η c decay constants; • Section VI the time-moments of c vector-vector correlators and the c quark hadronic vacuum polarisation contribution to a µ ; • Section VII gives our conclusions.
Each section of results includes a description of the pure QCD calculation followed by a determination of the impact of quenched QED on the result and then a discussion subsection including comparison to both experiment and previous lattice QCD calculations, where applicable. Finally Section VII collects up all of our results and summarises our conclusions.

II. LATTICE SETUP
We perform calculations on a total of 15 gluon field ensembles that include the effects of light (up/down with the same mass), strange and charm quarks in the sea (n f = 2 + 1 + 1) all using the HISQ action [8] and generated by the MILC collaboration [9,10]. These include lattices at 6 different β (the bare QCD coupling) values corresponding to 6 different sets of lattice spacing values with the finest lattice reaching a spacing of ∼ 0.03 fm. Although 'topology freezing' has been seen on the very finest lattices used here, we do not expect this to have significant impact on the charmonium quantities we study here because no valence light quarks are involved in the calculations [11]. We use ensembles with sea u/d masses at the physical point on all but the finest two lattice spacings. We also employ three ensembles with shared parameters except for their spatial extent in units of the lattice spacing L s . These ensembles allow us to investigate finite volume effects in our QED analysis. The gluon action on these ensembles is improved so that discretisation errors through O(α s a 2 ) are removed [12]. Parameters for the ensembles are given in Table I.
On these gluon field configurations we calculate propagators for valence c quarks by solving the Dirac equation for a source consisting of a set of Gaussian random numbers across a timeslice (a random wall source). We use multiple time sources per configuration to improve statistical accuracy. The number of configurations used and the number of time sources is given in Table I. The table also gives the valence c quark masses in lattice units, which may differ from those in the sea because they are tuned 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 from the Wilson flow parameter, w0 [13], with w0/a values given in column 2. The lattice spacing can be determined in fm by using w0 = 0.1715 (9) fm [14] (fixed from fπ). Ls and Lt are the lattice spatial and temporal extents in lattice units. am val c is the valence c quark mass with Naik the corresponding Naik parameter (see text). The column headed N cfg × Nt shows both the number 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 referred to as very coarse (a ≈ 0.15 fm), sets 4-8 as coarse (a ≈ 0.12fm), sets 9-11 as fine (a ≈ 0.09 fm), 12 and 13 as superfine (a ≈ 0.06 fm), 14 as ultrafine (a ≈ 0.045 fm) and 15  below. The HISQ action [8] includes an improved discretisation of the covariant derivative in the Dirac equation. This removes tree-level a 2 discretisation errors by the addition of a 3-link 'Naik' term to the symmetric 1-link difference. For heavy quarks the coefficient of the Naik term is adjusted from 1 to 1+ Naik to remove (am) 4 errors at tree-level [8]. A closed-form expression for in terms of the tree-level quark pole mass is given in [15] along with the formula for the tree-level quark pole mass in terms of the bare mass. Table I gives the values of that we use. We combine charm quark and antiquark propagators to calculate two types of quark-line connected two-point correlation functions: pseudoscalar and vector. The ground state of the pseudoscalar correlation function corresponds to the η c meson and the vector correlation function, to the J/ψ. When using staggered quarks, as here, the different spin structures are implemented using position dependent phases in the operators at source and sink. The two-point 'Goldstone' pseudoscalar (γ 5 ⊗ γ 5 in spin-taste notation) correlation functions are simply constructed from quark propagators S(x, 0) from the origin to x as where the factor of 4 accounts for the taste multiplicity with staggered quarks. For the vector correlation functions we use a local vector operator (spin-taste γ i ⊗ γ i ). The correlation functions then combine S(x, 0) with a propagator made from patterning the source with a phase (−1) xi and inserting (−1) xi at the sink timeslice as we tie the propagators together and sum over spatial sites. Our vector correlation functions average over all spatial polarisations, i, for improved statistical precision. Note that we do not calculate any quark-line disconnected correlation functions.
The HISQ local vector current is not conserved and requires renormalisation. For this purpose we use the RI-SMOM momentum subtraction scheme implemented on the lattice as discussed in [16]. In [16] it was shown that, because of the Ward-Takahashi identity, these renormalisation factors do not suffer any contamination by nonperturbative artefacts (condensates) and can therefore be safely used in calculations such as those presented here. The quenched QED correction to the RI-SMOM vector current renormalisation was also given in [16] and shown to be tiny (∼ 0.01%) for the HISQ action (as expected since the pure QCD Z V values only differ from 1 at the 1% level and quenched QED provides a small correction to this difference from 1). Here we use the Z V values from [16] at a scale µ of 2 GeV. We will also demonstrate (see Section V) that using µ = 3 GeV gives the same result as it must for a Z V that correctly matches the lattice to continuum physics.
Since we make use of an ensemble (set 14) with a finer lattice spacing than those studied in [16] we have directly calculated the value of Z V on set 14 at µ = 2 GeV in addition. We have, however, only used a small number of configurations (6) in that calculation due to the computational limitation of the stringent Landau gauge fixing required. We therefore double the statistical uncertainty for Z V on that ensemble. This has very little impact on our final results as the Z V uncertainty is small. See Ap-pendix A for a discussion of our Z V values, where we also derive a Z V value for set 15. In order to tune the mass of the valence charm quark we use bare charm mass values on each lattice that produce a J/ψ mass equal to the experimental value (both in pure QCD and in QCD+QED). We choose the J/ψ here rather than the η c because the relatively large width of the η c means that annihilation effects that we are not including could lead to small (order 0.1% ) uncertainties in the mass. This is mentioned in Section I and will be discussed 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.0969 GeV (with negligible uncertainty) [1]. The two panels of Fig. 1, where the horizontal line is the experimental value, 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 described above as a function of the time separation, t, between source and sink. The aim is to extract the energies (masses) and amplitudes (giving decay constants) of the ground-state mesons in each channel. However it is important to allow for the systematic effect of excited states that are present in the correlation functions and can affect the ground-state values if they are not taken into account. We do this by fitting the correlators to sums of exponentials associated with each energy eigenvalue and using Bayesian priors to constrain the (ordered) excited states in the standard way [17]. The pseudoscalar correlators are fit to f (E, t) = e −Et + e −E(T −Lt) .
The vector correlators require a more complicated form because of the presence of opposite parity states as a result of the use of staggered quarks: We cut out the correlator values at low values of t, below some value t min (5 − 10) where excited state contamination is most pronounced. We also use a standard procedure (see Appendix D of [18]) to avoid underestimating the low eigenvalues of the correlation matrix and hence the uncertainty.
The fit parameters that we need from Eq. 2 and 3 are the mass of the ground-state (E P 0 and E V 0 ) and the amplitude (A P 0 and A V 0 ). From the amplitude we determine the decay constant, see Section V.  Table I using the valence masses in that Table. Two error bars are shown. The error can be broken into parts that are uncorrelated between different sets and the contribution from fixing the lattice spacing from the physical value of w0 which is correlated. The outer error bar shows the full uncertainty and the inner bar the uncertainty without the contribution from w0. Bottom panel: The J/ψ masses on sets 2, 6, 10 and 12 with and without the inclusion of quenched QED. On each set the same valence mass (and lattice spacing) was used for both pure QCD and QCD+QED but the points have been separated on the x-axis for clarity. The error bars are the same as for the top panel, but note that here there is a correlation of the uncertainty from w0/a for the QCD and QCD+QED results on each set. The QCD+QED results are above the pure QCD results in every case.

B. Including quenched QED
We perform calculations in both lattice QCD and in lattice QCD with quenched QED. By quenched QED we mean that we include effects from the valence quarks having electric charge but we neglect effects from the electric charge of the sea quarks. To include quenched QED effects we generate a random momentum space photon field A µ (k) in Feynman gauge for each QCD gluon field configuration. This choice of gauge simplifies the generation of the photon field as the QED path integral weight takes the form of a Gaussian with variance 1/k 2 wherê k = 2sin(k µ /2). The results presented here do not depend on this gauge choice. Once the momentum space field is generated zero modes are set to zero using the QED L formulation [19]. A µ is then Fourier transformed into position space. We have checked that these Feynman gauge A µ fields produce the plaquette and average link expected from O(α QED ) perturbation theory (readily obtained from O(α s ) calculations in lattice QCD with Wilson glue [20,21]). These gauge fields are exponentiated as exp(ieQA µ ) to give a U(1) field which is then multiplied into the QCD gauge links before HISQ smearing. Q is the quark electric charge in units of the proton charge e.
Quenched QED affects only the valence quarks; the sea quarks remain uncharged. Hence there is no coupling of QED effects to purely gluonic quantities. This means that the Wilson flow parameter w 0 /a, measured on each ensemble, is unchanged. The physical value of w 0 that is used to determine the lattice spacing on each ensemble was determined in [14] by matching the decay constant of the π meson, f π , in lattice QCD to that obtained from experiment. The experimental value of f π is obtained from measurement of the rate for π → ν[γ] decay where [γ] indicates that the rate is fully inclusive of additional photons. The rate obtained is then adjusted to remove electromagnetic and electroweak corrections and to give a 'purely leptonic rate' corresponding to weak annihilation at lowest order in the absence of QED [22]. Combining this with a determination of |V ud | from nuclear β decay [1] gives an experimental value of f π ≡ f expt π which is a 'pure QCD' value, albeit that for a physical π + meson. The dominant uncertainty in f expt π is that from the remaining uncertainty in the electromagnetic corrections to the experimental rate, mainly from the hadronic-structure dependent contributions to the emission of additional photons. This is set at 0.1% in [22].
Because f expt π is a pure QCD quantity it can be used to set the lattice spacing in lattice QCD in a way that should be minimally different for lattice QCD+QED 2 . Small differences might still be expected between the lattice QCD f π and the experimental value from the way that the quark masses are tuned in a pure QCD scenario. The lattice QCD calculation in [14] used m u = m d and tuned the average mass, m l , to the experimental mass of the π 0 , which is the mass that both neutral and charged π mesons have in the absence of QED, up to quadratic corrections in the u − d mass difference. An uncertainty was included in the π 0 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 below 0.1%. They are at the same level as potential effects from the sea and would therefore be only possible to pin down with a calculation that included the impact of having electrically charged quarks in the sea.
These expectations are backed up by recent lattice QCD+QED results [6] that used the Ω baryon mass to determine w 0 . The impact of QED for the sea quarks was included to first order in α QED . No effects linear in 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 0 M Ω was found to be O(0.01%), whereas the effect of QED for the valence quarks (already allowed for in the f π analysis) was O(0.05%). The final value of w 0 using M Ω from [6] agrees well with the result using f π from [14], although the uncertainties in both cases are completely dominated by those from the pure QCD, isospin-symmetric part of the calculation.
From this we conclude that the largest effect from adding QED to lattice QCD will come from the interaction between the electric charges of the valence quarks and to study this effect we can compare lattice QCD plus quenched QED with pure lattice QCD using the same value of the lattice spacing, determined from f π , in both calculations.
An estimate of the size of corrections from quenched QED (often simply referred to as QED in what follows) in charmonium systems can be obtained by studying the effect on the J/ψ mass. The bottom panel of Fig. 1 shows the J/ψ mass for the same valence c mass for both QCD+QED and pure QCD calculations on sets 2, 6, 10 and 12. The QCD+QED and QCD results at the same lattice spacing are separated on the x-axis for clarity. All points share a correlated uncertainty (the outer error bar) from w 0 and this dominates the uncertainty. The uncorrelated error is shown by the smaller inner error bar. Note that the points at the same lattice spacing are also correlated through their w 0 /a value. The shift of the mass in QCD+QED compared to pure QCD is very small, at the level of 0.1%, and is upwards.
When discussing QCD+QED and pure QCD calculations of some quantity X we will use the notation X[QCD + QED] and X[QCD] respectively. We will often consider the ratio of the two for which we will use the shortened notation R (0) QED [X]. R 0 QED 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 includes the impact of retuning the c quark mass to give the experimental J/ψ mass in both the QCD+QED and pure QCD cases. So As shown in Fig. 1 the bare c quark mass has to readjusted downwards for QCD+QED relative to pure QCD.

C. Fitting strategy
We have results in pure QCD for all of the sets in Table 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 mistuning, which may be similarly sized, we choose to fit all of this data in a single fit for each quantity we consider. The generic form of the fit we use for a quantity X is Here Q is the valence quark electric charge (in units of e) used in the calculation and is therefore 0 in pure QCD. The pure QCD value of this fit at the physical point (the continuum limit with quark masses set to their physical values) is x. The value including quenched QED corrections is x[1 + α QED Q 2 c QED ]. Note that the factor of α QED multiplying the QED part of the fit function is there 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 pieces missing from our calculation and fit are α QED α s terms from 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 discretisation effects. Because we are dealing with heavy quarks here, the scale of discretisation effects can be set by m c and will typically be larger than those for light-quark quantities. Since m c > Λ QCD any discretisation effects set by scale Λ QCD will simply appear as m c -scale discretisation effects with a small coefficient.
The terms on the second line allow for mistuning of the sea u/d and s masses and discretisation effects in that mistuning (we shall see that those are important for the hyperfine splitting). The total of the mistuning of the sea masses is defined as in [25]: m phys s is taken from [25] or, where not available on the finest lattices, calculated from the tuned c quark mass and the m c /m s ratio given in [25]. The value of Λ in the discretisation effects multiplying the sea-quark mistuning is taken as 1 GeV (∼ m c ).
The effect of mistuning the charm mass in the sea is included in the third line of Eq.
The values of m phys c are taken from [25]. Although this used a slightly different tuning method the differences are negligible for this purpose. We have tested that including discretisation effects for this term has no effect on the fit.
Mistuning in the valence mass is accounted for through δ val,c m on the third line. We define this as The experimental value of the J/ψ mass is 3.0969 GeV [1] with negligible uncertainty. In order to make use of the correlations between our QCD+QED and pure QCD results on the same gluon field configurations we perform simultaneous fits to the correlators in each case. The fits then capture the correlations and we can propagate them to the fit of Eq. (5). At the same time it allows us to determine the ratio of QCD+QED to pure QCD for the quantities that we will study. We will give results for these ratios in the sections that follow.
The fit form of Eq. (5) has been constructed such that the coefficients (apart from x) are expected to be of order 1. We therefore use priors of 0 ± 1 for all fit parameters except x for which we take a prior width on its expected value of 20% (the prior mean for x depends on the quantity being fitted).

A. Pure QCD
The hyperfine splitting, ∆M hyp , is calculated on each ensemble as the difference of the vector and pseudoscalar ground state masses, in lattice units, divided by the lattice spacing. The results for aM ηc and aM J/ψ and their difference are given in Table II for the pure QCD case. Although the pseudoscalar and vector correlators on each configuration are correlated the fit outputs for the vector correlator dominate the uncertainties and so the correlations have very little effect as a result.
The pure QCD results are plotted in Fig. 2 along with the 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 our results because we have a highly-improved quark action with small discretisation errors. Since all tree-level a 2 errors have been removed, the shape of the curve reflects the fact that higher-order a 4 and a 6 errors are visible; discretisation errors of this kind are present in all formalisms of course, but often they are hidden below much larger a 2 effects and consequently overlooked. Note also the clear dependence on the light sea quark mass seen on the finest lattices. To pin down the value of the valence mass mistuning parameter, c val , we include results at deliberately mistuned c quark masses (see Table II    III. QCD+QED ηc and J/ψ masses and hyperfine splitting presented as the ratio of the QCD+QED result to the pure QCD one on that set. Correlations between the calculations in the QCD+QED and pure QCD cases are used in the determination of the ratio and result in the very high statistical accuracy obtained. Note that the ratio is calculated for the same amc value in the two cases i.e. the ratio given here does not include the impact of retuning the c quark mass. QCD case in the continuum limit and for physical quark masses is 118.6(1.1) MeV, which is higher than the experimental average value, as is clear in Fig. 2. In order to understand what this means, we need to quantify all possible sources of small systematic effects in our calculation, 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 given in Table III. The correlation between the QCD+QED and the pure QCD results enables very high statistical accuracy to be obtained in the ratio. The inclusion of quenched QED shifts both the η c and J/ψ masses up by O(0.1%), depending on lattice spacing, at a given am c value. Although these mass shifts are small, there is a difference between the shift for the J/ψ and that for η c and so the inclusion of quenched QED also changes the hyperfine splitting. The impact here is more substantial, 0.7%, because the hyperfine splitting is so much smaller. The size of the direct QED effect on the hyperfine splitting can be simply estimated by replacing C F α s by Q 2 α QED in a potential model estimate of the splitting. This gives a fractional effect of Fig. 3. This shows that the results are consistent across all lattice spacings and thus discretisation effects in this ratio are smaller than for the masses themselves.
Finite-volume effects are an issue in general for QED corrections to meson masses but we expect them to be small for the electrically neutral and spatially small charmonium mesons that we study here. In [26] it is shown that the finite volume expansion for electrically neutral mesons starts at O(1/L 4 s ). In Fig. 4 we compare results for the fractional effect of QED on the J/ψ and η c as a function of 1/L s . This calculation is done on sets 5, 6 and 7 (see Table I) which differ only in their spatial extent. We see no finite-volume effects to well within 0.01%, and we therefore ignore such effects.
Our results including both pure QCD and QCD+QED are shown in Fig. 5, plotted against (am c ) 2 . The fit curve from Eq. (5) at physical quark masses is also shown. The fit has a χ 2 /dof of 0.59 and gives a hyperfine splitting in the continuum limit at physical quark masses of 119.6(1.1) MeV. Taking the (correlated) ratio between the physical value of the full QCD+QED fit and the physical value from the fit at Q = 0 (i.e. the pure QCD result) we obtain R QED [∆M hyp ] = 1.00804 (43). This ratio now does include the effect of retuning the c quark mass to obtain the experimental J/ψ mass when quenched QED is included. This retuning requires a reduction of the bare c quark mass by O(0.1%) (see Table III) and this further increases the hyperfine splitting, but only by O(0.1%). The impact of QED here is therefore dominated by the 'direct' quenched QED effect.  There is an additional pure QED contribution to the J/ψ mass that has not been included yet since it is quarkline disconnected. This comes from a diagram in which the cc annihilate into a photon which converts back into cc. The contribution of this diagram is where ψ(0) is the nonrelativistic J/ψ wavefunction equal (in the normalisation being used here) to f J/ψ M J/ψ / √ 6 [27]. The contribution evaluates to +0.7 MeV. This is added to our hyperfine splitting result with a 30% uncertainty from possible QCD corrections to give a final result of The error budget for our hyperfine splitting result is given in Table IV. We follow Appendix A of [28] for the meaning of the uncertainties contributing to the error budget. The majority of the uncertainty is associated with the lattice spacing determination, either from the correlated w 0 uncertainty or the individual w 0 /a uncertainties. This is not surprising because the hyperfine splitting is sensitive to uncertainties in the determination of the lattice spacing for the reasons discussed in [29]. We have separated out the uncertainty arising from the pure QCD data and the R 0 QED [∆M hyp ] values from Table III which we label 'Pure QCD Statistics' and 'QCD+QED Statistics' in Table IV. The sea mistuning uncertainty comes from the c m coefficients in Eq. 5 and the valence mistuning uncertainty from the c val and c val,Q coefficients. The a 2 → 0 uncertainty is from the c a and c aQ coefficients.
Our final result is for the charmonium hyperfine splitting determined from quark-line connected correlation functions in QCD and including the impact of QED effects, through explicit calculation in quenched QED. We expect the effect of further QED effects in the sea to be negligible compared to our 1% uncertainty. The only significant Standard Model effect then missing is that of quark-line disconnected diagrams in which the cc annihilate to gluons. We expect this effect to be much larger for the η c than for the J/ψ so a comparison of our result for the hyperfine splitting to experiment can yield information on the size and sign of the annihilation contribution to the η c mass. This is discussed in the next subsection.

C. Discussion: Hyperfine Splitting
The experimental average value of the hyperfine splitting (113.0(5) MeV) from the Particle Data Group (PDG) [1] is calculated as the difference of the separate averages for the J/ψ and η c masses. The different experimental results contributing to the PDG average of the two masses are shown in Fig. 6. For the J/ψ mass the average is dominated by the most recent result from KEDR [41]. There are only three experimental results used in these analyses that can independently produce values for the hyperfine splitting. These are the KEDR experiment [34,41] and two LHCb analyses in different channels [31,32]. The LHCb result in [32] used the η c (2S) → pp decay while the analysis of [31] used η c (1S) → pp. In the comparison plot of Fig. 7 [31] is referred to as LHCb15 and [32] as LHCb17. Fig. 7 shows a comparison of lattice QCD results for the charmonium hyperfine splitting along with the PDG average value and separate experimental values that measured this splitting. Previous calculations on gluon field configurations that included n f = 2 + 1 flavours of sea quarks by HPQCD [29] and by Fermilab/MILC [38] both obtained values above the experimental average, although only by just over one standard deviation.
The result we present here is substantially more precise than these earlier studies and for the first time displays a significant, 6σ, difference from the experimental aver- The ηc results represent a recent subset of those used in the PDG average. The most recent result is from BELLE (denoted BELL18) [30]. There are three different determinations from LHCb [31][32][33], two of which also measured the hyperfine splitting. We include a KEDR measurement [34], two from different BaBar analyses [35] and two from BESIII [36,37].
age, clearly showing that the lattice result lies above the experimental one. We interpret this as the effect of ignoring annihilation to gluons in the calculation of the η c mass. From the comparison of our results to experiment we conclude that these annihilation effects increase the η c mass by 7.3(1.2) MeV, where the uncertainty is dominated by that from the lattice calculation.
Previous analyses of this issue have given mixed results. In NRQCD perturbation theory [8] we can relate the shift in the η c mass to its total (hadronic) width through the perturbative amplitude for cc → gg → cc at threshold [42]. Then The leading term here gives −3.1 MeV, but sub-leading corrections could easily change the sign. An alternative way to think about the effect is non-perturbatively and then the gluon annihilation allows mixing between the η c and other flavour-singlet pseudoscalar mesons. Since Comparison of different lattice results for the charmonium hyperfine splitting and separate experimental results that measure this difference, as well as the PDG average (pink band). The PDG average is obtained from taking the differences of the PDG J/ψ and ηc masses (see Fig. 6) rather than only from experiments that directly measure the splitting. The squares give lattice QCD results from calculations that include n f = 2+1 flavours of quarks in the sea. The hexagons give 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 neglecting ηc annihilation removed so that we might expect some difference between them and experiment (see text), but previous lattice QCD results have not been accurate enough to see this. The Fermilab/MILC result used Fermilab c quarks on gluon field configurations with asqtad sea quarks [38] and the previous HPQCD result [29] used HISQ quarks on the same asqtad-sea ensembles. Briceño et al [39] used a modification of the Fermilab approach known as relativistic heavy quarks on the n f = 2 + 1 + 1 HISQ sea quark ensembles that we use here. The χQCD result [40] used overlap quarks on gluon field configurations including n f = 2 + 1 domain-wall sea quarks. these are lighter than the η c this mixing could give a positive correction to the η c mass. Direct lattice QCD determination of the effect, by calculating the appropriate quark-line disconnected correlation functions, has so far not proved possible. This is because the lighter states that are introduced by the mixing make it very hard to pin down a small effect on the mass of a particle, the η c , which is so much further up the spectrum in this channel. An estimate of the mass shift of +1-4 MeV was obtained in the quenched approximation in which this mixing is not possible but where mixing with glueballs could happen instead [43].
Our result for the hyperfine splitting, by its accuracy, provides for the first time a clear indication of the size of the impact of the η c annihilation to gluons on its mass: IV. DETERMINATION OF mc

A. Pure QCD
In [44] we showed that it is possible to determine the strange and charm quark masses accurately in lattice QCD using an intermediate momentum-subtraction scheme. By intermediate we mean that the mass renormalisation factor to convert the tuned bare lattice quark mass to the momentum-subtraction scheme is calculated on the lattice. The conversion from the momentumsubtraction to the final preferred MS scheme is carried out using QCD perturbation theory in the continuum. To do this accurately it is important to use a momentumsubtraction scheme that has only one momentum scale, µ. This means that the squared 4-momentum on each leg of the vertex diagram, from which the mass renormalisation factor is calculated, is µ 2 . The RI-SMOM scheme [45] used in [44] is such a scheme. A further important point is that the mass renormalisation factor will be contaminated by nonperturbative (condensate) artefacts through its nonperturbative calculation on the lattice. To identify and remove these artefacts (that appear as inverse powers of µ) requires calculations at multiple values of µ and a fit to the results, as discussed in [44].
Below we briefly summarise the procedure followed in [44]: 1. Determine the tuned bare quark mass and lattice spacing at physical sea quark masses for each set of gluon field configurations at a fixed β value. We do 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 values of µ. We thereby obtain the quark mass in the RI-SMOM scheme at scale µ.
3. Convert the mass to the MS scheme at scale µ using a perturbative continuum matching calculation. We denote this conversion factor by Z 4. Run all the MS quark masses at a range of scales µ to a reference scale of 3 GeV using the four loop QCD MS β function. We denote these running factors r(3 GeV, µ).
5. Fit all of the results for the MS mass at 3 GeV to a function that allows for discretisation effects and condensate contamination, which begins at 1/µ 2 with the nongaugeinvariant A 2 condensate.
6. Obtain from the fit the physical value for the quark mass in the MS scheme at 3 GeV with condensate contamination removed.
Here we provide three small updates to [44]. We first list them and then discuss them in more detail below. The three updates are: Lattice spacing values and tuned c quark masses at physical sea quark masses for each group of ensembles at a fixed β value, denoted by their name in column 1 (see Table I). The lattice spacing value is given in units of w0 in column 2 and the c quark mass, fixed from the J/ψ meson mass, is given in GeV units in column 3. The first uncertainty on the mass is uncorrelated between lattice spacing values, and the second is correlated.
FIG. 8. mc(3 GeV) extrapolated to the continuum with a fit form that allows for condensate terms that behave like inverse powers of the renormalisation scale µ. This plot is an updated version of the upper section of Fig. 10 in [44], with added data points on ultrafine lattices (set 14) and a retuned bare c quark mass fixed from the J/ψ, rather than ηc, meson. The ultrafine points have slightly mistuned µ values compared to the corresponding lines (see text).
1. we improve the uncertainty in the tuning of the bare lattice c quark mass by using the J/ψ mass rather than the η c ; 2. we include results from a finer ensemble of lattices (set 14) to provide even better control of the continuum limit; 3. we use the new 3-loop-accurate SMOM to MS matching 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 bare charm quark mass is done. In [44] bare charm masses were 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 effects (from perturbation theory). A 100% uncertainty was included on the adjustment [25]. Now that we are explicitly including quenched QED it makes more sense to have a tuning process that uses an experimental meson mass with no adjustments. We also want to use the same tuning process for both the pure QCD case and the QCD+QED case to allow for a clear comparison and one that reflects the procedures that would be followed in a complete QCD+QED calculation. This means that we should use the J/ψ meson mass, as we have done in Section III. The J/ψ mass is more accurate experimentally than that of the η c and the J/ψ has a much smaller width, implying little effect on the mass from its 3-gluon annihilation mode. The impact of J/ψ annihilation to a single photon is a sub-1 MeV shift to the mass which is a 0.02% effect, so negligible.
Using our J/ψ meson masses and following the procedure of [25] we obtain tuned bare c masses for each β value. These are given in Table V along with the w 0 /a values corresponding to physical sea quark masses at that β value, which are also updated slightly from [44]. These slight changes in w 0 /a lead to small adjustments in the µ values relative to those given in Table IV of [44]. This is accounted for when we run the Z m to the correct reference scale in MS.
The second update is to include results from the ultrafine β = 7.0 ensemble (set 14). The appropriate tuned mass and w 0 /a value for physical sea quark masses is given 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 α 3 s correction to the SMOM to MS conversion factor, Z MS/SMOM m , for the mass renormalisation. This correction was recently calculated in continuum perturbative QCD [46,47]. For n f = 4, as here, the α 3 s correction is a small effect ( 0.2%), continuing the picture seen at O(α s ) and O(α 2 s ) and consistent with the uncertainty taken from missing it in [44].
Once we have determined results for m c (3 GeV) at a variety of lattice spacing values using the SMOM intermediate scheme at a variety of µ values, we need to fit the results to determine m c (3 GeV) in the continuum limit. We do this following our previous calculation [44], where the fit function is given in Eq. (26). The fit includes discretisation effects and condensate artefacts in the lattice calculation of Z SMOM m . In [44] we included a term in the fit, c α α 3 s (µ) (with α s in the MS scheme) to allow for the then-missing α 3 s term in the SMOM to MS conversion. Here we replace that term with c α α 4 s (µ) since the conversion is now calculated through α 3 s and included in our results. We take a prior value on c α of 0.0 ± 0.5. This allows for the coefficient of the α 4 s term in the conversion factor to be 4 times as large as the α 3 s coefficient. The updated fit to our results for m c in the MS scheme at 3 GeV is shown in Fig. 8. The fit has a χ 2 /dof of 0.71. The error budget for this calculation is shown in Table VI. Most of the entries are very similar to those in [44]. The contribution due to the continuum extrapolation has, unsurprisingly, dropped a little, as has the uncertainty from the missing higher order terms (here α 4 s ) in the SMOM to MS conversion. The correlated tuning uncertainty comes largely from the uncertainty in the TABLE VI. Error budget (in %) for the calculation of the charm quark mass in the MS scheme at a scale of 3 GeV using RI-SMOM as an intermediate scheme. The listed contributions have the same meaning as those in [44] except that we use r here for the running factor rather than R and we have an additional one labelled 'QED effects' which comes from the continuum extrapolation shown in Fig. 10. physical value of w 0 used to fix the lattice spacing. We will be able to reduce that uncertainty in future by improving the determination of the value of w 0 . Our updated pure QCD result is Running this down to a scale equal to the mass gives m c (m c ) QCD = 1.2723(78) GeV.
These results improve on and supersede the value in [44].

B. Impact of Quenched QED
To include quenched QED effects in the determination of the c quark mass we must determine both the bare quark mass and the mass renormalisation factor with quenched QED switched on.
We include quenched QED effects for Z m in the RI-SMOM scheme in the same way as that described for the vector current renormalisation in [16]. This involves the generation of U(1) fields in Landau gauge to remain consistent with the Landau gauge QCD configurations used in the pure QCD calculation. When converting from the RI-SMOM scheme for Z m to the MS scheme it is also necessary to include QED effects in the perturbative matching factor. We can evaluate the QED contribution to Z MS/SMOM m at order α QED by multiplying the known coefficient of α s [45] by (3/4)Q 2 to give the coefficient of α QED . The impact is very small (< 0.1%) and we therefore safely neglect higher order terms. The numerical values of the O(α QED ) term we do include for the RI-SMOM to MS matching are given in Table VII. Note that the inclusion of QED in the lattice correlation functions is fully nonperturbative and it is only the continuum matching and running that are simply done to O(α QED ).
To assess the QED impact on the tuned bare c mass we use the QCD+QED J/ψ masses given in Table III and shown in Fig. 1. As we have corrected the pure QCD determination of m c to be tuned to the experimental J/ψ mass this is the tuning we will use for the QCD+QED case as well. The fractional shift in am c required to obtain the correct J/ψ mass after QED has been included (which we denote R QED [am c ]) can be evaluated from the fractional change in the J/ψ mass. R QED [am c ] is the fractional change in am c required to return the J/ψ mass to the value it had in pure QCD (i.e. the experimental value) once QED is switched on. Thus the increase in J/ψ mass seen with QED must be compensated by a corresponding reduction in am c . From the deliberately mistuned am c values in Table II we see that the fractional change in the c mass is 1.5 times larger than the change seen in the meson mass. We can use this factor, along with R 0 QED [M J/ψ ] values from Table III and the required change of sign in the shift, to determine the retuning of 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 0 QED [M J/ψ ] by a factor of 5 to allow for the uncertainty in our conversion factor of 1.5 above.
We calculate the ratio of Z m in the SMOM scheme in QCD+QED to pure QCD (R QED [Z SMOM m ]) following the methods of [44]. The calculations are carried out at multiple values of the renormalisation scale µ and extrapolated to zero valence quark mass. Results are given in Table VII. Notice that these values are larger than 1 and so compensate to a large extent for the changes in the tuning of the bare lattice am c induced by QED. This reflects the fact that most of the shift is an unphysical ultraviolet self-energy effect. We then convert from the SMOM scheme to MS by making use of the ratio of the perturbative conversion factor for QCD+QED to pure QCD, i.e. the O(α QED ) piece of Z MS/SMOM m also given in Table VII. This is less than 1 but only by a very small amount.
Multiplying these two factors together gives the ratio of the lattice to MS mass renormalisation factors for QCD+QED to that for pure QCD, i.e. R QED [Z MS m ]. From Table VII, multiplying columns 3  The leading logarithm at each order can be derived from the anomalous dimensions of the mass, allowing us to write [48] Here C m is a constant, up to discretisation effects and higher order terms multiplying powers of log(aµ). Figure 9 plots our results for C m . These show very little variation with a and µ, confirming that the dependence on a and µ of R QED [Z MS m ] is almost entirely captured by Eq. (15).
Once the impact of QED on the c mass in the MS scheme at scale µ is obtained, as above, we then need to allow for QED effects in the running of the masses from µ to the reference scale of 3 GeV. This is done by multiplication by a factor r QED calculated in O(α QED ) perturbation theory and given in Table VII. These numbers are also very close to 1. The α s α QED term could in principle have some impact here but it is very small and we neglect it [48].
Multiplying R QED [am c ] and R QED [Z MS m ] together allows us to determine the ratio of the c quark mass in the MS scheme at 3 GeV from QCD+QED to that in pure QCD, i.e. R QED [m c (3 GeV)]. The values that we have for this ratio come from results at multiple values of µ and multiple values of the lattice spacing. To determine the physical ratio in the continuum limit with condensate contamination removed (in this case QED corrections to QCD condensates) we need to fit the results to a similar form to that used in [44]. We use The first term on the second line of Eq. (16) accounts for discretisation effects in R QED [am c ]; a scale of 1 GeV is chosen in these effects as this is close to the c quark mass. The term multiplying this (on the bottom two lines) models the a and µ dependence of the QED correction to Z m . This includes discretisation effects of the form (aµ) 2i and terms to model condensate contributions, starting at 1/µ 2 . The priors on all coefficients are taken as 0.0 ± 1.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 the fit output are shown in Figure 10. The fit has a χ 2 /dof of 0.87 and returns a physical value of R QED [m c (3 GeV)] of 0.99823 (17). We conclude that the impact of quenched QED is to lower the c quark mass, m c (3 GeV) by a tiny amount: 0.18(2)%.
very close to the pure QCD value at this scale. This is the first determination of the c quark mass to include QED effects explicitly, rather than estimate them phenomenologically as has been done in the past. The uncertainty achieved 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 continuum limit and from missing higher order terms in the SMOM to MS matching.
C. Discussion: mc Figure 11 gives a comparison of lattice QCD results for m c . We plot the masses at the scale of the mass, as is conventional even though this is a rather low scale. We restrict the comparison to results that were obtained on gluon field configurations including u, d, s and c quarks in the sea. The top result is our value from Eq. (17) that explicitly includes a calculation of the impact of quenched QED on the determination of the quark mass.
The top three results in the pure QCD section of the figure all include an estimate of, and correction for, QED effects. These corrections are made, however, by allowing for 'physical' QED effects such as those arising from the Coulomb interaction between quark and antiquark in a meson. They do not allow for the QED self-energy contribution which is substantial. Although a large part of this is cancelled by the impact of QED on the mass renormalisation, a consistent calculation has to include both effects, as we have done here.
An important point about Figure 11 is that the top three pure QCD results all have uncertainties of less than 1% and agree to better than 1%, using completely different 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 result 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 from the matrix element between the vacuum and a J/ψ meson at rest by where µ is the component of the polarisation of the J/ψ in the direction of the vector current. In terms of the ground state amplitude, A V 0 , and mass, M V 0 (≡ E V 0 ), obtained from the fit of Eq. (3) to the charmonium vector correlator it is (in lattice units) Z V is the renormalisation factor required to match the lattice vector current to that in continuum QCD if a nonconserved lattice vector current is used (as here). We discuss the renormalisation of vector currents using intermediate momentum-subtraction schemes in [16] and we will make use of the results based on the RI-SMOM scheme here (see Section II). Note that there is no additional renormalisation required to get from the RI-SMOM scheme to MS because the RI-SMOM scheme satisfies the Ward-Takahashi identity [16]. The partial decay width of the J/ψ to an + − pair ( = e, µ) is directly related to the decay constant. At leading order in α QED and ignoring (m /M J/ψ ) 4 correction terms, the relation is where Q c is the electric charge of the charm quark in units of the charge of the proton. Note that the formula contains the effective coupling, α QED,eff evaluated at the scale of M J/ψ but without including the effect of the J/ψ resonance in the running of α QED to avoid double-counting [51]. Experimental values of Γ(J/ψ → e + e − ) are obtained by mapping out the cross-section for e + e − → e + e − and e + e − → hadrons through the resonance region [52] or by using initial-state radiation to map out this region via e + e − → µ + µ − γ [53]. In either case initial-state radiation and non-resonant background must be taken care of [54,55]. A cross-section fully inclusive of final-state radiation is obtained; interference between initial and finalstate radiation is heavily suppressed [56]. The resonance parameter determined by the experiment is then the 'full' partial width [55,57], where Γ (0) is the partial width to lowest order in QED and Π 0 is the photon vacuum polarisation. The effect of the vacuum polarisation is simply to replace α QED in the lowest-order QED formula for the width with α QED,eff (M 2 ), as we have done in Eq. (21). The experimental determination of Γ is accurate to 2% for the J/ψ [1]. This allows us to infer a decay constant value from experiment, accurate to 1%, using Eq. (21).
The first uncertainty comes from the experimental uncertainty in Γ and the second is an O(α QED /π) uncertainty for higher-order in QED terms, for example from finalstate radiation, in the connection between Γ and f in Eq. (21). Note that using α QED of 1/137 would increase this number by 2.3% (9 MeV). This experimental value can then be compared to our lattice QCD results for a precision test of QCD. Here we improve on HPQCD's earlier calculation [29] by working on gluon field configurations that cover a wider range of lattice spacing values and with sea u/d quark masses now down to their physical values. In addition we now include c quarks in the sea and have a more accurate determination of the vector renormalisation factor Z V [16]. We will also test the impact on f J/ψ of the c quark's electric charge.  2 and 3 give the results in lattice units for the J/ψ and ηc decay constants respectively in pure QCD on the ensembles listed in Table I. The values of am val c are those given in Table I except for two cases with a deliberately mistuned c quark mass: set 6 denoted by a * where amc = 0.643 and set 14 denoted by a † where amc=0.188.
The results for f J/ψ do not include the multiplication by ZV needed to normalise them (Eq. (20)). Columns 4 and 5 give the electromagnetic corrections for the J/ψ and ηc decay constants, as the ratio of the QCD+QED result to that in pure QCD. Again, the electromagnetic corrections for f J/ψ do not include the corrections to ZV . ZV values are given in Table IX. The decay constant of the pseudoscalar η c meson is determined from our pseudoscalar correlators (of spintaste γ 5 ⊗γ 5 ) using the ground-state mass and amplitude parameters from the correlator fit, Eq. (2): Note that this is absolutely normalised and no Z factor is required. Because the η c does not annihilate to a single particle there is no experimental process from which we can directly determine f ηc . Nevertheless it is a useful quantity to calculate for comparison to f J/ψ and to fill out the picture of these hadronic parameters from lattice QCD [59]. Again we will improve on HPQCD's earlier calculation [60] as discussed above for the J/ψ.

A. Pure QCD
The second column of Table VIII gives our results for the (unnormalised) values of af J/ψ in pure QCD on 12 of the sets from Table I. We multiply af J/ψ /Z V by the value of Z V and convert to physical units using the inverse lattice spacing. Z V values are taken from [16] except for a new value calculated here for β = 7.0 (ultrafine) set 14. We collect these values in Table IX. See Appendix A for a discussion of the Z V results. The Z V values are very TABLE IX. Vector current renormalistion constants, ZV (µ), using the RI-SMOM scheme at µ = 2 GeV (column 2) and µ = 3 GeV (column 3) in pure QCD for each β value corresponding to a group of ensembles in Table I. Column 4 gives the QED correction to ZV at 2 GeV in the form of the ratio of the QCD+QED value to that of pure QCD. Most of these values are taken from [16] although the ZV value at β = 7 and the QED correction at β = 5.8 are new here. precise and so have little impact on the uncertainty in f J/ψ .
Our results for f J/ψ for the pure QCD case are shown in Fig. 12 plotted against the square of the lattice spacing (in units of the bare c quark mass). Clear dependence on the lattice spacing is seen. This dependence comes from the amplitudes of the two-point correlators; the lattice spacing 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 the prior on the physical value of f J/ψ (i.e. x in Eq. (5)) of 0.4(1). The χ 2 /dof of the fit is 0.43. The agreement with the result derived from experiment can clearly be seen.
We will discuss this result further in Section V C. 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 discretisation effects and results for the physical quantity, f J/ψ , using different renormalisation scales µ should agree in the continuum limit. Here we verify that this is the case using 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 QCD case. The difference between the two values of µ is barely visible. The values at µ = 3 GeV give a continuum limit result of f J/ψ = 408.7(1.8) MeV, in good agreement with that at µ = 2 GeV in Eq. (26) but slightly less accurate. The χ 2 /dof of the fit was 0.45.
Our results for af ηc , the η c decay constant in lattice units, are given in the third column of Table VIII for the pure QCD case. After conversion to physical units, they are plotted in Figure 14. The curve is similar to that for f J/ψ but with somewhat smaller discretisation effects. We also plot the results of performing the same fit as for f J/ψ using Eq. (5). The χ 2 /dof of the fit is 0.88 giving a result for the decay constant in pure QCD of f ηc,QCD = 397.5(1.0) MeV. (27) This agrees well with the earlier HPQCD value on n f = 2 + 1 gluon field configurations [60] of 0.3947(24) GeV but has half the uncertainty. In [60] the effects from neglecting the charm quark in the sea are estimated to be O(0.01%) which is negligible and means that the two calculations should give the same result. Figure 15 shows our results for the ratio of f J/ψ to f ηc in pure QCD. A lot of the discretisation effects cancel in the ratio, as is evident in comparing this figure to Figures 12 and 14. Systematic uncertainties, for example from the determination of the lattice spacing, are also reduced. The shape of the curve again, as in the hyperfine splitting case, reflects the fact that we have successfully reduced sources of a 2 error to the point where a 4 and a 6 are visible. We fit the ratio to the same fit as before (Eq. (5)) with a prior on the physical value of 1.0(1). The fit has a χ 2 /dof of 0.62 and returns a physical value for the decay constant ratio in pure QCD of f J/ψ,QCD f ηc,QCD = 1.0285 (18).  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 values of f J/ψ and f ηc , both for the pure QCD case and the QCD+QED case to be discussed in Section V B. The contributions from different sources are very similar between the two decay constants. It is clear from this that the dominant sources of error are related to the determination of the lattice spacing, as for the hyperfine splitting.
The error budget presented here for the J/ψ decay constant is markedly different from that of [29]. There the dominant contribution to the error was from the vector renormalisation constant, Z V , obtained using a matching between lattice time moments and high order perturbative QCD [61]. Here that error is substantially reduced by using the Z V values obtained in lattice QCD fully nonperturbatively in the RI-SMOM scheme [16].

B. Impact of Quenched QED
Including quenched QED effects into our calculations allows us to determine the effect on the J/ψ and η c decay constants of the electric charge of the valence c quarks. Because the J/ψ and η c are electrically neutral particles, there is no long-distance infrared component to cause problems (as there is for f π + ) and we can simply proceed to determine the decay constants after the addition of the QED field as we did in the pure QCD case. The fractional QED effect on the J/ψ and η c decay constants at fixed bare c quark mass in lattice units is given in Table VIII. We see a 0.3% increase, offset slightly by the change in Z V in the J/ψ case. The fractional QED effect on Z V is given in Table IX. The fractional QED effect at fixed bare mass is plotted in Figure 16. We see that the effect is similar for the J/ψ and η c in the continuum limit and shows very little dependence on the lattice spacing.
The volume dependence of the fractional QED effect is shown in Figure 17 on sets 5-7. We find that the effect is negligible well below the 0.1% level.  Table I in pure QCD (red points) and including also quenched QED (blue points) plotted against the square of the bare c quark mass in lattice units. The green curve marks our extrapolation to the physical point. The black cross gives the earlier HPQCD result on n f = 2 +1 gluon field configurations from [60].
We now combine our QCD+QED results with our pure QCD results and the full fit of Eq. (5), which takes into account the retuning of the quark mass needed when quenched QED is included. Figure 18 shows our pure QCD results, QCD+QED results and fit curve for the J/ψ decay constant. We obtain This is a 0.2% increase over the value in pure QCD (Eq. (26)) because retuning reduces the quark mass and offsets some of the impact of quenched QED seen in Table VIII. A more accurate statement is that the final fractional effect from quenched QED is .00188 (36). A very similar picture is seen for f ηc in Figure 19. We obtain 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 extrapolated to the continuum limit. This gives f J/ψ,QCD+QED f ηc,QCD+QED = 1.0284 (19). (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 and QCD+QED results to previous results including n f = 2 + 1 flavours of sea quarks for f ηc [60] and f J/ψ [29]. There is good agreement. These earlier calculations also used HISQ quarks but our new results are more accurate, particularly for f J/ψ because of the use of more accurate values of Z V [16].
There have also been calculations that use n f = 2 flavours of sea quarks. It is harder to make a comparison to these results because it is not clear what the systematic error is from not including at least the s quarks in the sea, and no uncertainty is included for this. The calculation of [62] uses twisted-mass quarks on n f = 2 gluon The σ here is that from the n f = 2 results since our uncertainty is much smaller. The n f = 2 results for f J/ψ are compatible with each other and with our result, again at 2σ. As discussed in Section V, the J/ψ decay constant is the hadronic quantity that is needed to determine the rate of J/ψ annihilation to leptons via Eq. (21). Our result for Γ using Eq.
The first uncertainty is from our lattice QCD+QED result for f J/ψ and the second uncertainty allows for a relative O(α QED /π) correction to Eq. (21) from higher-order effects. Figure 22 compares this width Γ(J/ψ → e + e − ) to results from experiment. Recent experimental results from KEDR [52] and BES III [53] are shown along with the Particle Data Group average [1] (as a grey band). With new results expected from the Fermilab g − 2 experiment soon there has been a concerted effort by the lattice community to understand and control systematic effects in the lattice QCD calculation of the hadronic vacuum polarisation (HVP) contribution to the anomalous magnetic moment of the muon. Since the most accurate values for the HVP currently come from experimental results on R(e + e − → hadrons), it is also important to compare lattice QCD results to these, disaggregated by flavour where possible.
The first calculation of the quark-line connected cquark contribution to the HVP, a c µ , was given in [64] using results for the time-moments of vector charmonium current-current correlators calculated in [29]. The time moments are defined by where C V (t) is the vector current-current correlator and Z V is the vector current renormalisation factor, discussed in Section V. Note that t ∈ {−L/2 + 1, −L/2 + 2, . . . , L/2}. The even-in-n time-moments for n ≥ 4 can be related to the derivatives at q 2 = 0 of the renormalised vacuum polarisation function [65],Π(q 2 ) ≡ Π(q 2 ) − Π(0), by G n = (−1) n/2 ∂ n ∂q n q 2Π (q 2 ) This means thatΠ(q 2 ) can be reconstructed, using Padé approximants, from the G n [64] and fed into the integral over q 2 that yields the quark-line connected HVP contribution to a µ [66]. Only time-moments of low moment number are needed to give an accurate result for a c µ because the integrand is dominated by small q 2 . We will give results for the four lowest moments, n = 4, 6, 8, 10, improving on the values given in [29]. The improvement comes mainly through the use of a more accurate vector current renormalisation as well as an improved method for reducing lattice spacing uncertainties but we also use second-generation gluon field configurations that include c quarks in the sea and calculate, rather than estimate, the impact of the leading QED effects. Π(q 2 ) and hence a c µ can also be determined from experimental results for R c ≡ R(e + e − → cc → hadrons) as a function of squared centre-of-mass energy, s. This can be done using inverse-s moments R c is obtained from the full e + e − rate from just below the c threshold upwards by subtracting the background contribution from u, d, and s quarks perturbatively, see e.g. [61]. FIG. 23. A study of the volume dependence of the electromagnetic correction to the first four time moments of the vector current-current correlator and to a c µ on sets 5-7. There is no observable dependence on the lattice spatial extent, Ls, as can be judged by comparison to the dashed horizontal lines at the weighted averages of the points.
The relationship between G 2k+2 and M k is then A comparison of our correlator time-moments calculated on the lattice and extrapolated to the continuum limit to the inverse-s moments determined from experiment is equivalent to a test of the agreement of the results for a c µ in the two cases.
A. Vector correlator moments: Pure QCD and QCD+QED results Table XI gives our raw results for the time-moments of the same vector current-current correlators from which we have determined the mass and decay constant of the J/ψ meson in Sections III and V. Notice that the statistical uncertainties are tiny. The correlators make use of a local vector current that must be renormalised as discussed in Section V. The results in Table XI are calculated before renormalisation and are given in lattice units. The quantity that is tabulated is We take the (n − 2)th root to reduce all results to the same dimensions [29]. To normalise the time-moments we use the Z V values at µ= 2 GeV given in Table IX that were used for f J/ψ in Section V. Table XI also gives the result of including quenched QED as the ratios R 0 QED for each rooted moment (at fixed am c ). These values are all very slightly less than 1,  Table I. The results tabulated are values of (Gn/Z 2 V ) 1/(n−2) in lattice units along with their statistical uncertainties, given for n = 4, 6 , 8 and 10 in columns 2, 3, 4 and 5. Uncertainties are statistical only. Note that the results for different moments are correlated because they are determined from the same correlation functions. The results marked with a * and † are for deliberately mistuned amc values as detailed in the caption to Table II. Also included are the ratios at fixed amc, denoted R 0 QED , of QCD+QED to QCD results for each rooted moment on a subset of ensembles. by up to 0.1% for n=4, and 0.2% for n=10. We can also test the finite-volume dependence of the quenched QED effect using sets 5, 7 and 8 and the results are shown in Fig. 23. There is no visible volume dependence in the QED effect on time-moments at the level of our statistical uncertainties. This is to be expected, as seen for the J/ψ mass and decay constant in Sections III and V, since the vector current being used here is electrically neutral. To fit the results as a function of lattice spacing it is convenient to work with the dimensionless combination: using our M J/ψ masses from Table II. This reduces the uncertainty coming from the value of the lattice spacing. At the same time it also reduces the sensitivity to mistuning of the valence c quark mass. We fit the quantity defined in Eq. (38) as a function of lattice spacing, allowing for quark-mass mistuning effects of both valence and sea quarks to derive results in the physical continuum limit. We do this as before using the fit function given in Eq. (5). Values in the physical continuum limit are then divided by M J/ψ from experiment to obtain our final results for G 1/(n−2) n in GeV −1 . In doing this we allow an uncertainty of 0.7 MeV in M J/ψ from annihilation to a photon (see Section III) since this effect is not included in our results. Fig. 24 plots the results for the rooted moments multiplied by M J/ψ as a function of (am c ) 2 . Also shown is the fit result from Eq. (5). Only the pure QCD lattice results are shown for clarity; those including the effect of quenched QED are very close to them. The fit result plotted 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 continuum limit for the first four time-moments (column 2) compared with 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. Column 4 gives the effect of quenched QED as a ratio of the physical results in QCD+QED to those in pure QCD.   (37) 1.00037 (10) 6th, 8th and 10th rooted moments in the physical continuum limit. These are obtained from a simultaneous fit to all of the moments, including the correlations between them (since they are derived from the same correlators). The fit has a χ 2 /dof of 0.62 for an svdcut of 1 × 10 −3 . Column 3 of Table XII gives the results derived from experimental data for R(e + e − → hadrons) and Γ and mass for the J/ψ and ψ by [67] for comparison (see Section VI). We have converted these results into the quantity that we calculate according to Eq. (36). We see agreement within 2σ for n ≥ 6 with the values derived from experiment, but a 2.4σ tension at n = 4. The results given from [67] correspond to their standard selection of experimental datasets. Results shift by ±1σ for other selections. Note that the lattice QCD results are now much more precise than those determined from experiment. The comparison between lattice QCD and experiment will be discussed further in Section VI B.
Column 4 of Table XII gives the final impact of ], are greater than 1, in contrast to the R 0 QED of Table XI that are less than 1. As discussed in Section II B, the inclusion of QED means that the c quark mass must be retuned downwards. The rooted time-moments are approximately inversely proportional to the quark mass and so they increase under this retuning, more than offsetting the direct effect of QED seen in R 0 QED . HPQCD's 2012 result on n f = 2 + 1 gluon field configurations with HISQ valence c quarks is marked 'HPQCD HISQ 2012' [29]. JLQCD's 2016 result using n f = 2 + 1 domainwall quarks is marked 'JLQCD DW 2016' for the 6th moment only. We also compare to results (open red squares) denoted 'pheno.' that are derived from experimental data for the cross-section for e + e − to hadrons as a function of centreof-mass energy, by determining the cc component. The points plotted come from [61] and [67]. Open orange circles show alternative selections of datasets from [67]; the upper value is for the 'maximal' set and the lower value for the 'minimal' set.

B. Discussion: Vector correlator moments
Our new results for the time-moments of vector current-current correlators improve significantly on earlier lattice QCD calculations and are now more accurate than results derived from experiment. Figure 25 compares our new results for the 4th and 6th moments to earlier lattice QCD results. Comparison for the 8th and 10th moments gives a very similar picture and so is not shown. The first lattice QCD calculation of the time-moments of vector charmonium current-current correlators was given by HPQCD in [29] using HISQ valence quarks on n f = 2 + 1 asqtad gluon field configurations. Our new results have an uncertainty almost ten times smaller than these. The error budget in the earlier results was dominated by the uncertainty in Z V from the use of continuum perturbation theory in the matching factors and there was also a sizeable uncertainty from the lattice spacing. These uncertainties have been enormously reduced here and in addition we no longer need an uncertainty from missing QED effects. We also show a comparison with the subsequent results from the JLQCD collaboration [68] using domain-wall quarks on n f = 2+1 gluon field configurations. JLQCD do not give a value for the 4th moment because of discretisation effects in their formalism (tree-level O(a 2 )). The dominant uncertainties in their results are from statistics and from the value of t 1/2 0 used to fix the lattice spacing, on which they have a 2% uncertainty. Good agreement is seen for all of the lattice QCD results.
Two of the most recent results from phenomenological determinations of the moments [61,67] are also compared in Fig. 25. The results from [67] include experimental datasets for the inclusive cross-section that are both older and newer than those used in [61]. Results from [67]'s 'standard' selection of datasets were given in Table XII and are shown in Fig. 25 in red. We also show, in orange, the results from the 'maximal' set (all experimental 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, keeping the most accurate results). Note that the resonance parameters are the same for all selections. We see that the variation with dataset selection covers almost 1σ for the 4th moment, but much less for the 6th moment. This is also 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 is defined to be the result from diagrams with a charm quark loop connected to a photon at both ends [61] i.e. the quark-line connected vector current-current correlator that we study on the lattice. The subtracted background includes QED effects for the non-charm and singlet (quark-line disconnected) contributions. The remainder, R c then includes the QED effects associated with the cc loop. The dominant source of uncertainty in R c comes from the charmonium resonance (J/ψ and ψ ) region and is set by the uncertainty in Γ ee for these states. The fractional uncertainty is approximately the same for all moments [61,67]. When the (n − 2)th root is taken the fractional uncertainty then falls with increasing n.
Good agreement is seen between the phenomenological results and our new lattice results for n = 6, 8 and 10, although the lattice results are systematically at the upper end of the phenomenological range. The largest discrepancy is a 2.8σ tension for the 4th moment between us and the results of [67] for their minimal selection of datasets. The tension is 2.4σ for the standard selection, and below 2σ for the maximal selection and for the results of [61]. The σ here is that for the phenomenological results since the lattice uncertainty is much smaller. Because the 4th moment dominates the determination of a c µ , this tension between lattice QCD+QED and some of the phenomenological results carries over to a c µ , to be discussed in the next section.
The time-moments can also be used to determine a value for m c by comparing to O(α 3 s ) continuum QCD perturbation theory and this was the focus of [61,67]. We do not do this here because the scale of α s is rather low in these determinations meaning that uncertainties from missing higher-order corrections can be substantial. We prefer instead the method of [25], which enables a higher scale to be used in the perturbation theory. We have checked, however that the m c value that would be obtained from the time-moments is consistent with both [25] and the value given in Section IV.
C. a c µ : Pure QCD and QCD+QED results To calculate the quark-line connected HVP contribution to a µ from c quarks, a c µ , we can either use the physical results for the vector current-current correlator timemoments discussed in the previous subsection or we can calculate a c µ on each lattice ensemble and perform a fit as a function of lattice spacing to extrapolate directly to the continuum limit. We will do the latter here.
The values of a c µ on each lattice are given in Table XIV. These are determined from the time-moments rescaled by the J/ψ mass on each lattice : As discussed in Section VI A rescaling by M J/ψ reduces the lattice spacing uncertainty and the impact of mistuning the c quark mass. This was used for lighter quark masses in [69] (see also [70]). The reduced effect of mistuning is clear from comparing the mistuned results in Table XIV to those in Table XI.  Table XIV also gives the direct effect of quenched QED through the ratio R 0 QED [a c µ ]. Because the rescaled moments have less sensitivity to the c quark mass these numbers are larger than 1 (unlike the results in Table XI) and reflect more closely the final impact of QED on a c µ . We observe no finite-volume dependence for the quenched QED corrections to a c µ , as for the correlator time moments. This is shown in Fig. 23.
The results from Table XIV a c µ = 14.638(47) × 10 −10 , QCD + QED along with R QED [a c µ ] = 1.00214 (19). These results agree well (within 1σ) with those obtained using the extrapolated values for the moments from Table XII and calculating a c µ in the continuum, as is seen in Fig. 26. The error budget for our final value of a c µ is given in Table XV. The largest uncertainties come from the deter-  [25] using values of time-moments determined in [29] using HISQ valence quarks on gluon field configurations including n f = 2 + 1 flavours of asqtad sea quarks. The other results all include n f = 2 + 1 + 1 flavours of sea quarks. The result labelled 'ETMC twisted mass 2017' uses the twisted mass formalism [71] and that labelled 'BMW stout stagg. 2017' a smeared staggered quark action [72]. Our new result (from Eq. (40)) labelled 'HPQCD HISQ' agrees with, but is much more accurate than, these earlier results. mination of the lattice spacing, although this uncertainty is much reduced by our rescaling with M J/ψ . D. Discussion: a c µ Our new result for the pure QCD case (Eq. (40)) is compared to earlier lattice QCD results with realistic sea quark content in Fig. 27. Our new result (the top point on the plot) is much more accurate than the earlier results. With respect to the first calculation of a c µ that also used HISQ quarks [25] we have greatly reduced the previous dominant sources of uncertainty from Z V and the determination of the lattice spacing.  Table XIV (red open symbols) and those from the BMW collaboration in [72] (filled blue circles) that use a less highly-improved staggered quark action. The filled blue circles only give estimates of the position of the BMW data points and do not indicate the statistical uncertainties which are much smaller than the size of the points.
With respect to results using other formalisms, we give one figure to demonstrate the control of discretisation effects that is possible with the HISQ formalism. Figure 28 compares the approach to the continuum limit of our results (from Table XIV) with results from BMWc [72], for which the continuum extrapolated result is shown in Fig. 27. The points plotted from [72] are estimates of the positions read from Figure S4, and do not include any indication of statistical uncertainties. The BMWc stout staggered quark action has O(a 2 ) discretisation errors at tree-level since it uses an unimproved derivative in its version of the Dirac equation. Figure 28 shows that the price to be paid for not improving the discretisation is a very large discretisation effect. This is particularly evident when working with heavier quarks such as charm, since it means that the dominant (Λa) 2 effects have Λ ≈ 1 GeV, as here. This means, for example, that the BMWc points at finest lattice spacing (∼0.06 fm) are about as far (20% away) from the continuum limit as our points at our coarsest lattice spacing (∼0.15 fm). Our results at a ∼ 0.06 fm are within 1% of the continuum limit, allowing us to achieve a sub-1% uncertainty in the final value.
We find that the impact of quenched QED on the result for a c µ is +0.214(19)% (Eq. 40). This is a shift in a c µ in the presence of the dominant QED effect of δa c µ = +0.0313(28) × 10 −10 .
We can compare this to the result obtained by ETMC in [5] following work on the QED effect on the renormalisation of quark bilinears in [23]. The ETMC value is +0.0182(36) × 10 −10 , 2.8σ smaller than ours. The two calculations of the quenched QED effect are different. Ours is a direct calculation of the quenched QED effect, retuning the valence quark mass through determination of a meson mass in the usual way. The The red open squares, denoted 'pheno.', use the charmonium moments from R e + e − given in [61,67] to determine a c µ . The orange open circles give the alternative maximal (upper) and minimal (lower) inclusive cross-section datasets from [67].
ETMC calculation is perturbative in quenched QED and fixes the valence quark masses so that they agree in the MS scheme at 2 GeV in QCD+QED and pure QCD. Our results in Section IV show that the quark mass in the MS scheme is lower in QCD+QED than in QCD (at scales above m c ) when the quark mass is tuned so that M J/ψ agrees with experiment in the two cases. This leads us to expect a larger result for the impact of quenched QED on a c µ with our tuning. Once full lattice QCD+QED calculations are underway tuning of the quark masses will be done 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 moments determined from J/ψ and ψ properties and the inclusive cross-section for e + e − → hadrons, removing contributions from quarks other than c. The results for these moments from [61,67] were discussed and compared to our results in Section VI B. In Fig. 29 we have converted the moments into a result for a c µ for comparison. As with the moments, we see a 2.5σ tension with the [67] result using the standard selection of datasets (a c µ = 14.03(24) × 10 −10 ), and a larger tension with the minimal selection of datasets. There is agreement within 2σ for the maximal selection of datasets from [67] and with [61]. Our result may then provide a pointer to the selection of R e + e − datasets for the phenomenological determination and/or indicate an issue with the perturbation theory used to subtract the u/d/s contribution to obtain R c (s).
A more recent determination of the complete HVP contribution to a µ using results from R e + e − is given in [51]. Although the c component is not separated out, the contribution from the J/ψ resonance is given as 6.26(19) × 10 −10 . We can also readily determine the contribution to a c µ from the J/ψ resonance alone since in that case [73] G n = n! f 2 Using our value for f J/ψ from Eq. (29) and the experimental value for M J/ψ [1] gives a c µ (J/ψ) = 6.345(53) × 10 −10 .
This is in good agreement with the result determined from the experimental J/ψ parameters, but more accurate, reflecting the fact that our result for Γ e + e − in Fig. 22 agrees with experiment but has smaller uncertainty. The J/ψ contribution provides almost half of a c µ -we conclude that it is the rest of the contribution, from the inclusive cross-section above the resonance region, that causes the tension between our results for a c µ and that from R e + e − for some selections of experimental datasets. The tension then amounts to 7(3)% of this non-resonant cross-section.

VII. CONCLUSIONS
We have performed the first n f = 2+1+1 lattice QCD computations of the properties of ground-state charmonium mesons. These have been done using the HISQ action to calculate quark-line connected two-point correlation functions on gluon field configurations that include u/d quark masses going down to the physical point. The small discretisation effects in the HISQ action and high statistics achievable have given us good control over both the continuum and chiral extrapolations and has enabled us to obtain smaller uncertainties than previous lattice QCD computations of these properties (including previous calculations by HPQCD). At the same time we have improved the tuning of the bare c quark mass to update the value of m c . From the same correlators that we use for the masses and decay constants we have also derived an improved result for the c quark HVP contribution to the anomalous magnetic moment of the muon.
The precision possible for c quark correlators with the HISQ action makes it possible to determine the impact of the c quark's electric charge. We do this directly and nonperturbatively in quenched QED by multiplying an appropriate U(1) field into our gluon field configurations. We tune the bare c quark mass so that the J/ψ mass agrees with experiment in both QCD+QED and QCD and we calculate mass and vector renormalisation constants in the RI-SMOM scheme in both cases, performing 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 (from Eqs. (10), (17), (29), (30), (40) Error budgets are given in Tables IV, VI, X and XV respectively. The precision of our result for the charmonium hyperfine splitting allows us to resolve, for the first time, the sign and magnitude of the anticipated difference between the lattice and experimental results arising from the fact that we do not include quark-line disconnected correlation functions. We take this difference to be the effect of the η c decay to two gluons which is prohibited in the lattice calculation, and conclude that ∆M annihln ηc = +7.3(1.2) MeV.
The effect of QED on the hyperfine splitting is fairly substantial (1.4%) and the largest effect that we observe here. This has 3 components that all act in the same direction: a direct effect of 0.7%, 0.1% from retuning the c quark mass and 0.6% from J/ψ annihilation to a photon that we add by hand.
Our updated value of the charm quark mass (in the MS scheme at a scale of 3 GeV) includes the effect of QED on the bare c quark mass, tuned so that M J/ψ matches experiment, 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 lattice spacing has improved the uncertainty slightly over HPQCD's earlier result [44].
The impact of QED on m c is small since QED effects tend to cancel between the retuning needed and changes in 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 most precise to date and acts as a subpercent test of QCD. The gain in precision from the 2012 HPQCD calculation of [29] is a result of the use of a more accurate renormalisation of the vector current [16] as well as gluon field configurations with a wider range of sea u/d quark masses and lattice spacing values. We can use our result for f J/ψ to determine a value for the width for J/ψ decay to a lepton-antilepton pair (repeating Eq. (32)): Γ(J/ψ → e + e − ) = 5.637(47)(13) keV. (45) This is more accurate than the current average of experimental results [1]. The impact of quenched QED on f J/ψ is +0.2%, since the retuning of the c quark mass offsets some of the direct effect. The effect on f ηc is almost the same so that the ratio of f J/ψ to f ηc remains the same. This ratio is determined here to be 1.0289 (19), so that it is definitely greater than 1; this was not completely clear from earlier calculations.
Our results for the time-moments of the charmonium vector current-current correlators also provide a new level of accuracy for these quantities, improving by a factor of 10 over the first such calculations in [29]. Our results are given in Table XII. We see some tension for the lowest (4th) moment with phenomenological results derived from R(e + e − → hadrons) in [67] when particular selections of experimental datasets are made.
We use the time-moments to derive the connected c quark HVP contribution to the anomalous magnetic moment of the muon, a c µ . Although this is not a large part of the total HVP contribution and so improving its uncertainty has little impact on the full HVP contribution, nevertheless it is a piece that can be calculated very accurately in lattice QCD and provides a test case for comparison of lattice calculations and a comparison with phenomenology.
Our result for a c µ improves the accuracy by a factor of 3 over earlier lattice QCD values [25,71,72]. Comparison of our result for a c µ to that determined from phenomenology can be divided into contributions from narrow resonances and from the continuum e + e − → cc. The contribution from the J/ψ agrees well between our result and phenomenology, with our result being more accurate (see Section VI D). This reflects the situation described above for f J/ψ and Γ . The total a c µ derived from the time-moments determined from experimental data on R(e + e − → hadrons) when the component from a cc loop is separated out shows some tension with our results, depending on which experimental datasets are used above 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 determinations 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 (repeating Eq. (41)) δa c µ = +0.0313(28) × 10 −10 .
This is a 0.2% effect and dominated by the impact of retuning 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 the full HVP contribution to a µ . We have shown that this is possible, for a small piece of the HVP, at least.   (11) there, we fix the α s coefficient to its known perturbative value of -0.1164(3). Fig. 30 shows the Z V data from [16] as hexagons, coloured according to their µ value. The fit lines for 2 and 3 GeV are also shown. The fit has a χ 2 /dof of 0.93. We note that the results from set 14 are included in this fit. The results from set 14, however, also agree with the fit result when they are not included in the fit. As discussed in Section IV µ was slightly mistuned on set 14; the true values are 2.04 and 2.98 GeV rather than 2 and 3 GeV. At this small lattice spacing the variation in Z V with µ (which is a discretisation effect) is small enough that this small mistuning can be neglected.
Note that the underlying perturbation theory of Eq. (A1) should agree with that obtained from the fit in [74] and it does. The discretisation effects are different in the two cases, of course. For RI-SMOM there is a momentum scale µ which we take as 2 GeV (although it can be taken as much smaller for Z V [16]) and this sets the size of discretisation effects as is clear from Figure 30.
Using the fit of Eq. (A1) we may also extract Z V values at finer lattice spacings where it may not be practical for us to perform direct calculations due to the computational cost of Landau gauge fixing. This includes the exafine lattices of set 15 in Table I with a lattice spacing adjusted for sea-mass mistuning of 0.032 fm. The value of 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 the lattice to the RI-SMOM scheme, Z SMOM m (µ), determined on ultrafine set 14 are given in Table XVI at µ = 2 GeV and 3 GeV. These are calculated in the same way as that discussed in [44].