QED interaction effects on heavy meson masses from lattice QCD+QED
QQED interaction effects on heavy meson masses from lattice QCD+QED
D. Hatton, ∗ C. T. H. Davies, † and G. P. Lepage (HPQCD collaboration) , ‡ SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK Laboratory for Elementary-Particle Physics, Cornell University, Ithaca, New York 14853, USA (Dated: September 17, 2020)Hadron masses are subject to few MeV corrections arising from QED interactions, almost entirelyarising from the electric charge of the valence quarks. The QED effects include both self-energycontributions and interactions between the valence quarks/anti-quarks. By combining results fromdifferent signs of the valence quark electric charge we are able to isolate the interaction term which isdominated by the Coulomb piece, (cid:104) α QED e q e q /r (cid:105) , in the nonrelativistic limit. We study this for D s , η c and J/ψ mesons, working in lattice QCD plus quenched QED. We use gluon field configurationsthat 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 toimprove phenomenological models for the QED contributions. The QED interaction term carriesinformation about meson structure; we derive effective sizes (cid:104) /r eff (cid:105) − for η c , J/ψ 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 highlevel of accuracy for ground-state meson masses. For ex-ample, a recent calculation of the mass splitting betweenthe
J/ψ and η c achieved an accuracy of 1 MeV [1]. Thisprecision requires that QED effects arising from the elec-tric charge of the quarks be included in the calculationand this is now being widely done, with a variety of ap-proaches [1–6].The QED effects arise almost entirely from the elec-tric charge of the valence quarks. To O ( α QED ) we thenexpect the impact of QED on a meson made of quark q and antiquark q to take the form∆ M q q = Ae q e q + Be q + Ce q , (1)ignoring the much smaller effects from the electric chargeof the sea quarks (suppressed by powers of α s and seaquark mass effects). Here e q and e q are the electriccharges of quarks q and q in units of e , the magnitude ofthe charge on the electron. The last two terms are domi-nated by ‘self-energy’ shifts in the valence quark masses.These are unphysical because they amount purely to arenormalisation of the quark mass by QED. The firstterm, with coefficient A , is physical, however. It is dom-inated, for nonrelativistic quarks, by the Coulomb inter-action between the valence quark and antiquark in themeson. This effect depends on the average separation ofthe quarks and so provides a measure for the size of themeson. Its accurate determination requires a calculationthat fully controls the QCD effects that bind the quarkand antiquark into the meson, i.e. the use of lattice QCD. ∗ [email protected] † [email protected] ‡ URL:
We will use lattice QCD calculations to which wealso add the effect of QED on the valence quarks in anapproach known as ‘quenched QED’ [7]. This is sim-ply achieved by generating a random photon field inmomentum-space and then packaging the field in posi-tion space into a compact U(1) variable that can be mul-tiplied into the gluon field as the Dirac equation is solvedfor each quark propagator. Since QCD is responsible forbinding the quark and antiquark into the meson and theeffect of QED is simply a perturbation to the meson mass,then the QED interaction term Ae q e q in Eq. (1) can beisolated by comparing results from lattice calculations inwhich we flip the sign of the electric charge for one of thequarks. We have M ( e q , e q ) − M ( − e q , e q ) = 2 Ae q e q . (2)We focus here on studying this QED effect for rela-tively heavy mesons ( η c , J/ψ and D s ) to test our un-derstanding of the impact of QED. The reason for thisis that the internal structure of these mesons is reason-ably well understood and in the past we have made useof estimates of the Coulomb interaction effects to assessthe impact of QED on these meson masses [8, 9]. In Sec-tion II we describe our lattice calculation and the resultsand in Section III we compare to these earlier estimatesfrom phenomenological models. Section IV gives our con-clusions. II. LATTICE QCD CALCULATION
We work on n f = 2 + 1 + 1 gluon field configurationsgenerated by the MILC collaboration [12, 13]. Theseconfigurations include the effect of u/d , s and c quarksin the sea using the Highly Improved Staggered Quark(HISQ) action [14]. Details of the parameters for theconfigurations are given in Table I.In [1] we analysed charmonium correlators calculatedin pure QCD and in QCD + quenched QED using these a r X i v : . [ h e p - l a t ] S e p TABLE I: The parameters of the ensembles used in our calculation, numbered in column 1. Column 2 gives the QCD gaugecoupling and column 3 the lattice spacing in units of the Wilson flow parameter, w [10]. The lattice spacing in fm is then givenby using w = 0.1715 fm, fixed from f π [11]. L s and L t are the lattice spatial and temporal extents in lattice units. Columns 6and 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 D s mesons. Columns 9 and 10 give the sea and valence c quark massesin lattice units, respectively. Not all sets are used for all calculations; * indicates that the set was used for charmonium, † thatthe set was used for D s and ‡ that the set was used for valence masses of 2 m c . Column 11 gives the corresponding number ofconfigurations used from the set.Set β w /a L s L t am sea l am sea s am val s am sea c am val c N cfgs ∗ † ∗† ∗ / 140 † ∗ †‡ † / 184 ‡ ∗ † ‡ ∗ / 87 † / 199 ‡ TABLE II: Results for the charmonium case. Column 2 givesthe ground-state η c (upper rows) and J/ψ (lower rows) mesonmasses in lattice units in the pure QCD case for the gluon fieldconfiguration sets given in column 1. Column 3 gives the ratioof the mass difference for the physical and unphysical QEDscenarios (see Eq. (3)) to the pure QCD mass. Column 4 givesthe finite-volume correction needed on that gluon configura-tion set for the unphysical QED scenario (Eq. (4) for mesoncharge 4 e /3). The uncertainty in ∆ FV comes mainly fromthe uncertainty in the lattice spacing and does not includethe systematic error from missing higher orders in 1 /L s (seetext). Finally column 5 gives the extracted coefficient, A η c or A J/ψ (Eq. (5)).Set aM QCD η c R η c ∆ FV [MeV] A η c [MeV]1 2.305364(39) -0.002080(39) -1.0308(54) 8.16(14)3 1.848041(35) -0.001806(25) -0.9600(51) 7.139(91)4 1.342455(21) -0.0017726(58) -0.8795(47) 6.944(42)6 0.896675(24) -0.001641(21) -1.3373(75) 7.020(80)Set aM QCD
J/ψ R J/ψ ∆ FV [MeV] A J/ψ [MeV]1 2.39308(14) -0.001342(23) -1.0295(54) 5.844(85)3 1.914749(67) -0.001144(12) -0.9589(51) 5.057(77)4 1.391390(43) -0.001063(17) -0.8785(47) 4.688(63)6 0.929860(54) -0.000883(25) -1.3352(75) 4.580(90)TABLE III: Results, as in Table II, but now for heavyoniummesons using quarks with mass 2 m c . Column 2 gives theground-state ‘ η c ’ (upper rows) and ‘ ψ c ’ (lower rows) mesonmasses in lattice units in the pure QCD case for the gluon fieldconfiguration sets given in column 1. Column 3 gives the ratioof the mass difference for the physical and unphysical QEDscenarios to the pure QCD mass. Column 4 gives the finite-volume correction needed on that gluon configuration set forthe unphysical QED (Q=4 e/
3) scenario. Finally column 5gives the extracted coefficient, A η c or A ψ c .Set aM QCD η c R η c ∆ FV [MeV] A η c [MeV]5 ‡ ‡ aM QCD ψ c R ψ c ∆ FV [MeV] A ψ c [MeV]5 ‡ ‡ (and further sets) of gluon field configurations. This en-abled us to determine accurately how the η c and J/ψ meson masses shift (for a fixed valence c quark mass)when the 2 e/ c quarksis included. The shifts are very small, upwards by ∼ c quark mass should be retuned when QED isswitched on. We chose the natural tuning procedure inwhich the c quark mass is adjusted in both QCD andQCD+QED until the J/ψ meson mass determined onthe lattice agrees with experiment. This led us to adetermination of the c quark mass in the MS schemeof m c (3 GeV) QCD+QED =0.9841(51) GeV. This value isthen 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 becausethe positive self-energy terms in Eq. (1) raise the me-son mass. The Coulomb interaction is attractive insidea charmonium meson, however, and so must lower themeson mass. Here we set out to isolate the Coulomb-dominated piece of the QED effect.As described in Section I we can do this by comparingtwo calculations in QCD+QED (which will be shorthandfor QCD + quenched QED in what follows). One calcu-lation is the normal QCD+QED charmonium calculationwith c and c quarks with opposite electric charge. Thesecond calculation is one in which the c quark electriccharge is flipped but not that of the c . The differencebetween the two results then gives twice the QED inter-action contribution to the meson mass (Eq. (2)).Note that the second calculation is for an unphysicalscenario as far as QED is concerned. The underlyingQCD physics is the same in both cases. We use the samevalence quark mass in the two calculations, i.e. a massclose to the tuned c quark mass in QCD+QED for thephysical scenario. Our valence c quark masses are givenin Table I. These are the same as masses used in [1], withtuning errors below 0.5%.For the two calculations we combine c and c propa-gators to generate two-point correlation functions thatwe average over the gluon field configurations. We fit .
00 3 .
25 3 .
50 3 .
75 4 .
00 4 . L s [fm] − . − . ∆ F V [ M e V ] LONLONNLO
FIG. 1: A plot to show the size of finite-volume shifts neededin the η c case (for the unphysical QED scenario with mesoncharge Q = 4 /
3) as a function of lattice spatial size. Theplot compares the leading-order 1 /L s calculation, which isindependent of meson mass, to the result of adding in higherorder terms in 1 /L s . these as a function of time separation between sourceand sink to determine the ground-state masses in latticeunits. The procedure for including quenched QED andfor fitting correlation functions is exactly the same asthat described in [1] and we do not repeat the discussionof either procedure here.In Table II we give our results for the η c and J/ψ mesons. We calculate the ground-state masses in thepure QCD case and also in the physical and unphysicalQCD + quenched QED scenarios. It is convenient togive the QCD+QED results for the masses as a ratio tothe value in pure QCD. In Table II we therefore givevalues for the difference of the ratios in the physical andunphysical QCD+QED cases: R = M ( e q , e q ) − M ( − e q , e q ) M (0 , . (3) e q and e q are the electric charges of the quark and an-tiquark in units of e ; for the charmonium case these are2 / − /
3. Multiplying R by the mass in the pureQCD case, M (0 ,
0) (column 2 of Table II), then gives themass difference needed for Eq. (2).One difference between the physical and unphysicalQED scenarios that we must take into account, however,is that of finite-volume effects from QED. In the physi-cal scenario the charmonium meson is electrically neutraland finite-volume effects are negligible, as demonstratedin [1]. In the unphysical scenario the meson has an elec-tric charge of 4 e/ /L s which takes the form∆ FV ( L s ) = M ( L s ) − M ( ∞ ) (4)= − Q α QED κ L s (cid:18) M L s (cid:19) . . . . . am h ) A η c [ M e V ] . . . . . am h ) A J / ψ [ M e V ] FIG. 2: The coefficient, A , of the QED interaction effect(Eq. (1)) in the η c (upper plot, blue symbols) and J/ψ (lowerplot, purple symbols) meson masses, shown as a function ofsquared lattice spacing (given in units of the quark mass, de-noted m h but here m c ). The fit is described in the text andshown by the curves in each plot. The pink symbols show thesame results, but for mesons made from a quark-antiquarkpair with quark mass m h twice that of the c quark. The sym-bol shape denotes the gluon field configurations used and isthe same as that for the matching charmonium calculation. with κ = 2.8373 and Q the meson electric charge in unitsof e . The leading term, which is independent of mesonmass, takes a value of Q × FV for Q = 4 / /L s term of Eq. (4)to the result of including both the 1 /L s and 1 /L s terms.We also show the impact of next-to-next-to-leading-order(NNLO) terms at 1 /L s from [16]. We take the value of (cid:104) r (cid:105) that appears in the 1 /L s terms from vector mesondominance as 6 /M J/ψ (since we have shown in [17] thatvector dominance works well for the electromagnetic formfactor of mesons at small momentum-transfer, includingfor the η c ). We estimate the systematic uncertainty frommissing out the 1 /L s terms at 0.005 MeV, which is neg-ligible compared to other sources of uncertainty.We then combine the mass differences and finite-volume shifts to isolate the QED interaction effect forthe η c and J/ψ (Eq. 2). The coefficient, A , is determinedas: A = 12 e q e q ( R × M (0 ,
0) + ∆ FV ) . (5)These values are given for each ensemble in Table II.We plot A η c and A J/ψ in Fig. 2 as a function of lat-tice spacing. We see that, as expected, the attractiveCoulomb interaction yields a negative contribution to themeson masses because A is positive and e q e q is nega-tive. The A values are not the same for the η c and J/ψ mesons because of the QED hyperfine interaction, whichacts in the same direction as the QCD hyperfine interac-tion raising the
J/ψ mass relative to the η c [1].In order to obtain a value for the coefficient A in thecontinuum limit we use a fit that allows for discretisationerrors as well as possible effects from the mistuning of thecharm quark valence mass and the mistunings of the seaquark masses from their physical values. The fit form weuse is similar to that in [1]: A ( a , δm ) = A (cid:34) (cid:88) i =1 c ( i ) a ( am c ) i + c m, sea δ sea ,udsm + c c, sea δ sea ,cm + c c, val δ val ,cm (cid:35) . (6)The mass mistuning terms here are defined as in [1]: δ sea ,udsm = 2 m sea l + m sea s − m phys l − m phys s m phys s , (7) δ sea ,cm = m sea c − m phys c m phys c ,δ val ,cm = M J/ψ − M expt J/ψ M expt J/ψ .M J/ψ is the lattice value in the QCD+QED case withthe physical QED scenario. For the experimental
J/ψ mass we use 3.0969 GeV [18]. We use priors of 0(1) forthe c ( i ) a , c m, sea and c c, val coefficients and a prior of 0 ± . c c, sea . The mistuning terms in the fit have verylittle effect but including them allows us to incorporateuncertainties from them in the final result.With χ / dof of 0.06 and 0.1 respectively we find A η c = 6 . A J/ψ = 4 . . The uncertainty is dominated by that from the extrap-olation to zero lattice spacing and is much larger thanthat from possible systematic errors in the finite-volumecorrection discussed above. The fit is able to pin downthe coefficient of the ( am c ) term ( c (1) a ) to be within 0.3 of zero. This is consistent with the expectation that thiscoefficient should be of size O ( α s ) [14].The Coulomb interaction effect probes the internalstructure of the meson at short distances between thequark-antiquark pair. It is therefore interesting to askhow the coefficient A changes for heavier quarks thanthe c quark. In Table III we give our results for a heavy-onium meson made from a quark-antiquark pair withquark mass twice that of the c quark (but the same elec-tric charge). Again we use these results to determine thecoefficient A (which is independent of electric charge) inthis case. These results are also plotted in Fig. 2. Thecoefficient A is substantially larger for the heavier masscase.We perform fits to the heavier mass points also usingthe fit form of Eq. (6), but with am c now replaced with2 am c and dropping the a terms because we have resultson fewer ensembles for this case. The functional formof the lattice spacing dependence should be the same inthe m c and 2 m c cases up to possible dependence on thesquared velocity of the heavy quark inside the bound-state in higher order coefficients in a [14]. We thereforeuse the results of the m c fit as prior information to con-strain the coefficients of the lattice spacing dependence( c (1) a and c (2) a ) in the fit for the 2 m c case. This amountsto choosing a prior width of 0.3 for the c (1) a coefficientand 0.7 for c (2) a . We find A η c = 8 . A ψ c = 6 . . The fits give χ / dof of 0.39 and 0.09 respectively.We will discuss a comparison of the values for A for charmonium and heavyonium with those determinedfrom static QCD potentials in Section III.We can contrast the heavyonium case with that of aheavy-light meson. The simplest meson to use for thiscase is the heavy-strange meson since this has no valencelight quark. We carry out the same analysis for the D s asfor 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 un-physical scenario (where the ‘ D s ’ has the smaller charge Q = 1 / FV the difference δ ∆ FV of the finite-volume effects for the Q = 1 and Q = 1 / s quark massesthat we use are given in Table I and are those obtainedfrom the m s tuning exercise in [19]. Our results for the D s meson mass in pure QCD, along with the ratio R ofEq. (3) and the finite-volume shifts discussed above, aregiven in Table IV.The results for the QED interaction coefficient A forthis case are shown in Fig. 3. A D s is positive, and com-bined with a positive product of electric charges gives, asexpected, a positive shift to the meson mass because theCoulomb interaction inside an electrically charged mesonis repulsive. We perform the same continuum extrapola-tion fit as for the charmonium case, with the same priors, TABLE IV: Results that we use to obtain the QED interac-tion effect for the D s meson. Column 2 gives the ground-state D s meson mass in lattice units in the pure QCD case for thegluon field configuration sets given in column 1. Column 3gives the ratio of the mass difference for the physical andunphysical QED scenarios (Eq. (3)) to the pure QCD mass.Column 4 gives the finite-volume correction needed on thatgluon configuration set for the difference between the physicalQED scenario (with meson charge 1) and the unphysical QEDscenario (with meson charge 1/3). Finally column 5 gives theextracted coefficient of the effect on the mass from the quarkelectric charge interaction term, A D s .Set aM QCD D s R D s δ ∆ FV [MeV] A D s [MeV]2 1.52428(16) 0.000921(39) 0.3917(24) 5.01(18)3 1.22386(17) 0.000820(29) 0.4880(30) 4.74(13)5 0.87740(10) 0.000891(47) 0.3345(20) 4.70(21)6 0.59203(22) 0.00075(12) 0.6839(46) 4.86(52)TABLE V: Results as for Table IV but now for a heavy quarkmass with value 2 m c , for which we denote the meson D s, c .Set aM QCD D s, c R D s, c δ ∆ FV [MeV] A D s, c [MeV]5 ‡ ‡ see Eq. (6). Our fit returns a value A D s = 4 . χ / dof of 0.015.We can also, as in the heavyonium case, work withheavy quarks with mass 2 m c . The results for this massare given in Table V and also plotted in Fig. 3. In con-trast to the heavyonium case, we find that the coefficient A hardly changes as we change the heavy quark mass in . . . . . am h ) A D s [ M e V ] FIG. 3: The coefficient, A , of the QED interaction effect(Eq. (1)) in the D s meson mass, shown with purple symbolsas a function of squared lattice spacing (in units of the heavyquark mass, here m c ). The fit is described in the text. Thepink symbols show the same results, but for mesons madefrom a heavy quark and strange antiquark with heavy quarkmass twice that of the c quark. Symbol shapes match thoseof the D s results on the same gluon field configurations. the heavy-strange meson to be 2 m c . From a fit to thiscase we obtain A D s, c = 4 . χ / dof of 0.002. This agrees very well with that forthe D s above, consistent with the fact that the points areon top of each other in Fig. 3. III. DISCUSSION
The coefficient A is a physical quantity, encoding in-formation about meson structure. The quantitative in-formation that lattice QCD results for A provide can beused to calibrate more qualitative model approaches forcomparable quantities. We discuss this below first forheavyonium and then heavy-light mesons.The language of potential models provides a reason-ably good approximation for heavyonium. A simple Cor-nell potential [20] of the form V ( r ) = − κr + rb (12)can readily be tuned to give the radial excitation energyof charmonium with an accuracy of ∼ κ is4 α s / b is the inverse string tension. The parame-ters used are: κ = 0 .
52 and b = 2 .
34 GeV − , along with a c quark mass in the kinetic energy term of Schr¨odinger’sequation of 1.84 GeV [20]. It is then straightforward toperturb the coefficient of 1 /r in Eq. (12) by α QED toinclude the Coulomb interaction effect and determine avalue for the ground-state energy shift which is the po-tential model value for A for charmonium, A potl c . Alter-natively 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 ashift of 2.6 MeV downwards in the meson mass when mul-tiplied by e q e q for charmonium [8]). This result is forthe leading spin-independent central potential of Eq. (12)and does not include any spin-dependent effects. Our lat-tice QCD results, on the other hand, are for the η c and J/ψ mesons separately. To compare our lattice results tothose from a spin-independent potential we need to takethe spin average: A c = A η c + 3 A J/ψ A c = 5 . A potl c , given above is15(3)% larger than our lattice QCD value. The un-certainty here comes from the lattice QCD calculationwhere it can be quantified. Clearly more sophisticatedpotential models, including potentials derived from lat-tice QCD [21], could be used to improve on the potentialmodel result. Our value for A c in Eq. (14) can also beused to tune the parameters of potential models. Fre-quently the tuning is done using quantities such as thewavefunction at the origin, along with the spectrum (see,for example, [22, 23]). The wavefunction at the origin isnot a physical quantity, however, and there are sizeableuncertainties associated with renormalising this to relateit to experimental decay rates. In contrast the quantity A is a physical, renormalisation-group invariant quantitythat can be compared much more precisely. A systematicuncertainty of order 10% on A potl c might be expected onthe potential model result from missing O ( v ) relativisticcorrections. However this could be ameliorated by tuningthe potential.We can also compare the lattice and Cornell potentialresults for the heavier quark mass of 2 m c . Then ourlattice spin-averaged result, using the values from Eq. (9)is A c = 6 . . (15)The result for A potl2 c from the same Cornell potential asfor the m c case is 9.1 MeV, now 30% too large.Our results show a variation of A with quark massthat behaves approximately as √ m . We can comparethis to what might be expected from scaling argumentsfor a potential of the form Cr N . Then, as we changethe quark’s reduced mass, µ ≡ m/
2, we obtain the samesolution for a rescaled distance λr where [24, 25] λ ∝ µ − / (2+ N ) . (16)A √ µ behaviour for A potl ≡ (cid:104) α QED /r (cid:105) would then corre-spond to N ≈
0. Such a form for the heavy quark poten-tial is in fact a standard one that has been successful inobtaining spectra, either taking N to be a small value ortaking V ( r ) to be logarithmic [26, 27]. These forms forthe potential give a wavefunction that does not grow sorapidly with mass at small distance as the Cornell poten-tial and might give results for A potl c and A potl2 c in betteragreement with our lattice QCD value. See [22, 28] fora comparison of spectrum and wavefunction results fordifferent potential forms.The difference of our results for A η c and A J/ψ gives the‘direct’ QED effect on the charmonium hyperfine split-ting, when multiplied by -4/9. Note that in [1] we also in-cluded a quark-line disconnected contribution from QEDto the hyperfine splitting that is not included here. Wehave a difference of A J/ψ and A η c of -2.5(3) MeV fromEq. (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 re-spect to the pseudoscalar by about 1 MeV in both cases(using electric charge 2 / m c and 2 m c cases here it is suf-ficient to take results from a single ensemble, set 6, as aguide to variation with heavy quark mass. The results inTables II and III then yield a pure QCD hyperfine split-ting of 111 MeV for the m c case and 75 MeV for the 2 m c case, i.e. a fall of 30% on doubling the quark mass. Thisis to be compared to a QED contribution that does notchange (at the level of our uncertainties). A key differ-ence between the QCD and QED hyperfine splittings isthe effective inclusion (implicit in our lattice QCD calcu-lation) of a running coupling constant in the QCD casewhich reduces the splitting as the mass increases.The QED hyperfine effect can also be compared to theexpectation from a potential model calculation, by deter-mining the impact of the perturbation from the Coulombterm on the wavefunction at the origin. The leading termin the hyperfine splitting from a potential model is pro-portional to the square of the wavefunction at the originand so the percentage change in the hyperfine splitting issimply twice the percentage shift in ψ (0). For the Cornellpotential discussed above we find the percentage changein the hyperfine splitting (for e q = 1) to be 1.92% for the m c case and 2.74% for the 2 m c case. This shows an in-crease in the percentage QED effect that grows with thequark mass, as we find from our lattice calculation. Tocompare more quantitatively to our results we multiplythese percentages by the pure QCD hyperfine splitting onset 6 given above. This gives a QED hyperfine effect (for e q = 1) of 2.1 MeV for both cases, in good agreementwith 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 finestructure from a potential model, so a semi-quantitativecomparison is the best we can do here.We now turn to the heavy-light meson case. In [9] weanalysed a model of QED and light-quark mass effectsin heavy-light pseudoscalar meson masses to isolate theQED interaction term phenomenologically. We used [29] M ( e q h , e q l , m q ) = M + Ae q h e q l + Be q l + C ( m q l − m l )(17)where e q h and e q l are the electric charges of the heavyquark and light antiquark respectively and m q l is thelight quark mass, m l being the average u/d quark mass.The coefficient A gives the QED interaction term thatwe are interested in here. The coefficient B is that of thelight-quark QED self-energy, assumed to be independentof light quark mass. No term is included for the heavyquark self-energy because that cancels in the differencesof heavy-light meson masses for the same heavy quarkthat we will use to fix the coefficients. The coefficient C allows for linear dependence on the light quark mass, in-dependent of QED effects. From heavy quark symmetrywe can expect that A , B and C will be constant up toΛ /m h corrections as the heavy quark mass m h → ∞ andindependent of m q l up to small chiral corrections. Thismodel was also used in [13].If we assume that the coefficient A ≡ A phen is indepen-dent of heavy quark mass (i.e. we ignore Λ /m h terms) wecan easily determine it from experimental information. Ifwe add the experimental mass difference of B + and B to the mass difference of D + and D [30] then the termswith coefficients C and B (if independent of heavy quarkmass) cancel out. We have23 A phen + 13 A phen = 4 . − . ,A phen = 4 . . (18)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 quarkmass twice that of charm (Eq. (11)) we are able to showthat indeed A is independent of heavy quark mass at thelevel of our uncertainties ( ∼ D meson fromthe D s to test for any dependence of A on the light quarkmass. IV. CONCLUSIONS
We have shown here how to separate out the QED in-teraction piece from the self-energy terms in the determi-nation of the effect of including QED for valence quarkson heavyonium and heavy-light meson masses. LatticeQCD calculations are now accurate enough that the ef-fect of QED, at least for the valence quarks, can havean impact. The full effect of QED needs to be includedin order to tune parameters, such as quark masses, bytuning meson masses until they take their experimentalvalue in the QCD+QED calculation (this was done, forexample, for QCD + quenched QED in [1]).There are multiple reasons for wanting to separate outthe QED interaction piece from the self-energy terms ofthe QED effect, however. One is to test our understand-ing of the physical contribution of QED by comparingto phenomenological model calculations. Another is touse the effect as a probe of meson, and more generally,hadron structure by using it to determine an effectiveaverage radial separation of the valence quarks.We have determined the coefficient A of the QED in-teraction piece for the η c , J/ψ and D s mesons as wellas for the corresponding mesons constructed by doublingthe c quark mass (see Eqs (8), (9), (10) and (11)). Theuncertainties we obtain at the physical point are 5% forthe heavyonium case and 10% for heavy-light.A simple potential model gives results for the Coulombinteraction effect in charmonium in reasonable agreementwith the lattice QCD numbers (spin-averaged to removespin effects). We suggest that the lattice results couldbe used to tune potential models more accurately. Thisin turn could improve results for calculations, for exam-ple 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 phenomeno-logical model based on heavy quark symmetry gives goodagreement with our D s results. We are able to demon-strate in that case that A is independent of heavy quarkmass.Since A is dominated by the Coulomb interaction effectfor heavy mesons we can define an effective size parame-ter (cid:104) /r eff (cid:105) by dividing our results for A by α QED . Thisgives values for η c , J/ψ and D s mesons of1 / (cid:104) /r eff (cid:105) = 0 . , η c (19)= 0 . , J/ψ = 0 . , D s . The η c result can be compared to the value for (cid:112) (cid:104) r (cid:105) using (cid:104) r (cid:105) = 6 /M J/ψ which is in reasonable agreementwith our results for the electromagnetic form factor ofthe η c at small squared momentum-transfer, q [17]. Thiswould give (cid:112) (cid:104) r (cid:105) = 0.156 fm. We also find that thesize parameter from Eq. (19) falls for heavier heavyoniummasses approximately as √ m but does not change at allfor heavy-light mesons as the mass is increased.We believe that this could be a useful approach to as-sessing the size of other hadrons because it requires onlythe calculation and fitting of correlated 2-point functions.The noncompact QED action is simply being used as aconvenient way to probe the r -dependence so a largervalue of α QED than the physical one can be used to in-crease the signal for the perturbation [7]. For these pur-poses it might also be easier to use a purely Coulomb pho-ton on each timeslice of the lattice as the direct Fouriertransfrom of 1 /r . By giving electric charge to pairs ofquarks in more complicated hadrons such as baryons,tetraquarks or pentaquarks it might be possible to distin-guish diquark-like configurations where they occur. Weplan to test this out. Acknowledgements
We are grateful to the MILC collaboration for the useof their configurations and QCD code. We adapted thisto include quenched QED. Computing was done on theCambridge service for Data Driven Discovery (CSD3),part of which is operated by the University of Cam-bridge Research Computing on behalf of the DIRACHPC Facility of the Science and Technology FacilitiesCouncil (STFC). The DIRAC component of CSD3 wasfunded by BEIS capital funding via STFC capital grantsST/P002307/1 and ST/R002452/1 and STFC operationsgrant ST/R00689X/1. DiRAC is part of the nationale-infrastructure. We are grateful to the CSD3 supportstaff for assistance. Funding for this work came fromthe UK Science and Technology Facilities Council grantsST/L000466/1 and ST/P000746/1 and from the NationalScience Foundation. [1] D. Hatton, C. Davies, B. Galloway, J. Koponen, G. Lep-age, and A. Lytle (HPQCD), (2020), arXiv:2005.01845[hep-lat] .[2] S. Borsanyi et al. , Science , 1452 (2015),arXiv:1406.4088 [hep-lat] .[3] D. Giusti, V. Lubicz, C. Tarantino, G. Martinelli, F. San-filippo, S. Simula, and N. Tantalo, Phys.Rev.D ,114504 (2017), arXiv:1704.06561 [hep-lat] .[4] P. Boyle, V. Glpers, J. Harrison, A. Jttner, C. Lehner,A. Portelli, and C. T. Sachrajda, JHEP , 153 (2017),arXiv:1706.05293 [hep-lat] .[5] S. Basak et al. (MILC), Phys. Rev. D , 034503 (2019),arXiv:1807.05556 [hep-lat] .[6] Z. Kordov, R. Horsley, Y. Nakamura, H. Perlt, P. Rakow,G. Schierholz, H. Stben, R. Young, and J. Zanotti(CSSM/QCDSF/UKQCD), Phys. Rev. D , 034517(2020), arXiv:1911.02186 [hep-lat] .[7] A. Duncan, E. Eichten, and H. Thacker, Phys. Rev. Lett. , 3894 (1996), arXiv:hep-lat/9602005 .[8] C. Davies, E. Follana, I. Kendall, G. Lepage, andC. McNeile (HPQCD), Phys. Rev. D , 034506 (2010),arXiv:0910.1229 [hep-lat] .[9] C. Davies, C. McNeile, E. Follana, G. Lepage, H. Na,and J. Shigemitsu, Phys. Rev. D , 114504 (2010),arXiv:1008.4018 [hep-lat] .[10] S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, et al. , JHEP , 010 (2012), arXiv:1203.4469 [hep-lat].[11] R. Dowdall, C. Davies, G. Lepage, and C. Mc-Neile (HPQCD), Phys.Rev. D88 , 074504 (2013),arXiv:1303.1670 [hep-lat] .[12] A. Bazavov et al. (MILC), Phys. Rev.
D87 , 054505(2013), arXiv:1212.4768 [hep-lat] .[13] A. Bazavov et al. , Phys. Rev.
D98 , 074512 (2018),arXiv:1712.09262 [hep-lat] .[14] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong(HPQCD, UKQCD), Phys. Rev.
D75 , 054502 (2007),arXiv:hep-lat/0610092 [hep-lat] .[15] M. Hayakawa and S. Uno, Prog. Theor. Phys. , 413(2008), arXiv:0804.2044 [hep-ph] .[16] Z. Davoudi and M. J. Savage, Phys. Rev.
D90 , 054503(2014), arXiv:1402.6741 [hep-lat] .[17] C. Davies, J. Koponen, P. G. Lepage, A. T. Lytle,and A. C. Zimermmane-Santos (HPQCD), PoS
LAT-TICE2018 , 298 (2018), arXiv:1902.03808 [hep-lat] .[18] M. Tanabashi et al. (Particle Data Group), Phys. Rev.D , 030001 (2018).[19] B. Chakraborty, C. T. H. Davies, B. Galloway, P. Knecht,J. Koponen, G. C. Donald, R. J. Dowdall, G. P. Lep-age, and C. McNeile, Phys. Rev. D91 , 054508 (2015),arXiv:1408.4169 [hep-lat] .[20] E. Eichten, K. Gottfried, T. Kinoshita, K. Lane, andT.-M. Yan, Phys. Rev. D , 203 (1980).[21] G. S. Bali, Phys. Rept. , 1 (2001), arXiv:hep-ph/0001312 .[22] E. J. Eichten and C. Quigg, Phys. Rev. D , 1726(1995), arXiv:hep-ph/9503356 .[23] E. J. Eichten and C. Quigg, (2019), arXiv:1904.11542[hep-ph] .[24] C. Quigg and J. L. Rosner, Phys. Rept. , 167 (1979).[25] C. Davies, Lect. Notes Phys. , 1 (1998), arXiv:hep-ph/9710394 .[26] A. Martin, Phys. Lett. B , 338 (1980).[27] C. Quigg and J. L. Rosner, Phys. Lett. B , 153 (1977).[28] E. J. Eichten and C. Quigg, Phys. Rev. D , 5845(1994), arXiv:hep-ph/9402210 .[29] J. L. Goity and C. P. Jayalath, Phys. Lett. B , 22(2007), arXiv:hep-ph/0701245 .[30] P. Zyla et al. (Particle Data Group), PTEP2020