Lattice calculation of $\chi_{c0} \rightarrow 2\gamma$ decay width

We perform a lattice QCD calculation of the $\chi_{c0} \rightarrow 2\gamma$ decay width using a model-independent method which does not require a momentum extrapolation of the corresponding off-shell form factors. The simulation is performed on ensembles of $N_f=2$ twisted mass lattice QCD gauge configurations with three different lattice spacings. After a continuum extrapolation, the decay width is obtained to be $\Gamma_{\gamma\gamma}(\chi_{c0})=3.65(83)_{\mathrm{stat}}(21)_{\mathrm{lat.syst}}(66)_{\mathrm{syst}}\, \textrm{keV}$. Albeit this large statistical error, our result is compatible with the experimental results within 1.3$\sigma$. Potential improvements of the lattice calculation in the future are also discussed.


I. INTRODUCTION
Charmonium physics lives in an energy regime where both perturbative and nonperturbative features of quantum chromodynamics (QCD) intertwine.Notably, Charmonium decay has played an important role in establishing the asymptotic freedom of QCD and served as a clean platform to probe the interplay between perturbative and nonperturbative dynamics.In particular, the two photon annihilation rates of charmonium are extremely helpful for the understanding of quark-antiquark interaction and the decay mechanisms [1,2].
In this paper, we study the two-photon decay width of χ c0 , which has been extensively studied from both experimental and theoretical sides.On the experimental side, using the decay of ψ(3686) → γχ c0 , χ c0 → γγ, both CLEO-c and BESIII collaborations reported results of the two-photon decay width Γ γγ (χ c0 ) [3,4]: (35) stat (22) syst keV Γ BESIII γγ (χ c0 ) = 2.03 (8) stat (14) syst keV Γ P DG γγ (χ c0 ) = 2.20 (22)keV (1) where the first line from CLEO-c, the second from BE-SIII and the last line being the PDG quoted value with combined errors.It is expected that more accurate results for these decay width will become available in the near future.
On the theoretical side, it is fair to say that the situation is far from satisfactory.Theoretical results for the decay rate have been obtained using a non-relativistic approximation [5,6], potential model [7], relativistic quark model [8][9][10][11], non-relativistic QCD (NRQCD) factorization [12][13][14][15][16][17][18], effective Lagrangian [19], Dyson-Schwinger equations (DSEs) [20], as well as quenched [21] and unquenched lattice calculations [22].These results are listed in Table I, which scatter quite a lot although all fall in the right ballpark.Note that within the framework of NRQCD, the leading-order (LO) prediction is close to the experimental measurements, but this process is extremely sensitive to high-order QCD radiative corrections and relativistic corrections.Therefore, only the LO predictions are listed in Table I .
In the last line of Table I, we list two existing lattice QCD results so far.The first one from Dudek et al is a quenched lattice computation on a single lattice spacing [21].The systematic error they quote mainly come from quenching.The second one from CLQCD is an unquenched study using N f = 2 twisted mass fermions at two distinct lattice spacings.The authors found that the lattice artifacts are substantial and only quoted results from a finite lattice spacing, without an error estimate of the finite lattice spacing errors.The number quoted in Table I is the result from the finer lattice spacing [22].Therefore, in both lattice studies, systematic effects such as finite lattice errors are not fully investigated which was found to be large in the second study [22].Obviously, in order to fully compare with the upcoming experiments, one needs to work in a theoretical framework that allows an improvable error control and in this respect, lattice computation obviously has an advantage over other phenomenological methods listed in Table I.
In this paper, we try to improve on the existing lattice computation of Γ γγ (χ c0 ) in two major aspects: First, in previous lattice studies, many systematic effects are not yet fully taken into account, the most important of which being the finite lattice spacing effect, which has been observed in Ref. [22].Second, one normally computed the off-shell form factors at various discrete photon virtualities.In order to obtain the physical decay width, an extrapolation of these results are required, introducing a model-dependent systematic error.
In this work, we have made the following improvements: First, to attack the lattice artifacts, we perform our calculation on ensembles with three different lattice spacings, allowing us to perform a reliable continuum extrapolation.Second, we adopt a novel method to extract the on-shell form factor directly, by-passing the conventional momentum extrapolation and therefore avoids the corresponding model-dependent extrapolation errors.We have also taken the excited-state contamination into consideration, further improving our results on the physical form factor. Similar procedures has been successfully utilized to two-photon decay of η c [26].We hope that these improvements could also shed some light on the two-photon decay of χ c0 .
This paper is organized as follows.In Sect.II, the methodology for extracting on-shell form factor is introduced.The Sect.III is divided into several parts: In Sect.III A, the information of the configurations and operators used in this work are introduced.In Sect.III B, the mass spectrum of χ c0 is presented.In Sect.III C, we give the renormalization factor and the spectrum weight factor.In Sect.III D, numerical results of the form factor in three different lattice spacings are presented.Then in Sect.III F, extrapolation of the results to continuum is performed, yielding our final result for the decay rate.We also compare our result with both experimental and theoretical results.The main sources of error in our work are discussed and possible solutions in the future are proposed.

II. METHODOLOGY
In this section, we outline the methodology for the calculation of two-photon decay width of χ c0 .In the traditional approach [27], using the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula and integrating out the QED part to O(α em ), the amplitude for twophoton decay of charmonium can be obtained as follows [21], where ϕ M ( x, t f ) is an appropriate composite operator which creates a desired meson M (in our case, the χ c0 meson) from the QCD vacuum; µ , ν are the polarization four-vectors for the two final photons; J µ = q e q qγ µ q (e q = 2/3, −1/3, −1/3, 2/3 for q = u, d, s, c) is the electromagnetic current operator due to the quarks, with e being the elementary charge unit.In this work, we only consider the connected contributions arising from the charm quark current.Disconnected contributions are neglected.These contributions are extremely costly for lattice computations and are assumed to be small in charmonium physics [21,28,29].Then, the matrix element in Eq. ( 2) relevant for χ c0 decay can be parameterized in terms of the form factor (3) where q 1 , q 2 are the two four-momenta of the final photons while 2 are the virtualities of the two photons.The mass of χ c0 is denoted as m χ and the polarization vectors of the two photons are given by 1 and 2 .The physical decay width is related to the on-shell form factor which is obtained by a momentum extrapolation towards the physical point: Thus, in this conventional approach, in order to have a better control on the extrapolation, one needs to compute the matrix element at various different non-physical virtuality combinations, thereby also introducing extra computational costs.The extrapolation itself also brings about model-dependent systematic errors.In the new approach introduced in this work, we adopt a method that requires no off-shell form factor calculations at all and therefore by-passing the model-dependent extrapolation in photon virtualities.The method has been successfully utilized in two-photon decays of η c [26].We now briefly outline the major steps for the case of χ c0 below.
One first relates the on-shell decay amplitude of χ c → 2γ to an infinite-volume hadronic tensor F µν (p) which is the Fourier transform of the real-space tensor H µν (t, x) in continuum Euclidean space, where we have chosen the rest-frame of the χ c0 meson so that k = (im χ , 0).Note that we have fixed the four-momentum for one of the final photons to be p = (im χ /2, p) with | p| = m χ /2, making it on-shell explicitly and energy-momentum conservation then guarantees the other photon with four-momentum p is also on-shell.With this choice, the on-shell decay amplitude may be written as, According to the quantum number of χ c0 , the hadronic tensor can be parameterized as (repeated indices are summed), The approach to extract the on-shell form factor F χc0γγ here is also slightly different from the conventional one.
By further multiplying the Lorentz structure factor in the above equation, the hadronic tensor can be contracted to a scalar including only the form factor F χc0γγ with a constant factor.Then the form factor can be derived by dividing the coefficient as follows, Until now, all derivations are in the continuum Euclidean space.We now utilize the spatial isotropy symmetry to average over the spatial direction of p, where j n (x) are the spherical Bessel functions.Finally the scalar from factor is expressed as where i = 1, 2, 3 take spatial indices and are assumed to be summed over.
To obtain the hadronic tensor H µν (t, x) in Eq. ( 9), we utilize the variational method to find the optimal interpolation operators to create the χ c0 meson state [30].The physical decay width of χ c0 is given by Therefore, one only needs to compute the Euclidean correlation functions H 0i and H ii that are directly relevant for the on-shell amplitude and substitute the results into Eq.( 9) to arrive at the physical decay width Γ γγ (χ c0 ) in Eq. (10).This completely avoids the on-shell extrapolation process in the conventional lattice approach.

A. Lattice setup
We utilize three N f = 2-flavor twisted mass gauge field ensembles generated by the Extended Twisted Mass Collaboration(ETMC) with lattice spacing a 0.0667, 0.085, 0.098 fm, respectively.The parameters of these ensembles are presented in Table .II.The valence charm quark mass parameter µ c is tuned so that the mass of the η c meson for each ensemble reproduces its correct Table II.From left to right, we list the ensemble name, the lattice spacing a, the spatial and temporal lattice size L and T , the number of the measurements of the correlation function for each ensemble N conf ×T , the light quark mass aµ l , the pion mass mπ and the range of the time separation t h between χc0 and photon.
physical value.For more details, we refer the reader to Ref. [31,32], Before getting into the simulation details, there remains one subtlety to clarify that is related to the twisted mass fermion.Since the twisted mass action breaks parity P by O(a 2 ) effects, the basis operator O 1 = cc for χ c0 would unfortunately mix with O 2 = cγ 5 c which has the opposite parity.This mixing implies that a specific combination of these operators will be relevant to create a physcial scalar charmonium in the twisted mass action [30] The two-point correlation function |0 can be derived by multiplying the corresponding coefficients with the basis correlation functions . Therefore, after choosing a time slice t 0 , one could disentangle the mixing of the two operators by solving a generalized eigenvalue problem (so-called GEVP procedure): where the generalized eigenvalues λ i behave like e −Ei(t−t0) at large time separation.In practice, we fix t 0 = 1 and solve Eq. ( 12) on each time-slice independently and use them to reconstruct the three-point correlation functions.

B. Mass spectrum for χc0
Since the generalized eigenvalues in Eq. ( 12) decay exponentially, the corresponding mass eigenvalues can be extracted easily from, Since we want to extrapolate the form factor to eliminate the excited state contamination, we therefore use the following two-state fit form for the χ c0 correlator, with V being the spatial volume, m 0 the ground state mass and m 1 the first excited state mass.The factors (with i = 0, 1) are the overlap amplitudes for the ground and the first excited state, respectively.The corresponding mass plateaus and the masses are illustrated in Fig. 1 for the three ensembles we utilize in this work.The left column of the panels show the effective mass on each time slice.The right panels denote the mass values fitted from two-point correlation functions, the upper one for the first excited state and the bottom one for the ground state.As the grey bands in the left panels indicate, the starting time slices are adjusted according to χ 2 /d.o.f of the fit while the ending time slices are fixed to be t max = 27, 20, 20 for ensemble I, II, III, respectively.Noting that the grey band of Ens.I is obviously different from the other two ensembles, the ground state mass m 0 might be underestimated.So we calculate the result for another plateau with green mark and take the difference of them as the major source of systematic uncertainty.
The results for the mass values are summarized in Table III.Note that we use the η c mass to fix the valence charm quark mass aµ c in this work.And the χ c0 experiment mass is 3414.7(3)MeVquoted by PDG [33].The hadronic tensor H µν contains the electromagnetic current operators J µ from all flavor of quarks.However, since we neglect the disconnected diagrams in this study, we only need to consider the charm quark current J (c) µ = cγ µ c( x, t).Since we adopt the local current form, there exists an extra multiplicative renormalization factor Z V that can be calculated by a ratio of the two-point function and the three-point function as in Eq. ( 15).In principle this renormalization factor does not depend on the particle state used to calculate it.For a better signal, we choose to use the η c correlators instead of χ c0 .Taking account of the around-of-world effect, we use the following relation to extract Z V .
The results for Z V are listed in Table IV.
When computing the scalar form factor F χc0γγ in Eq. ( 9) on the lattice, the integration over space-time are replaced by discrete summations.When two identical currents in H µν (t, x), meaning that they share the same Lorentz index, are at the same space-time point, an extra renormalization is needed to take into account of the contact term.This is due to a new type of composite operator that is not properly renormalized yet, even if each current is already properly renormalized by the factor Z V .In order to take this effect into account, one needs to impose another appropriate renormalization condition for this new composite operator.In this work, we choose not to sum the same space-time point contributions for identical currents and thereby avoiding this potential renormalization.To summarize, the above mentioned procedures taken on already O(a)-improved ensembles will at most introduce an extra O(a 2 ) lattice artifact on physical observables which will be taken care of in the final continuum extrapolation.

D. The scalar form factor Fχ c0 γγ
When computing the hadronic tensor, we evaluate the three-point correlation function J µ (x)J ν (0)O † χc0 (−t h ) .To produce the static meson state, we use the Z 4stochastic wall source placed at time-slice −t h .This cuts the uncertainty by nearly a half when compared with the simple point source for the meson mass.We also apply the APE [34] and Gaussian smearing [35] for the gauge field and χ c0 operator.We utilize the random point source propagator for the current to arrive at the threepoint correlation function.In practice, the hadronic tensor with current J ν (0) placed at zero point is actually an average of all the time slices and a random positions on each time slice.
Consequently, the scalar form factor we computed according to Eq. ( 9) on the lattice F χc0γγ actually suffers from excited state contamination due to higher excitation states of χ c0 .What we really need is the ground state χ c0 .This effect can be taken care of by considering t h dependence of the form factors. Therefore, we computed several different separations t h and perform the following fit, where F χc0γγ and ξ are the two free parameters.For the parameters m 0 and m 1 , we take the values presented in Table III.The form factor with different time separation t h together with the ground state extrapolation values for F χc0γγ for three set of ensembles are illustrated in Fig. 2.

E. Comparison of the form factor with previous lattice results
The most recent lattice computation on χ c0 → γγ decay in the literature is the one from CLQCD [22], which happened to use exactly the same set of ensembles as this work.This allows a more detailed comparison on the level of dimensionless form factors for each of the common lattice spacings.For this purpose, we decide to convert our results for F χc0γγ into dimensionless quantities which could be taken as either Γ γγ (χ c0 )/m χc0 or the dimensionless form factor G(0, 0) that is utilized in Ref. [22].The relation between these two dimensionless quantities is easily found to be, (17) Dimensionless quantities have the advantage that they are independent of the scale setting process for the lattice spacings, which is subject to its own errors depending on how the scale was set.Since the scale setting processes for lattice calculations also progress over the years, the information for the lattice spacing in physical units, both the central values and the errors, are also changing with time even for a given particular ensemble.It is therefore better to attach these errors due to scale-setting at the very end when comparing with the experiments.In the intermediate step when comparing with other lattice computations, it is easier to directly compare the dimensionless quantities if possible.In fact, this allows us to compare with previous lattice results in Ref. [22] at each individual lattice spacing, namely Ens.I and Ens.II that have also been utilized.Of course, when quoting the final physical decay width, the effect of scale setting will be taken into account together with its associated errors.
In Table V, the dimensionless form factor G(0, 0) obtained via Eq.( 17) from F χc0γγ for all three ensembles are listed together with the corresponding results for Ens.I and Ens.II from Ref. [22].Ens.III was not utilized in the study of Ref. [22].Two entries for Ens.I, labelled as Ens.I(a) and Ens.I(b) corresponds to the two different ways of extracting χ c0 masses as discusses in Fig. 1.The errors quoted for G(0, 0) in this work are obtained using the conventional jackknife method.As for the three errors for the results from Ref. [22], they stand for errors from statistical, from momentum extrapolations and estimates of the finite lattice spacing errors, respectively.V. Dimensionless form factors G(0, 0) obtained in this work and those obtained in Ref. [22] for each ensemble.Ens.I(a) and Ens.I(b) denotes two different results obtained by taking two different χc0 mass values as discussed in Fig. 1.Ensemble III was not available in Ref. [22].Errors quoted for G(0, 0) in this work are purely statistical that are obtained using the conventional jackknife method.Three errors for the results from Ref. [22] stands for errors from statistical, from momentum extrapolations and estimates for the finite lattice spacings.
We notice that the central values for dimensionless form factors G(0, 0) differ by almost a factor of two for Ens.I and Ens.II.The reason of this apparent discrepancy is still unknown to us.One possibility could be the under estimation of the lattice artifacts for each of the ensemble in Ref. [22].

F. Continuum extrapolation and the final result and discussions
After obtaining the dimensionless form factors for three different lattice spacings, we could investigate the continuum limit of this quantity.For this purpose, we decide to perform this extrapolation using the more physical quantity Γ γγ (χ c0 )/m χc0 , which is proportional to the norm-squared dimensionless form factor |G(0, 0)| 2 as indicated in Eq. ( 17).The continuum extrapolation is done by performing a linear fit in a 2 for the three ensembles and the result after the continuum extrapolation, together with the results for each ensemble, are illustrated in Fig. 3.Here the horizontal error bars for the data points indicate the errors in a 2 inferred from Ref. [31,32] It is seen that the three data points fit nicely on a straight line yielding a reasonable χ 2 /d.o.f .The two points near a 2 = 0 with larger error bars designate two different results obtained from the fit with and without considering the horizontal a 2 -errors for the lattice spacings.Below the two data points, we also plot the corresponding experimental value from PDG for this ratio.The two extrapolated results share almost identical central values.They differ only by their errors.The point with slightly larger error (the one slightly to the left) is the one that takes into account of the horizontal a 2errors while the other one is the one without considering a 2 -errors.Finally, there is another source of systematic errors arising from the different plateaus in the mass as discussed in Sec.III B. Therefore, we finally quote the result of the decay width in physical units as Γ γγ (χ c0 ) = 3.65(83) stat (21) lat.syst (66) syst keV, (18) where the first two errors represent the error obtained without/with the a 2 -errors.It should be interpreted as follows: the first error is the error without considering a 2 -errors.The second one with the subscript lat.syst indicate the extra amount of error if one would consider the a 2 -errors.In other words, one could add the first two errors in quadrature to obtain the error with a 2errors taken into consideration, which is shown by the left point near a 2 = 0 in Fig. 3.The last error with subscript syst reflects the systematic error from different mass plateaus in Ens.I(a) and Ens.I(b).In this manner, we separate different sources of systematic errors that have been studied in this paper.
It is evident that the central value for the decay width obtained in this paper is larger than the PDG value.But due to our large statistical and systematical uncertainty, it is still compatible with the experimental results within 1.3σ.
We have tried to estimate the systematic uncertainties that might influence our final result quoted above in Eq. ( 18).This includes choosing different plateaus for the mass, for the renormalization factor Z V , spectral weight factor Z i , and the number of time-slices we choose for the extrapolation of the ground state form fac-tor, etc.It turns out that only the two plateaus presented in Sec.III B contribute to a visible deviation in the central value, which we add in the third error in Eq. 18.
Needless to say, there are also other source of systematic errors that are more difficult to quantify, say neglecting the disconnected contributions, quenching of the the strange and charm quarks, etc.The disconnected diagrams contributions are believed to be suppressed in the charmonium system [28,29] due to the Okubo-Zweig-Iizuka(OZI) rule.Furthermore, the non-physical masses of up and down quarks usually only result in a small effect which is indicated in the previous lattice calculations [36].Therefore, the major direction in future improvements points to the deduction of the statistical noise in χ c0 correlation functions.Only after the large statistical uncertainty is fully under control, should we worry about other remaining systematic effects.
Part of the large statistical error in our study can be traced back to the mixing of χ c0 and η c in the twistedmass formulation of lattice QCD.To entangle this mixing, we have utilized a GEVP procedure that projects out the operators best overlapped with η c and χ c0 as discussed in Sec.III A. Although this procedure works perfectly for the ground state η c , the efficiency for χ c0 is not quite satisfactory, rendering the two-point and threepoint correlation functions of χ c0 much noisier than that of η c , resulting in a much larger error for the decay rate of χ c0 .Possibilities to get around this difficulty could be simply increasing the statistics of the ensembles, using more interpolating operators as the basis operators or simply using a formulation that does not suffer from this mixing effect at all, e.g.utilizing the clover-improved Wilson fermion configurations.

IV. CONCLUSION
In this paper, we report a new lattice QCD computation of the scalar charmonium χ c0 to two-photon decay width.We have performed our study using three ensembles of N f = 2 twisted mass gauge field configurations at three different lattice spacings.This allows us to perform a more reliable continuum extrapolation therefore eliminating the substantial finite lattice spacing errors observed in previous lattice studies.We also adopt a new method that directly extracts the relevant on-shell form factor, by-passing the extrapolation in the photon virtualities.We obtain the decay width of χ c0 meson to be Γ γγ (χ c0 ) = 3.65(83) stat (21) lat.syst (66) syst keV.Albeit the large errors in this computation, the result is compatible with the existing experimental values within 1.3σ.Further possible improvements are also discussed.This calculation and possible future more systematic studies will await the new experimental results that will become available soon.

Figure 1 .
Figure 1.The left panels show the effective mass at different time slices together with the corresponding fitting ranges (grey bands) and the right panels are the ground and excited state mass values fitted from two-point correlation functions using Eq.(14).The black symbols denote the chosen m0, that correspond to the grey band to its left.The green symbols in (a) denote another choice for the m0 and m1.
TableIII.Mass value m0 and spectral weight Z0 for ground state and the first excited state mass m1 on each ensemble respectively.Ens.I(a) and Ens.I(b) are corresponding to the black and green symbols respectively.

Figure 2 .
Figure 2. The left column represents the plateaus of the form factor with different t h while the right column shows the extrapolation to the ground state contribution.The label (a) and (b), (c) and (d), (e) and (f) are for Ens.I(a), Ens.II, Ens.III respectively.

Figure 3 .
Figure 3. Continuum extrapolation for the ratio Γγγ(χc0)/mχ c0 .The three data points with both horizontal and vertical error-bars are results from three ensembles.The extrapolated results are shown by two side-by-side points near a 2 = 0: The one with a smaller error bar (the right one) represents the extrapolation result without considering lattice spacing errors.The other one (left one) is the result with lattice spacing errors taken into consideration.The data point (blue) below these two with a smaller error is the PDG-fit value for this ratio.At the upper right corner, we have also indicated the result of the width in physical units.
Ensemble a (fm) L 3 × T N conf aµ l mπ (MeV) t h