Lattice study of two-photon decay widths for scalar and pseudo-scalar charmonium

This exploratory study computes two-photon decay widths of pseudo-scalar ( ) and scalar ( ) charmonium using two ensembles of twisted mass lattice QCD gauge configurations. The simulation is performed using two lattice ensembles with lattice spacings fm with size and fm with size . The decay widths for the two charmonia are obtained within the expected ballpark, but are however smaller than the experimental ones. Possible reasons for these discrepancies are discussed.


Introduction
Charmonium physics plays a key role in the foundation of quantum chromodynamics (QCD), which is believed to be a fundamental theory for strong interactions. Owing to its intermediate energy scale and special features of QCD, both perturbative and non-perturbative physics appear within charmonium physics [1]. To date, the optimal way to study non-perturbative QCD is the lattice QCD, which is a quantum field theory defined on the discrete Euclidean space-time. Within this formalism, physical quantities are encoded in various Euclidean correlation functions, which in turn can be measured by performing Monte Carlo simulations [2,3].
Recently, the two photon decay branching fractions for charmonium have been attracting considerable atten-tion. 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 sensitivity test for making corrections to the non-relativistic approximation via quark models or effective field theories, such as the non-relativistic QCD (NRQCD).
e + e − pp Furthermore, considerable progress has been made in recent years in the field of physics of charmonia via the investigations from various experimental collaborations, such as Belle, BaBar, CLEO-c, and BES [4][5][6][7]. Generically, there are two approaches to measure the two-photon branching fraction for charmonium: one by reconstructing the charmonium in light hadrons using two-photon fusion at machines, and the other from annihilation to charmonium.
In this study, we calculate the two-photon decay width of the pseudoscalar and scalar charmonium, i.e., and , in lattice QCD using two ensembles of dynamical twisted mass fermion configurations that were generated by the European Twisted Mass Collaboration (ETMC). This ensures the so-called automatic improvement for on-shell observables when the twisted mass fermions are at the maximal twist [8]. Lattice computation for the process was perfomred using the same set of configurations in Ref. [9]. In this work, we consider and simultaneously, since they mix with each other due to the lattice artifact of . This paper is organized as follows. In Sec. 2, 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. To improve the signals of the corresponding correlation functions, we apply the variation method to construct the relevant interpolate operators. In Sec. 3, the simulation details are divided into several parts: In Sec. 3.1, we provide the parameters relevant with the configurations used in this study. In Sec. 3.2, the mass spectra of and are obtained. In Sec. 3.3, we present the results of the renormalization factor of the electromagnetic current operators. In Sec. 3.4, numerical results of the form factors are presented, which are then converted to the two photon decay width of and . Final results for the decay widths are presented and compared with previous lattice computations and experiments. The discussion and conclusion are presented in Sec. 4.

Strategies for computation
In this section, we briefly review the methods for the calculation of the two-photon decay rate of charmonium, which was presented in Ref. [10]. According to the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula, we can express the amplitude for the two-photon decay of charmonium in the following form: Here, denotes the QCD vacuum state; M represents either the or meson state, depending on which one is calculated, and is the corresponding meson state with four-momentum , whereas ( ) is a single photon state that has the polarization vector with four-momentum and helicity . Then, treating the 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 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, a three-point function of the form needs to be com-puted, which is non-perturbative in nature and can be found using lattice QCD methods. Under certain conditions, Eq. (2) can be analytically continued from Minkowski space to Euclidean space, yielding where is the field operator for the meson M; is the spectral weight factor of the two-point function; is the energy for the photon at time-slice ; and is the energy for the corresponding meson with momentum . Then, the desired amplitude is obtained once the energies Chinese Physics C Vol. 44, No. 8 (2020) 083108 083108-2 and the corresponding overlap matrix element are known. These can be obtained from appropriate two-point functions. In this study, we use the variational method to find the optional interpolation operators to create/annihilate the and meson states [11].
Generally, an operator with definite can generate all the QCD eigenstates with the right quantum numbers To effectively generate the desired hadrons from vacuum, one could employ a basis of interpolators that share the same quantum numbers and construct a two-point correlation function matrix as follows Here, the operators 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: v n i Therefore, the best estimate for the weights must be obtained according to the solution of the following generalized eigenvalue problem Here, is the matrix whose elements are the correlation functions constructed from the basis of N operators, and is a generalized eigenvector. The generalized eigenvalues, or principal correlators, , behave like 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, such that for each state n, we obtain a time series of generalized eigenvectors . We use chosen on a single time-slice to construct the optimized operators in Eq. (6).
Apart from the two-point functions of and , we also need the three-point functions given by: We simulate on lattices across the temporal direction, while the sink of the meson is fixed. Then, we repeat this with a varying to integrate across the three point function using an exponential weight , and then to extract the matrix element in Eq. (3). In particular, we use the optimized interpolators in Eq. (6) to generate the or state for the field interpolating operator in previous formulas.
For the two photon decay of the meson, the matrix element in Eq. (3) can be parameterized using form factor as follows, where are the polarization vectors, are teh virtualities, and are the four-momenta for the two photons. The corresponding decay width can be expressed in terms of as follows, with being the fine structure constant. Similarly, for , we have another form factor , with the decay width given by 3 Simulation results In this study, we utilize two ensembles with (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 improvement [8]. The explicit parameters for these ensembles are presented in Ref. [12], and the two ensembles that we utilized are presented in Table 1. For the valence sector, we adopt the Osterwalder-Seiler setup, which introduces two extra twisted doublets for each nondegenerate quark flavor, namely, and with twisted masses and , respectively [13][14][15][16][17]. The explicit value of on Ens. is , whereas it is for Ens. . In this simulation, we use the physical mass of to set the value of , yielding explicit values of and for Ens. and Ens. , respectively. In each doublet, the Wilson parameters have opposite signs ( ). Performing an axial (or chiral) transformation, quark fields in the physical basis transform into the twisted basis [13]; i.e., Chinese Physics C Vol. 44, No. 8 (2020) 083108

Simulation setup
where is the twist angle, and represents the maximal twist. Then, the left-side of the abovementioned equations correspond to quark fields in the physical basis, and the right-side correspond to quark fields in the twisted basis.
Before writing the explicit form of the meson operators, one should exploit the symmetry properties of the twisted mass LQCD. We follow the discussion in reference [18] below. Isospin I and parity are broken by effects in the twisted mass LQCD. Meanwhile, a specific combination (i.e., light flavor exchange combined with parity) remains a symmetry of the 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 [16]. For the purpose of and , we use two basis operators , . According to Eq. (13), the two basic operators in the twisted basis are given by and , which appear to have opposite parity. However, because twisted mass lattice QCD breaks parity, they in fact mix with each other. Taking into account this mixing is crucial. The solution of the generalized eigenvalue problem in Eq. (7) must be found to obtain the optimized operator that will create the and meson from the vacuum. Without performing this generalized eigenvalue separation, the correct signal of cannot be observed in the two-point functions. The signal of can naturally be observed even without considering this mixing effect, as it is the lightest state under consideration.

Mass spectra for and
The eigenvalue in Eq. (7) corresponding with the corresponding meson state, i.e., represents meson, and represents the meson. Therefore, we use the anti-periodic boundary condition, In practice, we use eigenvalue with to fit the spectral weight , and the explicit value of this factor is and on Ens. and Ens. , respectively. Meanwhile, we use eigenvalue with to fit the spectral weight , and the explicit value of this factor is and on Ens.
and Ens. , respectively. Thus, the mass can also be extracted from with denoting the mass of meson and that of . The effective mass plateaus of these mesons for Ens. and Ens.
are illustrated in Fig. 1. From these mass plateaus, the masses of the mesons are determined, and the statistical errors are obtained using the jackknife method. The numerical results for the masses are summarized in Table 2. Notably, the mass values for are utilized to fix the valence charm quark mass parameter . Therefore, only the masses of are predicted from this lattice computation.
In principle, glueball states with the same quantum numbers are also present in a similar energy range [20]. However, in this lattice calculation, we have only utilized the quark bilinear operators for charmonium states and did not observed the sign of the glueballs.
The current operators in Eq. (8) are electromagnetic current operators. In principle, they contain all the 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. When considering the charm quark, we only need to consider the current . A subtlety in the lattice computation is that, with being the bare charm/anticharm quark field on the lattice, composite operators such as the current need an extra multiplicative renormalization factor , which can be extracted by the ratio of two-point function with respect to the three-point function for [21]: where is the Dirac index, which we assume to be zero. and are the two point correlation function and three point correlation function relevant to , respectively. The explicit forms are:    with a particular value for and .
The computation must cover the physical kinematic regions of interest. For this purpose, we must scan the corresponding parameter space of the two virtualities and . Basically, we follow the following strategy: we first fix the four-momentum of and , and place it on a given time-slice . In this simulation, we only compute the case of , and E is simply the mass of the or meson. Then, we judiciously choose several values of virtuality around the physical point . To be specific, we pick the range GeV on Ens. and GeV on Ens. , which satisfies the constraint [9]. For a given , a choice of completely specifies due to . Therefore, we take several choices of by changing three-dimensional integer . Then the energy of the first photon is also obtained    Chinese Physics C Vol. 44, No. 8 (2020) 083108 083108-6 using either the continuum or the lattice dispersion relations as follows: where is the lattice version for the virtuality. Apparently, we can also compute the virtuality of the second photon, since both and are constrained by the energy-momentum conservation. One has to make sure that the values of thus computed do satisfy the constraint , otherwise it is omitted. This procedure is summarized as follows: 1. Judiciously choose several values of in a suitable range. We picked seven values of ; Pick different values of , such that . As described above, this fixes both and using the energy-momentum conservation. This is done using either the continuum or the lattice dispersion relations. To be specific, for each , we picked four different ;  In this approach, we have obtained a total of 28 points on the plane around the origin. As an example, the distribution of these virtualities for the two mesons are shown in Fig. 5 for the case of lattice dispersion relations. One could also peform the same procedure using the conventional continuum dispersion relations. The difference of these two approaches will ultimately provide us with an estimate of the finite lattice spacing error of the calculation.
In our real lattice QCD computation, the integration (3) is replaced by discrete summation over using the trapezoid rule. The resulting values exhibit a plateau behavior with respect to t, which is then utilized to extract the corresponding form factor. In this approach, we have obtained numerical results for and at 28 different points in the plane of the two virtualities. As an example, the form factor plateaus for are illustrated in Fig. 6 for the case of . The corresponding case for is shown in Fig. 7.

Fitting of the form factor and physical decay widths
To obtain the physical decay width, we only need the values of the form factors at the physical photon point, namely . Fig. 5 shows the distribution of the 28 data points in the plane, and we could implement cuts in the plane. For a given value of , we select the points that satisfy the following inequality: Evidently, taking a large enough will allow all of the points be considered, while selecting a smaller value for will account for only those points, whose distances are closer than . In contrast, for a given value of , we could utilize a different fitting formula to obtain the corresponding values of the form factor. Because 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 and . Furthermore, due to boson symmetry, this function must be symmetric with respect to and . Therefore, by varying the cut value 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 and to the third power to fit the data of the form factor. For the meson, we use: Chinese Physics C Vol. 44, No. 8 (2020) 083108 083108-7 and we apply a similar form for the form factor . Notably, 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 and , have also been attempted. This implies that we are fitting the data points with , , and parameters, respectively, as terms at the same orders of and must be included or excluded on the same footing. Correlated fits are performed in all the cases. Depending on the num- ber of points considered, which is controlled by , and the number of terms kept in the fitting polynomial, we finally arrive at the values for the form factors at the origin, namely for both the ensembles. Similar procedures have been implemented as well for , resulting in the values for . 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 performed for either the pseudo-scalar or scalar meson on either of the two ensembles utilized in this calculation. Therefore, we perform the fitting procedure in eight different cases. The difference between the corresponding results obtained from different dispersion relations will then inform us   As an illustration, in Fig. 8, the fitting results for and on Ens. C1 using lattice dispersion relations are shown. Here the horizontal axis denotes the cut values , while the vertical axis represents the values for or together with the errors (data points with error-bars). For each fixed value of , we have performed three fits with , , and parameters. These three points are shifted horizontally by a small amount to make them recognizable. The corresponding values of for each fit are also shown as points without errorbars. By inspecting the plot, we obtain an estimate regarding the consistency and quality of these different fits, and the differences among the values of can also offer us an estimate of the systematics for the fitting procedure.
Having obtained these eight plots, we proceed as follows: • For each of these plots, taking the example above, we pick the case with lowest value of as the final result for 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 with comparable . This then yields the result for with a certain type of dispersion relations on a particular ensemble.
• By comparing the difference in 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.
• A similar procedure can be applied to as well.
In this way, we obtain the results of and on two ensembles, as 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 is the finite lattice error estimated using the two different dispersion relations. Evidently, in all the cases, the results are dominated by systematic errors, especially the finite lattice spacing errors. In fact, the results on the two different ensembles are consistent within this estimate of finite lattice spacing errors. We therefore decide not to make any continuum extrapolations. The computations with more values of lattice spacings will clearly be crucial to determine these large lattice spacing errors.
To compare with the previous lattice computations, we note that the result for is slightly larger than the previous result presented in Ref. [9] when the same set of configurations were used. This difference is attributed to the fact that the mixing of and in the twisted mass setup was not fully disentangled in the previous calculation in Ref. [9]. In the case where a properly chosen operator that mixes with both the -like and like interpolating operators is not used, it would not be possible to observe the correct signal as discussed at the end of subsection 3.1.
Finally, we convert the results in the form factors into their corresponding ones in the decay widths. We simply add all the errors in the form factors in the quadrature and neglect the ones in the mass of the mesons. This then leads to the following results for the decay These numbers are all in the same ballpark as the experimental ones, although they remain somewhat smaller. However, because 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 should be conducted at different lattice spacings to control the lattice artifacts in a systematic manner. Another source of systematic errors could come from the mixing with the nearby glueball states. Thus far, no lattice calculations have considered such an effect. Considerable efforts are needed in the future to take the aforementioned factors into account.

Conclusions
In this exploratory study, we calculate the two-photon decay width for and using unquenched twisted mass fermions. The computation is performed with two lattice ensembles (coarser and finer) at two different lattice spacings. The mass spectra for the and meson state are obtained by solving a generalized eigenvalue problem, which disentangles parity mixing between the two mesons.
Our results for the decay width and are summarized in Eq. (27) for the two ensembles utilized in this computation. With these 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 certainly required to control the finite lattice spacing errors, which appear to be a dominant source. Meanwhile, the disconnected contributions will be a good supplement to this study. Possible mixing with the gluonic excitations must also be addressed. It will also be helpful to research unquenched configurations with other ensembles or other methods. We also expect to obtain more precise experiments to be performed on double-photon decays of charmonium.