Two Photon Decays of $\eta_c$ from Lattice QCD

We present an exploratory lattice study for the two-photon decay of $\eta_c$ using $N_f=2$ twisted mass lattice QCD gauge configurations generated by the European Twisted Mass Collaboration. Two different lattice spacings of $a=0.067$fm and $a=0.085$fm are used in the study, both of which are of physical size of 2$fm$. The decay widths are found to be $1.025(5)$KeV for the coarser lattice and $1.062(5)$KeV for the finer lattice respectively where the errors are purely statistical. A naive extrapolation towards the continuum limit yields $\Gamma\simeq 1.122(14)$KeV which is smaller than the previous quenched result and most of the current experimental results. Possible reasons are discussed.


I. INTRODUCTION
Charmonium systems play a major role in the understanding of the foundation of quantum chromodynamics (QCD), 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, making it an ideal testing ground for our understanding of QCD from both sides.
Two-photon decay width of η c has been attracting considerable attention over the years from both theory and experiment sides. For example, it is related to the process gg → η c relevant for charmonia production at Large Hadron Collider (LHC) and the small-x gluon distribution function from the inclusive production of η c which describes the non-leptonic B mesons decays [1]. Furthermore, two-photon branching fraction for charmonium provides a probe for the strong coupling constant at the charmonium scale via the two-photon decay widths, which can be utilized as a sensitive test for the corrections for the non-relativistic approximation in the quark models or the effective field theories such as non-relativistic QCD (NRQCD).
On the experimental side, considerable progress has been made in recent years in the physics of charmonia via the investigations from Belle, BaBar, CLEO-c and BES [2][3][4][5]. Two methods can be utilized to measure the two-photon branching fraction for charmonium. One is reconstructing the charmonium in light hadrons with two-photon fusion at e + e − machines. The other one is to make pp pairs annihilated to charmonium with decay * Corresponding author. Email: liuchuan@pku.edu.cn and then to detect the real γγ pairs. Improvements of the measurement for two-photon branching fraction of charmonia will soon be reached in the future.
On the theoretical side, charmonium electromagnetic transitions have been investigated using various theoretical methods [6][7][8][9][10][11][12][13][14]. In principle, these processes involve both electromagnetic and strong interactions, the former being perturbative in nature while the latter being non-perturbative. Therefore, the study for charmonium transitions requires non-perturbative theoretical methods such as lattice QCD. Normal hadronic matrix element computations are standard in lattice QCD, however, processes involving initial or final photons are a bit more subtle. Since photons are not QCD eigenstates, one has to rely on perturbative methods to "replace" the photon states by the corresponding electromagnetic currents that they couple to. The details of this idea was illustrated in Ref. [15,16]. Using this technique, the first ab initio quenched lattice calculation of two photon decay of charmonia was reported in Ref. [17]. They found a reasonable agreement with the experimental world-average values for η c and χ c0 decay rates. However, an unquenched lattice study is still lacking. In this paper, we would like to fill this gap by exploring the two photon decay rates of η c meson in lattice QCD with N f = 2 flavors of light quarks in the sea. The gauge configurations utilized in this study are generated by the European Twisted Mass Collaboration (ETMC) [18][19][20][21][22][23][24][25][26][27][28][29], where the twisted mass fermion parameters are set at the maximal twist. This ensures the so-called automatic O(a) improvement for on-shell observables where a is the lattice spacing [30].
This paper is organized as follows. In Section II, we briefly review the calculation strategies for the matrix element for two-photon decay of η c . The matrix element is normally parameterized using a form factor, which in turn is directly related to the double photon decay rates.
The following section III are divided into three parts containing details of the simulation: In section III A, we introduce the twist mass fermion formulation and give the parameters of the lattices used in our simulation. In section III C, the continuum and lattice dispersion relations for η c are checked. In section III E, numerical results of the form factor are presented which are then converted to the decay width of η c meson. Our final number comes out to be smaller than the world-average experimental result and barely agrees with the previous quenched result. Possible reasons are discussed for this discrepancy. In Section IV, we discuss possible extensions of this cal-culation in the future and conclude.

II. STRATEGIES FOR THE COMPUTATION
In this section, we briefly recapitulate the methods for the calculation of two-photon decay rate of η c presented in Ref. [17]. The amplitude for two-photon decay of η c can be expressed in terms of a photon two-point function in Minkowski space by means of the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula, Here |Ω designates the QCD vacuum state, |η c (p) is the state with an η c meson of four-momentum p and |γ(q i , λ i ) for i = 1, 2 denotes a single photon state with corresponding polarization vector (q i , λ i ), with q i and λ i being the corresponding four-momentum and helic-ity, respectively. Then one utilizes the perturbative nature of the photon-quark coupling to approximately integrate out the photon fields and rewrites the corresponding path-integral as, The integration over the photon fields can be carried out by Wick contracting the fields into propagators. Neglect-ing the disconnected diagrams, one arrives at the following equation, In this equation, is the free photon propagator, which in momentum space will cancel out the inverse propagators outside the integral in Eq. (1) and Eq. (3) when the limit is taken. Effectively, each initial/final photon state in the problem is replaced by a corresponding electromagnetic current operator which couples to the photon and eventually one needs to compute a three-point function of the form Ω|T j ρ (z)j σ (w) |η c (p f ) . This quantity is nonperturbative in nature and should be computed using lattice QCD methods.
The current operators such as j ρ (x) appearing in Eq. (3) are electromagnetic current operators due to all flavors of quarks. However, we will only consider the charm quark in this preliminary study. Contributions due to other quark flavors, e.g. up, down or strange, only come in via disconnected diagrams which are neglected in this exploratory study. Another subtlety in the lattice computation is that, with c(x)/c(x) being the bare charm/anti-charm quark field on the lattice, composite operators such as the current j ρ (x) = Z V (g 2 0 )c(x)γ ρ c(x) needs an extra multiplicative renormalization factor Z V which we infer from Ref. [31]. To be specific, for the two set of lattices used in this study, the values of the renormalization factor Z V (g 2 0 ) are 0.6103(3) and 0.6451(3) for the lattice size 24 3 × 48 at β = 3.9 and 32 3 × 64 at β = 4.05, respectively. Annihilation diagrams of the charm quark itself are also neglected due to OZIsuppression. In fact, in our twisted mass lattice setup, we introduce two different charm quark fields with degenerate masses, so that this type of diagram is absent, see subsection III A.
The resulting expression (3) can then be analytically continued from Minkowski to Euclidean space. This continuation works as long as none of the q 2 i is too time-like. To be precise, the continuation is fine as long as the virtualities of the two photons where M V is the mass of the lightest vector meson in QCD [15,17]. For quenched lattice QCD, the lightest vector meson is J/ψ. However, for our unquenched study, it is safe to take M V = m ρ , i.e. the mass of the ρ meson. Using suitable interpolating operator (denoted by O ηc (x)) to create an η c meson from the vacuum and reversing the operator time-ordering for later convenience, we finally obtain, where O ηc (x) is an interpolating operator that will create an η c meson from the vacuum and ω 1 is the energy of the first photon. The kinematics in this equation is such that four-momentum conservation p f = q 1 + q 2 is valid. This equation serves as the starting point for our subsequent lattice computation. Basically, the current that couples to the first photon is placed at the source time-slice t i , the second current is at t while the final η c meson is at the sink time-slice t f and we are led to the computation of a three-point function of the form Of course, one has to compute the above three-point functions for each t i and perform an integration (summation) over t i . Apart from the above mentioned three-point functions, we also need information from η c two-point function. For example, in the above equation, Z ηc (p f ) is the spectral weight factor while E ηc (p f ) is the energy for η c with fourmomentum p f = (E ηc , p f ). These can be inferred from the corresponding two-point functions for η c . For this purpose, two-point correlation functions for the interpolating operators O ηc are computed in the simulation: where Z ηc (p f ) = Ω|O ηc |η c (p f ) is the corresponding overlap matrix element.
The three-point functions, denoted by G µν (t i , t), that need to be computed in our simulation are of the form, Keeping the sink of η c fixed at t f = T /2, we compute G µν (t i , t) across the temporal direction for all t i and t on our lattices. For a fixed t i , one has to use sequential source technique to obtain the t dependence of the three-point function. Then, the same calculation is repeated with a varying t i . Then, according to Eq. (5), the desired matrix element is obtained by using the results of G µν (t i , t) for different combinations of t i and t and integrate over t i with an exponential weight e −ω1|ti−t| . In practice, the integral is replaced by a summation over t i . To explore the validity of this replacement, we have checked the behavior of the integrand some of which are illustrated in Fig. 1. It is seen that these integrand as a function of t i indeed peak around the corresponding t values.
The matrix element η c |γ(q 1 , λ 1 )γ(q 2 , λ 2 ) can be parameterized using the form factor F (Q 2 1 , Q 2 2 ) as follows, where µ (q 1 , λ 1 ), ν (q 2 , λ 2 ) are the polarization vectors of the photons while q 1 and q 2 are the corresponding fourmomenta. The physical on-shell decay width Γ for η c to two photons is related to the form factor at Q 2 1 = Q 2 2 = 0, which will be referred to as the physical point in the following, via where α em (1/137) is the fine structure constant. Therefore, to extract the physical decay width, we simply compute the corresponding three-point functions in Eq. (5) and then extract the form factors F (Q 2 1 , Q 2 1 ) at various virtualities close to the physical point. Then, we can extract the information for F (0, 0) yielding the physical decay width. Although the physical decay width is only related to F (0, 0), the behavior of F (Q 2 1 , q 2 2 ) at non-zero virtualities are also of physical relevance when studying processes involving one or two virtual photons.

A. Simulation setup
In this study, we use twisted mass fermions at the maximal twist. The most important advantage of this setup is the so-called automatic O(a) improvement for the physical quantities. To be specific, we use N f = 2 (degenerate u and d quark) twisted mass gauge field configurations generated by the European Twisted Mass Collaboration (ETMC). The other quark flavors, namely strange and charm quarks, are quenched. These quenched flavors are introduced as valence quarks using the Osterwalder-Seiler (OS) type action [18,32]. Following the Refs. [18,22,23], in the valence sector we introduce three twisted doublets, (u, d), (s, s ) and (c, c ) with masses µ l , µ s and µ c , respectively. Within each doublet, the two valence quarks are regularized in the physical basis with Wilson parameters of opposite signs (r = −r = 1). The fermion action for the valence sector reads One can perform a chiral twist to transform the quark fields in physical basis to the so-called twisted basis as follows: where ω = π/2 implements the full twist. Two sets of gauge field ensembles are utilized in this work, each containing 200 gauge field configurations. We shall call them Ensemble I and II respectively. The explicit parameters are listed in Table I. The corresponding renormalization factor Z V (g 2 0 ) and the valence charm quark mass parameter µ c are taken from Ref. [31]. For the meson operators, in the physical basis, we use simple quark bi-linears such asqΓq and the corresponding form in twisted basis will be denoted asχ q Γ χ q which can be readily obtained from Eq. (11). For later convenience, these are tabulated in table II together with the possible J P C quantum numbers in the continuum and the names of the corresponding particle in the light and the charm sector. The current operators that appear in Eq. (7) are also listed. TABLE II. Local interpolating operators for vector and pseudo-scalar states and the current operators that appear in Eq. (7) in both physical and twisted basis,qΓq =χqΓ χq. The names of the corresponding particle and their J P C quantum numbers in the continuum are also listed. The index for i, µ and ν are 1, 2, 3.

B. Twisted boundary conditions
In order to increase the resolution in momentum space, particularly close to the physical point of Q 2 1 = Q 2 2 = 0, it is customary to implement the twisted boundary conditions (TBC) [31,[33][34][35] in recent lattice form factor computations, see e.g. [36]. We have also adopted the   twisted boundary conditions for the valence quark fields, also known as partially twisted boundary conditions. The quark field ψ θ (x, t), when it is transported by an amount of L along the spatial direction i(i = 1, 2, 3), will change by a phase factor e iθi , where θ = (θ 1 , θ 2 , θ 3 ) is the twisted angle for the quark field in spatial directions which can be tuned freely. In this calculations, we only twist one of the charm quark field in both vector currents, the other charm quark fields remain un-twisted. If we introduce the new quark fieldŝ it is easy to verify thatĉ (x, t) satisfy the conventional periodic boundary conditions along all spatial directions; satisfies the twisted boundary conditions (12). For Wilson-type fermions, this transformation is equivalent to the replacement of the gauge link; i.e, for µ = 0, 1, 2, 3 and θ µ = (0, θ). In other words, each spatial gauge link is modified by a U (1)-phase. Then the current vectors that appear in Eq. (7) are constructed using the hatted and the original charm quark field as, The allowed momenta on the lattice are thus modified to for i = 1, 2 where n i ∈ Z 3 is a three-dimensional integer. By choosing different values for θ, we could obtain more values of q 1 and q 2 than conventional periodic boundary conditions. In this paper, apart from the untwisted case of θ = (0, 0, 0), we have also computed the following cases: θ = (0, 0, π), (0, 0, π/2), (0, 0, π/4) and (0, 0, π/8).
These choices offer us many more data points in the vicinity of the physical kinematic region.

C. Meson spectrum and the dispersion relations
Before calculating the matrix element with two photon decay from η c , the mass for η c and ρ state and the energy dispersion relations for η c must be verified. This is particularly important for our study due to the following reasons. Firstly, we use the N f = 2 twisted mass configurations, the sea quarks contains u and d quark field. therefore virtual ρ state can enter the game. Thus, we should calculate the ρ mass so as to ensure the photon virtualities Q 2 1 , Q 2 2 > −m 2 ρ in this simulations. Secondly, we do need the information from η c correlation functions, the value of E ηc (p) and Z ηc (p) in order to extract the relevant matrix elements. Finally, we should also check the dispersion relation of η c which is quite heavy in lattice units (around 0.95) in our simulation and therefore some  kinematic factors (q ρ 1 and q σ 2 ) that enter Eq.(8) might need modifications accordingly.
Following Eq. (6), the energy E ηc (p f ) for η c state with three-momentum p f can be obtained from the corresponding two-point function via The two point function is symmetric about t = T /2. In real simulation we average the data from two halves about t = T /2 to improve statistics. We use the effective mass plateaus at zero three-momentum for the η c and ρ state to obtain the masses which are then listed in Table III. The mass of the η c comes out to be lighter than its physical value since these values are still finite lattice spacing values. When extrapolated towards the continuum limit, the mass will become compatible with the experimental value. The mass of the ρ here serves to restrict our kinematic regions where analytic continuation is justified. Similarly, we obtain the energies for η c at nonvanishing momenta via Eq. (17) which then can be utilized to verify the the following two dispersion relations: the conventional one in the continuum, and its lattice counterpart, For free particles, the constants Z cont and Z latt should be close to unity. In Fig. 2, we show this comparison for the two dispersion relations of the η c states in our simulation. In the left/right panel, the dispersion relations are illustrated using lattice/continuum dispersion relations, respectively. In both panels, points with errors are from simulations on 32 3 × 64 (open circles) or 24 3 × 48 (stars) lattices. Straight lines are the corresponding linear fits to the data. It is seen that, although both dispersion relations can be fitted nicely using linear fits, the slope for the naive continuum dispersion relation, i.e. Z cont is definitely different from unity, see e.g. right panel of Fig. 2, while its lattice counterpart Z latt is close. This suggests that, for the η c state, we should use the lattice dispersion relations instead of the naive continuum dispersion relation. This is not surprising since η c is quite heavy in lattice units. This modification of the dispersion relation does have consequences on our determination of the form factor. To illustrate this difference further, we plot the quantity E 2 (p)/(m 2 + p 2 ) as a function of p 2 in lattice units (i.e. a 2 p 2 , note that at these small values of p 2 , the difference between p 2 and the lattice versionp 2 is negligible) for two of our ensembles. This is shown in Fig. 3. It is seen that this quantity deviates from unity by as much as 10% even at rather small values of p 2 . This is actually caused by the difference between the rest mass and the kinetic mass of the η c meson.

D. Kinematics
In order to fully explore the form factor close to the physical point Q 2 1 = Q 2 2 = 0, we performed a parameter scan in the two virtualities. The following notations will be utilized. First of all, in the continuum, we will use q 1,2 to designate the four-momentum of the two photons. We will also use ω 1,2 to denote the temporal component of q 1,2 , i.e. ω 1,2 ≡ q 0 1,2 . When the photons are on-shell, we have ω 1,2 = |q| 1,2 with q 1,2 being the corresponding three-momentum. The so-called virtuality of the photons are defined as the corresponding four-momentum squared: Q 2 1,2 ≡ (−q 2 1,2 ). On the lattice, however, there are also lattice counterparts of the above notations, arising from the lattice dispersion relation (19). For that we simply add a hat on the corresponding variable. For example, we will usê ω 1 = 2 sinh(ω 1 /2) to denote the lattice version of ω 1 .
The computation has to cover the physical interesting kinematic region. For this purpose, we have to scan the corresponding parameter space. We basically follow the following strategy: We first fix the four-momentum of η c , p f = (E ηc , p f ), and place it on a given time-slice t f = T . Note that we just have to fix p f = n f (2π/L) and E ηc can be obtained from the dispersion relation (19). This effectively puts η c on-shell. Here we also have the freedom to pick a value for the twist angle θ. Then, we judiciously choose several values of virtuality Q 2 1 around the physical point Q 2 1 = 0. To be specific, we picked the range Q 2 1 ∈ [−0.5, +0.5]GeV 2 , which satisfies the con-straint Q 2 1 > −m 2 ρ . 1 Since p f = q 1 + q 2 , this means that, for a given p f , a choice of q 1 completely specifies q 2 and vice versa. We therefore take several choices of q 1 = n 1 (2π/L) by changing three-dimensional integer n 1 . At this stage, we can compute the energy of the first photon ω 1 , since ω 2 1 = q 2 1 − Q 2 1 . It turns out that we can also compute the virtuality of the second photon, Q 2 2 = |q 2 | 2 − ω 2 2 , since ω 2 = E ηc − ω 1 and q 2 is also known by the choice of q 1 . 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. Pick p f and θ. Obtain E ηc (p f ) from dispersion relation (19); 2. Judiciously choose several values of Q 2 1 in a suitable range, say Q 2 1 ∈ [−0.5, +0.5]GeV 2 ; 3. Pick values of n 1 such that q 1 = n 1 (2π/L). This fixes both ω 1 and Q 2 2 , using energy-momentum conservation; 4. Make sure all values of Q 2 1 , Q 2 2 > −m 2 ρ , otherwise the choice is simply ignored; 5. For each validated choice above, compute the threepoint functions (7), the two-point functions (6) and eventually obtain the hadronic matrix element using Eq. (5).

E. Form factors
In order to compute the desired hadronic matrix element η c (p f )|γ(q 1 , λ 1 )γ(q 2 , λ 2 ) in Eq. (5), we choose to place η c state at a fixed sink position t f = T /2. This sink position is then used as a sequential source for a backward charm propagator inversion. We compute this with all possible source positions t i and insertion point t. This method allows us to freely vary the value of ω 1 , Q 2 1 (as discussed in previous subsection) and to directly inspect the behavior of the integrand in Eq. (5).
Taking p f = 0 for the η c state as an example,, we show the behavior of the integrand in Fig. 1 for insertion positions t = 4, 8, 12, 16, 20 for ensemble I and t = 4, 8, 12, 16, 20, 24, 28 for ensemble II. It is seen that the integrand is peaked around t i = t, making the contributions close to this point the dominant part of the matrix element. For the lattice theory, the integration of t i in Eq. (5) is replaced by a summation.
When passing from the matrix element to the form factors, one should be careful about the form of the momenta to use. Recall that these momentum factors originate from derivatives in the continuum. On the lattice, they should be replaced by the corresponding finite differences, i.e. one should use the lattice version of the momentum: q 0 → 2 sinh(q 0 /2) and q i → 2 sin(q i /2). Since the spatial momenta that we are using are relatively small in lattice units, the effect of this replacement might be optional. However, for the 0-th component, since each of the photon is roughly half of the η c energy which is large in lattice units as we discussed in subsection III C, this replacement does make a difference.
According to Eq. (5), the matrix element and therefore also the form factor F (Q 2 1 , Q 2 2 ) should be independent of the insertion point t. We indeed observe this plateau behavior in our data which is illustrated in Fig. 4 for the case of Q 2 1 = 0 as an example. Other cases are similar. Fitting these plateaus then yields the corresponding values for the matrix element η c |γ(q 1 , λ 1 )γ(q 2 , λ 2 ) or equivalently the form factor F (Q 2 1 , Q 2 2 ). To describe the virtuality dependence of the form factor, we adopt a simple one-pole parametrization to fit our data.
where F (Q 2 1 , 0) and µ 2 (Q 2 1 ) are regarded as the fitting parameters at the given value of Q 2 1 . Since measurements at different values of Q 2 1 or Q 2 2 are all obtained on the same set of ensembles, we adopt the correlated fits, taking into account possible correlations among different Q 2 values. The covariance matrix among them are estimated using a bootstrap method.
As an example, taking Q 2 1 = −0.5GeV 2 , the fitting results are shown in Fig. 5. It is seen that this simple formula describes the data rather well even for quite large values of Q 2 2 . We therefore have taken all available values of Q 2 2 into the fitting process. Notice also that, by using the twisted boundary conditions together with different combinations of the lattice momenta, we are able to populate the physical region close to Q 2 1 = Q 2 2 = 0 rather effectively. We have tried both correlated and uncorrelated fits on our data. The central values for the fitted parameters are compatible, however, the error estimates are somewhat different. We adopt the correlated fits as our final results. Fits for other set of parameters are similar and the final results are summarized in Table IV for reference. Having obtained the results for F (Q 2 1 , 0), we can fit it again with another one-pole form, with F (0, 0) and ν 2 being the fitting parameters. This is illustrated in Fig. 6 for two of our ensembles. Again, correlated fits are adopted here. Apart from fitting the data in a two-step procedure as described above, we have also tried to fit the data in a one-step method. When we plug Eq. 21 into Eq. 20 and assuming that we are only interested in the value of the form factor close to the physical point, we may Taylor expand it assuming both Q 2 1 and Q 2 2 are small, Thus, we could fit the data in a region close to the origin with F (0, 0), a and b being the fitting parameters. This is illustrated in Fig. 7 for two of our ensembles. In each case, 35 data points of (Q 2 1 , Q 2 2 ) close to the origin are taken and the corresponding form factors F (Q 2 1 , Q 2 2 ) are obtained. Then using a linear fit in both Q 2 1 and Q 2 2 , c.f. Eq. (22), the form factors at the origin are obtained for both ensembles. Again, correlated fits are adopted here. The fitting results are summarized in Table VI. When computing the physical double photon decay width, according to Eq. (9), one has to plug in the mass of the η c meson. What we really compute on the lattice is the combination of correlation functions which is related to the matrix element η c |γ(q 1 , λ 1 )γ(q 2 , λ 2 ) via Eq. (5). When we parameterize this particular matrix element in terms of form factor in Eq. (8), the relation involves m ηc as well. Therefore, the decay width turns out to be proportional to m 3 ηc : Γ ∝ m 3 ηc | η c |γ(q 1 , λ 1 )γ(q 2 , λ 2 ) | 2 . Here it is then quite different if one substitutes in the value of m ηc obtained on the lattice, or the true physical value of m phys. ηc = 2.98GeV, the two differs by about 10% for the coarser lattice and about 5% for the finer lattice. Therefore, if one would substitute in the true physical mass, it will result in a 15% difference in the value of Γ for the finer lattice and about 30% for the coarser one.
The reason for the above mentioned difference is the following. We are taking the value of the valence charm quark mass parameter µ c from Ref. [31]. There, it is assumed that, when the continuum limit is taken, the value of m ηc will recover its physical value. However, being on a finite lattice, the computed value of m ηc comes out to be less than the corresponding physical value. The difference of the two is in fact an estimate of the finite lattice   spacing error. In fact, m ηc is not the only factor which affects the results. The renormalization factor Z V (g 2 0 ) that we quoted in Table I also depends on the lattice spacing. Therefore, we think it is more consistent to substitute in the values of m ηc computed on each ensembles. In the end, of course, one should try to take the continuum limit when the lattice computations are performed on a set of ensembles with different lattice spacings.
If using the two-step fitting procedure using Eq. (20) and using Eq. (21), with the values of m ηc obtained from each ensemble substituted in, we obtain for the decay width Γ = 1.019(3)KeV for the coarser and Γ = 1.043(3)KeV for the finer lattice ensembles. These results for the form factor F (0, 0) together with the corresponding results for the decay width are summarized in Table V.
As for the one-step fitting procedure using Eq. (22), we obtain Γ = 1.025(5)KeV for the coarser and Γ = , 0) and µ 2 (Q 2 1 ) using Eq. (20) for Ensemble I (left four columns) and Ensemble II (right four column). The total χ 2 value and the corresponding total degrees of freedom is also listed in the columns labelled by χ 2 /dof . Ensemble I Ensemble II  1.062(5)KeV for the finer lattice ensembles. The results for the form factor F (0, 0) and Γ(η c → γγ) are consistent with each other for Ensemble I using two different types of fitting procedure. However, for Ensemble II, a combined fitting using Eq. (22) gives a larger result for both F (0, 0) and Γ(η c → γγ). We think the value using the combined fit is more reliable since it gives a much less value of χ 2 /dof . In the combined fitting, the naive continuum extrapolated result for the decay width reads Γ = 1.122 (14)KeV. The fitted results for F (0, 0) together with the corresponding results for the decay width are summarized in Table VI. Let us now discuss the possible systematic errors. Although the mass of the pion in the two ensembles are relatively heavy, we do not expect the double photon decay width to be very sensitive to the pion mass. Also, since both of our ensembles have m π L ∼ 3.3, we do not expect very large finite volume errors as well. Since we have only two ensembles, it is not possible to make reliable extrapolation towards the continuum limit. However, if one would try a naive continuum limit extrapolation, assuming an O(a 2 ) error, we obtain Γ = 1.082(10)KeV which is also listed in Table V. There are of course other sources of systematic errors, e.g. the neglecting of the so-called disconnected contributions, the quenching of the strange quark, etc. Therefore, we decided not to quantify the systematic errors in this exploratory study. However, as we discussed above, the difference in the η c mass already indicates that there might be a finite lattice spacing error at the order of 15% for the finer and 30% for the coarser ensembles, respectively. Now let us briefly discuss the implications of our lattice results for Γ. First of all, our results are substantially smaller than the previously obtained quenched lattice result in Ref. [17], which we quote: Γ = 2.65(26) stat (80) scal. (53) quen. KeV. Second, our value is also much smaller than most of the experimental values. The current experimental value, according to PDG, is about 5.0 KeV with an error of 0.4 KeV [37]. However, one should note that the most recent determination of this quantity by Belle [38] with the result 5.8 ± 1.1 KeV is not a direct measurement of Γ γγ itself, but the product Γ γγ B(η c → η π + π − ) 50.5 eV. Therefore, the value of Γ γγ is usually extracted by inferring to earlier measurements in other channels, making the final results for Γ γγ differ quite a bit. For example, if we blindly use the PDG quoted value of the branching ratio, B(η c → η π + π − ) = 0.041 ± 0.017, we arrive at Γ γγ 1.25 KeV, which is comparable to our lattice result. However, if we would infer Γ γγ from the ratio of Γ γγ Γ(KKπ)/Γ tot = 0.407 ± 0.027KeV and Γ(KKπ)/Γ tot = (7.0 ± 1.2) × 10 −2 , we end up with Γ γγ = 5.8±1.1KeV as in Ref. [38]. Therefore, it is highly desirable to have a more precise and/or direct measurement of this quantity in future experiments.
The possible reasons for these apparent discrepancies can come from several sources to be discussed below. First, we have used different configurations from the quenched calculations. Our calculation takes into account the sea quark contributions from u and d quarks while in Ref. [17] these have been ignored. Although the quenching errors have been estimated in Ref. [17], it is well-known that this type of systematic is very difficult to quantify accurately. It is therefore quite possible that these effects have been under-estimated in Ref. [17].
Another possibility is that we have a rather large systematic errors which is not fully quantified in this ex-ploratory study. It is seen that our statistical errors seem to be small. However, as mentioned above, we do observe a large finite lattice spacing error of about 15-30% just from the mass of the η c . Since we have only two lattice spacings, the continuum limit extrapolation is also not well-controlled. In fact, if we blindly ascribe an error of about 15% for the numbers of the decay width for the two ensembles in Table V, it is possible that we could end up with a number that is close to the quenched result but with a rather large error coming from the continuum limit extrapolation. Of course, it is also possible that this disagreement is due to the combination of the above mentioned sources. In any case, a more systematic study with more lattice ensembles will definitely help to clarify these issues.

IV. CONCLUSIONS
In this exploratory study, we calculate the decay width for two-photon decay of η c using unquenched N f = 2 twisted mass fermion configurations. The computation is done with two lattice ensembles at two different lattice spacings. The mass spectrum and dispersion relations for the η c state are first examined. It is verified that lattice dispersion relations are better than the continuum ones. The implication of this is carried over to the computation of hadronic matrix element and the corresponding form factors.
By calculating various three-point functions, twophoton decays of η c matrix element are obtained at various of virtualities. It is particularly helpful to implement the so-called twisted boundary conditions which enable us to populate the physical region well. The matrix element is decomposed into kinematic factors and one form factor F (Q 2 1 , Q 2 2 ) which is obtained in a region close to the physical point. Then, we adopt a simple  , 0) is again fitted with a one-pole form: F (Q 2 1 , 0) = F (0, 0)/(1 + Q 2 1 /ν 2 ) are shown for the two ensembles (the first two lines). In the last column, we show the decay width obtained using Eq. (9). A naive continuum extrapolation are shown in the third line for reference.
Our result is significantly smaller than both the quenched result and the experimental values quoted by the PDG. However, taking into account of the possibly large systematic errors in the present lattice computations and the large uncertainties in the experimental result itself, it is still premature to say that there is a severe discrepancy here. Obviously, future more systematic lattice studies with various lattice spacings and more statistics are very much welcome here. It would also be helpful to estimate the disconnected contributions that has been neglected in this exploratory study. It will also be helpful to use other types of unquenched configurations, e.g. with 2 + 1 flavors or even 2 + 1 + 1 flavors in order to estimate the effects for the quenching of the other quark flavors. Last but not the least, more precise experimental results on double photon decays of charmonium are crucial in this area as well.