Charmonium properties from lattice $\mathrm{QCD}+\text{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 Highly Improved Staggered Quark (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}$ $(\text{connected})=120.3(1.1)\mathrm{MeV}$ and interpret the difference with experiment as the impact on $M_{{\eta }_c}$ of the ${\eta }_c$ decay to gluons, missing from the lattice calculation. This allows us to determine $\mathrm{\Delta }{M}_{{\eta }_c}^{\text{annihiln}}=+7.3(1.2)\mathrm{MeV}$, giving its value for the first time. Our result of $f_{J/\psi }=0.4104(17)\mathrm{GeV}$ gives $\mathrm{\Gamma }(J/\psi \to e^{+}{e}^{-})=5.637(49)\mathrm{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)\mathrm{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 hadronic vacuum polarization (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}^{-}\to \text{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}$.

Figure 1

A test of the impact of strong-isospin breaking in the sea through a comparison of results for the charmonium quantities given in Table 2 for sets 3A and 3B (Table 1). These sets differ only in the $u/d$ quark masses in the sea. Set 3A has $m_u=m_d$ and set 3B has $m_d/m_u=2.18$ with the same average value. The plot above gives the ratio of results in lattice units, for the same $a{m}_c^{\mathrm{val}}$, for the $J/\psi$ and ${\eta }_c$ masses and decay constants. The masses agree to within their 0.001% statistical errors and the decay constants agree to within their statistical errors (0.05% for the $J/\psi$).

Figure 2

Top panel: Pure QCD $J/\psi$ masses on all 15 sets from Table 1 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 $w_0$ which is correlated. The outer error bar shows the full uncertainty and the inner bar the uncertainty without the contribution from $w_0$. Bottom panel: The $J/\psi$ 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 $\mathrm{QCD}+\mathrm{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 $w_0/a$ for the QCD and $\mathrm{QCD}+\mathrm{QED}$ results on each set. The $\mathrm{QCD}+\mathrm{QED}$ results are above the pure QCD results in every case.

Figure 3

The charmonium hyperfine splitting as a function of lattice spacing on pure QCD ensembles. Our results are from ensembles including $u$, $d$, $s$ and $c$ quarks in the sea with varying $u/d$ average quark mass. The raw lattice results, for well-tuned $c$-quark masses, are given by symbols with error bars. The error bars include both statistical uncertainties and those from determining the lattice spacing from $w_0$, which are correlated between values. The results on each group of ensembles with approximately the same lattice spacing are given the same symbol. Within these groups, the results go from right to left as the $u/d$ quark mass changes from $m_s/5$ to the physical value. Notice the small range on the $y$ axis; this is the result of discretization effects being so small for the HISQ action. Results at mistuned $c$ masses are not plotted but are included in the fit. The fit line is the output of the fit from Eq. (5) at physical quark masses and with $Q=0$. The orange cross gives our result in the continuum limit for physical quark masses. The black cross gives the experimental average result [1].

Figure 4

The fractional effect of quenched QED on the charmonium hyperfine splitting plotted against $(a{m}_c{)}^2$. The fractional effect is determined at the same $a{m}_c$ value, i.e., it does not include $c$-quark mass retuning effects. The dashed line is horizontal and shows the weighted average value. The results show the precision that can be obtained by capitalizing on the correlations between $\mathrm{QCD}+\mathrm{QED}$ and QCD. This enables a clear demonstration that the impact of quenched QED here is not dependent on the lattice spacing.

Figure 5

The fractional shift in the $J/\psi$ and ${\eta }_c$ masses from the inclusion of quenched QED plotted as a function of $1/L_s$ on sets 5, 6 and 7. The dashed lines are horizontal lines at the weighted average values.

Figure 6

The charmonium hyperfine splitting as a function of lattice spacing, including both $\mathrm{QCD}+\mathrm{QED}$ and pure QCD points. The red open triangles are the same lattice results as in Fig. 3. The additional $\mathrm{QCD}+\mathrm{QED}$ points are given as blue open hexagons. The green fit band is the output of the fit from Eq. (5), but now with the $c$ quark electric charge, $Q=2/3$. The orange cross and orange band gives our result in the continuum limit for physical quark masses. The black cross and black band gives the experimental average result [1].

Figure 7

A stability analysis for our fit (see text) to the charmonium hyperfine splitting. Our full fit result is shown at the top. The lines with error bars below this show the result for different modifications of this fit. From the top, these include missing out results on particular sets of configurations, then missing out the results using mistuned $c$-quark masses and then changing the priors by a factor of 10 on different sets of parameters in the fit.

Figure 8

Comparison of different experimental results for the $J/\psi$ and ${\eta }_c$ masses along with the PDG average values. The ${\eta }_c$ results represent a recent subset of those used in the PDG average. The most recent result is from BELLE (denoted BELL18) [36]. There are three different determinations from LHCb [34, 35, 37], two of which also measured the hyperfine splitting. We include a KEDR measurement [33], two from different BABAR analyses [38] and two from BESIII [39, 40].

Figure 9

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/\psi$ and ${\eta }_c$ masses (see Fig. 8) rather than only from experiments that directly measure the splitting. The squares give lattice QCD results from calculations that include $n_f=2+1$ flavors of quarks in the sea. The hexagons give results that include $n_f=2+1+1$ flavors of sea quarks, including the results we present here at the top of the plot. All lattice QCD results have had uncertainties from neglecting ${\eta }_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 [41] and the previous HPQCD result [31] used HISQ quarks on the same asqtad-sea ensembles. Briceño et al. [42] 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 $\chi \mathrm{QCD}$ result [43] used overlap quarks on gluon field configurations including $n_f=2+1$ domain-wall sea quarks.

Figure 10

${\overline{m}}_c(3\mathrm{GeV})$ extrapolated to the continuum with a fit form that allows for condensate terms that behave like inverse powers of the renormalization scale $\mu$. This plot is an updated version of the upper section of Fig. 10 in [46], with added data points on ultrafine lattices (set 14) and a retuned bare $c$-quark mass fixed from the $J/\psi$, rather than ${\eta }_c$, meson. The ultrafine points have slightly mistuned $\mu$ values compared to the corresponding lines (see text).

Figure 11

The factor $C_m$ that forms part of the QED effect in the lattice to $\overline{\mathrm{MS}}$ mass renormalization, as defined in Eq. (15), plotted against the square of the lattice spacing. The fact that $C_m$, as shown here, has almost no $\mu$ or $a$ dependence demonstrates that the $\mu a$ dependence seen in $R_{\mathrm{QED}}[Z_m^{\overline{\mathrm{MS}}}]$ from columns 3 and 4 of Table 8 is simply that expected from perturbation theory.

Figure 12

QED correction to the charm quark mass in the $\overline{\mathrm{MS}}$ scheme at a scale of 3 GeV. The different $\mu$ values, shown as different colors and shapes, have all been run to 3 GeV and only differ by discretization and condensate effects. The red point on the left is the result for $R_{\mathrm{QED}}[{\overline{m}}_c(3\mathrm{GeV})]$ returned by the fit of Eq. (16).

Figure 13

Comparison of lattice QCD results for $m_c$ that include $u$, $d$, $s$ and $c$ quarks in the sea. The top two results are the ones from this paper. Our $\mathrm{QCD}+\text{quenched}$ QED result is given in Eq. (17). Our pure QCD result, Eq. (13), supersedes our earlier result in [46]. The Fermilab/MILC/TUM result is from [51] and uses a method based on charm-light meson masses. The “HPQCD HISQ JJc” result is from [26] and uses current-current correlator techniques. These three results agree to better than 1%. The ETMC result is from [52] and uses the RI-MOM intermediate scheme. The grey band gives the $\pm 1\sigma$ uncertainty band from the Particle Data Group [1].

Figure 14

The $J/\psi$ decay constant calculated on the ensembles of Table 1 in pure QCD and plotted against the square of the bare $c$-quark mass in lattice units. The different red shapes correspond to different groups of ensembles with similar lattice spacing and the error bars shown include the full uncertainty on the points. Points at mistuned $m_c$ are not plotted but are included in the fit. The green curve marks our extrapolation to the physical point, where the black cross shows the result determined from the experimental average for ${\mathrm{\Gamma }}_{e^{+}{e}^{-}}$ from Eq. (24).

Figure 15

The continuum extrapolation of $f_{J/\psi }$ using $Z_V(2\mathrm{GeV})$ (results and fit curve in red, same values as in Fig. 14) and $Z_V(3\mathrm{GeV})$ (results and fit curve in blue). The continuum extrapolated results agree in the two cases as they should. A black cross shows the experimental average result from Eq. (24).

Figure 16

The ${\eta }_c$ decay constant calculated on the ensembles of Table 1 in pure QCD and plotted against the square of the bare $c$-quark mass in lattice units. The different red shapes correspond to different groups of gluon field ensembles with similar lattice spacing. The error bars on each point are the full uncertainty, including correlated uncertainties from, for example, the determination of the lattice spacing. The green curve shows our fit and extrapolation to the physical point. The black cross gives the earlier HPQCD result on $n_f=2+1$ gluon field configurations from [62].

Figure 17

The ratio of $J/\psi$ to ${\eta }_c$ decay constant determined in pure QCD, plotted against the square of the bare $c$-quark mass in lattice units. The different red shapes correspond to different groups of gluon field ensembles with similar lattice spacing. The error bars on each point are the full uncertainty, including correlated uncertainties from, for example, the determination of the lattice spacing. The green curve shows our fit and extrapolation to the physical point.

Figure 18

The fractional QED correction to the $J/\psi$ and ${\eta }_c$ decay constants as a function of lattice spacing. The horizontal dashed lines mark the weighted average of the points.

Figure 19

The volume dependence of the fractional QED effect on the $J/\psi$ and ${\eta }_c$ decay constants measured on sets 5–7. The dashed lines are horizontal and indicate the weighted average of the points. There is no observable finite-volume effect at the level of our statistical uncertainties.

Figure 20

The $J/\psi$ decay constant calculated on the ensembles of Table 1 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, where the black cross shows the experimental average result from Eq. (24).

Figure 21

The ${\eta }_c$ decay constant calculated on the ensembles of Table 1 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 [62].

Figure 22

The ratio of $J/\psi$ to ${\eta }_c$ decay constants calculated on the ensembles of Table 1 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.

Figure 23

A comparison of our new results for $f_{{\eta }_c}$ and $f_{J/\psi }$ with earlier lattice QCD results, also by HPQCD, on gluon field configurations that include $n_f=2+1$ flavors of quarks in the sea. The results labeled “HPQCD HISQ” are from this paper and the results labeled “HPQCD HISQ $2010/12$” are from [31, 62]. The grey bands are the $\pm 1\sigma$ bands from our new $\mathrm{QCD}+\mathrm{QED}$ results.

Figure 24

A comparison of the width for $J/\psi$ decay to $e^{+}{e}^{-}$ implied by our new results for $f_{J/\psi }$ to that obtained from two recent experiments. The point labeled “KEDR 18” is from [54] and the point labeled “BES III 16” is from [55]. The grey bands are the $\pm 1\sigma$ bands from the Particle Data Group average [1].

Figure 25

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_{\mu }^c$ on sets 5–7. There is no observable dependence on the lattice spatial extent, $L_s$, as can be judged by comparison to the dashed horizontal lines at the weighted averages of the points.

Figure 26

The four lowest time moments and their extrapolation to $a=0$. The symbols give results for the rooted moment multiplied by $M_{J/\psi }$ with different symbols denoting different groups of ensembles with similar lattice spacing. The different colors pick out the different moments, from $n=4$ at the bottom to $n=10$ at the top of the plot. Only the results from the pure QCD calculation at well-tuned $c$-quark masses are shown for clarity. Uncertainties are too small to be visible. The extrapolation for all the moments are performed simultaneously including correlations between moments determined from the same current-current correlators. The dashed lines give the $\mathrm{QCD}+\mathrm{QED}$ fit curve using Eq. (5) and the colored crosses mark the result in the continuum limit at physical quark masses.

Figure 27

Comparison of our new result to those of previous lattice QCD calculations for the 4th and 6th time moments (appropriately rooted) of the charmonium vector current-current correlator. Our new result obtained here on $n_f=2+1+1$ gluon field configurations and including the effect of quenched QED is given at the top (blue cross). HPQCD’s 2012 result on $n_f=2+1$ gluon field configurations with HISQ valence $c$ quarks is marked “HPQCD HISQ 2012” [31]. JLQCD’s 2016 result using $n_f=2+1$ domain-wall 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 center-of-mass energy, by determining the $c\overline{c}$ component. The points plotted come from [63, 69]. Open orange circles show alternative selections of datasets from [69]; the upper value is for the “maximal” set and the lower value for the “minimal” set.

Figure 28

Extrapolation to the continuum physical point of the connected charm HVP contribution to the anomalous magnetic moment of the muon. Different symbols denote results on groups of ensembles with similar lattice spacing. Results at deliberately mistuned $c$-quark masses are not plotted but are included in the fit. The red points correspond to pure QCD, the light blue points to $\mathrm{QCD}+\mathrm{QED}$ and the dashed green fit curve plotted is that for $\mathrm{QCD}+\mathrm{QED}$. The continuum result (red cross) is compared to the result (open black square) obtained by calculating $a_{\mu }^c$ from the individually extrapolated time moments in Sec. 6a.

Figure 29

Comparison of lattice QCD results (not including QED) for the connected $c$ quark HVP contribution to $a_{\mu }$, $a_{\mu }^c$. Results are divided according to the number of sea quark flavors included in the gluon field configurations on which the calculation was done. The first result, labeled “HPQCD HISQ 2014” is from [26] using values of time moments determined in [31] using HISQ valence quarks on gluon field configurations including $n_f=2+1$ flavors of asqtad sea quarks. The other results all include $n_f=2+1+1$ flavors of sea quarks. The result labeled “ETMC twisted mass 2017” uses the twisted mass formalism [73] and that labeled “BMW stout stagg. 2017” a smeared staggered quark action [74]. Our new result [from Eq. (40)] labeled “HPQCD HISQ” agrees with, but is much more accurate than, these earlier results.

Figure 30

Comparison of the discretization effects in $a_{\mu }^c$ in our results from Table 15 (red open symbols) and those from the BMW Collaboration in [74] (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.

Figure 31

Comparison of our lattice $\mathrm{QCD}+\mathrm{QED}$ result (blue cross) for the connected $c$ quark HVP contribution to $a_{\mu }$, $a_{\mu }^c$, with results determined from experimental information. The red open squares, denoted “pheno.,” use the charmonium moments from $R_{e^{+}{e}^{-}}$ given in [63, 69] to determine $a_{\mu }^c$. The orange open circles give the alternative maximal (upper) and minimal (lower) inclusive cross-section datasets from [69].

Figure 32

Results for the local vector current renormalization factor in the RI-SMOM scheme compared to the fit form given in Eq. (a1). All of the data points shown are included in the fit but only the 2 and 3 GeV fit lines are drawn. The data points not on the fit lines correspond to $\mu =2.5\mathrm{GeV}$ (light orange) and 4 GeV (light green). The dashed line shows the value of $a$ on the exafine lattices that allows us to determine a $Z_V$ value there.