Precise determination of decay rates for <math display='inline'><mrow><msub><mrow><mi>η</mi></mrow><mrow><mi>c</mi></mrow></msub><mo stretchy='false'>→</mo><mi>γ</mi><mi>γ</mi></mrow></math> , <math display='inline'><mrow><mi>J</mi><mo>/</mo><mi>ψ</mi><mo stretchy='false'>→</mo><mi>γ</mi><msub><mrow><mi>η</mi></mrow><mrow><mi>c</mi></mrow></msub></mrow></math> , and <math display='inline'><mrow><mi>J</mi><mo>/</mo><mi>ψ</mi><mo stretchy='false'>→</mo><msub><mrow><mi>η</mi></mrow><mrow><mi>c</mi></mrow></msub><msup><mrow><mi>e</mi></mrow><mrow><mo>+</mo></mrow></msup><msup><mrow><mi>e</mi></mrow><mrow><mo>−</mo></mrow></msup></mrow></math> from lattice QCD

I. INTRODUCTION

Decay rates of mesons via annihilation to photons, or radiative transitions with emission of a photon, can in principle provide stringent tests of our understanding of the internal structure of these mesons from strong interaction physics. The strong interaction effects are parametrized by a decay constant or a form factor and these can be calculated from first principles using the techniques of lattice QCD. The decay rates are free from the uncertainties that arise in weak decay processes from sometimes poorly known Cabibbo-Kobayashi-Maskawa matrix elements. This means that the combination of accurate lattice QCD and experimental results can directly test both QCD and the Standard Model. An example is that of the leptonic decay rate of the J=ψ meson via a photon. The decay constant of the J=ψ was recently calculated with an uncertainty of 0.4% (including effects from the electric charges of the valence c quarks) giving a value for ΓðJ=ψ → e þ e -Þ accurate to 0.9% [1]. The lattice QCD result agrees well with the experimental average which has an uncertainty of 1.8%. Here we study two further processes of this kind for ground-state charmonium mesons, η c → γγ and J=ψ → γη c .

There are a few experimental results for the decay width for J=ψ → γη c [2][3][4][5]; the Particle Data Group (PDG) [6] gives a branching fraction of 1.7(4)% as an average of results from the Crystal Ball [2] and CLEO [3], with the uncertainty increased by a factor of 1.5 to allow for the tension between them. The average corresponds to a partial decay width, ΓðJ=ψ → γη c Þ, of 1.57 (37) keV.

Following early work in lattice QCD in the quenched approximation [7] and including u=d quarks in the sea [8,9], a result with a more realistic n f ¼ 2 þ 1 quark sea was obtained [10]. Despite the varying numbers of sea quarks the lattice QCD calculations consistently give values for the partial decay width that are higher than the PDG average [6] of the experimental results. Here we aim to shed further light on this issue by improving the accuracy from lattice QCD and including now u=d, s and c quarks in the sea.

For η c → γγ the experimental and theoretical picture is less clear. The PDG [6] combines multiple sets of products of branching fractions involving η c → γγ to obtain a fit value for that branching fraction with a 7% uncertainty [1.61ð12Þ × 10 -4 ]. Individual experimental results are typically much less accurate than this, however. The fitted branching fraction corresponds to a partial decay width of 5.15 (35) keV.

Lattice QCD calculations for η c → γγ also show an uncertain picture on the theoretical side. Early results in the quenched approximation [11] and subsequent results including u=d quarks in the sea [12,13] gave results for the decay rate with a central value much less than the PDG fit value above. Further recent results from lattice QCD including u=d sea quarks [14,15] give larger values for the decay rate in agreement with the PDG fit value [14] or exceeding it [15]. Here we significantly improve the theoretical understanding of this decay rate by performing the first 1%-accurate lattice QCD calculation of it, and we include a realistic sea quark content (u=d, s and c).

The accurate determination of Γðη c → γγÞ here, along with HPQCD's previous accurate determination of ΓðJ=ψ → e þ e -Þ [1], allows us to test the relationship between these two quantities expected at leading order (LO) in nonrelativistic QCD (NRQCD). The ratio of rates in NRQCD is [16] ΓðJ=ψ

Here Q c is the electric charge of the c quark in units of e. Such a simple formula is possible because the hadronic parameters, here the "wave function at the origin" cancel out at leading order. Sizeable radiative and relativistic corrections could be expected to this ratio but there is evidence in [16], calculating through Oðα 2 s Þ, that there is some cancellation between these corrections. Here we can determine the ratio of these two decay rates accurately and fully nonperturbatively, including the complete relativistic dynamics of the c quarks inside the mesons, using lattice QCD. This allows us to assess how closely the relationship of Eq. ( 1) is followed in full QCD.

For J=ψ → η c decay our lattice QCD calculation involves calculating a form factor as a function of the 4-momentum transfer, q, between parent and daughter mesons. The value of the form factor at q 2 ¼ 0 is the appropriate one for the radiative decay of a J=ψ to η c accompanied by a real photon. Here we calculate the form factor across the full q 2 range and so can also provide predictions for the case with an off-shell photon, i.e. J=ψ → η c e þ e -. The rate for the equivalent process for the D s meson, D  s → D s e þ e -has been measured experimentally [17]; these decays provide an additional test of our understanding of meson structure in QCD.

A further simple leading-order relationship between ΓðJ=ψ → γη c Þ, ΓðJ=ψ → e þ e -Þ and Γðη c → γγÞ was suggested many years ago by Shifman in [18], based on the approximation of J=ψ dominance of the vector cc current. He gives

By combining the J=ψ → η c γ and η c → γγ results we obtain here with HPQCD's earlier values for ΓðJ=ψ → e þ e -Þ [1], we can also test how well this works in full QCD.

We are able to perform these calculations accurately in lattice QCD because we use the highly improved staggered quark (HISQ) discretization of the Dirac equation [19]. The HISQ action has particularly small discretization effects and this means that c quark physics can be handled accurately in lattice QCD on lattices with moderate values of the lattice spacing [1,10]. This in turn means that a wide range of lattice spacings can be covered for accurate extrapolation to the continuum a → 0 limit. We use gluon field configurations generated by the MILC collaboration that include u=d, s and c quarks in the sea with lattice spacing values ranging from 0.15 to 0.06 fm.

The paper is laid out as follows. In Sec. II we discuss the calculation of the rate for η c → γγ using our lattice QCD determination of the amplitude. This includes first a discussion of the method, followed by a description of the lattice calculation with HISQ quarks and then a discussion of our results, including comparison to earlier lattice calculations and to experiment along with tests of expectations in the nonrelativistic limit. In Sec. III we follow the same path through the calculation of the rate for J=ψ → γη c and J=ψ → η c e þ e -. Section IV summarizes our results and gives our conclusions.

A. Method

To determine the decay rate of the η c to two photons, we need to calculate the matrix element between the η c and two on-shell photons (each with squared 4-momentum, q 2 ¼ 0). The procedure for the calculation is similar to that for a meson-to-meson transition form factor except that we must take a weighted integral over the time insertion point for one of the vector currents to fix the energy of the final state to that of the required photon.

We start from a standard lattice QCD 3-point function constructed from c quark propagators on each gluon field configuration and averaged over the ensemble, see Fig. 1,

project out 3-momenta q 1 and 0, respectively, and a is the lattice spacing. j μ and j ν are lattice vector currents, cγ μ c and cγ ν c, and O η c is a pseudoscalar or temporal axial current operator that couples to pseudoscalar charmonium states. We can couple an external photon field to jμ by multiplying by a photon propagator [20,21] so that

is the photon propagator in Euclidean time for a photon with 3-momentum q 1 , where ω 1 ¼ jq 1 j. The sum in Eq. ( 5) is over all t γ 1 on the lattice. The construction of Eq. ( 5) yields a 3-point function between a tower of η c states (at rest) created by O η c at t η c and a photon at t 0 induced by a vector current at t γ 1 . The 4-momentum transferred by the current is constrained by energy-momentum conservation. If we choose

where M η c is the η c mass, then for the ground-state η c this 3-point function encapsulates an η c → γγ transition with two real photons in the final state. The ground-state contribution to the 3-point function is

F μν is the matrix element between η c and γγ states that will allow us to determine the decay rate.

To obtain F μν from the 3-point function it is convenient to first peel off the final-state photon by dividing by the left-most factor in Eq. (8). Instead of Eq. ( 5) we construct in practice

At the same time we construct the standard 2-point function

By fitting Cμν and C η c simultaneously we can determine the contribution of the ground state η c to both 2-point functions and use this to obtain F μν ðη c → γγÞ.

The fit form for C η c can be written as

and the ground-state energy, corresponding to n ¼ 0, is

with [see Eq. ( 8)]

For n > 0, the b n correspond to matrix elements between excited η c states and one on-shell and one off-shell photon. Using parity and Lorentz invariance we can define a transition form factor Fðq 2 1 ; q 2 2 Þ byfoot_0

This equation makes clear that a nonzero result will only be obtained (for an η c at rest) if there is a component of the (equal-and-opposite) spatial momentum of the two photons that is orthogonal to the polarization vectors of both photons (which must themselves be orthogonal). The specific configurations of momenta and polarizations that we use will be described below. Here we are interested in determining Fð0; 0Þ (i.e. with two on-shell photons) and so the kinematic factors in Eq. ( 16) mean that F is obtained from b 0 via

The decay amplitude is given by

allowing for interchange of the two photons and inserting electric charge factors (Q c ¼ 2=3 for the c quark) and polarization vectors. The decay rate is then

The factor of 1=2 above avoids double counting the identical photons and there are two spin combinations for the final state, both with jε 1 × ε 2 j ¼ 1.

In the next section we give more details of how we set up our lattice calculation to determine Fð0; 0Þ.

Ensembles and parameters

We use ensembles with N f ¼ 2 þ 1 þ 1 flavors of dynamical HISQ sea quarks from the MILC collaboration [22,23]. Details of the ensembles are tabulated in Table I. The six ensembles that we use give us four different lattices spacings, a ≈ 0.06, 0.09, 0.12 and 0.15 fm. The sea u and d quark masses are taken to be the same and denoted m l (l for light). Two ensembles, sets 2 and 4, have l quarks in the sea with the physical mass, m l ¼ ðm u þ m d Þ=2, and the other four have heavier-than-physical light quarks. This allows us to test the dependence of our results on the mass of the sea light quarks.

On these ensembles we calculate correlation functions constructed from valence c quark propagators, also using the HISQ formalism. The valence c quark mass values that we use here are the same as those given in Table I of [1]. These masses are more accurately tuned than those of the c quark in the sea. The tuning is done by comparing the result for the mass of the J=ψ meson obtained on each ensemble to its physical mass from experiment, as discussed in [1].

TABLE I. Parameters for the MILC ensembles of gluon field configurations. The sets are numbered but also given a "handle" to distinguish them in the second column. The lattice spacing is determined from the Wilson flow parameter, w 0 [24], and values of w 0 =a are given in column 4, following the gauge coupling, β, in column 3. The physical value w 0 ¼ 0.1715ð9Þ fm was fixed from f π in [25]. Sets 1 and 2, 3 and 4, 5, and 6 have a ≈ 0.15; 0.12; 0.09; 0.06 fm, respectively. The number of lattice points in space, N x , and time, N t , are given in column 5. Sea quark masses in lattice units are given in columns 6, 7 and 8. All the configuration sets have equal-mass u and d quarks with m u ¼ m d ¼ m l . Sets 1, 3, 5 and 6 have heavier-than-physical mass m l such that m s =m l ¼ 5 and sets 2 and 4 have m l close to the physical average of u and d quark masses. As described in the text, we use valence c quark masses, given in column 9, that differ from the sea c quark masses, being more closely tuned to the physical value. The ϵ Naik parameter that accompanies am val c in the HISQ action [19,26] are given in the next column. On set 3 we also include results from a deliberately mistuned valence c quark mass of 0.654, and denote this calculation as "3A." We use 1000 gluon field configurations from each set, with two time sources on each configuration to increase statistics, except for set 6 where we include only one time-source per configuration.

Set

Handle β w 0 =a N

3 x × N t am sea l am sea s am sea c am val c ϵ Naik 1 Very-coarse 5.80 1.1119(10) 16 3 × 48 0.013 0.065 0.838 0.888 -0.3820 2 Very-coarse-physical 5.80 1.1367(5) 32 3 × 48 0.00235 0.0647 0.831 0.863 -0.3670 3 Coarse 6.00 1.3826(11) 24 3 × 64 0.0102 0.0509 0.635 0.664 -0.2460 3A Coarse 6.00 1.3826(11) 24 3 × 64 0.0102 0.0509 0.635 0.654 -0.2402 4 Coarse-physical 6.00 1.4149(6) 48 3 × 64 0.00184 0.0507 0.628 0.643 -0.2336 5 Fine 6.30 1.9006(20) 32 3 × 96 0.0074 0.037 0.440 0.450 -0.1250 6 Superfine 6.72 2.8941(48) 48 3 × 144 0.0048 0.0240 0.286 0.274 -0.0491 2. Correlation functions C μν and C η c

On the ensembles of Table I we calculate 2-point and 3-point correlation functions, as described in Sec. II A. The 3-point correlation function is that between a source current that couples to the η c and its excitations and two vector currents, depicted in Fig. 1. Note that we construct only one combination of quark propagators on the lattice and determine F from that; the factor of 2 for interchanging the photons is allowed for in Eq. (18). We include only the quark-line connected diagram of Fig. 1. There are quarkline disconnected diagrams (quark loops connected only by gluons) that contribute but we expect their contribution to be small because the c quark mass is large. We will estimate the systematic uncertainty from neglecting quark-line disconnected diagrams in Sec. II B 7.

The correlation function C μν ðt γ 1 ; t γ 2 ; t η c Þ [Eq. ( 3)] is calculated using the standard sequential source technique. On timeslice t η c , we set up a random wall source [27] from which two c quark propagators are calculated. The first propagator, evaluated at timeslice t γ 1 and multiplied by γ μ , is used as the source vector for another c quark propagator calculation. The result of this calculation, the extended/ sequential propagator, is contracted with the second c quark propagator from t η c , inserting a γ ν before summing over indices, to obtain the 3-point function. The second c propagator from t η c is calculated with an additional γ t γ 5 multiplying the random wall source to achieve the quantum numbers of the η c when contracted. The 3-point correlation functions obtained are averaged over all gluon fields in the ensemble. We use two time sources for t η c on each configuration (except for set 6) to improve statistical accuracy. We take μ; ν ¼ z, x and take the photon momentum (discussed below) to be in the y direction to satisfy the requirements for a nonzero result from Eq. ( 16).

Because we are using HISQ valence c quarks, the γ matrices in the paragraph above become spatial-positiondependent phases with which the sources and sinks are patterned to achieve the required spin. We also have to consider operator "taste" [19], also represented by γ matrices, and we must choose tastes for the different operators such that the product of the taste γ matrices gives 1, for a nonzero correlation function. We want the spin-taste representation of the two vector currents to be the same (except for z ↔ x) for symmetry, and we want to avoid point splitting of the vector currents along the y direction in which the spatial momentum flows. Our preferred setup uses a temporal axial current for O η c with spin taste in the standard notation (see, for example, [19]) γ 5 γ t ⊗ γ x γ z and vector currents with spin-taste γ x ⊗ γ x and γ z ⊗ γ z . In this case the two vector currents are local and this has the advantage that they have no tree-level discretization errors. The temporal axial current is one-link point split in the y direction. Use of the temporal axial current (rather than the pseudoscalar current) avoids any temporal point splitting.

We will call this the "LOCAL" setup because of the nature of the vector currents used. An alternative, that we will use as a test on a subset of ensembles (sets 1, 3 and 5) is to take all of the operators to be "tasteless" i.e. with spin-tastes γ 5 γ t ⊗ 1, γ x ⊗ 1 and γ z ⊗ 1. Now the vector currents have a one-link point splitting along the x and z directions, respectively, and the temporal axial current has a 3-link point splitting (from one corner to the opposite of a cube). This is the "ONE-LINK" setup.

Calculating HISQ c quark propagators is numerically relatively inexpensive and so we cover the full range of values of t γ 1 by calculating C μν [Eq. ( 3)] for t γ 1 from t η c to t η c -N t =2, obtaining the other values by periodicity. This allows us to test the behavior of the integral/sum over t γ 1 in Eq. (9). Results at all values of t γ 2 are obtained in the final contraction of the sequential propagator with that from t η c . This amount of calculation is not in fact necessary, as we discuss below.

In η c → γγ decay each photon carries away spatial momentum with magnitude M η c =2 in the η c rest frame. In our calculation spatial momentum is inserted into the c quark propagator that connects the operators at timeslices t γ 1 and t γ 2 (the sequential propagator, see Fig. 1) using twisted boundary conditions [28,29]. This enables us to choose any value of the spatial momentum, q 1 , using a twist angle θ. One photon will have momentum q 1 and the other -q 1 . Given our vector current polarizations, q 1 must have a y component and we simply take q 1 to be in the y direction. Its value is then related to θ by

where N x is the number of lattice points in a spatial direction. The kinematics that are specified in Eq. ( 7) then require a twist angle of

Since ω 1 ¼ jq 1 j and both of these quantities are set in lattice units, photon 1 will be exactly on shell. It is harder to arrange for photon 2 to be exactly on shell and there will inevitably be a slight mistuning of the on shell condition for this photon. Its spatial momentum has magnitude ajq 1 j and its energy is aM latt η caω 1 in lattice units, where M latt η c

is the mass of the η c obtained on that ensemble from the lattice calculation. Photon 2 will only be on shell if aω 1 is exactly aM latt η c =2. The mistuning depends on what value of M η c is used in determining the twist angle for q 1 . Various approaches that are equivalent in the continuum limit are possible. Here, for the LOCAL setup, we choose to fix the M η c value used in q 1 [Eq. (21)] to the value we obtain from connected correlation functions in lattice QCD in the physical continuum limit, since this corresponds to the value to which our η c masses will converge in that limit. This value was obtained by calculating the charmonium hyperfine splitting in [1] for both pure QCD and QCD plus quenched QED. Here we use the pure QCD result,

This mass differs slightly (by 5.6 MeV or 0.2% of the mass) from the experimental value [6]. The most likely explanation for this is that it represents the impact of quark-line disconnected correlation functions (not included in the lattice calculation) that allow the η c to mix with lighter flavor-singlet mesons such as the η and η 0 . The twists from Eqs. (20) and (22) used on each set are given in Table II (top section). For the ONE-LINK setup we chose to test a different tuning, equivalent in the continuum limit. In that case we determined aM latt η c for that taste of η c on that ensemble in a separate calculation and then used that value in our choice for θ [Eq. ( 21)], so that both photons 1 and 2 should be exactly on shell for each ensemble. The values for the twists used in that case are given in the lower section of Table II.

The accuracy with which photon 2 is tuned to the on-shell point will be discussed further in Sec. II B 5. The mis-tuning is small in both cases here (very small for the ONE-LINK case) but we will take account of it in our continuum/chiral extrapolation for Fð0; 0Þ. We will also estimate and include an error associated with the fact that M phys;latt η c is not equal to the experimental value. As well as C μν [Eq. (3)] we also calculate the 2-point correlation function C η c [Eq. (10)] constructed from the same O η c operator used in C μν . Fitting C η c and Cμν [derived from C μν , see Eq. ( 9)] simultaneously allows us to determine the form factor for η c → γγ.

Determining Cμν

The 2-point function Cμν ðt γ 2 ; t η c Þ is constructed from a weighted sum over all lattice time slices of C μν as given in Eq. ( 9), using aω 1 ¼ aq y 1 as discussed in Sec. II B 2. It is normalized so that the lattice vector currents match those in the continuum as discussed in Sec. II B 3. In Fig. 2 we plot a quantity proportional to the summand of Eq. ( 9) to illustrate how the sum works. The quantity plotted is

TABLE II. Twist θ used for the sequential propagator in C μν [Eq. (3) and Fig. 1] from which η c → γγ form factor is extracted. The corresponding 3-momentum component is given by aq y 1 ¼ θπ=N x and given in column 4. For the LOCAL setup (local vector currents, top table) the twist is chosen to achieve q y 1 ¼ M phys;latt η c =2 with M phys;latt η c given in Eq. (22). For the ONE-LINK setup (onelink vector currents, lower table) the twist was chosen so that aq y 1 ¼ aM latt η c =2, i.e. half the η c meson mass on that ensemble. LOCAL Set θ aq y 1 1 5.928 1.1640 2 11.598 1.1386 3 7.151 0.9361 4 13.976 0.9147 5 6.936 0.6809 6 6.833 0.4472 ONE-LINK Set θ aq y 1 1 6.0759 1.1930 3 7.2399 0.9477 5 6.9866 0.6859 TABLE III. Vector current renormalization constants, Z V ðμÞ, using the RI-SMOM scheme. Values are taken from Tables III and VI in [30] where they were calculated in pure QCD for each β value corresponding to a group of ensembles in Table I. We use the values given at μ ¼ 2 GeV but note that the small μ dependence in this quantity is purely a lattice artifact, vanishing in the continuum limit. Here we give values for each set (listed in column 1), but those with the same β value are the same. Z γ i ⊗γ i V (column 3) is the renormalization constant for the local current used in the LOCAL setup and Z γ i ⊗1 V (column 4) is the 1-link current renormalization needed for the ONE-LINK setup. Note, in the latter case, that this is for the 1-link tadpole-improved current constructed with a "thin-link" included. The tadpoleimprovement factor u 0 , given by the mean value of the gluon field U μ in Landau gauge, is listed in column 5. Since we use a 1-link vector current that is not tadpole-improved, the renormalization factor for the ONE-LINK case is

5.80 0.95932(18) 0.93516(16) 0.81960 2 5.80 0.95932(18) 0.93516(16) 0.82042 3 6.00 0.97255(22) 0.94966(20) 0.83461 4 6.00 0.97255(22) 0.94966(20) 0.83505 5 6.30 0.98445(11) 0.96695(11) 0.85248 6 6.72 0.99090(36) 0.97996(34) 0.87094

This is plotted as a function of the time separation between vector current operators for a fixed value of t η ct γ 2 of N t =4. We divide by expð-M η c ðt η ct γ 2 ÞÞ to remove the time dependence related to the η c mass expected from Eq. ( 14). The η c mass used here is the one from the fit to C η c [Eq. ( 11)] (i.e. M latt η c ). We also divided by a 0 , which is the ground-state amplitude from this fit. This means that the integral of R μν is equal to b 0 [Eq. ( 14)] up to excited-state contamination, which is very small at this large value for t η ct γ 2 . Indeed the figure is unchanged over a wide range of t η ct γ 2 values since excited state contamination falls off on distance scales of Oð0.5 fmÞ and N t =4 ≈ 2 fm.

We see from Fig. 2 that R μν is strongly dominated by the region of t γ 1 very close to t γ 2 . This is because, once a photon is emitted from the η c , the resulting system is far off shell and decays exponentially fast in t. The region in jt γ 1t γ 2 j for which the integrand is nonzero is less than about 0.5 fm.

Because we use staggered quarks the integrand has a component that oscillates in time. On summing/integrating, however, this component reduces to an afoot_1 discretization effect, because the oscillations get closer together on finer lattices, as we show in Appendix A. Discretization effects of this kind are allowed for in our extrapolation of our results for F to the a → 0 continuum limit.

We construct Cμν ðt γ 2 ; t η c Þ for all values of t η ct γ 2 and then fit it as described in the next section, in conjunction with C η c , to determine b 0 . Our fits take full account of excited-state contamination.

Fitting the correlators Cμν and C η c

We fit Cμν and C η c simultaneously to standard staggeredquark fit forms for two-point correlation functions. These are the same as those given in Eqs. ( 11) and ( 14) except that there are additional terms that oscillate in time. The fit form for C η c becomes

with fðE; tÞ given in Eq. ( 12). The oscillating terms arise from opposite parity states, all of which are heavier than the ground-state η c and decay faster in the large-time limit.

The fit-form for Cμν given in Eq. ( 14) is likewise extended to include oscillating terms from the same opposite parity states.

We use the CORRFITTER package [31] to fit our two-point correlation functions, taking as parameters the logarithms of the energy differences (to keep the energy levels ordered) and the logarithms of the amplitudes. We take the prior on the energy differences to be 0.5 (2) GeV and the prior width on the amplitudes a n to be 20% and b n [see Eq. ( 14)] 50%. All prior widths are orders of magnitude larger than the uncertainties returned by the fit for ground-state values. We find χ 2 =dof ≲ 1 for all fits. 2 Our final fits use 4 exponentials, i.e. N n ¼ N o ¼ 4 in Eq. ( 24). We do not include results for values of t below t min in our fits where t min =a varies between 6 and 9 depending on lattice spacing for Cμν and is N t =8 for C η c .

The aim of the fits is to determine the ground-state parameters corresponding to n ¼ 0: E 0 , a 0 and b 0 . From b 0 we can determine the form factor for η c → γγ using Eq. ( 17). This equation is derived from Eqs. ( 15) and ( 16) and assumes that both photons are exactly on shell. As discussed in Sec. II B 2, photon 1 is exactly on shell but photon 2 is not. We must therefore modify our determination of F from b 0 to take account of this, following Eq. ( 16). Instead of Eq. ( 17) we must use

where q 2 2 is the q 2 value for photon 2 and b 0 is in lattice units. aM latt η c is the value of the η c mass in lattice units given by E 0 .

Table IV gives our results for aM latt η c and F latt ð0; q 2 2 Þ=a on each set of gluon field configurations and for both the FIG. 2. R μν , defined in Eq. ( 23), as a function of t γ 1t γ 2 . We use a fixed value for the separation between t η c and t γ 2 of N t =4. Results are given for the LOCAL setup and plotted in physical units for two values of the lattice spacing, corresponding to very coarse (set 1: a ≈ 0.15 fm) and superfine (set 6: a ≈ 0.06 fm). Integrating R μν gives b 0 and hence F [see Eq. ( 17)]; note the narrow width of the region of support for the integral. Oscillations are a result of using staggered quarks. Their impact on the integral is a discretization effect-see Appendix A.

LOCAL and ONE-LINK setups. The results also include those for a mistuned valence c mass on set 3, denoted set 3A. Notice the small statistical uncertainties, well below 1% in F=a and aM η c , typical of lattice QCD calculations with heavy quarks. We include values for the off shellness of photon 2, q 2 2 , given by

on the coarsest lattices for the LOCAL setup. We will discuss how we fit the results for F (allowing for the off shellness q 2 2 ) to obtain a value for Fð0; 0Þ in the physical continuum limit in Sec. II B 6.

First we demonstrate that results can be obtained with a subset of the correlation functions that we have calculated here. Figure 3 shows the results for F if we restrict the range of integration over t γ 1 to a time distance of t width either side of t γ 2 . In keeping with Fig. 2, we see that the result for F reaches its final value very quickly as a function of t width (in less than 1 fm). The sum over t γ 1 could be truncated in this case with no loss of accuracy. We perform the full sum here, however.

In Appendix B we further show that we can restrict the fit of the two-point function Cμν ðtÞ to a set of specific t ≡ t γ 2t η c values rather than fitting the full t range. Here, however, we use the results from our full fit.

As discussed in Sec. II B 2, we use different staggered spin-taste representations for the mesons in different parts of our calculations. Here we test that the different representations agree in mass in the continuum limit. We use an interpolating operator with spin taste γ 5 γ t ⊗ γ x γ z for the pseudoscalar η c meson in our LOCAL setup and γ 5 γ t ⊗ 1 for our ONE-LINK setup. In the study of J=ψ → γη c , we instead take a γ 5 ⊗ γ 5 (Goldstone) interpolator for the η c meson. This latter spin taste corresponds to the lightest η c in the taste multiplet and the one used, for example, in [1]. In Fig. 4, we compare the masses of the different η c tastes determined from our fits to the suite of correlation functions

TABLE IV. Results for F latt ð0; q 2

2 Þ and M η c in lattice units obtained from our correlator fits [Eqs. (24) and (25)]. The top table shows results from our LOCAL setup (with local vector current), the lower table those from the ONE-LINK setup (with 1-link vector current). Column 5 gives the value for q 2 2 (q 2 for photon 2) in lattice units from Eq. ( 26) (these values are close to zero for the ONE-LINK setup from the way that the momentum twists were chosen). Set 3A for the LOCAL case corresponds to a deliberately mistuned valence c quark mass (see Table I). We fit the correlators for sets 3 and 3A simultaneously so that correlations between them are fully taken into account.

Taking the physical-continuum limit

To obtain a physical result for the form factor Fð0; 0Þ, we fit our lattice data for Fð0; q 2 2 Þ from Eq. ( 25) (given in Table IV and plotted in Fig. 5) to a function that accounts for discretization effects and mistunings of the sea and valence quark masses, as well as allowing for the small amount by which photon 2 is off shell. The data is fit using the Lsqfit package [33] and our fit form is

where F ðtÞ latt ð0; q 2 2 Þ is the lattice data with superscript t denoting the taste, i.e. either the LOCAL (local vector) or ONE-LINK (1-link vector) cases. Fð0; 0Þ is the form factor in the limit of vanishing lattice spacing and physical masses that we wish to determine. The factor of ð1q 2 2 =M 2 pole Þ allows us to adjust for the amount by which photon 2 is offshell with a simple one-pole parametrization which was tested in [7,12]. The pole form was found to work well, with the pole mass M pole taking a value around the J=ψ mass. Here we will take M pole as a fit parameter, with prior 3.0(3) GeV. As discussed in Sec. II B 4, our q 2 2 values are close to zero here (see Table IV) and the pole factor only has a small effect, at most 2% on our coarsest lattices in the LOCAL setup and less than 1% in other cases.

Equation ( 27) takes discretization effects into account, mainly through the first nontrivial term in the square brackets. These discretization effects arise from the HISQ action but also through the trapezoidal integration used to determine Cμν and the oscillating term contribution to that integral (Appendix A). We allow for the size of these discretization effects to be set by scale Λ and include terms up ðaΛ 2i max Þ. Since we want to fit both the LOCAL and ONE-LINK cases with the same fit form, although they differ significantly in their discretization effects (see Fig. 5), we give them different Λ parameters and choose Λ ðtÞ using the Empirical Bayes criterion. This means varying Λ ðtÞ in fits to the two cases (with independent fit parameters, κ ði;tÞ aΛ ) and taking the values that maximize the Bayes factor [34]. We also allow for discretization effects coming from the sea in the final term, including terms up to ða ΛÞ 4 . We take Λ ¼ 1 GeV in both cases but use independent fit parameters, κ ð1=2;tÞ sea;uds . We allow for mistuning of sea and valence quark masses in Eq. ( 27) with the terms containing δ. The mistunings of the c valence and sea quark masses and the u, d and s sea quark masses are expressed as

respectively. To tune the valence and sea c quark mass, we take am tuned c to be

where M expt J=ψ ¼ 3.0969 GeV from [6], and lattice values for aM J=ψ are obtained from Table III in [1]. The power of 1.5 FIG. 5. Data points show our lattice QCD results for the η c → γγ form factor, in red for the LOCAL setup (with local vector current) and green for the ONE-LINK setup (with 1-link vector current). The LOCAL points have been corrected for the slight off shellness of one photon (see text). The LOCAL points include one at a deliberately mistuned c quark mass (open red square). The blue and pink bands show our chiral/continuum fit to the points, applying Eq. ( 27) simultaneously for the LOCAL and ONE-LINK points. The fit bands are plotted as a function of a at the physical quark mass point. The physical result in the continuum limit is then shown by the black star.

is empirically chosen, based on the results from [1] (the power is not 1, because of the binding energy inside the J=ψ).

The tuned s quark mass is given by

from leading-order chiral perturbation theory. am tuned l is found by dividing the value for m tuned s in Eq. ( 32) by the ratio [35]

For all except set 5 we use δm sea uds =m s values from Table I of [36]. For set 5 in Table I, values in lattice units for M η s are taken from [37] and the "physical" value of the η s meson, 688.5(2.2) MeV, is taken from [25]. The η s is an unphysical pseudoscalar ss meson whose mass can nevertheless be determined in terms of π and K masses in lattice QCD. This gives a δ sea;uds value of 0.0297 (17) for set 5.

Since the dependence of Fð0; 0Þ on quark masses is a physical effect it should be the same (up to discretization errors) for the LOCAL and ONE-LINK cases. We therefore take the κ sea;c , κ val;c and κ ð0Þ sea;uds parameters to be the same in the two cases when fitting them simultaneously.

The parameters to be determined by the fit are M pole , κ ðiÞ aΛ for i ¼ 1 to i max , κ sea;c , κ val;c and κ ðjÞ sea;uds for j ¼ 0, 1, 2. As discussed above, we take the parameters to be independent for our two setups when they correspond to discretization effects, but otherwise take them to be the same. For priors, we use Fð0; 0Þ ¼ 0.1ð1Þ, M pole ¼ 3.0ð3Þ GeV, κ ðiÞ aΛ ¼ 0ð1Þ [except for i ¼ 1 for the ONE-LINK case where we take the prior to be 0(2) from inspection of Fig. 5], κ val;c ¼ 0ð2Þ (from comparison of results from sets 3 and 3A), κ ðkÞ sea;uds ¼ 0ð1Þ and κ sea;c ¼ 0.0ð1Þ (since we expect the effect of c in the sea to be minor). For our final fit we take the largest coefficient for discretization effects, i max ¼ 3. Our preferred fit is a joint fit to the LOCAL and ONE-LINK data, but we obtain almost identical results from fitting simply the LOCAL results.

Figure 5 shows our lattice results (now in physical units) for both the LOCAL and ONE-LINK setups as a function of squared lattice spacing. The lattice results have been adjusted to correspond to the Fð0; 0Þ on-shell point, i.e. F latt ð0; q 2 2 Þ has been multiplied by the pole term ð1q 2 2 =M 2 pole Þ. The uncertainty on the lattice results is dominated by the correlated uncertainty in the value of the lattice spacing. Also shown is our continuum fit using Eq. ( 27) to both sets of data simultaneously. Notice that the discretization effects are much larger in the ONE-LINK case than in the LOCAL case. Using the Empirical Bayes criterion [34] we find that the optimal Λ is 0.10 GeV for LOCAL and 0.49 GeV for ONE-LINK. The larger discretization effects for ONE-LINK are not surprising because the vector currents in that case are 1-link operators with tree-level discretization errors. The LOCAL case, in contrast, uses local vector current operators that have no tree-level errors at any order in a. The bands plotted on Fig. 5 correspond to the fit at tuned sea masses as a function of lattice spacing, i.e. Fð0; 0Þ½1 þ P i κ ði;tÞ aΛ ðaΛ ðtÞ Þ 2i  [see Eq. ( 27)], with i max ¼ 3.

The result we obtain for Fð0; 0Þ in the continuum limit from the joint fit to the LOCAL and ONE-LINK results is Fð0; 0Þ ¼ 0.08793ð29Þ GeV -1 with a χ 2 =d:o:f: of 0.8. This confirms that the LOCAL and ONE-LINK results can readily be fit to a common continuum value, providing a test that HISQ taste effects are purely lattice artifacts. The result from fitting LOCAL alone is very similar, not surprisingly because we have the best coverage of lattice spacing and sea masses in that case and discretization effects are smaller than for the ONE-LINK case. In Sec. II C we will discuss additional sources of systematic error that must be accounted for in our final result.

First we discuss the stability of our fitted value for Fð0; 0Þ as we change details of our fits. This is shown in Fig. 6, where we plot the value of Fð0; 0Þ in the physical/ continuum limit on changing some aspect of either the correlator fits or the chiral/continuum fit. The preferred (base) fit described above is given on the left. Variations include making all the priors a factor of 2 smaller or larger and dropping datasets at either end of the lattice spacing FIG. 6. The value of Fð0; 0Þ in the limit of vanishing lattice spacing and physical quark masses obtained from variations to our base fit. These include (from left to right) dropping the coarsest and finest datasets, changing all the prior widths in our correlator fits, changing all the prior widths in our chiral/ continuum fits and adding an additional normal and oscillating exponential to our correlator fits. Note that the values for Λ in Eq. ( 27) are fixed (see text) under these fit variations.

range. We see very little variation in the final answer under any of these variations, showing that our result is robust.

Additional systematic uncertainties

As discussed in Sec. II A our lattice QCD calculation does not include quark-line disconnected contributions or QED effects. Here we estimate the size of these and include an additional systematic uncertainty to allow for them.

In [1] the small difference between the mass of the η c obtained from a lattice QCD calculation including connected correlation functions only and the mass found in experiment was interpreted as the effect on the mass of missing quark-line disconnected correlation functions. The Wick contractions for the disconnected correlation functions include the annihilation of the η c into gluons and allow mixing between the η c and other flavor-singlet pseudoscalar states. The difference found in [1] amounted to 0.2% of the η c mass. A reasonable estimate of the impact of these disconnected diagrams on the η c wave function (which can be related to the decay amplitude to γγ in the nonrelativistic limit) would then also be 0.2%, giving a 0.4% uncertainty in the decay width. That this is reasonable is confirmed by our fit, which includes the dependence of Fð0; 0Þ on the η c mass and returns a coefficient for κ val;c close to 1. This dependence is visible in Fig. 5 and Table IV, comparing values for tuned and mistuned m c .

We can also consider the impact of another class of quark-line disconnected diagrams in which the η c radiates a photon before annihilating to gluons. The quark loop generated from the gluons then also annihilates to a photon. The impact of this diagram should be very small, because of suppression both by powers of α s and by quark mass effects since the sum of the electric charges of the light quarks in the sea is zero. We therefore expect a contribution smaller than a relative size of α 2 s m 2 s =m 2 c ≈ 0.2%. This does not then modify our estimate of the uncertainty from missing disconnected diagrams as 0.2% from above.

The impact of the c quark's electric charge on the decay amplitude (also missing in our calculation) can be estimated from the impact of QED on the η c decay constant, determined in [1]. This effect was 0.17%, so we allow an additional uncertainty of 0.2% in Fð0; 0Þ for this. Further QED corrections to the decay width coming from additional radiation are tiny since there are no electrically charged particles in either the initial or final states. By charge conjugation, the η c cannot decay to γγγ, so QED corrections to the decay rate would come from η c to γγγγ. This would be suppressed by a further 2 powers of α, and thus negligible.

Adding these two 0.2% systematic uncertainties in quadrature gives an additional uncertainty of 0.3% in Fð0; 0Þ and 0.6% to the decay rate.

C. Results

We take our final result for the form factor for η c → γγ from the joint fit to the LOCAL and ONE-LINK cases, giving Fð0; 0Þ ¼ 0.08793ð29Þ fit ð26Þ syst GeV -1 ; ð34Þ which has a total uncertainty of 0.4%. The second uncertainty in Eq. ( 34) comes from the additional systematic uncertainties that we estimate in Sec. II B 7.

The error budget for this result is given in Table V. We see that the dominant uncertainty is that from fixing the lattice spacing using w 0 . The next most important uncertainties come from statistics and from systematic errors from missing quark-line disconnected diagrams and QED effects.

We use our value for Fð0; 0Þ in Eq. ( 34) and the formula in Eq. ( 19) to find the decay width Γðη c → γγÞ. For the η c mass in Eq. ( 19), we use the experimental value M exp η c ¼ 2.9839ð4Þ GeV (the average from [6]) since this is a purely kinematic factor. We also take 1=α ¼ 137.036 from [6], noting that the momentum scale for this decay is a relatively low one. We obtain the decay width with an uncertainty of 0.9% as Γðη c → γγÞ ¼ 6.788ð45Þ fit ð41Þ syst keV:

Using the experimental value for the η c total width of 32.0 (7) MeV [6], this corresponds to a branching fraction of Bðη c → γγÞ ¼ 2.121ð14Þ fit ð13Þ syst ð46Þ expt × 10 -4 : ð36Þ

TABLE V. Error budgets for our results for Fð0; 0Þ and Vð0Þ (discussed in Sec. III C). The top seven entries come from our fits to Eqs. ( 27) and ( 49) respectively. The uncertainty labeled "q 2 dependence" arises from the tuning to q 2 ¼ 0 in the Fð0; 0Þ case and from capturing the q 2 dependence in the Vð0Þ case.

The lowest two entries are additional systematic uncertainties from missing quark-line disconnected correlations functions and QED effects. These are discussed in Sec. II B 7 for Fð0; 0Þ and in Sec. III B 3 for Vð0Þ, where the two sources of uncertainty are combined together. Fð0; 0Þ Vð0Þ Statistics 0.17 0.29 w 0 =a 0.05 0.01 w 0 0.25 0.07 a 2 → 0 0.08 0.23 Valence mistuning 0.01 0.03 Sea mistuning 0.05 0.05 q 2 mistuning 0.01 0.10 Missing "disconnected" correlators 0.2 0.4 Missing QED 0.2 -Total 0.43 0.56

The third uncertainty here is from the experimental total width and it dominates over the lattice QCD uncertainties.

D. Discussion

Figure 7 compares our result for Γðη c → γγÞ from Eq. ( 35) with earlier lattice QCD results obtained on gluon field configurations that include sea quarks. It is more natural to compare results for Fð0; 0Þ from lattice calculations but some of the earlier results do not include this information. The values plotted come from [15] (MFLWZ21), [14] (LMZ20) and [13] (CLQCD20) that all work with twisted mass quarks on n f ¼ 2 gluon field configurations (i.e. including u and d quarks in the sea but no s quarks) at two or three values of the lattice spacing. We do not include the earlier value from CLQCD in [12] that is superseded by [13]. CLQCD20 does not include a continuum extrapolation; we take the value as that quoted for their finest lattice. The uncertainty quoted does not include systematic errors. We plot the continuum results for the η c → γγ decay rate quoted by MFLWZ21 and LMZ20, with systematic uncertainties combined in quadrature. Because the sea content is unphysical for the n f ¼ 2 case, those results will not necessarily agree with ours; no uncertainty is included in CLQCD20, LMZ20 or MFLW21 for missing s quarks in the sea. Our result is larger than the two earlier values with a tension exceeding 3σ, but it agrees well with that from MFLWZ21. We note that the three n f ¼ 2 results do not agree well with each other, however. Our value is obtained on n f ¼ 2 þ 1 þ 1 gluon field configurations (i.e. with a realistic sea quark content) at four values of the lattice spacing and includes systematic errors for missing quark-line disconnected diagrams and QED. Our total uncertainty is smaller than that of the earlier values. Further lattice calculations with a full sea quark complement and uncertainties comparable to ours are needed.

Figure 8 compares our result for Γðη c → γγÞ to results from experiment. There is a lot of experimental information that has a bearing on Γðη c → γγÞ, generally obtained as products of branching fractions depending on the method of η c production and the decay mode observed. The PDG [6] combines this information to yield a fit value for Γðη c → γγÞ of 5.15 (35) keV. This is shown by the blue band in Fig. 8. It is lower than our lattice QCD result by 4.6σ, where we have combined lattice and PDG fit uncertainties in quadrature (the PDG fit uncertainty dominates) to obtain σ. This denotes very significant tension between the Standard Model and experiment, which could be taken as an indication of new physics. The PDG fit returns a large χ 2 of 118 for 81 degrees of freedom [6], however, and this calls into question the reliability of both their central fit value and its uncertainty. We note that an alternative global fit of η c data which does not include any pre-1995 results, has a better χ 2 and gives a larger partial width for η c → γγ decay of 5.43 þ0. 41 -0.38 keV [38]. This shows a reduced, but still sizeable, tension of 3.3σ with our lattice QCD result. FIG. 7. A comparison of full lattice QCD results for the width for η c → γγ decay. The result obtained here is denoted "HPQCD23" (red asterisk, error bar same size as symbol) and uses gluon field configurations that include n f ¼ 2 þ 1 þ 1 sea quark flavors at four values of the lattice spacing to determine a physical result. Earlier results use n f ¼ 2 gluon field configurations at two values of the lattice spacing (orange filled circles) or three values (orange filled triangle). The points denoted "MFLWZ21" from [15] and "LMZ20" from [14] are from determinations of the rate in the continuum limit including an estimate of systematic errors, although no error for missing s quarks in the sea is included. The point denoted "CLQCD20" corresponds to the value quoted at the finest lattice spacing used in [13] with no extrapolation to the continuum limit. The red band carries our result down the plot for comparison. FIG. 8. A comparison of our lattice QCD result (HPQCD23, red asterisk and red band) for the decay width Γðη c → γγÞ to experimental values taken from [6]. Blue circles denote determinations of the product Γðη c → iÞΓðη c → γγÞ=Γ total ðη c Þ, where the decay channel i is listed on the right. The values plotted are derived from either single experimental results or PDG average values from several experiments for the product above, divided by the branching fraction BðiÞ, again using either a single result or the PDG average. Uncertainties are combined in quadrature. The blue square denotes a determination from the product of branching fractions for J=ψ → γη c and η c → γγ from BESIII [39] combined with the PDG average for the branching fraction BðJ=ψ → γη c Þ. The blue band shows the result of the PDG fit to the experimental data shown here along with other results (e.g. for ratios of branching fractions). The fit has a χ 2 of 118 for 81 degrees of freedom [6], reflecting the inconsistencies in the experimental data seen above.

Individual experimental results have much larger uncertainties and a large spread of central values. Figure 8 shows values derived from results listed in [6]. The blue circles use values quoted in the section headed "ΓðiÞΓðγγÞ=ΓðtotalÞ," where the η c is produced via two-photon fusion in e þ e - collisions and detected through its decay to channel i listed on the right of Fig. 8. These values are determined without reference to the PDG fit. They use either the PDG average value (where one is given) or the single experimental value (if there is no average) for the product above. We then divide by the branching fraction for channel i again using either the PDG average or the single experimental value quoted, if there is no average. Uncertainties are combined in quadrature. We see a large spread of experimental values in Fig. 8, several of which are in significant (4σ) tension with the PDG fit result. This is not surprising given the χ 2 value for the fit. Some of the low values for Γðη c → γγÞ seen are in disagreement with our result; the values from the ϕϕ and K þ K -π þ π -channels differ by an amount exceeding 6σ. On the other hand, the experimental result using the K Kπ channel is in good agreement with our value, within 2σ. The K Kπ channel has been studied by several experiments because it has a relatively large branching fraction. The value plotted comes from an average of 8 different experimental results for the product of rates (the average is dominated by results from CLEO [40] and BABAR [41]) and 10 for the branching fraction (where the average is dominated by results from BESIII [42]). This gives a final result for Γðη c → γγÞ of 5.90 (58) keV. The 10% uncertainty is the smallest relative uncertainty for any channel.

We conclude that the experimental picture of η c decay is not yet a very coherent one. Further experimental results, with small uncertainties, will be needed to resolve the issue of whether or not there is tension between experiment and lattice QCD/the Standard Model for η c → γγ decay.

The filled blue square in Fig. 8 comes from a determination of the product of branching fractions for J=ψ → γη c and η c → γγ by BESIII [39] combined with the PDG average for the branching fraction BðJ=ψ → γη c Þ [6]. The value agrees with our result, but has a large uncertainty. The PDG gives an average value for the branching ratio of η c → γγ by combining this result with an earlier one for the same product of branching fractions from CLEO [43]. 3 This gives an average branching fraction of 2.2 þ0.9 -0.6 × 10 -4 [44]. This also agrees with our value within its large uncertainties. We will discuss these results further in Sec. IV.

We can also explore what our results imply about the nonrelativistic nature of the c quarks inside the η c . As discussed in Sec. I there is a very simple relationship between Γðη c → γγÞ and ΓðJ=ψ → e þ e -Þ in LO NRQCD [Eq. ( 1)]. In Fig. 9 we compare this LO result, shown as a green dashed line, to our lattice QCD value (red asterisk). The central value for the LO result (7.5

keV) uses

ΓðJ=ψ → e þ e -Þ ¼ 5.637ð49Þ keV from lattice QCD þ QED [1] [this lattice QCD þ QED result agrees well with the experimental average of 5.53 (10) keV [6] but has a smaller uncertainty]. The LO NRQCD central value is then 10% above our lattice QCD result. Figure 9 shows a 30% error band in green around the LO NRQCD central value to allow for subleading corrections. The lattice QCD result, incorporating the full relativistic dynamics of the c quarks, falls well within this 30% band showing that the LO nonrelativistic approximation works well here.

The green points with error bars show results from two different calculations in continuum NRQCD, going beyond LO. In CM01 [16] higher-order QCD and relativistic corrections to Eq. (1) were added. The authors find substantial Oð30%Þ corrections from the two sources but also see a large amount of cancellation between them. They conclude with an estimate of a 10% upward shift of Γðη c → γγÞ compared to the LO expression, with a 5% uncertainty. They warn, however, that missing higherorder corrections might be substantial. The comparison with our result in Fig. 9 shows this to be the case because their 10% shift has taken them in the wrong direction from the LO result and their 5% uncertainty is insufficient to cover the gap. References [48,49] extend the analysis of FIG. 9. A comparison of our lattice QCD result (HPQCD23, red asterisk) for the decay width Γðη c → γγÞ to values obtained from theory calculations using NRQCD. The green dashed line is the leading order NRQCD result from Eq. ( 1), i.e. 4ΓðJ=ψ → e þ e -Þ=3 where we have taken the J=ψ leptonic width from lattice QCD [1] (which agrees well with experiment). The uncertainty in the LO NRQCD value from higher order corrections is 30% and denoted by the green band. The green points give results from higher-order calculations. CM01 is from [16], determining QCD and velocity-expansion corrections to the ratio of ΓðJ=ψ → e þ e -Þ=Γðη c → γγÞ. FJS17 [45] calculates NNLO QCD corrections to the η c → γγ branching fraction through v 2 order in the NRQCD velocity expansion. BC01 [46] and BCK18 [47] calculate the inverse branching fraction resumming QCD corrections in the large n f limit. They give results from two resummation methods, naïve non-Abelianization (NNA) and background-field gauge (BFG). We use the PDG average [6] for the η c total width to convert their branching fractions into a width for η c → γγ. The blue band repeats the PDG fit result shown in Fig. 8. this approach but do not quote final values, concluding that the theoretical uncertainties are very large for the charmonium case.

FJS17 [45] calculates next-to-next-to-leading-order (NNLO) QCD corrections to the η c → γγ branching fraction through v 2 order in the NRQCD velocity expansion and concludes that NRQCD factorization does not work well for this case, given the large value that they obtain for the branching fraction. BC01 [46] and BCK18 [47] resum QCD corrections to the inverse branching fraction in the large n f limit each using two different approaches, and allowing uncertainties for missing color-octet contributions. We use the PDG average [6] for the η c total width to convert their branching fractions into a width for η c → γγ. Given the large uncertainties they have, all of their results agree and both are in reasonable agreement with our value. It is clear, however, that accurate results for Γðη c → γγÞ are only currently available using lattice QCD.

We can extend our comparison to the LO NRQCD expectation by converting Eq. ( 1) into a relationship between the hadronic parameters Fð0; 0Þ and the J=ψ decay constant, f J=ψ . Using

along with Eqs. ( 1) and (19) we see that the LO NRQCD result implies

independent of Q c . Since there is no distinction between M η c and M J=ψ at this order in NRQCD, we can rewrite this as

R fF is now a simple combination of hadronic parameters that we can calculate directly in lattice QCD.

Figure 10 shows our results for R fF as a function of lattice spacing using Fð0; 0Þ values for the LOCAL setup (which corresponds to our largest set of results in Fig. 5). We take values of M J=ψ and f J=ψ from [1], not including the QED effects that are calculated there. We also determine values for the mistuned m c case (set 3A) from those results.

We obtain a value for the ratio in the continuum limit at physical quark masses by performing a fit of the same form as that in Eq. ( 27) to our lattice results for f J=ψ =½ð1q 2 2 =M 2 pole ÞFð0; q 2 2 ÞM 2 J=ψ . The values show strong lattice-spacing dependence coming from the decay constant results. This is because the annihilation of a J=ψ meson is a very short-distance process (shorter distance than η c → γγ where the energy is shared between two photons). Using the Empirical Bayes criterion we find a value for Λ for this fit of 0.86 GeV, larger than that seen in either of the earlier fits shown in Fig. 5. Because of the larger discretization effects we take i max ¼ 5 [see Eq. ( 27)] for this fit. The result in the continuum limit does not change between i max ¼ 4 or 5. Notice also that the results vary little between the tuned and mistuned c quark mass in contrast to what was seen in Fig. 5. From NRQCD we expect the result for the ratio in Eq. ( 39) to depend on the (heavy) quark mass only through subleading terms in the velocity expansion.

The result for the ratio of Eq. ( 39) that we obtain in the continuum limit and at physical quark masses is

Here we have added a second uncertainty of 0.3% to allow for additional systematic errors in Fð0; 0Þ as discussed in Sec. II B 7, giving a total uncertainty of 1.2%. The χ 2 =d:o:f: for the fit is 0.27. The result of Eq. ( 40) for R fF is (only) 4.5(1.2)% below the LO NRQCD value of 0.5 given in Eq. ( 39). This test of LO NRQCD is slightly different, although equally valid, to that given in Fig. 9. This is why they both result in a LO NRQCD value that is above the lattice QCD result, even though the quantities tested are inversely related to each other. The differences arise from small effects that are ignored in LO NRQCD. These include the FIG. 10. The ratio of J=ψ decay constant to Fð0; 0Þ multiplied by the square of the J=ψ mass (red points) as a function of squared lattice spacing using our Fð0; 0Þ values for the LOCAL setup (with a local vector current) and results from [1] for f J=ψ and M J=ψ . The values for F have been corrected for the slight off shellness of one photon (see text). The blue band gives our continuum fit (see text) and the black dotted line gives the LO NRQCD result of 0.5 [Eq. ( 39)].

differences in the scale of α used in ΓðJ=ψ → e þ e -Þfoot_3 and Γðη c → γγÞ for the comparison in Fig. 9 and the simplification of the combination of meson masses in going from Eqs. (38) and (39). Once both of these effects are taken into account, the two comparisons in Figs. 9 and 10 are consistent. Looking at either figure, we conclude that LO NRQCD is better than might have been expected as an approximation here. This only becomes clear with accurate lattice QCD results.

The J=ψ → γη c decay is an M1 electromagnetic mesonto-meson transition. The rate is straightforward to determine in lattice QCD. Using 3-point and 2-point correlation functions, to be discussed below, we can calculate the matrix element of the c electromagnetic current, j μ c ¼ cγ μ c, between the initial J=ψ and final η c states. Note that, in keeping with Sec. II, we do not include a factor of the electric charge in the current. The matrix element can be parametrized by the form factor Vðq 2 Þ as

where ϵ J=ψ σ is the polarization of the J=ψ meson and q ¼ pp 0 is the 4-momentum transfer.

The decay width is then given by [7] ΓðJ=ψ

where V is evaluated at q 2 ¼ 0 for an on-shell photon. jkj takes the value

corresponding to the spatial momentum of the η c (or photon) in the J=ψ rest frame at q 2 ¼ 0.

It is convenient for us to use Vðq 2 Þ rather than the more conventional Vðq 2 Þ ≡ 2 Vðq 2 Þ because in our lattice QCD calculation we will compute only one of the two diagrams that contribute identically. A photon can be emitted by either the c or c constituent quarks of the initial J=ψ meson and we will evaluate one of these two cases.

Whilst the decay width to a real photon is a function of Vðq 2 ¼ 0Þ, we can also map out the q 2 dependence of V up to q 2 max ¼ ðM J=ψ -M η c Þ 2 . This is needed to compute the width of the Dalitz decay ΓðJ=ψ → η c e þ e -Þ. We do this by defining the ratio

The derivative of R eeγ with respect to q 2 , the squared 4momentum of the virtual photon in J=ψ → η c e þ e -, can be written in terms of the form factor Vðq 2 Þ as [50]

By multiplying R eeγ by ΓðJ=ψ → γη c Þ [determined through Eq. ( 42)], we can then also find the decay width ΓðJ=ψ → η c e þ e -Þ.

Correlation functions

The matrix element given in Eq. ( 41) can be obtained by computing an appropriate 3-point correlation function on the lattice. This correlation function, summarized diagrammatically in Fig. 11, is similar in appearance to the correlator given in Fig. 1 used in the calculation for η c → γγ; the difference is that one of the vector operators now couples to a cc vector meson. We consider only the connected contribution to the decay J=ψ → η c γ ðÞ which comes from an insertion of a cγ μ c current onto one of the two quark lines between the initial and final states. There are two such diagrams which contribute equally; we compute only one of them. We expect quark-line disconnected contributions to be very small here because of the large mass of the c quark [7]. FIG. 11. Schematic diagram for the connected 3-point correlation function for the J=ψ transition to η c via an electromagnetic current. The lines between the operators represent quark propagators at the mass of the charm quark.

We will discuss the uncertainty associated with neglecting these diagrams in Sec. III B 3.

Our 3-point functions have three operator insertions, two vector and one pseudoscalar, just as in the 3-point function in Sec. II B 2. We use different spin-taste combinations here, however. We take the γ 5 ⊗ γ 5 interpolator for the η c meson, γ x γ t ⊗ γ 5 γ z to interpolate the J=ψ meson, and the local γ z ⊗ γ z operator for the electromagnetic vector current insertion. 5 As discussed in Sec. II B 3, this allows us to use the very precise values from [30] for the multiplicative renormalization Z V for the local HISQ-HISQ vector current.

Spatial momentum is applied to the charm propagator between the η c operator and the vector current insertion to give momentum to the η c meson while the J=ψ meson is always at rest on the lattice. As well as the three-momentum necessary to achieve q 2 ¼ 0, the case of an on-shell photon, we calculate correlation functions for various other 3momenta of the η c meson so that we obtain results for a range of q 2 values around q 2 ¼ 0. In addition we calculate 2-point correlation functions corresponding to the operator used for the J=ψ (at rest) and for the η c (at the spatial momentum values used).

As for η c → γγ (see Sec. II B 2), we insert spatial momentum by using a twist angle. The twist values used for each set are tabulated in Table VI. The twists are chosen to broadly cover the full range 0 ≤ q 2 ≤ q 2 max ¼ ðM J=ψ -M η c Þ 2 . Given that we have chosen x and z polarizations for the J=ψ and the vector current, respectively, then Eq. ( 41) requires the spatial momentum to have a component in the y direction. In fact we take the momentum to be parallel to the y axis.

We calculate 3-point correlation functions for several different values of the source-sink separation, T ¼ t J=ψt η c (see Fig. 11). This improves the determination of the ground-state to ground-state matrix elements that we are interested in by giving a better handle on the excited states present in the correlation function. The values of T that we use are also listed in Table VI.

As in Sec. II B 2, we fit the 2-point correlation functions to the form in Eq. ( 24). The 3-point functions (for all T and momenta) are simultaneously fit to the standard form

with the addition of terms that oscillate in time and involve opposite parity states [10] that are not shown here. The energies of the states with overlap onto γ x γ t ⊗ γ 5 γ z (the J=ψ and its excitations) are indexed by i, and the energies of the states with overlap onto γ 5 ⊗ γ 5 (the η c and its excitations) are indexed by j. The energies and the amplitudes a i=j are also parameters for the 2-point correlator fits and this enables us to extract the parameters V ij . Our fits account for four normal vector and pseudoscalar charmonium states in Eq. ( 46) (i.e. N n ¼ 4) and four additional oscillating states for each. As for the fit to C η c described in Sec. II B 5, we drop 2-point correlator data below t min =a ¼ N t =8 when fitting. We also drop 3-point correlator data below a t min =a [and ðT -tÞ min =a] value of between 2 and 4, increasing as the lattice spacing falls.

The ground-state to ground-state matrix element that characterizes the J=ψ → η c transition is related to the parameter V 00 of Eq. ( 46) via

where M latt J=ψ ≡ E i¼0 and E latt η c ≡ E j¼0 and the p factors are the relativistic normalizations of the states jJ=ψi and jη c i respectively. Z V is the multiplicative renormalization factor for the local γ z ⊗ γ z vector current. The factor of 1=2 present on the left-hand side above accounts for the second Wick contraction (that contributes equally) where the vector current insertion is placed onto the other quark line connecting the initial and final states.

where q y is the component of the η c momentum in the y direction, orthogonal to both the vector current and J=ψ polarizations (see Table VI for the values used).

TABLE VI. Parameters for the correlation functions for the decay J=ψ → γη c . The second column gives the set of twists θ used for each lattice implementing twisted boundary conditions in the ð0 1 0Þ direction. The corresponding 3-momentum component is given by q y ¼ θπ=aN x [see Eq. ( 20)]. Given the x and z polarizations chosen for the J=ψ and vector currents, respectively, this gives a nonzero matrix element in Eq. ( 41). The twists are chosen to cover the kinematic range from zero recoil to q 2 ¼ 0. The third column gives the different time separation between source and sink (T ¼ t J=ψt η c , see Fig. 11) of the 3-point correlation function.

Taking the physical-continuum limit

We fit our lattice results for V, which we denote Vlatt (see Table VIII), to the following function:

sea;c δ sea;c þ κ ð0;kÞ sea;uds δ sea;uds × n 1 þ κ ð1;kÞ sea;uds ð ΛaÞ 2 þ κ ð2;kÞ sea;uds ð ΛaÞ 4 o  : ð49Þ

TABLE VII. Results for the masses of the J=ψ and η c mesons from our combined 2-and 3-point correlator fits aimed at determining the form factor for J=ψ → γη c decay. The η c meson here has spin-taste γ 5 ⊗ γ 5 and so has a different mass to the values given in Table IV. Set aM latt J=ψ aM latt η c 1 2.43418(38) 2.331959(86) 2 2.38671(15) 2.287729(36) 3 1.94525(21) 1.876352(68) 3A 1.92625(21) 1.857007(69) 4 1.90225(7) 1.833910(22) 5 1.41616(13) 1.366980(46) 6 0.92973(12) 0.896714(42) TABLE VIII. Our lattice results for Vlatt ðq 2 Þ on each ensemble from Table I and at each value of the momentum used.

The corresponding values of q 2 are also listed in lattice units. q 2 is calculated from

, where the ground-state J=ψ mass and η c energy, M latt

J=ψ and E latt η c , are obtained from the correlator fits and q 2 is determined from the twist angles imposed (see

Table VI). M latt J=ψ and M latt η c values are given in Table VII. Set 1 a 2 q 2 0.009220(74) 0.006762(74) 0.004301(74) 0.001843(73) Vlatt ðq 2 Þ 1.834(38) 1.834(23) 1.833(18) 1.831(15) Set 2 a 2 q 2 0.005680(28) 0.001564(28) -0.002551ð27Þ    Vlatt ðq 2 Þ 1.858(12) 1.8520(91) 1.8478(76)    Set 3 a 2 q 2 0.003956(27) 0.003165(26) 0.002374(26) 0.001583(26) Vlatt ðq 2 Þ 1.887(53) 1.881(38) 1.878(31) 1.875(27) a 2 q 2 0.000793(26) 0.000001(26) -0.000789ð26Þ -0.001580ð26Þ Vlatt ðq 2 Þ 1.873(25) 1.872(23) 1.870(21) 1.869(20) Set 3A a 2 q 2 0.004002(27) 0.003211(27) 0.002419(27) 0.001628(27) Vlatt ðq 2 Þ 1.888(54) 1.880(38) 1.877(32) 1.874(28) a 2 q 2 0.000837(27) 0.000045(27) -0.000745ð27Þ -0.001536ð26Þ Vlatt ðq 2 Þ 1.872(25) 1.870(23) 1.869(21) 1.867(20) Set 4 a 2 q 2 0.0020259(86) -0.0006190ð85Þ -0.0032636ð84Þ    Vlatt ðq 2 Þ 1.8558(47) 1.8530(34) 1.8501(29)    Set 5 a 2 q 2 0.001833(12) 0.001247(12) 0.000661(12) 0.000075(12) Vlatt ðq 2 Þ 1.854(14) 1.856(10) 1.8565(83) 1.8562(73) a 2 q 2 -0.000510ð12Þ -0.001097ð12Þ       Vlatt ðq 2 Þ 1.8555(66) 1.8547(62)       Set 6 a 2 q 2 0004578(76) -0.0001746ð75Þ -0.0008068ð74Þ    Vlatt ðq 2 Þ 1.873(13) 1.8723(93) 1.8704(79)   

This takes the same form as that for Fð0; 0Þ in Eq. ( 27), except that we must allow here for dependence on q 2 . We do this through a truncated Taylor series in q 2 =M 2 J=ψ , M J=ψ being the appropriate mass for a form factor induced by a cc vector current. The coefficients A ðkÞ for each term allow the fit to adjust the mass away from M J=ψ if required. The q 2 range is very small here (relative to M 2 J=ψ ), so we expect the form factor to be very flat as a function of q 2 and do not need many terms in the polynomial. We include terms up to q 4 . The J=ψ mass used in the q 2 =M 2

J=ψ term is that obtained from the fit (see Table VII). The fit form allows for independent discretization effects and quark mass mistuning terms for each power of q 2 =M 2 J=ψ (denoted by k). The form of these terms is the same as in the fit function for Fð0; 0Þ in Eq. ( 27) and definitions for δ can be found in Eqs. ( 28)- (30).

We take priors of A ðkÞ ¼ 2ð1Þ. The choice is informed by leading-order NRQCD where the rate for the radiative transition depends on the wave function overlap between J=ψ and η c modulated by a Bessel function. For small momentum transfer and ignoring relativistic corrections that generate spin-dependent differences in the wave function, this wave function overlap is 1. Then we expect Vð0Þ ≈ 2 from comparing Eq. ( 5) of [51] to Eq. ( 42) and taking m c ≈ M η c =2 ≈ M J=ψ =2. Our lattice QCD calculation includes relativistic effects fully so it will give a much more accurate result for Vð0Þ than this leading-order nonrelativistic argument. How different the result is, we will see below. The priors for A ðk>0Þ encompass q 2 dependence following a pole form ð1q 2 =M 2 J=ψ Þ -1 . As with the fit of Fð0; 0Þ in Sec. II B 6, we test the size of discretization effects using the Empirical Bayes approach. We do not expect large discretization effects here based on the nonrelativistic arguments above which tell us that discretization effects can only enter through the small momentum transfer between J=ψ and η c and through small spin-dependent effects on the wave function overlap. We find that the Bayes factor is maximized by Λ ¼ 0.12 GeV. This is close to the largest jqj value in the kinematic range of the decay and so seems a reasonable value to set the scale for discretization effects.

Following Sec. II B 6 we take the priors κ ði;kÞ aΛ ¼ 0ð1Þ, κ ðkÞ val;c ¼ 0ð1Þ (little difference is seen between our results on sets 3 and 3A), κ ðkÞ sea;c ¼ 0.0ð1Þ and κ ðj;kÞ sea;uds ¼ 0ð1Þ for j ¼ 0, 1, 2. Our preferred fit takes the number of discretization terms included, i max ¼ 3.

In the limit of vanishing lattice spacing and physical quark masses, the form factor is then given by

The value of M J=ψ used here is that from experiment [6]. Note that our c quark mass is tuned so that our J=ψ masses match that value [see Eq. ( 31)]. In particular the form factor at q 2 ¼ 0 needed to determine ΓðJ=ψ → γη c Þ via Eq. ( 42), Vð0Þ, is simply given by the parameter A ð0Þ . In Fig. 12, we plot our lattice results for Vlatt ðq 2 Þ from Table VIII against q 2 for each ensemble (upper plot). The results show little dependence on q 2 , as expected. Note that the statistical uncertainties on the points are smaller for values at low q 2 in comparison to those at zero recoil. The statistical noise of the η c correlators increases with spatial momentum, but this is offset by the increased value of V 00 so that its relative uncertainty falls. We also show the fit band in gray that corresponds to the fit of Eq. ( 49) evaluated in the continuum limit with physical quark masses [see Eq. ( 50)].

The lower plot of Fig. 12 shows the lattice results interpolated to q 2 ¼ 0 on each ensemble and plotted against the square of the lattice spacing. The interpolation in q 2 is FIG. 12. Upper plot: lattice results for Vðq 2 Þ plotted against q 2 . The different colors denote the different ensembles used (see Table I). The gray band gives the result, with 1σ error bars, from the fit to Eq. (49) evaluated in the continuum limit at physical quark masses. The dashed line corresponds to the maximum physical q 2 value, ðM J=ψ -M η c Þ 2 . Lower plot: lattice results on each ensemble interpolated to q 2 ¼ 0 using the continuum fit result. The gray band gives the result from the fit to Eq. ( 49) as a function of lattice spacing at physical quark masses. done using the continuum result for the q 2 dependence of V.

The figure shows, again as expected, very little dependence on the lattice spacing.

The value that we obtain from our fit in the continuum limit and at physical quark masses is

with an uncertainty of 0.4%. The fit has a χ 2 =dof value of 0.19. Our result is clearly distinguishable from the leadingorder NRQCD expectation of 2, indeed it differs by 7.2 (4)%. Before discussing additional systematic uncertainties that need to be included, we first discuss the stability of our result for Vð0Þ under changes to the parameters of our fits. We plot the impact on the value of Vð0Þ of changes to either our correlator and continuum/chiral fits in Fig. 13. Our preferred (base) fit described above is given on the left. Variations include dropping datasets and changing the priors by factors of 0.5 or 2.0. We see very little variation in the final answer under any of these variations, confirming that our result is robust.

Additional systematic uncertainties

In Sec. II B 7 we estimated the additional systematic uncertainty on Fð0; 0Þ for η c → γγ from missing quark-line disconnected correlation functions and from missing QED effects. Here we do the same for Vðq 2 Þ.

As discussed in Sec. II B 7, the missing disconnected correlation functions mean that there is a 7.3 MeV mismatch between the η c mass determined on the lattice in the continuum limit (tuning the c quark mass so that the J=ψ mass is correct) and that determined in experiment [1]. In Sec. III B 2 we described how leading-order NRQCD gives a result of 2 for Vð0Þ because it ignores the spin-dependent differences between the J=ψ and η c wave functions. Our results, using a fully relativistic approach, show that Vð0Þ differs from 2 by 7.2(4)%. However, the missing disconnected correlation functions mean that we are missing a small part of the effects that generate a difference between the J=ψ and η c (and result in their "wave function overlap" differing from 1). The 7.3 MeV shift is 6% of the mass difference between the J=ψ and η c mesons (the hyperfine splitting). We might therefore expect that the missing disconnected correlation functions could generate a shift in V of 6% of the 7.2% difference from 2, i.e. a 0.4% shift of V.

An additional effect to be considered is the identification of the value at q 2 ¼ 0. Because the lattice η c mass does not exactly match that in experiment, the q 2 ¼ 0 point will correspond to a slightly (6%) incorrect value for the η c spatial momentum, jqj. The q 2 dependence of V is so small, however, that shifting the q 2 value at which we determine V has negligible effect. Note also that our results on sets 3 and 3A for V (see Table VIII) show almost no difference when we change the η c mass by 1% (see Table VII). A 1% shift in η c mass is much larger than the 0.2% shift coming from missing quark-line disconnected correlation functions.

We conclude that a 0.4% systematic uncertainty in V from missing quark-line disconnected correlation functions is reasonable. This systematic uncertainty includes the effect of the fact that the c quarks in charmonium carry an electric charge. The hyperfine splitting was calculated in [1] in lattice QCD þ QED and the 7.3 MeV shift between the lattice and experimental hyperfine splittings quoted above includes this QED effect. The impact of final-state QED interactions should be negligible for ΓðJ=ψ → γη c Þ because there is no electric charge in the final state. We allow an additional uncertainty of Oðα=πÞ ¼ 0.2% for higher-order QED corrections when calculating ΓðJ=ψ → γη c Þ.

The systematic uncertainty from missing quark-line disconnected correlation functions discussed above will be largely independent of q 2 , since Vðq 2 Þ is such a flat function over the kinematic range of the decay. This means that it will cancel almost entirely in the ratio R eeγ of Eq. ( 44). When we come to consider ΓðJ=ψ → η c e þ e -Þ, however, we must allow a systematic uncertainty for finalstate QED interactions because of the charged particles produced. We will take an OðαÞ ≈ 1% systematic uncertainty for this.

C. Results

Combining our fit result of Eq. ( 51) with the additional systematic error discussed in Sec. III B 3 we obtain a final result of FIG. 13. The value of Vð0Þ in the limit of vanishing lattice spacing and physical quark masses obtained from variations to our base fit. These include (from left to right) dropping the coarsest and finest datasets, changing all the prior widths in our correlator fits, changing all the prior widths in our chiral/ continuum fits and adding an additional normal and oscillating exponential to our correlator fits. Note that the values for Λ in Eq. ( 49) are fixed (see text) under these fit variations.

Vð0Þ ¼ 1.8649ð73Þ fit ð75Þ syst :

The total uncertainty here is 0.56%. The error budget is given in Table V and can be compared to that previously discussed for Fð0; 0Þ. The main sources of uncertainty are somewhat different reflecting the fact that this is a dimensionless quantity, less sensitive to w 0 , but with larger statistical errors from fitting three-point correlation functions and larger uncertainties from a 2 and q 2 dependence. The relative total uncertainty is similar in the two cases.

We can use this value in Eq. ( 42) to determine the decay width. For the kinematic factors of masses on the righthand side of the equation we use experimental average values from [6]. These give jkj ¼ 110.9ð4Þ MeV. Taking α ¼ 1=137.036 [6], appropriate to the low moment-transfer here, we obtain ΓðJ=ψ → γη c Þ ¼ 2.219ð17Þ fit ð18Þ syst ð24Þ expt ð4Þ QED keV:

The third error here is from the experimental hyperfine splitting that appears in jkj. Since jkj is raised to the third power in Γ this gives the largest single uncertainty in the final answer. The fourth uncertainty is from additional QED effects in the rate, as discussed in Sec. III B 3. Our total uncertainty on Γ is 1.6%, adding the four contributions in quadrature. Taking the total J=ψ width as 92.6(1.7) keV from [6] gives a branching fraction of BrðJ=ψ → γη c Þ ¼ 2.40ð3Þ latt ð5Þ expt %:

We have combined the two lattice uncertainties (from the fit and from the additional systematics) into the first uncertainty here. The second comes from the experimental hyperfine splitting, the additional QED uncertainty and from the J=ψ total width, combined in quadrature. It is dominated by the error in the J=ψ total width. Turning now to J=ψ → η c e þ e -, we plot the expression in Eq. ( 45) for dR eeγ =dq 2 as a function of q 2 in Fig. 14, taking out the normalizing factor of 3π=α. The integrand is very singular at low q 2 values and requires some care to integrate. It is cut off at the lower kinematic end point, 4m 2 e . We find the integrated value to be (using α ¼ 1=137.036)

From Fig. 12, we know that the form factor V is very flat with respect to q 2 and we find that the integrand is virtually indistinguishable if we replace the factor of Vðq 2 Þ= Vð0Þ with 1.0. If we integrate Eq. ( 55) setting V to a constant we obtain a value of 0.006078421 (43) which differs from that above in Eq. ( 55) by 0.02%. Combining our results for ΓðJ=ψ → γη c Þ and R ee in Eqs. (53) and (55), respectively, we find the decay width

The first uncertainty comes from the combined "fit þ syst" lattice uncertainty on ΓðJ=ψ → γη c Þ from Eq. ( 53) and the second uncertainty from the experimental contribution to that. The third uncertainty is the additional systematic error from QED effects such as final-state interactions discussed in Sec. III B 3. Combining the value from Eq. ( 56) with the J=ψ total width gives a branching fraction BrðJ=ψ → η c e þ e -Þ ¼ 1.457ð16Þ latt ð15Þ QED ð31Þ expt × 10 -4 : ð57Þ

The third, experimental, uncertainty is dominated by that from the J=ψ total width.

D. Discussion

Figure 15 compares our final result from Eq. ( 52) for the form factor Vð0Þ needed to determine the rate for J=ψ → γη c decay to earlier lattice QCD calculations including different numbers of flavors of sea quarks. Our result is a big improvement in accuracy over the earlier calculations, as well as having the most realistic sea quark content with u, d, s and c quarks in the sea. We also show, as a shaded blue band, the value of Vð0Þ of 1.57 (18) inferred from the average branching fraction for J=ψ → γη c of 1.7(4)% and J=ψ width of 92.6(1.7) keV [6] using Eq. (42). We see good agreement between the lattice results but they are all higher than the value for Vð0Þ inferred from the experimental rate. The experimental average branching fraction has a large uncertainty, inflated by a scale factor of 1.5 because of poor agreement between experiments. This means that the FIG. 14. 3π=α × dR eeγ =dq 2 from Eq. ( 45) plotted against q 2 . Note the log scale on the y axis. The vertical dashed lines mark the kinematic end points of the integral at 4m 2 e and ðM J=ψ -M η c Þ 2 . The integrand for the case where Vðq 2 Þ is set equal to Vð0Þ, i.e. the form factor is taken to be completely flat in q 2 , is visually indistinguishable from this. tension with our lattice QCD result is 1.6σ where σ comes from the experimental value.

We show more detail of the experimental picture in Fig. 16. There we plot our result for ΓðJ=ψ → γη c Þ [from Eq. ( 53)] along with the three most recent experimental results, from KEDR [4], CLEO [3] and Crystal Ball [2]. Uncertainties quoted on the experimental values are combined in quadrature. We also include the PDG average value [6] which is an average of the branching fractions from CLEO and Crystal Ball (that we multiply by the average total J=ψ width). We see that the Crystal Ball result is the lowest, over 3σ below our value from lattice QCD. The KEDR result is much higher, 2σ above our value. The KEDR analysis includes some model dependence and they quote their result as Γ 0 γη c , which is the value we plot in Fig. 16. The CLEO result is in between and in good agreement (within 1.5σ) with our result.

Figure 17 shows a comparison between our lattice result for ΓðJ=ψ → γη c Þ and a selection (not intended to be exhaustive) of earlier theoretical results using different techniques. It is clear that lattice QCD is able to provide a much more accurate result for this decay width than previous approaches. A comparison is nevertheless useful to allow an assessment of these other approaches for use in cases that are not as amenable to lattice QCD calculations. Indeed here the experimental picture is not at all clear and lattice QCD results can substitute for an accurate experimental value in this assessment, as for Γðη c → γγÞ in Sec. II D.

We see in Fig. 17 that theory results for ΓðJ=ψ → γη c Þ have covered a wide range (compared to the lattice QCD results of Fig. 15). A traditional approach has been to use a nonrelativistic potential, but often these results are quoted with no error bars. In Fig. 17 we give two results from nonrelativistic potentials. One, from [51], takes the wave function overlap between J=ψ and η c to be 1, and the FIG. 15. A comparison of values of Vð0Þ from lattice QCD. The result from the work here is labeled "HPQCD23" (red asterisk) and includes u, d, s and c quarks in the sea with results at multiple lattice spacing values. The results with u, d and s quarks are denoted with purple filled triangles. "HPQCD12" used HISQ c quarks on gluon field configurations including sea asqtad staggered quarks and two values of the lattice spacing [10]. "Hadspec23" used clover quarks on anisotropic lattices at one value of the lattice spacing but include an estimate of systematic errors in their quoted uncertainty. "ETM12" (filled orange circle) used the twisted mass formalism with gluon field configurations including u and d sea quarks only and four values of the lattice spacing [9]. The blue band shows the value for Vð0Þ inferred from the average experimental branching fraction [6]. The red band carries our result down the plot for comparison. FIG. 16. A comparison of our result for ΓðJ=ψ → γη c Þ from lattice QCD to values from experiment. The result from the work here is labeled "HPQCD23" (red asterisk). The filled blue circles are results from individual experiments: "CBALL85" is from the Crystal Ball [2], "CLEO08" is from CLEO [3] and "KEDR14" is from KEDR [4] (plotting the quantity denoted Γ 0 γη c ). The blue band shows the average of the CLEO and Crystal Ball results [6]. The red band carries our result down the plot for comparison. FIG. 17. A comparison of our result for ΓðJ=ψ → γη c Þ from lattice QCD to a selection of values from other theoretical approaches. The result from the work here is labeled "HPQCD23" (red asterisk). The left-pointing green triangle labeled sum rules is an early result with this method from [52]. Three values are given in the row labeled "potentials" to encompass the range of results. The rightmost point (green triangle) is from [51], the middle point (green circle) is the starting point for a pNRQCD analysis from [53]. The numbers plotted for these two cases have been updated to use the current value for jkj (see text). The leftmost point (green square) is from a relativistic quark model [54]. Two pNRQCD analyses are shown; the lower value (green diamond) is from [55] using the weak-coupling limit and the upper one (green circle) from [53] including corrections through α 2 s and v 2 (this value has also been corrected to use the current value of jkj). The right-pointing green triangle labeled "LCSR" is a recent result using light-cone sum rules from [56]. The blue band shows the experimental average of the CLEO and Crystal Ball results [6]. The red band carries our result down the plot for comparison. other, from [53], is the leading-order result for a potential Non-Relativistic QCD (pNRQCD) analysis. To make a fair comparison with our value, we have corrected both of these results to use the value for jkj (see Eq. ( 43) obtained from current average experimental masses [6] rather than the values they used from earlier PDG reports. Since the experimental average hyperfine splitting has changed by several percent over time and jkj appears cubed in Γ, this has some impact. These two nonrelativistic potential model results differ by 20%, bracketing our value, and this reflects reasonably the range of results (compare, for example, more recent values in [57]). The pNRQCD approach systematically adds correction terms to this [53,55]. We can see how this works most clearly in the result of [53] since the size of corrections through Oðα 2 s Þ and Oðv 2 Þ are tabulated. The v 2 corrections are large but of opposite sign to those at Oðα s Þ and Oðα 2 s Þ. Disappointingly we see, by comparing the two green circles in Fig. 17, that the net effect of these corrections is to move the result further from our value rather than towards it. Higher-order corrections are still likely to be sizeable and the good news is that the uncertainty estimates attached to the pNRQCD results means that they are in good agreement with our value. In contrast a recent result from light-cone sum rules has significant tension, over 3σ, with our value. The result from a relativistic quark model [54] also looks in disagreement.

We now turn to a test of the relationship between ΓðJ=ψ → γη c Þ, Γðη c → γγÞ and ΓðJ=ψ → e þ e -Þ suggested by Shifman many years ago and given in Eq. ( 2). The accurate results that we now have from lattice QCD for these decay widths enables us to see how well this approximate relationship works. It is easiest to do this by converting the expression of Eq. ( 2) into a connection between the hadronic parameters, Vð0Þ, Fð0; 0Þ and f J=ψ . Equation (2) becomes

In terms of the ratio R fF defined in Eq. ( 39) this reduces to

Our results for Vð0Þ from Eq. ( 52) and R fF from Eq. ( 40) yield

The expectation from the ratio of masses on the right-hand side of Eq. ( 59) gives 0.982 using experimental averages from [6]. This differs from our result for Vð0ÞR fF by 10%, which is well within the leeway provided on the 0.982 by possible Oðα s Þ corrections. In the nonrelativistic limit, the mass ratio would simply be 1.0, and 0.982 is a slight improvement on this, in terms of being closer to the full lattice QCD value that we obtain in Eq. ( 60). Going further in the nonrelativistic direction, Eq. ( 59) would reduce to the expectation that Vð0Þ ¼ 2 (when R fF is assumed to take the value 1=2), which in fact works just as well in comparison to our results (since 2.0 differs from our value for Vð0Þ by 7%). We conclude from the comparison of our results to other theory values that nonrelativistic approaches to J=ψ → γη c work at the 10-20% level at leading order, but it is hard to improve on this by adding corrections.

Finally, we note that our result for R eeγ in Eq. ( 55) agrees with that given in [58] using a simple ψ 0 pole model for Vðq 2 Þ. This is not surprising since, as discussed in Sec. III C, the result is insensitive to details of Vðq 2 Þ over the short q 2 range of the decay.

IV. CONCLUSIONS

We give improved lattice QCD results for the hadronic matrix elements needed to determine Γðη c → γγÞ, ΓðJ=ψ → γη c Þ and ΓðJ=ψ → η c e þ e -Þ, These results include u, d, s and c quarks in the sea for the first time. The HISQ action gives us small discretization errors and we use a wide range of lattice spacing values and sea u=d quark mass values for good control of the physical-continuum limit. Our lattice QCD results now have smaller uncertainties than the corresponding experimental values. This means that new experimental results with improved uncertainties could have considerable impact as stringent tests of QCD.

Our results for η c → γγ transform the theoretical picture for this decay. Our determination of the hadronic form factor Fð0; 0Þ, defined in Eq. ( 16), is [repeating Eq. ( 34)] Fð0; 0Þ ¼ 0.08793ð29Þ fit ð26Þ syst GeV -1 :

This gives a decay width of Γðη c → γγÞ ¼ 6.788ð45Þ fit ð41Þ syst keV;

repeating Eq. ( 35). The first uncertainty comes from the lattice calculation and the second is from remaining systematic errors (see Sec. II B 7) from missing quark-line disconnected diagrams and QED effects.

As discussed in Sec. II D, our result has 4.6σ tension with the PDG fit result of 5.15 (35) keV [6]. The PDG fit has a poor χ 2 and the difficulty of determining a reliable fit value and uncertainty from the wide range of indirect experimental results that exist is clear in Fig. 8. We believe that this fit needs to be revisited. Instead our result for the decay width agrees within 2σ with the value for the width of 5.90 (58) keV obtained from the PDG average [6] of a more restricted set of experimental results, those for η c production via 2-photon fusion using the η c decay mode to K Kπ. This channel also has the advantage of having the smallest relative uncertainty.

Our results, along with earlier 0.4%-accurate HPQCD calculations of the J=ψ decay constant [1], allow us to test how well earlier expectations using NRQCD work when confronted with the results from a fully relativistic QCD calculation. Leading-order NRQCD gives a simple relationship between Fð0; 0Þ for η c → γγ decay and the J=ψ decay constant given in Eq. (38). We determine the ratio [repeating Eq. ( 40)]

The leading-order NRQCD expectation of 0.5 is not far from this number, suggesting that there is significant cancellation of the NRQCD higher-order corrections, expected to be individually of order 30% [16,48,49]. The comparison of NRQCD calculations to our result is shown in Fig. 9.

The relative success of the leading-order NRQCD expectation for the ratio above suggests that it may also be used to predict η b → γγ. We therefore expect that

Using HPQCD's lattice QCD results for Γðϒ → e þ e -Þ of 1.292 (37)(3) keV [59] we obtain the prediction

The second uncertainty allows for 20% variation from missing higher-order radiative and relativistic corrections.

For the b case we expect the Oðα s ≈ 0.2Þ corrections to be larger than the relativistic corrections [Oðv 2 ≈ 0.1Þ] and the cancellation of these may then not work as well as in the c case. Indeed the authors of CM01 [16] are quoted in [60] as determining a width Γðη b → γγÞ ¼ 570ð50Þ eV including corrections through α 2 s in continuum NRQCD to Eq. (64). A more recent analysis [49] gives a similar central value but larger uncertainty, with a decay width of 540(150) eV. This would represent a larger correction to the leading-order NRQCD result than seen in the c case but going in the same (positive) direction. We saw from Fig. 9, however, that the full lattice QCD result for Γðη c → γγÞ is below the LO NRQCD result rather than above.

Ultimately the Standard Model (SM) value for Γðη b → γγÞ will be resolved by an accurate calculation in lattice QCD. Now that we have demonstrated the accuracy possible for the η c calculation using HISQ quarks we envisage increasing the mass up towards the b and applying the techniques used in [59] to do this. A prediction ahead of possible experimental results from Belle II would be timely.

Our results for the form factor for the M1 radiative transition J=ψ → γη c also represent a step up in accuracy over previous lattice QCD results. We find a value for the form factor at q 2 ¼ 0, repeating Eq. ( 52),

This gives a decay width, repeating Eq. ( 53)

ΓðJ=ψ → γη c Þ ¼ 2.219ð17Þ fit ð18Þ syst ð24Þ expt ð4Þ QED keV:

The first uncertainty here comes from the lattice fit, the second from additional systematic errors from missing quark-line disconnected correlation functions, the third uncertainty is from the experimental average hyperfine splitting that enters the kinematic factors converting the squared form factor into a decay width [see Eq. ( 42)] and the fourth from additional QED effects in the rate. For the related Dalitz decay, J=ψ → η c e þ e -we predict ΓðJ=ψ → η c e þ e -Þ ¼ 0.01349ð15Þ latt ð15Þ expt ð13Þ QED keV; ð68Þ repeating Eq. ( 56).

Our result for ΓðJ=ψ → γη c Þ joins earlier lattice QCD values in being higher than the experimental average value [6] obtained from averaging results from Crystal Ball [2] and CLEO [3]. The small uncertainty of our result makes the tension more compelling than before. Our result is in 3σ tension with that from Crystal Ball but agrees within 1.5σ with that from CLEO. See Fig. 16 for the comparison.

We can use our results to calibrate other theoretical approaches (see Fig. 17) and also to test the suggested [18] simple relationship between Vð0Þ, Fð0; 0Þ and f J=ψ , by determining [repeating Eq. ( 60)]

This is to be compared with 0.982 from Eq. ( 59), or 1.0 in the nonrelativistic limit. Once again we find that the simple expectations work fairly well, at the ∼10% level, but this does not of course provide the kind of accuracy available now from lattice QCD, as we have shown here. Finally, for an alternative perspective on the experimental situation, we multiply together the branching fractions we have determined for J=ψ → γη c and η c → γγ [in Eqs. (36) and (54)] and compare to direct experimental determination of this product by CLEO [43] and BESIII [39]. This comparison is shown in Fig. 18. We see good agreement between the lattice QCD result and the experimental values, particularly that from BESIII. The experimental determinations have large uncertainties at present but improvements here to provide more stringent tests against lattice QCD would be very useful, as we have stressed throughout this paper.

APPENDIX A: TRAPEZOIDAL INTEGRATION AND OSCILLATING CONTRIBUTIONS

Here we discuss the accuracy of representing the integral over t γ 1 in the continuum with a sum over lattice times in Eq. ( 9) and the impact on this of the presence of oscillating terms arising from the use of a staggered quark formalism.

We want to approximate the integral

with a sum

that includes a term oscillating in time through the factor ð-1Þ n . f and g are smooth functions, continuous at t ¼ 0, that vanish in a well-behaved way as t → ∞.

For the sum over f Eq. (A2) is using the standard trapezoidal rule, which has a 2 errors. For the sum over g we combine three adjacent terms, reducing the sum to even values of n only, a X n¼þ∞;n even n¼-∞;n even

Splitting the sum into two pieces, for positive and negative n, we have for positive n a X n¼þ∞ n¼0;2;4

if g and g 0 vanish at t → ∞. Negative n gives a similar result, with opposite sign, so that we have a total for the sum over g in Eq. (A2) of

The result is a discretization effect, proportional to the discontinuity in the derivative of g at t ¼ 0 but vanishing as a → 0 as a 2 . We conclude that summing over lattice time slices to obtain Cμν [Eq. ( 9)] as an approximation to the time integral that sets the photon on shell introduces discretization errors proportional to a 2 at leading order. These come both from the trapezoidal integration implied by the sum and from the oscillating terms in the correlation function that arise from the use of staggered quarks. Such discretization errors are taken into account by our fit to the results for the form factor as a function of lattice spacing and removed in our result for Fð0; 0Þ in the continuum limit.

A toy model illustrates this further. We take f and g to be single exponentials, f ¼ expð-M n jtjÞ and g ¼ expð-M o jtjÞ. For f we have FIG. 18. A comparison of our result for the product of branching fractions BrðJ=ψ → γη c Þ × Brðη c → γγÞ from lattice QCD to values from experiment. The result from the work here is labeled "HPQCD23" (red asterisk). The filled blue circles are from CLEO [43] (labeled "CLEO08") and BESIII [39] (labeled "BESIII13"). The blue band shows the average of the CLEO and BESIII results [6]. The red band carries our result down the plot for comparison.

Z þ∞ -∞ dt e -M n jtj ≈ a X þ∞ j¼-∞ e -aM n j ;

on the lattice. The left-hand side of Eq. (A6) evaluates to 2=M n . The sum on the right-hand side is a geometric series, so we have

As expected the integral obtained in this manner is accurate up to ðaM n Þ 2 errors. We repeat this exercise for the oscillating exponential to obtain

As expected the impact of the oscillatory contributions vanish in the continuum limit as a 2 and the result matches that from Eq. (A5). We reach the same conclusions about discretization error from the toy model as from the more general f and g functions of Eq. (A2).

APPENDIX B: FITTING A SUBSET

OF η c → γγ DATA As remarked in Sec. II B 5, it is possible to obtain accurate results for the η c → γγ form factor with a lot less numerical work than we have expended here. We showed, in Fig. 3, that our results converge very rapidly to their final value a function of the t width region over with the 3-point function is integrated to obtain Cμν . Here we provide another test that significantly reduces the amount of computation needed.

Figure 19 shows the results obtained if we restrict the fit of the 2-point function Cμν ðtÞ to a set of specific t ≡ t η ct γ 2 values rather than fitting the full t range. The left-hand point shows the full fit that we use here and the right-hand points show the results from selecting 4, 5 or 6 specific t values. We need a mix of even and odd t values for an optimal fit because of the oscillating terms from opposite parity states [see Eq. ( 24)]. We must also adjust the separation in lattice units between t values as the lattice spacing changes. On the very-coarse lattices a separation of 1 is appropriate, whereas for coarse and fine lattices a separation of 3 or 5 gives a better range of t to reproduce the full fit. We conclude from Fig. 19 that we could obtain similar uncertainties to our full fit by calculating correlation functions at only 4 or 5 time separations t η ct γ 2 if these were well chosen. This is possible because lattice data is correlated as a function of t. Usually lattice 2-point functions are calculated at all t values because there is no significant time saving in making a t selection. Here, because we actually calculate a 3-point function, working with selected t values reduces the computational cost considerably and we will make use of this in future. Here, however, we use the results from our full fit. FIG. 19. Fitted results for F latt ð0; q 2 2 Þ, comparing values obtained from our full fit with those that use a subset of time separations, t, between source and sink in Cμν ðtÞ. Results are given for the LOCAL setup on set 1 (very coarse, a ≈ 0.15 fm), set 3 (coarse, a ≈ 0.12 fm, fit simultaneously with set 3A) and 5 (fine, a ≈ 0.09 fm); all values are in lattice units. The values on the left give results from the full fit to all t values except those discarded at early times (below t min , see text in Sec. II B 5). The values on the right give fit results for a selected 4, 5 or 6t values. We see that uncertainties close to that of the full fit are possible even with 4t values. Even and odd t values are needed and the spacing between t values and range covered is adjusted as a function of lattice spacing.