QED interaction effects on heavy meson masses from lattice QCD+QED

Hadron masses are subject to few MeV corrections arising from QED interactions, almost entirely arising from the electric charge of the valence quarks. The QED effects include both self-energy contributions and interactions between the valence quarks/anti-quarks. By combining results from different signs of the valence quark electric charge we are able to isolate the interaction term which is dominated by the Coulomb piece, $\langle \alpha_{\mathrm{QED}}e_{q_1}e_{\overline{q}_2}/r \rangle$, in the nonrelativistic limit. We study this for $D_s$, $\eta_c$ and $J/\psi$ mesons, working in lattice QCD plus quenched QED. We use gluon field configurations that include up, down, strange and charm quarks in the sea at multiple values of the lattice spacing. Our results, including also values for mesons with quarks heavier than charm, can be used to improve phenomenological models for the QED contributions. The QED interaction term carries information about meson structure; we derive effective sizes $\langle 1/r_{\mathrm{eff}} \rangle^{-1}$ for $\eta_c$, $J/\psi$ and $D_s$ of 0.206(8) fm, 0.321(14) fm and 0.307(31) fm respectively.


I. INTRODUCTION
Lattice QCD calculations can now achieve a very high level of accuracy for ground-state meson masses. For example, a recent calculation of the mass splitting between the J/ψ and η c achieved an accuracy of 1 MeV [1]. This precision requires that QED effects arising from the electric charge of the quarks be included in the calculation and this is now being widely done, with a variety of approaches [1][2][3][4][5][6].
The QED effects arise almost entirely from the electric charge of the valence quarks. To O(α QED ) we then expect the impact of QED on a meson made of quark q 1 and antiquark q 2 to take the form ∆M q1q 2 = Ae q1 e q 2 + Be 2 q2 + Ce 2 q1 , ignoring the much smaller effects from the electric charge of the sea quarks (suppressed by powers of α s and sea quark mass effects). Here e q1 and e q2 are the electric charges of quarks q 1 and q 2 in units of e, the magnitude of the charge on the electron. The last two terms are dominated by 'self-energy' shifts in the valence quark masses. These are unphysical because they amount purely to a renormalisation of the quark mass by QED. The first term, with coefficient A, is physical, however. It is dominated, for nonrelativistic quarks, by the Coulomb interaction between the valence quark and antiquark in the meson. This effect depends on the average separation of the quarks and so provides a measure for the size of the meson. Its accurate determination requires a calculation that fully controls the QCD effects that bind the quark and antiquark into the meson, i.e. the use of lattice QCD. * daniel.hatton@glasgow.ac.uk † christine.davies@glasgow.ac.uk ‡ URL: http://www.physics.gla.ac.uk/HPQCD We will use lattice QCD calculations to which we also add the effect of QED on the valence quarks in an approach known as 'quenched QED' [7]. This is simply achieved by generating a random photon field in momentum-space and then packaging the field in position space into a compact U(1) variable that can be multiplied into the gluon field as the Dirac equation is solved for each quark propagator. Since QCD is responsible for binding the quark and antiquark into the meson and the effect of QED is simply a perturbation to the meson mass, then the QED interaction term Ae q1 e q 2 in Eq. (1) can be isolated by comparing results from lattice calculations in which we flip the sign of the electric charge for one of the quarks. We have M (e q1 , e q 2 ) − M (−e q1 , e q 2 ) = 2Ae q1 e q 2 . (2) We focus here on studying this QED effect for relatively heavy mesons (η c , J/ψ and D s ) to test our understanding of the impact of QED. The reason for this is that the internal structure of these mesons is reasonably well understood and in the past we have made use of estimates of the Coulomb interaction effects to assess the impact of QED on these meson masses [8,9]. In Section II we describe our lattice calculation and the results and in Section III we compare to these earlier estimates from phenomenological models. Section IV gives our conclusions.

II. LATTICE QCD CALCULATION
We work on n f = 2 + 1 + 1 gluon field configurations generated by the MILC collaboration [12,13]. These configurations include the effect of u/d, s and c quarks in the sea using the Highly Improved Staggered Quark (HISQ) action [14]. Details of the parameters for the configurations are given in Table I.
In [1] we analysed charmonium correlators calculated in pure QCD and in QCD + quenched QED using these I: The parameters of the ensembles used in our calculation, numbered in column 1. Column 2 gives the QCD gauge coupling and column 3 the lattice spacing in units of the Wilson flow parameter, w0 [10]. The lattice spacing in fm is then given by using w0 = 0.1715 fm, fixed from fπ [11]. Ls and Lt are the lattice spatial and temporal extents in lattice units. Columns 6 and 7 give the sea light quark masses in lattice units, with the sea u and d quark masses taken to be the same and denoted l. Column 8 gives the valence s quark mass used in the Ds mesons. Columns 9 and 10 give the sea and valence c quark masses in lattice units, respectively. Not all sets are used for all calculations; * indicates that the set was used for charmonium, † that the set was used for Ds and ‡ that the set was used for valence masses of 2mc. Column 11 gives the corresponding number of configurations used from the set.   (3)) to the pure QCD mass. Column 4 gives the finite-volume correction needed on that gluon configuration set for the unphysical QED scenario (Eq. (4) for meson charge 4e/3). The uncertainty in ∆FV comes mainly from the uncertainty in the lattice spacing and does not include the systematic error from missing higher orders in 1/Ls (see text). Finally column 5 gives the extracted coefficient, Aη c or A J/ψ (Eq. (5)).   Table II, but now for heavyonium mesons using quarks with mass 2mc. Column 2 gives the ground-state 'η2c' (upper rows) and 'ψ2c' (lower rows) meson masses in lattice units in the pure QCD case for the gluon field configuration sets given in column 1. Column 3 gives the ratio of the mass difference for the physical and unphysical QED scenarios to the pure QCD mass. Column 4 gives the finitevolume correction needed on that gluon configuration set for the unphysical QED (Q=4e/3) scenario. Finally column 5 gives the extracted coefficient, Aη 2c or A ψ 2c . (and further sets) of gluon field configurations. This enabled us to determine accurately how the η c and J/ψ meson masses shift (for a fixed valence c quark mass) when the 2e/3 electric charge of the valence c quarks is included. The shifts are very small, upwards by ∼0.1%, but clearly visible. From this we could work out how the c quark mass should be retuned when QED is switched on. We chose the natural tuning procedure in which the c quark mass is adjusted in both QCD and QCD+QED until the J/ψ meson mass determined on the lattice agrees with experiment. This led us to a determination of the c quark mass in the MS scheme of m c (3 GeV) QCD+QED =0.9841(51) GeV. This value is then 0.2% lower than in pure QCD [1].
The reason that the inclusion of QED lowers the c quark mass (tuning to a fixed meson mass) is because the positive self-energy terms in Eq. (1) raise the meson mass. The Coulomb interaction is attractive inside a charmonium meson, however, and so must lower the meson mass. Here we set out to isolate the Coulombdominated piece of the QED effect.
As described in Section I we can do this by comparing two calculations in QCD+QED (which will be shorthand for QCD + quenched QED in what follows). One calculation is the normal QCD+QED charmonium calculation with c and c quarks with opposite electric charge. The second calculation is one in which the c quark electric charge is flipped but not that of the c. The difference between the two results then gives twice the QED interaction contribution to the meson mass (Eq. (2)).
Note that the second calculation is for an unphysical scenario as far as QED is concerned. The underlying QCD physics is the same in both cases. We use the same valence quark mass in the two calculations, i.e. a mass close to the tuned c quark mass in QCD+QED for the physical scenario. Our valence c quark masses are given in Table I. These are the same as masses used in [1], with tuning errors below 0.5%.
For the two calculations we combine c and c propagators to generate two-point correlation functions that we average over the gluon field configurations. We fit these as a function of time separation between source and sink to determine the ground-state masses in lattice units. The procedure for including quenched QED and for fitting correlation functions is exactly the same as that described in [1] and we do not repeat the discussion of either procedure here.
In Table II we give our results for the η c and J/ψ mesons. We calculate the ground-state masses in the pure QCD case and also in the physical and unphysical QCD + quenched QED scenarios. It is convenient to give the QCD+QED results for the masses as a ratio to the value in pure QCD. In Table II we therefore give values for the difference of the ratios in the physical and unphysical QCD+QED cases: e q1 and e q 2 are the electric charges of the quark and antiquark in units of e; for the charmonium case these are 2/3 and −2/3. Multiplying R by the mass in the pure QCD case, M (0, 0) (column 2 of Table II), then gives the mass difference needed for Eq. (2). One difference between the physical and unphysical QED scenarios that we must take into account, however, is that of finite-volume effects from QED. In the physical scenario the charmonium meson is electrically neutral and finite-volume effects are negligible, as demonstrated in [1]. In the unphysical scenario the meson has an electric charge of 4e/3 and QED finite-volume effects are much larger. The finite-volume effects have been calculated analytically as an inverse power series in the spatial extent of the lattice [2,15,16]. We need only the (universal) result up to 1/L 2 s which takes the form  1)) in the ηc (upper plot, blue symbols) and J/ψ (lower plot, purple symbols) meson masses, shown as a function of squared lattice spacing (given in units of the quark mass, denoted m h but here mc). The fit is described in the text and shown by the curves in each plot. The pink symbols show the same results, but for mesons made from a quark-antiquark pair with quark mass m h twice that of the c quark. The symbol shape denotes the gluon field configurations used and is the same as that for the matching charmonium calculation.
with κ = 2.8373 and Q the meson electric charge in units of e. The leading term, which is independent of meson mass, takes a value of Q 2 × 0.5 MeV on a 4 fm lattice. We see then that the finite-volume effects are small here, but not negligible compared to our QED shifts. We handle them by correcting our finite-volume masses using the formula above for the cases where we have an (unphysical) electrically charged charmonium meson. The finite-volume shifts in each case are given in Table II. Figure 1 plots ∆ FV for Q = 4/3 as a function of spatial lattice size for the range of lattice sizes that we use here. The plot compares the leading 1/L s term of Eq. (4) to the result of including both the 1/L s and 1/L 2 s terms. We also show the impact of next-to-next-to-leading-order (NNLO) terms at 1/L 3 s from [16]. We take the value of r 2 that appears in the 1/L 3 s terms from vector meson dominance as 6/M 2 J/ψ (since we have shown in [17] that vector dominance works well for the electromagnetic form factor of mesons at small momentum-transfer, including for the η c ). We estimate the systematic uncertainty from missing out the 1/L 3 s terms at 0.005 MeV, which is neg-ligible compared to other sources of uncertainty. We then combine the mass differences and finitevolume shifts to isolate the QED interaction effect for the η c and J/ψ (Eq. 2). The coefficient, A, is determined as: These values are given for each ensemble in Table II. We plot A ηc and A J/ψ in Fig. 2 as a function of lattice spacing. We see that, as expected, the attractive Coulomb interaction yields a negative contribution to the meson masses because A is positive and e q1 e q 2 is negative. The A values are not the same for the η c and J/ψ mesons because of the QED hyperfine interaction, which acts in the same direction as the QCD hyperfine interaction raising the J/ψ mass relative to the η c [1].
In order to obtain a value for the coefficient A in the continuum limit we use a fit that allows for discretisation errors as well as possible effects from the mistuning of the charm quark valence mass and the mistunings of the sea quark masses from their physical values. The fit form we use is similar to that in [1]: The mass mistuning terms here are defined as in [1]: The uncertainty is dominated by that from the extrapolation to zero lattice spacing and is much larger than that from possible systematic errors in the finite-volume correction discussed above. The fit is able to pin down the coefficient of the (am c ) 2 term (c (1) a ) to be within 0.3 of zero. This is consistent with the expectation that this coefficient should be of size O(α s ) [14].
The Coulomb interaction effect probes the internal structure of the meson at short distances between the quark-antiquark pair. It is therefore interesting to ask how the coefficient A changes for heavier quarks than the c quark. In Table III we give our results for a heavyonium meson made from a quark-antiquark pair with quark mass twice that of the c quark (but the same electric charge). Again we use these results to determine the coefficient A (which is independent of electric charge) in this case. These results are also plotted in Fig. 2. The coefficient A is substantially larger for the heavier mass case.
We perform fits to the heavier mass points also using the fit form of Eq. (6), but with am c now replaced with 2am c and dropping the a 6 terms because we have results on fewer ensembles for this case. The functional form of the lattice spacing dependence should be the same in the m c and 2m c cases up to possible dependence on the squared velocity of the heavy quark inside the boundstate in higher order coefficients in a 2 [14]. We therefore use the results of the m c fit as prior information to constrain the coefficients of the lattice spacing dependence (c The fits give χ 2 /dof of 0.39 and 0.09 respectively.
We will discuss a comparison of the values for A for charmonium and heavyonium with those determined from static QCD potentials in Section III.
We can contrast the heavyonium case with that of a heavy-light meson. The simplest meson to use for this case is the heavy-strange meson since this has no valence light quark. We carry out the same analysis for the D s as for charmonium, but now the QED finite-volume effects (Eq. (4)) apply to both the physical scenario (since the D s meson is electrically charged with Q=1) and the unphysical scenario (where the 'D s ' has the smaller charge Q = 1/3). In Eq. (5) we therefore substitute for ∆ FV the difference δ∆ FV of the finite-volume effects for the Q = 1 and Q = 1/3 cases. The valence s quark masses that we use are given in Table I and are those obtained from the m s tuning exercise in [19]. Our results for the D s meson mass in pure QCD, along with the ratio R of Eq. (3) and the finite-volume shifts discussed above, are given in Table IV. The results for the QED interaction coefficient A for this case are shown in Fig. 3. A Ds is positive, and combined with a positive product of electric charges gives, as expected, a positive shift to the meson mass because the Coulomb interaction inside an electrically charged meson is repulsive. We perform the same continuum extrapolation fit as for the charmonium case, with the same priors, IV: Results that we use to obtain the QED interaction effect for the Ds meson. Column 2 gives the ground-state Ds meson mass in lattice units in the pure QCD case for the gluon field configuration sets given in column 1. Column 3 gives the ratio of the mass difference for the physical and unphysical QED scenarios (Eq. (3)) to the pure QCD mass. Column 4 gives the finite-volume correction needed on that gluon configuration set for the difference between the physical QED scenario (with meson charge 1) and the unphysical QED scenario (with meson charge 1/3). Finally column 5 gives the extracted coefficient of the effect on the mass from the quark electric charge interaction term, AD s .  with a χ 2 /dof of 0.015. We can also, as in the heavyonium case, work with heavy quarks with mass 2m c . The results for this mass are given in Table V and also plotted in Fig. 3. In contrast to the heavyonium case, we find that the coefficient A hardly changes as we change the heavy quark mass in the heavy-strange meson to be 2m c . From a fit to this case we obtain A Ds,2c = 4.68(66) MeV (11) with χ 2 /dof of 0.002. This agrees very well with that for the D s above, consistent with the fact that the points are on top of each other in Fig. 3.

III. DISCUSSION
The coefficient A is a physical quantity, encoding information about meson structure. The quantitative information that lattice QCD results for A provide can be used to calibrate more qualitative model approaches for comparable quantities. We discuss this below first for heavyonium and then heavy-light mesons.
The language of potential models provides a reasonably good approximation for heavyonium. A simple Cornell potential [20] of the form can readily be tuned to give the radial excitation energy of charmonium with an accuracy of ∼10%. Here κ is 4α s /3 and b 2 is the inverse string tension. The parameters used are: κ = 0.52 and b = 2.34 GeV −1 , along with a c quark mass in the kinetic energy term of Schrödinger's equation of 1.84 GeV [20]. It is then straightforward to perturb the coefficient of 1/r in Eq. (12) by α QED to include the Coulomb interaction effect and determine a value for the ground-state energy shift which is the potential model value for A for charmonium, A potl c . Alternatively this can be obtained by integrating over α QED /r weighted by the square of the ground-state wavefunction.
Doing this gives a value for A potl c of 5.9 MeV (this is a shift of 2.6 MeV downwards in the meson mass when multiplied by e q1 e q 2 for charmonium [8]). This result is for the leading spin-independent central potential of Eq. (12) and does not include any spin-dependent effects. Our lattice QCD results, on the other hand, are for the η c and J/ψ mesons separately. To compare our lattice results to those from a spin-independent potential we need to take the spin average: Our results from Section II (Eq. (8)) yield The potential model result, A potl c , given above is 15(3)% larger than our lattice QCD value. The uncertainty here comes from the lattice QCD calculation where it can be quantified. Clearly more sophisticated potential models, including potentials derived from lattice QCD [21], could be used to improve on the potential model result. Our value for A c in Eq. (14) can also be used to tune the parameters of potential models. Frequently the tuning is done using quantities such as the wavefunction at the origin, along with the spectrum (see, for example, [22,23]). The wavefunction at the origin is not a physical quantity, however, and there are sizeable uncertainties associated with renormalising this to relate it to experimental decay rates. In contrast the quantity A is a physical, renormalisation-group invariant quantity that can be compared much more precisely. A systematic uncertainty of order 10% on A potl c might be expected on the potential model result from missing O(v 4 ) relativistic corrections. However this could be ameliorated by tuning the potential.
We can also compare the lattice and Cornell potential results for the heavier quark mass of 2m c . Then our lattice spin-averaged result, using the values from Eq. (9) is The result for A potl 2c from the same Cornell potential as for the m c case is 9.1 MeV, now 30% too large.
Our results show a variation of A with quark mass that behaves approximately as √ m. We can compare this to what might be expected from scaling arguments for a potential of the form Cr N . Then, as we change the quark's reduced mass, µ ≡ m/2, we obtain the same solution for a rescaled distance λr where [24,25] A √ µ behaviour for A potl ≡ α QED /r would then correspond to N ≈ 0. Such a form for the heavy quark potential is in fact a standard one that has been successful in obtaining spectra, either taking N to be a small value or taking V (r) to be logarithmic [26,27]. These forms for the potential give a wavefunction that does not grow so rapidly with mass at small distance as the Cornell potential and might give results for A potl c and A potl 2c in better agreement with our lattice QCD value. See [22,28] for a comparison of spectrum and wavefunction results for different potential forms.
The difference of our results for A ηc and A J/ψ gives the 'direct' QED effect on the charmonium hyperfine splitting, when multiplied by -4/9. Note that in [1] we also included a quark-line disconnected contribution from QED to the hyperfine splitting that is not included here. We have a difference of A J/ψ and A ηc of -2.5(3) MeV from Eq. (8) for the m c case and -2.4(8) MeV from Eq. (9). The QED effect is then to raise the vector mass with respect to the pseudoscalar by about 1 MeV in both cases (using electric charge 2/3).
The hyperfine splitting itself falls with increasing quark mass, so the relative QED effect (for the same electric charge) is growing. In [1] we did a complete analysis of the charmonium hyperfine splitting, including QED effects, as a function of lattice spacing and sea quark mass. To compare the m c and 2m c cases here it is sufficient to take results from a single ensemble, set 6, as a guide to variation with heavy quark mass. The results in Tables II and III then yield a pure QCD hyperfine splitting of 111 MeV for the m c case and 75 MeV for the 2m c case, i.e. a fall of 30% on doubling the quark mass. This is to be compared to a QED contribution that does not change (at the level of our uncertainties). A key difference between the QCD and QED hyperfine splittings is the effective inclusion (implicit in our lattice QCD calculation) of a running coupling constant in the QCD case which reduces the splitting as the mass increases.
The QED hyperfine effect can also be compared to the expectation from a potential model calculation, by determining the impact of the perturbation from the Coulomb term on the wavefunction at the origin. The leading term in the hyperfine splitting from a potential model is proportional to the square of the wavefunction at the origin and so the percentage change in the hyperfine splitting is simply twice the percentage shift in ψ(0). For the Cornell potential discussed above we find the percentage change in the hyperfine splitting (for e q = 1) to be 1.92% for the m c case and 2.74% for the 2m c case. This shows an increase in the percentage QED effect that grows with the quark mass, as we find from our lattice calculation. To compare more quantitatively to our results we multiply these percentages by the pure QCD hyperfine splitting on set 6 given above. This gives a QED hyperfine effect (for e q = 1) of 2.1 MeV for both cases, in good agreement with our results from the difference of A ηc and A J/ψ . Note, however, that sizeable (O(30%) for charmonium) systematic errors are to be expected in analyses of fine structure from a potential model, so a semi-quantitative comparison is the best we can do here.
We now turn to the heavy-light meson case. In [9] we analysed a model of QED and light-quark mass effects in heavy-light pseudoscalar meson masses to isolate the QED interaction term phenomenologically. We used [29] M (e q h , e q l , m q ) = M 0 + Ae q h e q l + Be 2 q l + C(m q l − m l ) (17) where e q h and e q l are the electric charges of the heavy quark and light antiquark respectively and m q l is the light quark mass, m l being the average u/d quark mass. The coefficient A gives the QED interaction term that we are interested in here. The coefficient B is that of the light-quark QED self-energy, assumed to be independent of light quark mass. No term is included for the heavy quark self-energy because that cancels in the differences of heavy-light meson masses for the same heavy quark that we will use to fix the coefficients. The coefficient C allows for linear dependence on the light quark mass, independent of QED effects. From heavy quark symmetry we can expect that A, B and C will be constant up to Λ/m h corrections as the heavy quark mass m h → ∞ and independent of m q l up to small chiral corrections. This model was also used in [13].
If we assume that the coefficient A ≡ A phen is indepen- This result agrees well with our lattice determination of A for the D s meson in Eq. (10). From our calculation of results with a heavy quark mass twice that of charm (Eq. (11)) we are able to show that indeed A is independent of heavy quark mass at the level of our uncertainties (∼10%). It would be straightforward to extend our calculation to the D meson from the D s to test for any dependence of A on the light quark mass.

IV. CONCLUSIONS
We have shown here how to separate out the QED interaction piece from the self-energy terms in the determination of the effect of including QED for valence quarks on heavyonium and heavy-light meson masses. Lattice QCD calculations are now accurate enough that the effect of QED, at least for the valence quarks, can have an impact. The full effect of QED needs to be included in order to tune parameters, such as quark masses, by tuning meson masses until they take their experimental value in the QCD+QED calculation (this was done, for example, for QCD + quenched QED in [1]).
There are multiple reasons for wanting to separate out the QED interaction piece from the self-energy terms of the QED effect, however. One is to test our understanding of the physical contribution of QED by comparing to phenomenological model calculations. Another is to use the effect as a probe of meson, and more generally, hadron structure by using it to determine an effective average radial separation of the valence quarks.
We have determined the coefficient A of the QED interaction piece for the η c , J/ψ and D s mesons as well as for the corresponding mesons constructed by doubling the c quark mass (see Eqs (8), (9), (10) and (11)). The uncertainties we obtain at the physical point are 5% for the heavyonium case and 10% for heavy-light. A simple potential model gives results for the Coulomb interaction effect in charmonium in reasonable agreement with the lattice QCD numbers (spin-averaged to remove spin effects). We suggest that the lattice results could be used to tune potential models more accurately. This in turn could improve results for calculations, for example involving excited states and hadronic decay channels, that are currently more readily done in a potential model than using lattice QCD. We also find that a phenomenological model based on heavy quark symmetry gives good agreement with our D s results. We are able to demonstrate in that case that A is independent of heavy quark mass.
Since A is dominated by the Coulomb interaction effect for heavy mesons we can define an effective size parameter 1/r eff by dividing our results for A by α QED . This gives values for η c , J/ψ and D s mesons of 1/ 1/r eff = 0.206(8) fm, η c (19) = 0.321 (14) fm, J/ψ = 0.307(31) fm, D s .
The η c result can be compared to the value for r 2 using r 2 = 6/M 2 J/ψ which is in reasonable agreement with our results for the electromagnetic form factor of the η c at small squared momentum-transfer, q 2 [17]. This would give r 2 = 0.156 fm. We also find that the size parameter from Eq. (19) falls for heavier heavyonium masses approximately as √ m but does not change at all for heavy-light mesons as the mass is increased.
We believe that this could be a useful approach to assessing the size of other hadrons because it requires only the calculation and fitting of correlated 2-point functions. The noncompact QED action is simply being used as a convenient way to probe the r-dependence so a larger value of α QED than the physical one can be used to increase the signal for the perturbation [7]. For these purposes it might also be easier to use a purely Coulomb photon on each timeslice of the lattice as the direct Fourier transfrom of 1/r. By giving electric charge to pairs of quarks in more complicated hadrons such as baryons, tetraquarks or pentaquarks it might be possible to distinguish diquark-like configurations where they occur. We plan to test this out.