A Lattice Study of the Two-photon Decay Widths for Scalar and Pseudo-scalar Charmonium
Ying Chen, Ming Gong, Ning Li, Chuan Liu, Yu-Bin Liu, Zhaofeng Liu, Jian-Ping Ma, Yu Meng, Chao Xiong, Ke-Long Zhang
AA Lattice Study of the Two-photon Decay Widths for Scalar and Pseudo-scalarCharmonium
Ying Chen, Ming Gong, Ning Li, ∗ Chuan Liu, Yu-Bin Liu, ZhaofengLiu, Jian-Ping Ma, Yu Meng, Chao Xiong, and Ke-Long Zhang (CLQCD Collaboration) Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, ChinaSchool of Physics, University of Chinese Academy of Sciences, Beijing 100049, China School of Science, Xi’an Technological University, Xi’an 710032, China School of Physics and Center for High Energy Physics, Peking University, Beijing 100871, ChinaCollaborative Innovation Center of Quantum Matter, Beijing 100871, China School of Physics, Nankai University, Tianjin 300071, China Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China School of Physics, Peking University, Beijing 100871, China
In this exploratory study, two photon decay widths of pseudo-scalar ( η c ) and scalar ( χ c ) charmo-nium 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 size32 ×
64 and a = 0 .
085 fm with size 24 ×
48, respectively. The results for the decay widths for thetwo charmonia are obtained which are in the right ballpark however smaller than the experimentalones. Possible reasons for these discrepancies are discussed.PACS numbers: 12.38.Gc, 11.15.HaKeywords: charmonium, decay width, lattice QCD.
I. INTRODUCTION
Charmonium physics plays an important role in thefoundation of quantum chromodynamics (QCD) whichis believed to be the fundamental theory for the stronginteraction. Due to its intermediate energy scale andthe special features of QCD, both perturbative andnon-perturbative physics show up within charmoniumphysics [1]. Until now, the best way to study non-perturbative QCD is lattice QCD, a quantum field the-ory defined on discrete Euclidean space-time. Withinthis formalism, physical quantities are encoded in vari-ous Euclidean correlation functions which in turn can bemeasured using Monte Carlo simulations [2, 3].Recently, two photon decay branching fraction forcharmonium has been attracting considerable attentions.Theoretically, this quantity is considered to provide aprobe for the strong coupling constant at the charmo-nium scale. It has also been proposed as a sensitive testof the corrections to the non-relativistic approximationvia quark models or the effective field theories such asnon-relativistic QCD (NRQCD).In addition, considerable progress has been made in re-cent years in the physics of charmonia via the investiga-tions from various experimental collaborations, such as,Belle, BaBar, CLEO-c, and BES [4–7]. Generically, thereare two ways to measure two-photon branching fractionfor charmonium: one is by reconstructing the charmo-nium in light hadrons with two-photon fusion at e + e − machines, the other one os from the p ¯ p annihilation to ∗ Corresponding author. Email: [email protected] charmonium.In this paper, we calculate the two-photon decay widthof the pseudoscalar and scalar charmonium, i.e. Γ( η c → γγ ) and Γ( χ c → γγ ), in lattice QCD using two ensem-bles of N f = 2 dynamical twisted mass fermion configura-tions that were generated by the European Twisted MassCollaboration(ETMC). This ensures the so-called auto-matic O ( a ) improvement for on-shell observables whenthe twisted mass fermions are at the maximal twist [8].Lattice computation for the process η c → γγ has beendone using the same set of configurations in Ref. [9]. Inthis work, we consider η c and χ c simultaneously sincethey mix with each other due to lattice artifact of O ( a ).This paper is organized as follows. In Sect. II, the cal-culation strategy for the matrix element for two-photondecay of charmonium is reviewed which is then relatedto the double-photon decay rates. We also outline thestrategy for the mass spectrum and form factor of char-monium. In order to improve the signals of the corre-sponding correlation functions, we apply the variationmethod to construct the relevant interpolate operators.In Sect. III, the simulation details is divided into severalparts: In Sect. III A, we give the parameters relevantwith the configurations used in this work. In Sect. III B,mass spectrum of η c and χ c are obtained. In Sect. III C,we presents the results of renormalization factor of elec-tromagnetic current operators. In Sect. III D, numericalresults of the form factors are presented which are thenconverted to the two photon decay width of η c and χ c .Final results for the decay widths are presented and com-pared with previous lattice computations and the exper-iments. The discussion and the conclusion can be foundin Sect. IV. a r X i v : . [ h e p - l a t ] M a y II. STRATEGIES FOR THE COMPUTATION
In this section, we briefly review the methods for thecalculation of two-photon decay rate of a charmonium which was presented in Ref. [10]. According to Lehmann-Symanzik-Zimmermann (LSZ) reduction formula, we canexpress the amplitude for two-photon decay of charmo-nium in the following form: (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) = − lim q (cid:48) → q q (cid:48) → q (cid:15) ∗ µ ( q , ι ) (cid:15) ∗ ν ( q , ι ) q (cid:48) q (cid:48) (cid:90) d xd y e iq (cid:48) .y + iq (cid:48) .x (cid:104) Ω | T (cid:8) A µ ( y ) A ν ( x ) (cid:9) | M ( p f ) (cid:105) . (1)Here | Ω (cid:105) designates the QCD vacuum state; M representseither η c or χ c meson state, depending on which one iscalculated, and | M ( p f ) (cid:105) is the corresponding meson statewith four-momentum p f , and | γ ( q i , ι i ) (cid:105) ( i = 1 ,
2) is a sin-gle photon state which has the polarization vector (cid:15) ( q i , ι i ) with four-momentum q i and the helicity ι i . Then, treat-ing QED part perturbatively, we can replace the photonfield operators by their corresponding current operatorsin QCD. We finally arrive at the following equation atthe lowest order of QED, (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) = ( − e ) lim q (cid:48) → q q (cid:48) → q (cid:15) ∗ µ ( q , ι ) (cid:15) ∗ ν ( q , ι ) q (cid:48) q (cid:48) (cid:90) d xd yd w d ze iq (cid:48) .y + iq (cid:48) .x D µρ ( y, z ) D νσ ( x, w ) ×(cid:104) Ω | T (cid:8) j ρ ( z ) j σ ( w ) (cid:9) | M ( p f ) (cid:105) (2)with D µν ( y, z ) being the free photon propagator. Ba-sically, each initial/final photon state in the problemis replaced by the corresponding electromagnetic cur-rent operator that couples to the photon. Finally, oneneeds to compute a three-point function of the form (cid:104) Ω | T (cid:8) j ρ ( z ) j σ ( w ) (cid:9) | M ( p f ) (cid:105) which is non-perturbative in nature and can be computed using lattice QCD meth-ods.Under certain conditions, Eq. (2) can be analyticallycontinued from Minkowski space to Euclidean space,yielding (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) = lim t f − t →∞ e (cid:15) µ ( q , ι ) (cid:15) ν ( q , ι ) Z M ( p f )2 E M ( p f ) e − E M ( p )( t f − t ) (cid:90) dt i e − ω ( t i − t ) (cid:104) Ω | T (cid:110) (cid:90) d (cid:126)x e − i (cid:126)p f .(cid:126)x ϕ M ( (cid:126)x, t f ) (cid:90) d (cid:126)y e i (cid:126)q .(cid:126)y j ν ( (cid:126)y, t ) j µ ( (cid:126) , t i ) (cid:111) | Ω (cid:105) , (3)where ϕ M ( (cid:126)x, t f ) is the field operator for the mesonM, Z M ( p f ) is the spectral weight factor of two-pointfunction, ω is the energy for the photon at time-slice t i , and E M ( p f ) is the energy for corresponding mesonwith the momentum p f . Then, the desired amplitude (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) can be obtained once the ener-gies E M ( p f ) and the corresponding overlap matrix ele-ment Z M ( p f ) are known. These could be obtained fromappropriate two-point functions. In this study, we usethe variational method to find the optional interpola- tion operators to create/annihinite the η c and χ c mesonstate [11].Generally, an operator O † i having definite J P C can pro-duce all QCD eigenstates with the right quantum num-bers O † i | Ω (cid:105) = (cid:88) n | n (cid:105)(cid:104) n |O † i | Ω (cid:105) . (4)In order to create the desired hadrons from vacuum ef-fectively, one could employ a basis of interpolators O i that share the same quantum numbers and constructtwo-point correlation function matrix as follows C ij = (cid:104) Ω |O i ( t ) O † j (0) | Ω (cid:105) . (5)Here, the operators O i are color-singlet constructionsbuilt from the basic quark and gluon fields of QCD. Then,we can express the correlation functions in the followingform: C ij = (cid:88) n E n (cid:104) Ω |O i (0) | n (cid:105)(cid:104) n |O † j (0) | Ω (cid:105) e − E n t . and the optimized interpolators are:Ω † n = (cid:88) i v ni O † i . (6)Therefore, one should obtain the best estimate for theweights v ni according to the solution of the following gen-eralized eigenvalue problem C ( t ) v n = λ n ( t ) C ( t ) v n . (7)Here, C ( t ) is the N × N matrix whose elements are thecorrelation functions C ij ( t ) constructed from the basis of N operators, v n is a generalized eigenvector. The gener-alized eigenvalues, or principal correlators, λ n ( t ) behavelike e − E n ( t − t ) at large times, and can be used to deter-mine the spectrum of the states. In practice, we solveEq. (7) independently on each time-slice t , so that foreach state n , we obtain a time series of generalized eigen-vectors v n ( t ). We use v ni chosen on a single time-slice toconstruct the optimized operators in Eq. (6).Apart from the two-point functions of η c and χ c , wealso need the three-point functions G µν ( t i , t ) given by: G µν ( t i , t ) = (cid:104) | T (cid:110) (cid:90) d (cid:126)x e − i (cid:126)p f .(cid:126)x ϕ M ( (cid:126)x, t f ) × (cid:90) d (cid:126)y e i (cid:126)q .(cid:126)y j ν ( (cid:126)y, t ) j µ ( (cid:126) , t i ) (cid:111) | (cid:105) . (8)We simulate G µν ( t i , t ) on lattices across the temporal di-rection while the sink of the meson is fixed. Then, werepeat this with a varying t i to integrate over the threepoint function with an exponential weight e − ω t i , andthen to extract the matrix element in Eq. (3). In partic-ular, we use the optimized interpolators Ω n in Eq. (6) tocreate the η c , ( n = 1) or χ c , ( n = 2) state for the fieldinterpolating operator ϕ M in previous formulas.For the two photon decay of η c meson, the matrix ele-ment (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) in Eq. (3) can be param-eterized using form factor F ( Q , Q ) as, (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) = 2( e ) m − η c F ( Q , Q ) (cid:15) µνρσ (cid:15) µ ( q , ι ) (cid:15) ν ( q , ι ) q ρ q σ (9)where (cid:15) , (cid:15) are polarization vectors, Q , Q are virtuali-ties and q , q are the four-momenta for the two photons. The corresponding decay width can be expressed in termsof F (0 , η c → γγ ) = πα em m η c | F (0 , | (10)with α em (cid:39) /
137 being the fine structure constant. Sim-ilarly, for χ c we have another form factor G ( Q , Q ), (cid:104) γ ( q , ι ) γ ( q , ι ) | M ( p f ) (cid:105) = 2( e ) m − χ c G ( Q , Q )( (cid:15) · (cid:15) q · q − (cid:15) · q (cid:15) · q )(11)with the decay width given byΓ( χ c → γγ ) = πα em m χ c | G (0 , | . (12) III. SIMULATION RESULTSA. 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 themaximal twist to implement the so-called automatic O ( a )improvement [8]. The explicit parameters for these en-sembles are presented in Ref. [12] and the two ensem-bles that we utilized are tabulated in Table I. For the TABLE I: Configuration parameters. β a [fm]
V /a aµ sea m π [MeV] m π L N cfg
Ens. B . ×
48 0.004 315 3.3 200Ens. C .
05 0.067 32 ×
64 0.003 300 3.3 199 valence sector, we adopt the Osterwalder-Seiler setupwhich amounts to introducing two extra twisted doubletsfor each non-degenerate quark flavors, namely, ( u, d ) and( c, c (cid:48) ) with twisted masses aµ l and aµ c , respectively [13–17]. The explicit values of aµ l on Ens. B is 0 . .
003 for the Ens. C , respectively. In this simulation, weuse the physical mass of η c to set the value of the aµ c ,and the explicit values are 0 . . B and Ens. C respectively. In each doublet, the Wilsonparameters have opposite signs ( r = − r (cid:48) = 1). Perform-ing an axial (or chiral) transformation, quark fields in thephysical basis transform into the twisted basis [13]; i.e., (cid:32) ud (cid:33) = exp( iωγ τ / (cid:32) χ u χ d (cid:33)(cid:32) cc (cid:48) (cid:33) = exp( iωγ τ / (cid:32) χ c χ c (cid:48) (cid:33) , (13)where ω is the twist angle, and ω = π/ I and parity P are broken by O ( a ) effects in twisted mass LQCD. While, a specificcombination (i.e. light flavor exchange combined withparity) is still a symmetry of twisted mass LQCD. Wefirst write down the interpolating-field operators in thetwisted basis, and build the interpolating operators withthe same Wilson parameters [16]. For the purpose of η c and χ c , we use two basis operators O ( x ) = ¯ c ( x ) γ c ( x ), O = ¯ c ( x ) c ( x ). According to Eq. (13), the two basicoperators in twisted basis are given by O = ¯ χ c χ c and O = ¯ χ c γ χ c which appear to have opposite parity. How-ever, since twisted mass lattice QCD breaks parity, theyin fact mix with each other. Taking into account of thismixing is crucial. One needs to go through the solutionof the generalized eigenvalue problem in Eq. (7) to ob-tain the optimized operator that will create the η c and χ c meson from the vacuum. Without performing thisgeneralized eigenvalue separation, it is found that thecorrect signal of χ c cannot be observed in the two-pointfunctions. The signal of η c can of course be observedeven without considering this mixing effect since it is thelightest state under consideration. B. Mass spectra for η c and χ c The eigenvalue λ n in Eq. (7) corresponding with thecorresponding meson state, i.e. n = 1 represents η c me-son, and n = 2 represents χ c meson . Since, we use theanti-periodic boundary condition, λ n ( t, p f ) t (cid:29) −→ | Z M | E M ( p f ) e − E M ( p f ) · T × cosh (cid:20) E M ( p f ) · (cid:18) T − t (cid:19)(cid:21) . (14)In practical, we use eigenvalue λ with p f = (0 , , Z η c , and the explicit value ofthis factor is 0 . . B andEns. C , respectively. While, we use eigenvalue λ with p f = (0 , ,
0) to fit the spectral weight Z χ c , and theexplicit value of this factor is 0 . . B and Ens. C , respectively. It is easily seenthat the mass can also be extracted from,cosh( m n ) = λ n ( t −
1) + λ n ( t + 1)2 λ n ( t ) , (15)with m being the mass of η c meson and m being thatof χ c , respectively. The effective mass plateaus of thesemesons for Ens. B and Ens. C are illustrated in Fig. 1.From these mass plateaus, the masses of the mesons aredetermined and the statistical errors are obtained usingjackknife method. The numerical results for the massesare summarized in Tab. II. Note that the mass values for η c are utilized to fix the valence charm quark massparameter aµ c . Therefore, only the mass of χ c are pre-dictions from this lattice computation. TABLE II: Mass values for η c and χ c on Ens. B and Ens. C respectively. The last line cites the corresponding result fromPDG [19]. η c χ c Ens. B : Mass[MeV] 2978(1) 3454(6)Ens. C : Mass[MeV] 2970(1) 3416(13)PDG: Mass[MeV] 2983.9(5) 3414.71(30) In principle, glueball states with the same quantumnumbers are also present in the similar energy range [20].However, in this lattice calculation, we have only utilizedthe quark bilinear operators for the charmonium statesand have not observed the sign of the glueballs.
C. Renormalization factor Z V of electromagneticcurrent operators The current operators in Eq. (8) are electromagneticcurrent operators. In principle, they contain all flavorsof quarks weighted by the corresponding charges. Lightquark flavors will only enter the question via discon-nected diagrams which are neglected in this study. Onlyconsidering the charm quark, we only need to considerthe current ¯ cγ ρ c ( x ). A subtlety in the lattice computa-tion is that, with c ( x ) / ¯ c ( x ) being the bare charm/anti-charm quark field on the lattice, composite operatorssuch as the current j ρ ( x ) = Z V ¯ c ( x ) γ ρ c ( x ) needs an extramultiplicative renormalization factor Z V which can beextracted by the ratio of two-point function with respectto the three-point function for η c [21]: Z µV = p µ E ( p ) (1 / (2) η c Γ (3) η c , (16)where µ is the Dirac index and we take it to be zero.Γ (2) η c and Γ (3) η c are two point correlation function and threepoint correlation function relevant with η c . The explicitforms are:Γ (2) η c = (cid:88) x e − i p · x (cid:104)O η c ( x , t ) O † η c (0 , (cid:105) (17)Γ (3) η c = (cid:88) x , y e − i p f · x e i q · y (cid:104)O η c ( x , t f )¯ cγ c ( y , t ) O † η c (0 , (cid:105) (18)with O † η c and O η c creating and annihilating a state withthe quantum number of η c meson, respectively. Actually,we use the simple local operator,i.e.¯ cγ c . According toEq. (16), we can obtain the multiplicative renormaliza-tion factor Z V , and show it in Fig. 2. The values of the Ey( c0 )=3454(6)MevEy( c )=2978(1)Mev (a) Ens. B Ey( c0 )=3416(13)MevEy( c )=2970(1)Mev (b) Ens. C FIG. 1: Mass plateaus of η c (blue point) and χ c (red point) for Ens. B and Ens. C . The horizontal line segments denote thecorresponding mass plateaus. renormalization factor Z V are 0 . . B and Ens. C , respectively. D. Form factors for two photon decay of η c and χ c To compute the relevant matrix element in Eq. (3), weplace the meson at a fixed sink position t f which is cho-sen to be 24 for Ens. B and 32 for Ens. C , respectively.These sinks are then used as a sequential source for abackward propagator inversion. This allows us to inves-tigate all possible source positions t i . We can then freelyvary the values of ω , Q and inspect the integrand asa 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-tegrand for the insertion positions t = 4 , , , ,
20 forEns. B and t = 4 , , , , , ,
28 for Ens. C with aparticular value (cid:126)p f = (000) for η c and χ c .The computation has to cover the physical interest-ing kinematic regions. For this purpose, we have toscan the corresponding parameter space of the two vir-tualities Q and Q . Basically, we follow the follow-ing strategy: we first fix the four-momentum of η c and χ c , p f = ( E, p f ) and place it on a given time-slice t f = T /
2. In this simulation, we only compute the caseof p f = (0 , , E is simply the mass of η c or χ c meson. Then, we judiciously choose several values of vir-tuality Q around the physical point Q = 0. To bespecific, we picked the range Q ∈ [ − . , +0 .
1] GeV onEns. B and Q ∈ [ − . , +0 .
1] GeV on Ens. C , whichsatisfies the constraint Q > − m ρ [9]. For a given p f , achoice of q completely specifies q due to p f = q + q .Therefore, we take several choices of q = n (2 π/L ) bychanging three-dimensional integer n . Then the energyof the first photon is also obtained using either the con- tinuum or the lattice dispersion relations: ω = q − Q , (19)4 sinh ( ω /
2) = 4 (cid:88) i sin ( q i / − ˆ Q , (20)where ˆ Q = 4 sinh ( Q /
2) is the lattice version for thevirtuality. It turns out that we can also compute thevirtuality of the second photon, since both ω and q areconstrained by the energy-momentum conservation. Onehas to make sure that, the values of Q thus computed dosatisfy the constraint Q > − m ρ otherwise it is omitted.This procedure is summarized as follows:1. Judiciously choose several values of Q in a suitablerange. We picked 7 values of Q ;2. Pick different values of n such that q = n (2 π/L ). As described above, this fixes both ω and Q using energy-momentum conservation.This is done using either the continuum or the lat-tice dispersion relations. To be specific, for each Q , we picked 4 different q ;3. Make sure all values of Q , Q > − m ρ , otherwisethe choice is simply ignored;4. For each validated choice above, compute the three-point functions (8) and obtain the hadronic matrixelement using Eq. (3).In such a way, we have obtained altogether 28 points onthe ( Q , Q ) plane around the origin. As an example, thedistribution of these virtualities for the two mesons areshown in Fig. 3 for the case of lattice dispersion relations.One could also do the same thing using the conventionalcontinuum dispersion relations. The difference of these (a) Ens. B (b) Ens. C FIG. 2: Renormalization factors Z V for Ens. B and Ens. C . The horizontal line segments denote the corresponding valuesfor Z V . -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 Q -0.15-0.1-0.0500.050.10.150.20.25 Q Ens. C -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 Q -0.1-0.0500.050.10.150.20.25 Q Ens. C FIG. 3: The distribution of virtualities ( Q , Q ) (lattice units) for Ens. C computed for η c (left panel) and χ c (right panel). two treatments finally will provide us with an estimatefor the finite lattice spacing error of the calculation.In our real lattice QCD computation, the integrationof t i in Eq. (3) are replaced by discrete summation over t i using trapezoid rule. The resulting values exhibit aplateau behavior with respect to t which is then utilizedto extract the corresponding form factor. In such a way,we have obtained numerical results for F ( Q , Q ) and G ( Q , Q ) at 28 different points in the plane of the twovirtualities. As an example, the form factor plateaus for η c are illustrated in Fig. 6 for the case of Q = 0 GeV .The corresponding case for χ c are shown in Fig. 7. E. Fitting of the form factor and the physicaldecay widths
To obtain the physical decay width, we only need thevalues of the form factors at the physical photon point,namely Q = Q = 0. As we have seen in Fig. 3 thedistribution of 28 data points in the ( Q , Q ) plane, wecould implement cuts in the plane. For a given value of Q >
0, we select the points that satisfy the followinginequality: (cid:113) ( Q ) + ( Q ) ≤ Q . (21) t/a=4t/a=8t/a=12t/a=16t/a=20 Q = 0 GeV (a) Ens. B t/a=4t/a=8t/a=12t/a=16t/a=20t/a=24t/a=28 Q = 0 GeV (b) Ens. C FIG. 4: The integrand for insertion positions obtained from simulations on Ens. B (left figure), and Ens. C (right figure)respectively. We take n = ( − − − n f = (0 0 0) as an example. The insertion positions for lattice size: 24 ×
48 andlattice size: 32 ×
64 are t = 4 , , , ,
20 and t = 4 , , , , , ,
28 respectively. t/a=4t/a=8t/a=12t/a=16t/a=20 Q = 0 GeV (a) Ens. B t/a=4t/a=8t/a=12t/a=16t/a=20t/a=24t/a=28 Q = 0 GeV (b) Ens. C FIG. 5: The integrand for insertion positions obtained from simulations on Ens. B (left figure), and Ens. C (right figure)respectively. We take n = (00 − n f = (0 0 0) as an example. The insertion positions for lattice size: 24 ×
48 and latticesize: 32 ×
64 are t = 4 , , , ,
20 and t = 4 , , , , , ,
28 respectively.
Obviously, taking a large enough Q will allow all ofthe points be considered while selecting a smaller valuefor Q will take into account only those points whosedistances are closer than Q . On the other hand, fora given value of Q , we could utilize different fittingformula to obtain the corresponding values of the form factor. Since it is the value at the origin that is directlyrelated to the decay width, it is natural to use a polyno-mial type of fit in both Q and Q . What is more, dueto boson symmetry, this function needs to be symmetricwith respect to Q and Q . Therefore, by varying thecut value Q and various orders of polynomials in the n =(0 0 -3)n =(0 -1 -3)n =(-2 -2 -2)n =(-1 -1 -2) Ens. B n =(0 0 -3)n =(0 -1 -3)n =(-2 -2 -2)n =(-1 -1 -2) Ens. C FIG. 6: The plateaus of the pseudoscalar form factors obtained from integration of t i for three point function G µν ( t i , t ) withEns. B (left figure), and Ens. C (right figure) respectively. We take Q = 0 GeV ; n f = (0 0 0) as an example. The errors inthese figures are estimated using the conventional jack-knife method. n =(0 0 -3)n =(0 -1 -3)n =(-1 -1 -3)n =(0 -2 -2) Ens. B n =(0 0 -3)n =(0 -1 -3)n =(-1 -1 -3)n =(0 -2 -2) Ens. C FIG. 7: The plateaus of scalar form factor obtained from integration of varies t i for three point function G µν ( t i , t ) with Ens. B (left figure) and Ens. C (right figure) respectively. We take Q = 0 GeV ; n f = (0 0 0) as an example. The errors in thesefigures are estimated using the conventional jack-knife method. virtualities, we could investigate the values of the formfactors at the physical point.To be specific, we adopt a polynomial ansatz up to Q and Q to the third power to fit the data of the formfactor. For η c meson, we use: F ( Q , Q ) = a + a ( Q + Q )+ a ( Q + Q ) + a Q Q + a ( Q + Q ) + a ( Q Q + Q Q ) . (22)and a similar form for the χ c form factor G ( Q , Q ).Note that a ≡ F (0 ,
0) is the form factor at the phys-ical 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 and Q .have also been attempted. Note that this implies that weare fitting the data points with 2, 4 and 6 parameters,respectively since terms at the same orders of Q and Q should be included or excluded on the same footing.In all cases, correlated fits are performed. Depending onthe number of points taken into account which is con-trolled by Q , and the number of terms kept in thefitting polynomial, we finally arrive at the values for theform factors at the origin, namely F (0 ,
0) = a for bothensembles. Similar procedures have been implemented aswell for χ c , resulting in the values for G (0 , η c and χ c on Ens.C1 using lattice dispersion relations areshown. Here the horizontal axis denotes the cut val-ues Q while the vertical axis indicates the values for F (0 ,
0) or G (0 ,
0) together with the errors (data pointswith error-bars). For each fixed value of Q , we haveperformed three fits with 2, 4 and 6 parameters. Thesethree points are slightly shifted horizontally so that theyare recognizable. The corresponding values of χ / dof foreach fit are also shown as points without error-bars. Byinspecting the plot, we get a feeling about the consis-tency and quality of these different fits and the differ-ences among the values of F (0 ,
0) can also offer us anestimate 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, wepick the case with lowest value of χ / dof as thefinal result for F (0 ,
0) together its statistical errorin this particular case. • We further attribute a systematic error arising fromthe fitting procedure by taking the largest differ-ence in the central values of F (0 ,
0) with compara-ble χ / dof. This then yields the result for F (0 , • By comparing the difference in F (0 ,
0) between thetwo different dispersion relations, we further assigna systematic error, which is taken to be the dif-ference between the two values, arising from thelattice spacing. • Similar procedure can be applied to χ c as well.In such a way, we have obtained the results of F (0 , G (0 ,
0) on two ensembles which are explicitly listedbelow: F (0 , B = 0 . , (23) F (0 , C = 0 . , (24) G (0 , B = 0 . , (25) G (0 , C = 0 . , (26)In these expressions, the first error is statistical, the sec-ond is the error from the fitting procedure describedabove and the third one is the finite lattice error esti-mated from using two different dispersion relations. It isclearly seen that, in all cases, the results are dominated by systematic errors, especially the finite lattice spacingerrors. In fact, two results on two different ensembles areconsistent within this estimate of finite lattice spacingerrors. We therefore decide not to make any continuumextrapolations. Obviously, computations with more val-ues of lattice spacings will be crucial to nail down theselarge lattice spacing errors.To compare with previous lattice computations, we no-tice that the result for Γ( η c → γγ ) is slightly larger thanprevious result presented in Ref. [9] using the same setof configurations. This difference might come from thefact that the mixing of the η c and χ c in the twistedmass setup was not fully disentangled in the previouscalculation in Ref. [9]. It is found that, if we were notusing properly chosen operator that mix with both the η c -like and χ c -like interpolating operators, we would notbe able to observe the correct χ c signal as we have dis-cussed at the end of subsec III A.Finally, let us convert the results in the form factorsinto corresponding ones in the decay widths. We simplyadd all the errors in the form factors in quadrature andneglect the errors in the mass of the mesons. This thenleads to the following results for the decay widths:Γ( η c → γγ ) B = 1 . , Γ( η c → γγ ) C = 1 . , Γ( χ c → γγ ) B = 1 . , Γ( χ c → γγ ) C = 0 . , (27)These are to be compared with the following values givenby PDG: Γ( η c → γγ ) P DG = 5 . , Γ( χ c → γγ ) P DG = 2 . , (28)These numbers are all in the same ballpark as the exper-imental ones though still somewhat smaller. However,since no controlled continuum extrapolations have beenperformed yet, it is still premature to draw any conclu-sions for the discrepancy. Our large estimated finite lat-tice errors offer some hint that this might be the majorsource of errors. In the future, more studies at differ-ent lattice spacing s are needed in order to control thelattice artifacts in a systematic fashion. Another sourceof systematic error could come from the mixing with thenearby glueball states. So far, no lattice calculations haveconsidered such an effect. Considerable efforts are neededin the future in order to bring these into account. IV. CONCLUSIONS
In this exploratory study, we calculate the two-photondecay width for η c and χ c using unquenched N f = 2twisted mass fermions. The computation is done withtwo lattice ensembles (coarser and finer) at two differentlattice spacings. The mass spectrum for the η c and χ c parameters:6parameters:4parameters:2 (a) Ens. C η c parameters:6parameters:4parameters:2 (b) Ens. C χ c FIG. 8: The fiting results for η c and χ c on Ens.C1 using lattice dispersion relations. The horizontal axis denotes the cut values Q while the vertical axis indicates the values for F (0 ,
0) or G (0 ,
0) together with the errors (data points with error-bars).The integers in brackets along the horizontal axis indicate the number of data points below Q . Points without errors showthe corresponding values of χ /dof, the values of which is obtained to the right edge of each box. meson state are obtained by solving a generalized eigen-value problem which disentangles parity mixing betweenthe two mesons.Our results for the decay width Γ( η c → γγ ) andΓ( χ c → γγ ) are summarized in Eq. (27) for two ensem-bles utilized in this computation. With only two ensem-bles we only estimate the finite lattice spacing errors foreach ensemble and no continuum extrapolations are per-formed. Albeit without the continuum extrapolations,our results are in the right ballpark as the PDG valuesshown in Eq. (28).In the future, lattice calculations with more values oflattice spacings are definitely needed so as to control thefinite lattice spacing errors which appear to be a domi-nant source. Meanwhile, the disconnected contributionsis a good supplement to this study. Possible mixing withthe gluonic excitations should also be studied. It’s alsohelpful to research on unquenched configurations withother ensembles or other methods. We also expect moreprecise experiments on double-photon decays of charmo-nium. Acknowledgments
The authors would like to thank the European TwistedMass Collaboration (ETMC) to allow us to use their gauge field configurations. Our thanks also go to Na-tional Supercomputing Center in Tianjin (NSCC). Thiswork is supported in part by the National Science Foun-dation of China (NSFC) under the Project No. 11505132,No. 11335001, No. 11275169, No. 11405178, No.11575197, No. 11935017, No. 11575196, No. 11875169,No. 11775229. It is also supported in part by the DFGand the NSFC (No. 11261130311) through funds pro-vided to the Sino-Germen CRC 110 symmetries and theEmergence of Structure in QCD. This work is also fundedin part by National Basic Research Program of China(973 Program) under code number 2015CB856700. M.Gong and Z. Liu are partially supported by the YouthInnovation Promotion Association of CAS (2013013,2011013). This work is also supported by the ScientificResearch Program Funded by Shaanxi Provincial Edu-cation Department under the grant No. 19JK0391, andNatural Science Basic Research Plan in Shaanxi Provinceof China (Program No. 2019JM-001). [1] W. I. Eshraim and C. S. Fischer, Eur. Phys. J. A , 139(2018), arXiv:1802.05855 [hep-ph] [2] X. Garcia i Tormo, Mod. Phys. Lett. A28 , 1330028 (2013), arXiv:1307.2238 [hep-ph][3] J. M. Cornwall, Mod. Phys. Lett. A28 , 1330035 (2013),arXiv:1310.7897 [hep-ph][4] J. P. Lees et al. (BaBar), Phys. Rev.
D81 , 052010 (2010),arXiv:1002.3000 [hep-ex][5] T. N. Pham, Nucl. Phys. Proc. Suppl. , 291 (2013)[6] V. Savinov (Belle), Nucl. Phys. Proc. Suppl. , 287(2013)[7] M. Ablikim et al. (BESIII), Phys. Rev.
D85 , 112008(2012), arXiv:1205.4284 [hep-ex][8] R. Frezzotti, S. Sint, and P. Weisz (ALPHA collabora-tion), JHEP , 048 (2001), arXiv:hep-lat/0104014[hep-lat][9] T. Chen et al. (CLQCD), Eur. Phys. J.
C76 , 358 (2016),arXiv:1602.00076 [hep-lat][10] J. J. Dudek and R. G. Edwards, Phys. Rev. Lett. ,172001 (2006), arXiv:hep-ph/0607140 [hep-ph][11] C. J. Shultz, J. J. Dudek, and R. G. Edwards, Phys. Rev. D91 , 114501 (2015), arXiv:1501.07457 [hep-lat][12] B. Blossier et al. (ETM Collaboration), Phys.Rev.
D82 ,114513 (2010), arXiv:1010.3659 [hep-lat][13] Y. Chen, D.-C. Du, B.-Z. Guo, N. Li, C. Liu, et al. , Phys.Rev.
D84 , 034503 (2011), arXiv:1104.2655 [hep-lat][14] A. M. Abdel-Rehim, R. Lewis, R. Woloshyn, andJ. M. Wu, Phys.Rev.
D74 , 014507 (2006), arXiv:hep-lat/0601036 [hep-lat][15] B. Blossier et al. (European Twisted Mass Collabora-tion), JHEP , 020 (2008), arXiv:0709.4574 [hep-lat][16] B. Blossier et al. (ETM Collaboration), JHEP , 043(2009), arXiv:0904.0954 [hep-lat][17] R. Frezzotti and G. Rossi, JHEP , 070 (2004),arXiv:hep-lat/0407002 [hep-lat][18] M. Kalinowski and M. Wagner, Phys. Rev.