Ward Identity of the Vector Current and the Decay Rate of $\eta_c\rightarrow\gamma\gamma$ in Lattice QCD

Using a recently proposed method arXiv:1910.11597 (Yu Meng et al.), we study the two-photon decay rate of $\eta_c$ using two $N_f=2$ twisted mass gauge ensembles with lattice spacings $0.067$fm and $0.085$fm. The results obtained from these two ensembles can be extrapolated in a naive fashion to the continuum limit, yielding a result that is consistent with the experimental one within two standard deviations. To be specific, we obtain the results for two-photon decay of $\eta_c$ as $\mathcal{B}(\eta_c\rightarrow 2\gamma)= 1.29(3)(18)\times 10^{-4}$ where the first error is statistical and the second is our estimate for the systematic error caused by the finite lattice spacing. It turns out that Ward identity for the vector current is of vital importance within this new method. We find that the Ward identity is violated for local current with a finite lattice spacing, however it will be restored after the continuum limit is taken.


I. INTRODUCTION
The charmonium two-photon decay process η c → 2γ has long been an ideal testing ground for the understanding of non-perturbative nature of quantum chromodynamics (QCD) [2] due to the medium energy scale of charmonium systems in strong interactions [3]. On one hand, it offers an access to the strong coupling constant at the charmonium scale within the framework of perturbative QCD. On the other hand, it also provides a sensitive test for the application of effective field theories such as non-relativistic QCD (NRQCD) [4], which has become a useful tool in the treatment of quarkonium spectrum, decay and production. Therefore, the study of this particular process has been attracting considerable attentions over the years and discussed extensively in the literature from both experiments [5][6][7][8] and various theoretical methods, notably NRQCD and lattice QCD studies [9][10][11][12].
Combining the experimental results in recent years, the latest Particle Data Group (PDG) lists the branching fraction for the process as B(η c → 2γ) = (1.57 ± 0.12) × 10 −4 [13]. On the theoretical side, however, there have been no theoretical results so far that come even close to the experimental values. For example, within the framework of NRQCD factorization, the authors in Ref. [12] have computed the next-to-next-to-leading order QCD corrections to this process, yielding a value for B(η c → 2γ) that is about twice the value quoted by PDG. This indicates that the NRQCD might even break down for such processes due to non-perturbative effects.
It is then natural to turn to genuine nonperturbative methods such as lattice QCD (LQCD). With the proposal and realization of photon hadronic structure on the lattice in Refs. [14], such an idea has been widely used in photon structure functions [15], radiative transi-  [13] tion [16] and two-photon decays in charmonia [9][10][11] and neutral pion two-photon decay [17]. The first quenched LQCD calculation of η c → 2γ was presented in 2006 [9] and unquenched results have also showed up in recent years [10,11]. In Table. I, we have summarized the available lattice results so far. As a comparison, also listed in this table are the values from NRQCD and PDG. It is evident that none of these theoretical results can explain the PDG value satisfactorily.
In previous lattice calculations mentioned so far, the hadronic matrix elements relevant for charmonia double gamma decays are decomposed into form factors which are functions of the photon virtualities Q 2 i , i = 1, 2. Via an appropriate fitting of matrix element at different Q 2 i with a specific functional form, one obtains the complete off-shell form factors. Then, the physical decay width can be obtained by setting all virtualities to the on-shell values, namely Q 2 i = 0, yielding the final decay rate. However, the large deviations between the experiments and lattice results in Ref. [9][10][11] indicate that such methods might suffer from rather severe lattice artifacts and the decomposition itself might also be troublesome on the lattice with finite lattice spacings. This is understandable in a way because, for such processes, the photons in the final state are rather energetic (typically 1.5GeV in physical unit) in lattice units for commonly used lattice spacings.
Therefore, it is of great significance to explore new methods. In a recent work [1], we have proposed a new method to compute the three photon decay rate of J/ψ on the lattice directly with all polarizations of the initial and final states summed over. Such a method is originally put forward to avoid the complicated decomposition for the matrix element M (J/ψ → 3γ). In this paper, we plan to use it on the two-photon decay of η c , which is the simplest case that one could imagine. If we are only interested in the physical decay width, i.e. on-shell matrix elements, we can perform the summation over polarizations of the initial and final particles first, and then only the on-shell matrix elements needs to be computed on the lattice. In order to perform the summation over polarisations, Ward identities associated with the vector currents are crucial. It is well-known that, in the continuum Minkowski space, due to Ward identities, the summation over polarizations of the photons yields the Minkowski metric, e.g. λi ( λi µ (q i ) λi, * µ (q i ) ⇒ −g µµ . Generally speaking, Ward identity is broken for finite lattice spacing a. In other words, Ward identity breaking (WIB) corrections have to be considered when we sum over the photon polarizations on the lattice. However, as we will see, we can still make this substitution once the summation over all polarizations of initial and final particles are performed. This comes from the fact that the WIB effects for on-shell matrix element eventually vanish when the continuum limit is taken.
The rest of this paper is organized as follows. In Sec. II, we give a detailed derivation of the matrix element for the two-photon decay of η c . In Sec. III, we compare the new method that has been proposed in Ref. [1] with the conventional approaches and explain how the decay width can be obtained directly without the decomposition of the relevant form factor. In Sec. IV, details of simulations are given and the main results are presented. This section is divided into three parts: in Sec. IV A, the lattice dispersion relation for η c is checked; in Sec. IV B, the current renormalization constant is calculated; in Sec. IV C, numerical results of the matrix element squared and the corresponding WIB corrections are presented. These results are eventually converted into the two-photon decay width of η c . A naive continuum extrapolation is also performed and the final results are compared with the PDG value. It is found that our result is consistent with the PDG value within two standard deviations. Finally, we conclude in Sec. V.

II. APPROACH TO DECAY AMPLITUDE ON LATTICE
In this section, we recapitulate on the general method adopted in the calculation of the two-photon decay width of η c which have been utilized in previous lattice studies [9][10][11]. We start by expressing the decay matrix element of η c → 2γ in terms of the appropriate three-point function using Lehmann-Symanzik-Zimmermann reduc-tion formula in Minkoswki space and integrating out the photon fields perturbatively. It then follows that the relevant matrix element reads where the two functions H µν (x, y) and Q µν (x, y), which will be called the hadronic and the non-hadronic part respectively, will be defined shortly. For later convenience, we reverse the operator time ordering and the decay amplitude M on a finite lattice can be written as where V = L 3 , L is the space length and T is time length of the lattice. The factor V · T arises from the fourmomentum conservation δ-function on a finite lattice. In the following, we will fix the meson at time slice t f , first photon with four-momentum q 1 = (ω 1 , q 1 ) at time slice t i and the other at time slice t with four-momentum q 2 = (ω 2 , q 2 ).

A. The hadronic part Hµν
The hadronic part H µν (x, y) is given by To produce the meson η c with three-momentum p from the QCD vacuum state |0 , we introduce the interpolating field operatorÔ ηc (z, t f ) in coordinate space that carries the quantum number of η c , and the state |η c (p) may be obtained via, Taking the conjugate of this equation and substitute into Eq. (3) and inserting the completeness relation where |n, p is the eigenstate of QCD Hamiltonian. Here p indicates the meson momentum and n the corresponding energy level. For large eneugh t f , only the ground state n = 0 dominates. For simplicity, we denote E 0 (p) as E p and finally obtain the following expression for the hadronic part function H µν (x, y), with Z ηc (p) being the ground state amplitude η c (p)|Ô ηc (0)|0 .

B. The non-hadronic part Q µν
The non-hadronic, or to be more precise, the photonic part is given by where λi µ (q i ) denotes the photon polarization vector with arbitrary four-momentum q i and helicity λ i . It can be obtained by an appropriate Lorentz transformation from the standard basis 1 µ = (0, 1, 0, 0) and 2 µ = (0, 0, 1, 0). The free photon propagator D µµ (x, w) is given by which cancel out the inverse propagator outside the integral in Eq. (7) in momentum space. As explained in Ref. [14], the resulting expression of Q µν can be analytically continued from Minkowski to Euclidean space. This process introduces the photon virtualities Q 2 i = |q i | 2 −ω 2 i , which are not too time-like to produce any on-shell vector hadrons. More specifically, one needs V where M V is mass of the lightest vector meson. The expression in Eq. (7) would blow up with too timelike photon because of the possible presence of intermediate state with the energy M V and larger ones. Plug the expression of free photon propagator into Eq. (7), we have (9) where the standard Wick rotation t i → −it i , t → −it has been performed. Combining the Eq. (2), (6) and (9) together, the final result is obtained with form as M ∼ 1 T dt dt i (· · · ). It can be viewed as the summation average of various time slices. For the usual lattice simulation, the equivalent treatment is to replace the summation average of timeslice t by its corresponding platform, then statistical error would be greatly reduced without regarding to the too large time-slice t where the signal is noisy. Then the final result is written as The correlation function appearing in the above equation can be evaluated in lattice QCD in terms of quark propagators. For simplicity, we denote the matrix element in Eq. (10) as M = µ ν M µν . Each M µν can be computed on the lattice by searching a plateau behavior in t, as long as t f − t is large enough. For simplicity, we have used the local current j µ (x) =c(x)γ µ c(x) for the charm quark which can be renormalized by a multiplicative factor Z V . In real simulations, the integrals in Eq. (10) are also replaced by corresponding trapezoidal summation. We also mention here that it is impossible to put both photons with momenta q i on-shell on the lattice. Thus the matrix element M µν calculated on the lattice is always off-shell with some non-vanishing photon virtualities Q 2 i , i = 1, 2.

III. NEW APPROACH TO THE DECAY WIDTH ON THE LATTICE
In this section, we first discuss the relationship between amplitude M and decay width Γ with conventional method [9][10][11], and then introduce the new method that has been put forward in Ref. [1].
In conventional simulations, the hadronic matrix element M µν is parameterized using a form factor F (Q 2 1 , Q 2 2 ): The physical on-shell decay width Γ for η c decaying to two physical photons is related to the form factor at Q 2 where α em (1/137) is the fine structure constant in quantum electrodynamics (QED). We remind the reader that the decomposition (11) is obtained under the assumption of Lorentz invariance and Bose symmetry. However, when evaluated on the lattice, the matrix element M µν has only hypercubic symmetry. Therefore, strictly speaking, this decomposition can only be utilized when the relevant momenta, namely the components of q 1 and q 2 , are small in lattice units. This might become problematic since the typical momentum of each photon in the final state is roughly m ηc /2.
In this paper we proceed in another way as advocated in Ref. [1]. We introduce which will be called T -function in the following. In the above equation, we have considered the Ward identity of the currents. Then the summation over polarizations of the photons yields the Minkowski metrix, e.g.
In our actual simulations, we sum over all possible |M µν | 2 's. The physical decay width of η c → 2γ in η c in the center of mass frame can be expressed as where in the last line, the momenta in the T -function needs to be on-shell due to the energy-momentum conservation.
Due to the discreteness of the momenta on the lattice, however, it is impossible to exactly impose the onshell condition for all final particles, making the on-shell quantity T not directly accessible, but with these nonvanishing small virtualities. These matrix elements can be computed directly on the lattice, the norm of which we denote as T (Q 2 1 , Q 2 2 ). This differs from the on-shell Tfunction only because of the fact that some of the photons are still off-shell. We then try to estimate the on-shell quantity T by the following fitting formula, where everything is measured in lattice units. We expect such behavior since the final two photons are identical.

A. Ward identity breaking corrections
As we have pointed above, Lorentz invariance is broken on lattice, which leads to a breakdown of Ward identity. In principle, the replacement in Eq. (14) is not valid. This results in a correction in the T -function which will be called Ward identity breaking (WIB) correction. With WIB correction included, the summation over polarizations of the photons is modified as [19] λi λi µ (q i ) λi, * µ (q i ) ⇒ −g µµ + ∆ where we have introduced ∆ where is the WIB correction term. In principle, one expects that δT (∆) approaches to zero in the continuum limit. This could be verified numerically in the simulations.

IV. SIMULATIONS AND RESULTS
Our lattice simulations are performed using two N f = 2 flavour twisted mass gauge field ensembles generated by the Extended Twisted Mass Collaboration (ETMC) with lattice spacing a 0.067 fm and 0.085 fm, respectively. The corresponding physical pion masses are 300 MeV and 315 MeV. The most important advantage of these setups is so-called automatic O(a) improvement for the physical quantities with twisted mass quark action at maximal twist [20]. In Table. II, we list all ensembles used in this study with the relevant parameters. For more details about these ensembles, see Ref. [21,22]. We utilized the Osterwalder-Seiler setup for the valence sector of the charm quark [23]. Firstly, we determine the heavy quark mass by the physical η c mass with the corresponding meson operator O ηc (z) =c(z)γ 5 c(z) in physical basis, and verify the discrete dispersion relation in Eq. (20) by calculating the energies of η c at a series of three-momenta. The valence charm quark masses are tuned such that the mass of η c reproduces its true physical value in physical units. Such a verification is necessary because this particular discrete dispersion relation is to be utilized to obtain the photon energy ω i with given virtuality Q 2 i and three momentum q i . The discrete dispersion relation for the meson η c is, In Fig. 1, our results for the dispersion of η c are shown for two ensembles. It is found that the constant Z latt is almost 1. In what follows, we calculate the energies of the photons ω i with virtuality Q 2 i and momentum q i using this discrete dispersion relation (basically replacing m ηc by iQ i ). In this study, we have chosen 4 sets of suitable momenta and virtualities Q 2 i with the purpose of putting the T-function almost on-shell.

B. ZV and Zη c (p)
Secondly, we need to determinate current renormalization factor Z V which is needed to renormalize the photon current operator j µ (x) =cγ µ c(x). We extract the current renormalization factor Z V using a ratio of the three-point function over a two-point function [16] as given by where the factor 1/2 accounts for the equal contribution to the two-point function of the source at time slice 0 and the image of the source at time slice T . For decay particles in the rest frame, only µ = 0 is needed. For simplicity, this index µ will be omitted in the following. The plateau behavior of Z (0) V (t) across different time slice t then yields the value of the renormalization factor Z V . As an illustration, this is shown in Fig. 2 where the data points with errors are from our simulation and the horizontal bars indicate the intervals from which Z V are extracted. The final extracted Z V is given by Z V = 0.6237(2), 0.6523(1) for L = 24 and L = 32, respectively.
The value of Z ηc (p), which is also needed in the calculation of the decay width as indicated in Eq. (10), can be obtained from from the two-point function given by Eq. (22) below, where we have denoted Z ηc = Z ηc (0), E ηc = E ηc (0) for simplicity. Note that the meson η c is fixed at the source t f = T /2 and wall-source is adopted in our simulations.
C. The decay width of ηc → 2γ The conventional sequential method has been adopted to calculate the three-point function in Eq. (10). The meson is fixed at the source position with time slice t f = T /2. We put the sequential source on one current with time slice t i , and the contraction is performed on the other current at sink t. After the integration (summation) of time slice t i , the matrix element M µν , being a function of time slice t, can be obtained on the lattice.
The input parameters include photon momenta q i = 2π L n i , virtualities Q 2 i and energies ω i . It is found that, for each set of photon momenta, only some of the M µν values can be reached. We follow the strategies outlined in Ref. [9][10][11]. In fact, Q 2 2 is uniquely dependent on Q 2 1 due to energy-momentum conservation. In actual simulations, we choose a convention in which the photons share the same virtualities Q 2 1 = Q 2 2 = Q 2 m , which is then determined by After summations of all 16 matrix elements M µν , the T-function T (Q 2 m , t) can be obtained immediately and the results are shown in Fig. 4 in case of n q = (0, 2, 2). Finally, the on-shell T-function can be fitted by Eq. (16) where two variables are utilized. In the case of one Q 2 m variable, the on-shell fitting formula is reduced to with T (0) and a, b being the fitting parameters. One can also take the WIB correction terms in Eq. (19) into ac-count and estimate the WIB effect on the two-photon decay width of η c . We point out that the WIB effects result from the non-conservation of the local current arising from the finite lattice spacing. In fact, the difference between the two treatment can be viewed as an estimate of the lattice artifacts.
In the following, we denote T W as the T-function with WIB corrections included while T being the one without the corrections. Similar symbols will be utilized for the decay widths Γ W and Γ. Both the on-shell results of T W and T under two spacings are shown in the left panel of Fig. 5 and the corresponding values are also summarized in Table. III The errors here only include the statistical errors estimated using bootstrap method. Also included are the errors from the current renormalization factor Z V , ground state amplitude Z ηc and on-shell fitting process as suggested in Eq. (24). As seen from the results in Eq. (25) , there exist differences between Γ and Γ W in both ensembles at finite lattice spacings. These differences can be viewed as an estimate of the finite lattice spacing errors on these two ensembles, We therefore regard the difference between Γ and Γ W the systematic error and take the average value as the chosen decay width which is denoted asΓ. The values are then obtained as, Here the first error is statistical and the second represent the systematic error as described above.
We now turn to the naive continuum extrapolations. For the study of charmonium with N f = 2 configurations, one can assume an O(a 2 ) errors for the lattice results for the decay widths obtained above. This allows us to connect the the two results for Γ, Γ W andΓ at two lattice spacings and obtain the corresponding results at a = 0. We call this the naive continuum extrapolation. Admittedly, this is not a well-controlled continuum extrapolation. For that purpose, one needs results for at least three or more different lattice spacings. Taking the average of Γ and Γ W , namelyΓ as our final result, the decay width for the η c → 2γ is found to be, Γ(η c → 2γ) = 4.11(9)(58)KeV (27) Here the first error is statistical and the second is the estimate of the systematic error due to lattice spacing. These quantities are illustrated in Fig. 5. In the left panel, the on-shell fitting for T and T W under two different spacings are performed. It is found that the difference caused by the WIB effect is dependent on the lattice spacing. For the finer lattice spacing, the difference between the two is smaller. This is expected since the breaking of the Ward identity is due to a finite lattice spacing. In the right right panel of Fig. 5, we illustrate the naive continuum limit extrapolations for the decay width Γ(η c → 2γ) and Γ W (η c → 2γ), respectively. In this limit, Γ and Γ W are well consistent with each other as expected. Also shown in the right panel is the average of the Γ(η c → 2γ) and Γ W (η c → 2γ), namelyΓ for the two ensembles and together with its naive continuum extrapolation. It is seen that, when the finite lattice spacing errors are taking into account, our naive continuum extrapolated result in Eq. (27) is consistent with the experimental one within two standard deviation.
We emphasize that, all the continuum extrapolations, whether for Γ and Γ W , or Γ, are just naive due to the limited number of lattice spacings. Still, our final result for the decay width of η c → 2γ is encouraging. These are the first lattice results that is consistent with the experiments within 2σ level. There are also other source of systematic errors: finite volume effects, pion mass still far from physical and the contribution of disconnected diagrams. However, we think that finite lattice spacing errors are by far the most relevant error at present. Future lattice studies should aim to improve on this by utilizing more lattice ensembles which will substantially reduce this error.
The branching fraction, if the uncertainty of η c total width ignored, is given by B(η c → 2γ) = 1.29(3)(18) × 10 −4 , where the first error is statistical and the second is our estimates for the systematics due to finite spacing. The result is reliably consistent with the experiment result B exp (η c → 2γ) = 1.57(12) × 10 −4 [13]. Compared to the previous much smaller ones obtained with traditional method of form factor parameterizations, our results seem to indicate that the continuum form of parameterizations utilized might fail drastically for the calculation of the hadronic decays on the lattice.

V. CONCLUSIONS
In this paper, we calculate the two-photon decay rate of η c with all polarizations of the final photon states summed over, which is first proposed in Ref. [1]. Using two N f = 2 twisted mass gauge ensembles with two different spacings, we have obtained the branching fraction B(η c → 2γ) = 1.29(3)(18) × 10 −4 where the first error is statistical and the second is our estimated systematic error due to finite lattice spacing. This result is consistent with the experimental result quoted by PDG within two standard deviation. An improved result would be expected in the future if more lattice spacings are utilized.
Further more, we have demonstrated that Ward identity for the current, which is essential for our method to work, is in fact violated with a finite lattice spacing a for a local current. After a detailed comparison between the decay width of η c → 2γ with Ward identity breaking (WIB) effects included and excluded, we have shown that such a discrepancy vanishes in the continuum limit. This indicates that we can always replace the summation of photon polarizations safely by the Minkowski metric when we calculate the decay width of multi-photon final states as long as the continuum limit is taken in the end.