A Lattice Study of the Two-photon Decay Widths for Scalar and Pseudo-scalar Charmonium

In this exploratory study, two photon decay widths of pseudo-scalar ($\eta_c$) and scalar ($\chi_{c0}$) charmonium are computed using two ensembles of $N_f=2$ twisted mass lattice QCD gauge configurations. The simulation is performed two lattice ensembles with lattice spacings $a=0.067$ fm with size $32^3\times{64}$ and $a=0.085$ fm with size $24^3\times{48}$, respectively. The results for the decay widths for the two charmonia are obtained which are in the right ballpark however smaller than the experimental ones. Possible reasons for these discrepancies are discussed.


I. INTRODUCTION
Charmonium physics plays an important role in the foundation of quantum chromodynamics (QCD) which is believed to be the fundamental theory for the strong interaction. Due to its intermediate energy scale and the special features of QCD, both perturbative and non-perturbative physics show up within charmonium physics.
Until now, the best way to study nonperturbative QCD is lattice QCD, a quantum field theory defined on discrete Euclidean space-time. Within this formalism, physical quantities are encoded in various Euclidean correlation functions which in turn can be measured using Monte Carlo simulations [1,2].
Recently, two photon decay branching fraction for charmonium has been attracting considerable attentions. Theoretically, this quantity is considered to provide a probe for the strong coupling constant at the charmonium scale. It has also been proposed as a sensitive test of the corrections to the non-relativistic approximation via quark models or the effective field theories such as non-relativistic QCD (NRQCD).
In addition, considerable progress has been made in recent years in the physics of charmonia via the investigations from various experimental collaborations, such as, Belle, BaBar, CLEO-c, and BES [3][4][5][6]. Generically, there are two ways to measure two-photon branching fraction for charmonium: one is by reconstructing the charmonium in light hadrons with two-photon fusion at e + e − machines, the other one os from the pp annihilation to charmonium.
In this paper, we calculate the two-photon decay width of the pseudoscalar and scalar charmonium, i.e. Γ(η c → γγ) and Γ(χ c0 → γγ), in lattice QCD using two ensembles of N f = 2 dynamical twisted mass fermion configurations that were generated by the European Twisted Mass Collaboration(ETMC). This ensures the so-called automatic O(a) improvement for on-shell observables when the twisted mass fermions are at the maximal twist [7]. Lattice computation for the process η c → γγ has been done using the same set of configurations in Ref. [8]. In this work, we consider η c and χ c0 simultaneously since they mix with each other due to lattice artifact of O(a 2 ). This paper is organized as follows. In Sect. II, the calculation strategy for the matrix element for two-photon decay of charmonium is reviewed which is then related to the double-photon decay rates. We also outline the strategy for the mass spectrum and form factor of charmonium. In order to improve the signals of the corresponding correlation functions, we apply the variation method to construct the relevant interpolate operators. In Sect. III, the simulation details is divided into several parts: In Sect. III A, we give the parameters relevant with the configurations used in this work. In Sect. III B, mass spectrum of η c and χ c0 are obtained. In Sect. III C, we presents the results of renormalization factor of electromagnetic current operators. In Sect. III D, numerical results of the form factors are presented which are then converted to the two photon decay width of η c and χ c0 . Final results for the decay widths are presented and compared with previous lattice computations and the experiments. The discussion and the conclusion can be found in Sect. IV.

II. STRATEGIES FOR THE COMPUTATION
In this section, we briefly review the methods for the calculation of two-photon decay rate of a charmonium which was presented in Ref. [9]. According to Lehmann-Symanzik-Zimmermann (LSZ) reduction formula, we can express the amplitude for two-photon decay of charmonium in the following form: Here |Ω designates the QCD vacuum state; M represents either η c or χ c0 meson state, depending on which one is calculated, and |M (p f ) is the corresponding meson state with four-momentum p f , and |γ(q i , ι i ) (i = 1, 2) is a single photon state which has the polarization vector (q i , ι i ) with four-momentum q i and the helicity ι i . Then, treating QED part perturbatively, we can replace the photon field operators by their corresponding current operators in QCD. We finally arrive at the following equation at the lowest order of QED, with D µν (y, z) being the free photon propagator. Basically, each initial/final photon state in the problem is replaced by the corresponding electromagnetic current operator that couples to the photon. Finally, one needs to compute a three-point function of the form Ω|T j ρ (z)j σ (w) |M (p f ) which is non-perturbative in nature and can be computed using lattice QCD methods.
Under certain conditions, Eq. (2) can be analytically continued from Minkowski space to Euclidean space, yielding where ϕ M ( x, t f ) is the field operator for the meson M, Z M (p f ) is the spectral weight factor of two-point function, ω 1 is the energy for the photon at time-slice t i , and E M (p f ) is the energy for corresponding meson with the momentum p f . Then, the desired amplitude γ(q 1 , ι 1 )γ(q 2 , ι 2 )|M (p f ) can be obtained once the energies E M (p f ) and the corresponding overlap matrix element Z M (p f ) are known. These could be obtained from appropriate two-point functions. In this study, we use the variational method to find the optional interpola-tion operators to create/annihinite the η c and χ c0 meson state [10]. Generally, an operator O † i having definite J P C can produce all QCD eigenstates with the right quantum numbers In order to create the desired hadrons from vacuum effectively, one could employ a basis of interpolators O i that share the same quantum numbers and construct two-point correlation function matrix as follows Here, the operators O i are color-singlet constructions built from the basic quark and gluon fields of QCD. Then, we can express the correlation functions in the following form: and the optimized interpolators are: Therefore, one should obtain the best estimate for the weights v n i according to the solution of the following generalized eigenvalue problem Here, C(t) is the N × N matrix whose elements are the correlation functions C ij (t) constructed from the basis of N operators, v n is a generalized eigenvector. The generalized eigenvalues, or principal correlators, λ n (t) behave like e −En(t−t0) at large times, and can be used to determine the spectrum of the states. In practice, we solve Eq. (7) independently on each time-slice t, so that for each state n, we obtain a time series of generalized eigenvectors v n (t). We use v n i chosen on a single time-slice to construct the optimized operators in Eq. (6).
Apart from the two-point functions of η c and χ c0 , we also need the three-point functions G µν (t i , t) given by: We simulate G µν (t i , t) on lattices across the temporal direction while the sink of the meson is fixed. Then, we repeat this with a varying t i to integrate over the three point function with an exponential weight e −ω1ti , and then to extract the matrix element in Eq. (3). In particular, we use the optimized interpolators Ω n in Eq. (6) to create the η c , (n = 1) or χ c0 , (n = 2) state for the field interpolating operator ϕ M in previous formulas. For the two photon decay of η c meson, the matrix element γ(q 1 , ι 1 )γ(q 2 , ι 2 )|M (p f ) in Eq. (3) can be parameterized using form factor F (Q 2 1 , Q 2 2 ) as, where 1 , 2 are polarization vectors, Q 2 1 , Q 2 2 are virtualities and q 1 , q 2 are the four-momenta for the two photons.
The corresponding decay width can be expressed in terms of F (0, 0), with α em 1/137 being the fine structure constant. Similarly, for χ c0 we have another form factor G( with the decay width given by

A. Simulation setup
In this work, we utilize two ensembles with N f = 2 (degenerate u and d quarks) twisted mass configurations. These configurations are generated by the ETMC at the maximal twist to implement the so-called automatic O(a) improvement [7]. The explicit parameters for these ensembles are presented in Ref. [11] and the two ensembles that we utilized are tabulated in Table I. For the valence sector, we adopt the Osterwalder-Seiler setup which amounts to introducing two extra twisted doublets for each non-degenerate quark flavors, namely, (u, d) and (c, c ) with twisted masses aµ l and aµ c , respectively [12][13][14][15][16]. The explicit values of aµ l on Ens. B 1 is 0.004, and 0.003 for the Ens. C 1 , respectively. In this simulation, we use the physical mass of η c to set the value of the aµ c , and the explicit values are 0.2542 and 0.2018 for Ens. B 1 and Ens. C 1 respectively. In each doublet, the Wilson parameters have opposite signs (r = −r = 1). Performing an axial (or chiral) transformation, quark fields in the physical basis transform into the twisted basis [12]; i.e., where ω is the twist angle, and ω = π/2 represents maximal twist. Then, the left of the above equations correspond to quark fields in the physical basis, and the right correspond to quark fields in the twisted basis.
Before writing out the explicit form of meson operators, one should exploit the symmetry properties of twisted mass LQCD. We will follow the discussion in reference [17] below. Isospin I and parity P are broken by O(a 2 ) effects in twisted mass LQCD. While, a specific combination (i.e. light flavor exchange combined with parity) is still a symmetry of twisted mass LQCD. We first write down the interpolating-field operators in the twisted basis, and build the interpolating operators with the same Wilson parameters [15]. For the purpose of η c and χ c0 , we use two basis operators . According to Eq. (13), the two basic operators in twisted basis are given by O 1 =χ c χ c and O 2 =χ c γ 5 χ c which appear to have opposite parity. However, since twisted mass lattice QCD breaks parity, they in fact mix with each other. Taking into account of this mixing is crucial. One needs to go through the solution of the generalized eigenvalue problem in Eq. (7) to obtain the optimized operator that will create the η c and χ c0 meson from the vacuum. Without performing this generalized eigenvalue separation, it is found that the correct signal of χ c0 cannot be observed in the two-point functions. The signal of η c can of course be observed even without considering this mixing effect since it is the lightest state under consideration.

B. Mass spectra for ηc and χc0
The eigenvalue λ n in Eq. (7) corresponding with the corresponding meson state, i.e. n = 1 represents η c meson, and n = 2 represents χ c0 meson . Since, we use the anti-periodic boundary condition, In practical, we use eigenvalue λ 1 with p f = (0, 0, 0) to fit the spectral weight Z ηc , and the explicit value of this factor is 0.4416(8) and 0.2675(3) on Ens. B 1 and Ens. C 1 , respectively. While, we use eigenvalue λ 2 with p f = (0, 0, 0) to fit the spectral weight Z χc0 , and the explicit value of this factor is 0.6699(71) and 0.2983(33) on Ens. B 1 and Ens. C 1 , respectively. It is easily seen that the mass can also be extracted from, with m 1 being the mass of η c meson and m 2 being that of χ c0 , respectively. The effective mass plateaus of these mesons for Ens. B 1 and Ens. C 1 are illustrated in Fig. 1. From these mass plateaus, the masses of the mesons are determined and the statistical errors are obtained using jackknife method. The numerical results for the masses are summarized in Tab. II. Note that the mass values for η c are utilized to fix the valence charm quark mass parameter aµ c . Therefore, only the mass of χ c0 are predictions from this lattice computation. The current operators in Eq. (8) are electromagnetic current operators. In principle, they contain all flavors of quarks weighted by the corresponding charges. Light quark flavors will only enter the question via disconnected diagrams which are neglected in this study. Only considering the charm quark, we only need to consider the currentcγ ρ c(x). A subtlety in the lattice computation is that, with c(x)/c(x) being the bare charm/anticharm quark field on the lattice, composite operators such as the current j ρ (x) = Z Vc (x)γ ρ c(x) needs an extra multiplicative renormalization factor Z V which can be extracted by the ratio of two-point function with respect to the three-point function for η c [19]: where µ is the Dirac index and we take it to be zero. Γ (2) ηc and Γ ηc are two point correlation function and three point correlation function relevant with η c . The explicit forms are: with O † ηc and O ηc creating and annihilating a state with the quantum number of η c meson, respectively. Actually, we use the simple local operator,i.e.cγ 5 c. According to Eq. (16), we can obtain the multiplicative renormalization factor Z V , and show it in Fig. 2. The values of the renormalization factor Z V are 0.6296 (18)   To compute the relevant matrix element in Eq. (3), we place the meson at a fixed sink position t f which is chosen to be 24 for Ens. B 1 and 32 for Ens. C 1 , respectively. These sinks are then used as a sequential source for a backward propagator inversion. This allows us to investigate all possible source positions t i . We can then freely vary the values of ω 1 , Q 2 1 and inspect the integrand as a function of t i in Eq. (3) for a given insertion position t. As an example, in Fig. 4 and Fig. 5, we show the in- The computation has to cover the physical interesting kinematic regions. For this purpose, we have to scan the corresponding parameter space of the two virtualities Q 2 1 and Q 2 2 . Basically, we follow the following strategy: we first fix the four-momentum of η c and χ c0 , p f = (E, p f ) and place it on a given time-slice t f = T /2. In this simulation, we only compute the case of p f = (0, 0, 0), and E is simply the mass of η c or χ c0 . For a given p f , a choice of q 1 completely specifies q 2 due to p f = q 1 + q 2 . Therefore, we take several choices of q 1 = n 1 (2π/L) by changing three-dimensional integer n 1 . Then the energy of the first photon is also obtained using either the continuum or the lattice dispersion relations: whereQ 2 1 = 4 sinh 2 (Q 1 /2) is the lattice version for the virtuality. It turns out that we can also compute the virtuality of the second photon, since both ω 2 and q 2 are constrained by the energy-momentum conservation. One has to make sure that, the values of Q 2 2 thus computed do satisfy the constraint Q 2 2 > −m 2 ρ otherwise it is omitted. This procedure is summarized as follows: 1. Judiciously choose several values of Q 2 1 in a suitable range. We picked 7 values of Q 2 1 ; 2. Pick different values of n 1 such that q 1 = n 1 (2π/L). As described above, this fixes both ω 1 and Q 2 2 using energy-momentum conservation. This is done using either the continuum or the lattice dispersion relations. To be specific, for each Q 2 1 , we picked 4 different q 1 ; 3. Make sure all values of Q 2 1 , Q 2 2 > −m 2 ρ , otherwise the choice is simply ignored; 4. For each validated choice above, compute the threepoint functions (8) and obtain the hadronic matrix element using Eq. (3).
In such a way, we have obtained altogether 28 points on the (Q 2 1 , Q 2 2 ) plane around the origin. As an example, the distribution of these virtualities for the two mesons are shown in Fig. 3 for the case of lattice dispersion relations. One could also do the same thing using the conventional continuum dispersion relations. The difference of these two treatments finally will provide us with an estimate for the finite lattice spacing error of the calculation.
In our real lattice QCD computation, the integration of t i in Eq. (3) are replaced by discrete summation over t i using trapezoid rule. The resulting values exhibit a plateau behavior with respect to t which is then utilized to extract the corresponding form factor. In such a way, we have obtained numerical results for F (Q 2 1 , Q 2 2 ) and G(Q 2 1 , Q 2 2 ) at 28 different points in the plane of the two virtualities. As an example, the form factor plateaus for η c are illustrated in Fig. 6 for the case of Q 2 1 = 0 GeV 2 . The corresponding case for χ c0 are shown in Fig. 7.

E. Fitting of the form factor and the physical decay widths
To obtain the physical decay width, we only need the values of the form factors at the physical photon point, namely Q 2 1 = Q 2 2 = 0. As we have seen in Fig. 3 the distribution of 28 data points in the (Q 2 1 , Q 2 2 ) plane, we could implement cuts in the plane. For a given value of Q 2 cut > 0, we select the points that satisfy the following inequality: Obviously, taking a large enough Q 2 cut will allow all of the points be considered while selecting a smaller value for Q 2 cut will take into account only those points whose distances are closer than Q 2 cut . On the other hand, for   0.08 t/a=4 t/a=8 t/a=12 t/a=16 t/a=20 t/a=24 t/a=28 a given value of Q 2 cut , we could utilize different fitting formula to obtain the corresponding values of the form factor. Since it is the value at the origin that is directly related to the decay width, it is natural to use a polynomial type of fit in both Q 2 1 and Q 2 2 . What is more, due to boson symmetry, this function needs to be symmetric with respect to Q 2 1 and Q 2 2 . Therefore, by varying the cut value Q 2 cut and various orders of polynomials in the virtualities, we could investigate the values of the form factors at the physical point.
To be specific, we adopt a polynomial ansatz up to Q 2 1 and Q 2 2 to the third power to fit the data of the form  factor. For η c meson, we use: and a similar form for the χ c0 form factor G(Q 2 1 , Q 2 2 ). Note that a 0 ≡ F (0, 0) is the form factor at the physical photon point which is directly related to the decay width of the meson. Polynomial forms with less terms, i.e. with only up to first or second powers in Q 2 1 and Q 2 2 . have also been attempted. Note that this implies that we are fitting the data points with 2, 4 and 6 parameters, respectively since terms at the same orders of Q 2 1 and Q 2 2 should be included or excluded on the same footing. In all cases, correlated fits are performed. Depending on the number of points taken into account which is controlled by Q 2 cut , and the number of terms kept in the fitting polynomial, we finally arrive at the values for the form factors at the origin, namely F (0, 0) = a 0 for both ensembles. Similar procedures have been implemented as well for χ c0 , resulting in the values for G(0, 0).
The fitting procedures described above can be implemented using either the continuum or lattice version of the dispersion relations as indicated in Eq. (19) or Eq. (20). The procedure can be carried out for either pseudo-scalar or scalar meson on either of the two ensembles utilized in this calculation. Therefore, we carry out the fitting procedure in 8 different cases. The difference between the corresponding results obtained from different dispersion relations will then inform us about lattice artifacts of the calculation.
As an illustration, in Fig. 8, the fitting results for η c and χ c0 on Ens.C1 using lattice dispersion relations are shown. Here the horizontal axis denotes the cut values Q 2 cut while the vertical axis indicates the values for F (0, 0) or G(0, 0) together with the errors (data points with error-bars). For each fixed value of Q 2 cut , we have performed three fits with 2, 4 and 6 parameters. These three points are slightly shifted horizontally so that they are recognizable. The corresponding values of χ 2 /dof for each fit are also shown as points without error-bars. By inspecting the plot, we get a feeling about the consistency and quality of these different fits and the differences among the values of F (0, 0) can also offer us an estimate of the systematics for the fitting procedure.
Having obtained these 8 plots, we proceed as follows: • For each of these plots, say the example above, we pick the case with lowest value of χ 2 /dof as the final result for F (0, 0) together its statistical error in this particular case.
• We further attribute a systematic error arising from the fitting procedure by taking the largest difference in the central values of F (0, 0) with comparable χ 2 /dof. This then yields the result for F (0, 0) with a certain type of dispersion relations on a particular ensemble.
• By comparing the difference in F (0, 0) between the two different dispersion relations, we further assign a systematic error, which is taken to be the difference between the two values, arising from the lattice spacing.
• Similar procedure can be applied to χ c0 as well.
In such a way, we have obtained the results of F (0, 0) and G(0, 0) on two ensembles which are explicitly listed below: In these expressions, the first error is statistical, the second is the error from the fitting procedure described above and the third one is the finite lattice error estimated from using two different dispersion relations. It is clearly seen that, in all cases, the results are dominated by systematic errors, especially the finite lattice spacing errors. In fact, two results on two different ensembles are consistent within this estimate of finite lattice spacing errors. We therefore decide not to make any continuum extrapolations. Obviously, computations with more values of lattice spacings will be crucial to nail down these large lattice spacing errors.
To compare with previous lattice computations, we notice that the result for Γ(η c → γγ) is slightly larger than previous result presented in Ref. [8] using the same set of configurations. This difference might come from the fact that the mixing of the η c and χ c0 in the twisted mass setup was not fully disentangled in the previous calculation in Ref. [8]. It is found that, if we were not using properly chosen operator that mix with both the η c -like and χ c0 -like interpolating operators, we would not be able to observe the correct χ c0 signal as we have discussed at the end of subsec III A.
Finally, let us convert the results in the form factors into corresponding ones in the decay widths. We simply add all the errors in the form factors in quadrature and neglect the errors in the mass of the mesons. This then leads to the following results for the decay widths: These are to be compared with the following values given by PDG: Γ(η c → γγ) P DG = 5.02(51)KeV, Γ(χ c0 → γγ) P DG = 2.20(22)KeV, These numbers are all in the same ballpark as the experimental ones though still somewhat smaller. However, since no controlled continuum extrapolations have been performed yet, it is still premature to draw any conclusions for the discrepancy. Our large estimated finite lattice errors offer some hint that this might be the major source of errors. In the future, more studies at different lattice spacing s are needed in order to control the lattice artifacts in a systematic fashion.

IV. CONCLUSIONS
In this exploratory study, we calculate the two-photon decay width for η c and χ c0 using unquenched N f = 2 twisted mass fermions. The computation is done with two lattice ensembles (coarser and finer) at two different lattice spacings. The mass spectrum for the η c and χ c0 meson state are obtained by solving a generalized eigenvalue problem which disentangles parity mixing between the two mesons.
Our results for the decay width Γ(η c → γγ) and Γ(χ c0 → γγ) are summarized in Eq. (27) for two ensembles utilized in this computation. With only two ensembles we only estimate the finite lattice spacing errors for each ensemble and no continuum extrapolations are performed. Albeit without the continuum extrapolations, our results are in the right ballpark as the PDG values shown in Eq. (28).
In the future, lattice calculations with more values of lattice spacings are definitely needed so as to control the finite lattice spacing errors which appear to be a dominant source. Meanwhile, the disconnected contributions is a good supplement to this study. It's also helpful to research on unquenched configurations with other ensembles or other methods. We also expect more precise experiments on double-photon decays of charmonium.