Lattice QCD and QCD Sum Rule determination of the decay constants of eta_c, J/psi and hc states

We compute the decay constants of the lowest ccbar-states with quantum numbers J(PC)=0(-+) [eta_c], 1(--) [J/psi], and 1(+-) [hc] by using lattice QCD and QCD sum rules. We consider the coupling of J/psi to both the vector and tensor currents. Lattice QCD results are obtained from the unquenched (Nf=2) simulations using twisted mass QCD at four lattice spacings, allowing us to take the continuum limit. On the QCD sum rule side we use the moment sum rules. The results are then used to discuss the rate of eta_c -->gamma gamma decay, and to comment on the factorization in B -->X K decays, with X being either eta_c or J/psi.


Introduction
Charmonium systems provide us with a playground for understanding the features of quark confinement, for testing the validity of various quark models, and for describing processes that are interesting for the weak interaction phenomenology as well as for the search of physics beyond the Standard Model (BSM) [1]. Most of the quark models aim at describing the spectrum of charmonium states, including the orbital and radial excitations. Not many of these models, however, are used to describe the hadronic matrix elements as that requires a more detailed knowledge about the non-perturbative QCD dynamics of hadronic confinement.
The method of QCD sum rules (QCDSR) was first tested on charmonium systems, and the fact that a number of J P C = 1 −− charmonium excitations have been detected and their electronic widths measured, was actually used to fix some of the QCDSR parameters relevant to the non-perturbative QCD effects expressed in terms of power corrections (QCD vacuum condensates) [2,3,4]. In this paper we report on our results concerning the simplest matrix elements related to three charmonium states (η c , J/ψ, h c ) and focus on four decay constants (f ηc , f J/ψ , f T J/ψ , f hc ) defined via, where the µ-dependence of the couplings to the tensor current indicates the renormalization scale and scheme dependence. Of the above couplings only f J/ψ can be directly extracted from experiment via Γ(J/ψ → e + e − ) = 4πα em 3m J/ψ The other couplings are not as directly related to experiment but they are still very relevant for phenomenology. For example, f ηc enters decisively in the theoretical description of the γ * γ * → η c decay form factor, and of Γ(η c → γγ) in particular [5,6,7,8,9]. Similarly, the phenomenological studies of the small-x gluon distribution function from the inclusive production of η c requires the knowledge of f ηc [10]. Furthermore, such couplings can be helpful in describing the non-leptonic B-decays and to check for deviations between the measured and the results obtained by using the factorization approximation. For example, by combining the following decay modes [11], one can use f ηc /f J/ψ and the known information about the B → K form factors to check for validity of the factorization approximation. Otherwise, by imposing the factorization, one can get a useful information about the ratio of B → K form factors. On the other hand, a measurement of a non-zero B(B + → h c K + ), which is currently only bounded from above (B(B + → h c K + ) < 3.8 × 10 −6 [11]), could be interpreted as either a measurement of the deviation with respect to the factorization approximation, or a signal of the presence of coupling to the tensor operator that might appear only in the case of physics BSM. Finally, the coupling f T J/ψ may be interesting when checking for the presence of the New Physics operators in various processes.
In the first part of this paper we will discuss the computation of f ηc , f J/ψ , f T J/ψ , f hc by using the QCDSR. As we shall see the approximation of 'one resonance plus continuum', that we employ on the phenomenological side of sum rules, results in sizable error bars on the decay constants. In the second part, we compute the same quantities by means of numerical simulations of QCD on the lattice. Finally we compare our results and make a brief discussion of the impact of our results on two topics in phenomenology.

Two-point QCD sum rules
To estimate the hadronic properties of charmonium systems (masses and decay constants) by means of QCDSR one needs to compute the two-point correlation functions. Here we will focus to the following three: where V µ =cγ µ c, P = 2m c icγ 5 c, and T µν =cσ µν c, with σ µν = i/2 × [γ µ , γ ν ]. In terms of the Lorentz scalars the vector and tensor correlation functions can be written as: where the projectors separate the even and odd parity parts, and therefore Π + (q 2 ) will be used to discuss the h c (1 +− ) channel, while Π − (q 2 ) the ordinary J/ψ(1 −− ) state. Note that P i P j = 12q 4 δ ij in d = 4 dimensions. For the perturbative part, each of the invariant functions Π i (q 2 ) (i = P, V, +, −) satisfies the dispersion relation, with a suitable number of subtractions. Each spectral function, ρ i (s), is then computed in perturbation theory and can be written as, where the scale dependence is kept implicit. Besides the perturbative contribution, to the above Π i (q 2 ) one also needs to add the non-perturbative terms. The leading nonperturbative contributions to the correlation functions involving charmonia are power corrections proportional to the gluon condensate, αs π G a µν G µν a ≡ αs π G 2 , namely, where the Wilson coefficients C G i (Q 2 ) ∝ 1/Q 2n i are also computable perturbatively, with n i > 0 depending on the operators used.
Complete expressions for the spectral functions ρ pert i (s) as well as for gluon condensate contributions C G i (Q 2 ) are collected in Appendix, where a brief discussion about the calculation can be found. 1 When studying charmonia it is convenient to use the so-called moment sum rules [2,3,4]. One starts by defining the moments of eq. (7), at some spacelike Q 2 0 , far from the resonance region. In practice Q 2 0 is a parameter that is to be adjusted in order to improve the convergence of the integral on the right hand side (r.h.s.) of the above equation. Since the mass of the charm quark is large with respect to Λ QCD it is customary to use Q 2 0 = 4m 2 c ξ, and by changing the integration variable in the dispersion relation, s → v 2 = 1 − 4m 2 c /s, the theory part of the n th moment can be written as, On the other hand, the same moments (10) can be expressed in terms of hadronic quantities. By inserting all possible hadronic states H in the correlators (4) that can couple to each of the above operators, one can write where the sum runs over all possible single or multiparticle hadronic states, and J i stands for a generic bilinear quark operator. The situations in which the masses and couplings of the higher excited states in the sum (12) are experimentally established are extremely rare. A notable example is that of the first few J P C = 1 −− states for which both the masses and electronic widths, Γ(ψ(nS) → e + e − ), have been measured. This information was used to fix the value of the gluon condensate in ref. [2], and then further refined in ref. [12]. 2 In the most phenomenologically relevant situations, however, only the position of the first pole in the sum (12) is known, whereas for the rest of the sum one invokes the quark-hadron duality and replaces them by the spectral function ρ pert i (s) in the dispersion relation, starting from some threshold, . After using the definitions (1), we have where the renormalization scale is chosen to be µ 2 = m 2 c + Q 2 0 , with m c ≡ m MS c (m c ). After equating eqs. (11) and (13) we can define where v[s 0 ] = 1 − 4m 2 c /s 0 , so that For a recent review concerning the various estimates of the gluon condensate, please see ref. [13].
In other words the masses are obtained from the ratios of moments, while the decay constants are computed from one or several moments separately. Before discussing the practical procedure we use to get the results for the decay constants we need to stress that: (1) O(α s ) corrections to the functions Π ± (Q 2 0 ) are new. In ref. [3] the authors computed ρ   (2) Our results for C i G (Q 2 ) agree with those presented in refs. [3,14,15]. Here again, the result for C − G (Q 2 ) is new.

Evaluation of Sum Rules
In this section we discuss the evaluation of QCDSRs given in eq. (15). Our strategy for all sum rules, except the one for h c , consists in requiring that the mass of the lowest lying hadron obtained from the ratios of moments coincides with the experimentally established value to less then or equal to 1%. Only in the case of the sum rule for h c we allow that agreement to be within 5%. For the reader's convenience we quote the masses of the three lowest lying states we discuss in this paper [11]: In evaluating the left hand side (l.h.s.) of eq. (15) we take into account the charm quark mass and the value of the gluon condensate from ref. [12], that are found to be highly correlated (cf. fig. 5 in ref. [12]). 3 We take that correlation into account and also vary the threshold parameter s 0 above the square of the mass of the lowest state and its first radial excitation. More specifically, While m exp.   [11], the first radial excitation of the h c state could be extracted from lattice QCD study of ref. [17], m latt. h c = 3.639(1) GeV. With the sum rule parameters [m MS c (m c ), αs π G 2 , s i 0 ] varied in the intervals indicated above, we then look for the moments n such that δm QCDSR ηc,J/ψ /m exp. ηc,J/ψ ≤ 1% and δm QCDSR hc /m exp.
hc ≤ 5%. Furthermore we impose the standard QCDSR requirements, namely that the nextto-leading order correction to the moments represents less than 30% with respect to the leading order term, and that the contribution coming from the gluon condensate does not exceed 50% of the perturbative part. For the former requirement it is important to work with ξ = 0. We actually checked that for two values, ξ = 1 and ξ = 2, the range of values for the moments is such that the above criteria are fulfilled and the resulting values for the decay constants remain unchanged. The only exception is the sum rule for h c , for which the mass of h c , as obtained from the ratios of moments, is always larger than the physical one. This excess is less than 5% only for lower moments and for lower values of the threshold parameter s 0 . For that reason in the discussion of our results for m hc and f hc we will vary s hc 0 ∈ [3.6 2 , 3.8 2 ] GeV 2 . We attempted to enlarge the window in s 0 , but the impact on our final results was only marginal. We should stress that the value of n is not fixed to be common to all s 0 , but they were found for each s 0 separately. Therefore this somewhat implicitly corresponds to a strategy adopted in ref. [18] in which the threshold parameter was considered to be a function of the Borel parameter (or equivalently of n in the case of the moment sum rules). We do not introduce any extra parameter but verify that s 0 and n are indeed strongly correlated.
For example, and by using the central value of the charm quark mass and of the gluon condensate (17), and by varying s i 0 , we find that all the above criteria are satisfied for The above ranges of n are found for ξ = 2. They change with the value of ξ and for larger ξ the values of n satisfying our criteria become larger. With these values of n we then compute the decay constants. Illustration of the stability of the sum rule results is provided in fig. 1. We should note that each decay constant is highly sensitive to the mass of the hadron. To make the procedure fully self-consistent, in the evaluation of the sum rule for each decay constant we use the corresponding hadron mass obtained by the same sum rule. Had we used the physical mass of the hadron state instead of the one obtained from the sum rule, the resulting curves in fig. 1 would be considerably flatter. It turns out that the variation of the threshold parameter s 0 already covers most of the allowed values for the decay constants that are shown by the shaded areas in fig. 1. Note that these shaded intervals in fig. 1 are obtained by varying all of the QCD sum rule parameters: s 0 , n, m MS c (m c ), αs π G 2 , and for ξ ∈ {1, 2}. Another important comment is that we take into account the correlation between the charm quark mass and of the gluon condensate found in ref. [12]. In that latter paper the values of the charm quark mass and of the gluon condensate have been obtained from the vector-vector sum rule by using the three-loop perturbative expressions, and by including the loop corrections to the Wilson coefficient that multiplies the gluon condensate. More importantly, as far as the stability of the results is concerned, a rich experimental information about the spectral function in the M phen. V n (Q 2 0 ) has been included. We do not aim at that level of accuracy. Instead we content ourselves by working with the two-loop QCD expressions on the perturbative side, and only one resonance has been included to the hadronic side of the sum rules before evoking the quark-hadron duality. For that reason the expected accuracy of the sum rules on the decay constant will be relatively modest. Indeed we get therefore with about 10% uncertainty, which is a typical accuracy of the sum rule computation of the hadronic decay constants [4]. Note again that 335 MeV and 447 MeV correspond to the minimal and maximal value of f J/ψ obtained from the QCD sum rule after varying all the parameters in a way described above.
As for f ηc we find which is somewhat lower than f ηc = 356 (16) MeV found in ref. [19] in which the authors made additional assumptions about the contributions to the phenomenological side of the spectral function coming from the radially excited η c , or f ηc = 346(33) MeV found in ref. [?, 20] where the threshold parameter was assumed to be much larger than the lowest radially excited states. 4 In the earlier sum rule estimates another definition has been used, related to ours via g ηc = f ηc /2m c [21]. That definition is renormalization scale dependent but since the authors of ref. [21] used only the leading order expressions for the perturbative part of the spectral function, the choice of the scale and scheme could not be specified. Note also that in the past the computations were often done by using the pole charm quark mass, so that the approximation 2M c ≈ m ηc was justified. In that way the resulting value for f ηc was larger. Finally, the recent computation of this decay constant by using the Borel sum rule and somewhat different criteria for the choice of the sum rule parameters, lead to a much larger value [22].
Concerning the coupling f T J/ψ (µ), the discussion of this sum rule is qualitatively and quantitatively very similar to that of f J/ψ obtained from the vector-vector correlation function, and the moments for which the criteria discussed above are satisfied is essentially the same. We obtain, in the MS scheme. Our result agrees with f T J/ψ = (408 ± 26) MeV, presented in ref. [23]. In contrast to the latter paper we could also specify the renormalization scale at which f T J/ψ (µ) is defined because we included in our calculation the next-to-leading order QCD correction to the spectral function (cf. discussion in the appendix of the present paper).
Since the behavior of f T J/ψ (µ) with respect to the variation of the QCD sum rule parameters is very similar to that of f J/ψ , it is more judicious to compute the ratio of the two, which is illustrated in fig. 1 where we see that the ratio R T J/ψ is much more accurately estimated. By following the same criteria discussed above, and by using µ = 2 GeV, we get We emphasize that this is the first determination of f T J/ψ and R T J/ψ . Finally, our result for the decay constant of the h c state, is We stress that this result is obtained for low moments and that for larger moments the sum rules progressively deteriorates in the sense that the mass of the lowest lying state becomes much larger and the decay constant much smaller. To our knowledge, up to now, the only QCDSR analysis of the h c state has been made in ref. [24] in which the author reported f hc = 490(60) MeV, in clear disagreement with our value (25). We compared our expressions with those presented in ref. [24], and realized that the author of ref. [24] has calculated only a part of Π + (q 2 ) obtained using only the first part of P + µνρσ in (6), and therefore his result does not correspond to any physical state.

Lattice QCD results
We now compute the same quantities discussed above but by means of numerical simulations of QCD on the lattice. To that end we use the gauge field configurations generated by European Twisted Mass Collaboration (ETMC), in which the effect of N f = 2 dynamical ("sea") light quarks have been included by using the Wilson regularization of QCD on the lattice with the maximally twisted mass term, namely [25] where ∇ µ (∇ * µ ) stands for the forward (backward) covariant derivative, m cr is the critical mass term tuned to restore the chiral symmetry of the massless action, otherwise broken by the Wilson term (also in the brackets), and µ c is the bare charm quark mass. In the above action ψ(x) = [c(x) c (x)] T is a doublet of the charm quark field and its replica. The factor iγ 5 τ 3 r cures the pathology of the standard Wilson quark action by rotating the Wilson term to the imaginary axis which is why one can simulate with sea quark masses considerably closer to the chiral limit. The quark propagators S c (0, 0; x, t) and S c (0, 0; x, t) are then obtained by inverting the above Wilson-Dirac operator with r and −r, respectively. In practice r = 1. Finally, we should mention that the action (26) refers to the valence charm quarks, but the same one is used to generate the gauge field configurations but with µ c → µ q , mass of the light sea quark. Detailed information about the lattices used in this work are given in tab. 1.
Hadron masses and decay constants are extracted from the study of the two-point correlation functions with operators chosen with desired quantum numbers, namely: In the above expressions the dependence of the renormalization constants on the bare lattice coupling is implicit, namely Z A ≡ Z A (g 2 0 ), and Z T (µ) ≡ Z T (g 2 0 , µ). Notice also that the above definition of the pseudoscalar operator P is renormalization scale and scheme invariant both in the continuum and on the lattice with twisted mass QCD. To extract masses and decay constants one studies the large time separation between the operators in the two-point correlation functions. More specifically, Note that the action is written in the "physical basis" and not in the twisted one.   Table 1: Summary of the lattice ensembles used in this work (more information can be found in ref. [28]).
Data obtained at different β's are rescaled by r 0 /a, and the overall lattice spacing is fixed by matching f π computed on the lattice with its physical value, leading to r 0 = 0.440 (12) fm (c.f. ref. [27]). All quark masses are given in lattice units.
where i, j ∈ {1, 2, 3}, and T stands for the size of the periodic lattice in the time direction. Since our charmonia are taken to be at rest the matrix elements (1) that appear in (28) (β, µ sea , L)

read:
0|P |η c ( 0) = f ηc m 2 ηc , In eq. (28) we assumed the local source operators, which are needed for extraction of the decay constants. In practice, however, we implement the Gaussian smearing procedure in order to increase the overlap between the interpolating operator and the lowest state coupling to a given operator. The smearing procedure and the parameters used in actual computations have been discussed in refs. [29,30].
The above matrix elements are then extracted by dividing the local-smeared and smeared-smeared correlation functions, where the coupling to the smeared correlation functions can be studied from the smeared-smeared correlation functions in a way similar to eq. (28). Hadron masses am H (H = η c , J/ψ, h c ) are extracted from the fit to a constant on the plateau of the effective mass m eff with J = P, V, T (+) respectively. The results for the masses have been presented in our previous paper [29] and, for the reader's convenience, are presented in tab. 2 of the present paper. The novelty is that we could also check that the results for the mass of J/ψ state obtained from the correlation function C T (t) coincide with those we obtain from the study of C V (t) except that the errors are about 2 ÷ 3 times larger. Notice that only the mass of m ηc is given in the lattice units while the other masses are obtained from the fit to a constant of the ratios in which the statistical uncertainties cancel to a large extent. As for the decay constants, they are extracted in a way indicated in eq. (28) and by using the definitions (29). Their values are converted to physical units by using the lattice spacings quoted in tab. 1, and listed in tab. 2.
To reach a physically interesting results we need to extrapolate our decay constants, obtained at four lattice spacings, to the continuum limit. Since the physical quantities computed with maximally twisted mass QCD on the lattice are automatically O(a) improved, the leading terms are those proportional to a 2 . Furthermore, at each lattice spacing we computed the charmonium couplings for several values of the dynamical light quark masses, which is useful in order to check on their dependence on the sea quark mass, like we did in our previous paper where we showed that the masses of charmonia are completely insensitive to the sea quark mass [29]. To get the physically relevant result in the continuum limit, the decay constants from tab. 2 are therefore fit to the following form, where the parameter b H measures the dependence on the sea quark mass, denoted by m q ≡ m MS q (2 GeV), while the parameter c H measures the leading discretization effects. Division by a β=3.9 = 0.086 fm is made for convenience. The linear fit (32) in a 2 provides an adequate description of all our results if we leave out the data obtained at β = 3.80, as it can be appreciated from the plots provided in fig. 2. Note however that the extrapolation curve shown in fig. 2 takes into account the fact that at each lattice spacing the results are obtained for several values of the sea quark mass. The dependence on the sea quark  mass is shown in fig. 3. We therefore report the results of the fit to eq. (32) separately for the case in which the data at β = 3.8 are left out, and the results of the continuum extrapolation with all the lattice data included, cf. tab. 3. Although the quality of the fit deteriorates when all the lattice data are included, its χ 2 /d.o.f. is still acceptable and we prefer to use both results and include the difference in the estimate of the systematic uncertainty. That leads us to our final estimates: Two of the couplings discussed in this paper (f J/ψ and f ηc ) have been already computed on the lattice in an unquenched setup but with the different lattice regularization. By using the staggered quark action and by including N f = 2+1 dynamical light flavors, the authors of ref. [31] obtained f ηc = 395(2) MeV, in the continuum limit. With a similar setup, the same collaboration later reported f J/ψ = 405(6)(2) MeV [32]. Knowing that the lattice actions are very different, the fact that these results agree quite well in the continuum limit with our values (33), is a good indication of the robustness of the lattice QCD predictions. Our results indicate that there is no dependence of the charmonium quantities (masses decay constants and the form factors discussed in ref. [29]) on the light sea quark mass. The results presented in refs. [31,32] also suggest that the decay constants f ηc,J/ψ do not depend on the strange sea quark mass. Finally, we remark that the values for f T J/ψ and f hc are new.

Phenomenology
In this section we comment on two topics of phenomenological interest, already mentioned in introduction of the present paper, namely the η c → γγ ( * ) decay, and the factorization of the non-leptonic B-decays to two mesons, one of which is a charmonium.

η c → γγ ( * )
For a theoretical estimate of Γ(η c → γγ) the non-perturbative information is essential and is related to f ηc . In the standard derivation of the expression for Γ(η c → γγ) one starts from the η c → γ * γ * decay amplitude, where k 1,2 and e 1,2 stand for the momenta and polarization vectors of the two off-shell photons. One then assumes the validity of factorization of the soft QCD dynamics of η c and the hard rescattering of cc into photons. By taking one photon on-shell the other is expected to be energetic enough for factorization to be applicable. The resulting process η c → γγ * is then described by the form factor F γηc (q 2 ) ≡ F (q 2 , 0) which enters directly the expression for the η c → γγ decay rate as, The form factor F γηc (Q 2 ) can be studied experimentally through dσ(e + e − → e + e − η c )/dQ 2 (Q 2 = −q 2 > 0), a process driven by γγ * → η c . In this way, after a detailed measurement of such a process, the BaBar Collaboration was able to determine the shape of F γηc (Q 2 ) in a large energy window corresponding to Q 2 ∈ (0, 50) GeV 2 [33]. They found that the data are very well described by a single pole form, with the pole being at m pole = 2.9(1)(1) GeV. Such a pole-like behavior was predicted by (quenched) QCD on the lattice [34], and is compatible with the vector meson dominance. As for the intercept of the form factor, F γηc (0), different models give different answers [5]. For example, the perturbative QCD approach of ref. [6], results in where k 2 ⊥ is the mean transverse momentum of the c-quark with respect to the momentum of η c . After assuming k 2 ⊥ = 0, one gets the expression usually employed in the literature [5]. Similarly, the authors of ref. [7] used the heavy quark spin symmetry to estimate Γ(η c → γγ), and their expression for the form factor coincides with eq. (36) after replacing 2 k 2 ⊥ → b ηc m ηc , with b ηc = 2m c − m ηc . The latter quantity is clearly ambiguous as the quark mass is a renormalization scale and scheme dependent quantity. If one assumes m c to be the pole mass, the value of b ηc can be fixed if one knows f ηc and Γ exp. (η c → γγ). Taking b ηc = 0 (→ k 2 ⊥ = 0) reduces eq. (36) to the formula most frequently used in the literature.
In ref. [9], by imposing the local quark-hadron duality on the decay amplitudes, the authors derived a yet another expression for F ηcγ (0), namely where m V = 3.75(25) GeV has been fixed from the single pole fit to the BaBar data at large values of Q 2 . In all these expressions f ηc enters decisively and its impact on eq. (35) should be checked against the experimental data. Another possibility is to rely on the nearest vector meson dominance (VMD) hypothesis, namely [8], which can be tested since the value of the J/ψ → η c γ form factor is nowadays known from the lattice QCD studies of refs. [29,32]. By using eq. (35) or (36) we can write With the experimentally established B(η c → γγ) = (1.57 ± 0.12) × 10 −4 , and Γ(η c ) = 32.0(9) MeV, we have Γ exp (η c → γγ) = 5.0(4) keV, which together with our f ηc = 0.387 (8) GeV, allows us to deduce the value of δ = 0.15(5) GeV 2 . That value is too large to be interpreted as k 2 ⊥ = 0.81 (14) GeV, 6 and also too large to be identified as b ηc = δm ηc = 0.46 (16) GeV. Finally, we should say that the VMD is actually quite good an approximation. By using V J/ψ→ηc (0) = 1.92(3)(2) computed in ref. [29], together with our result for f J/ψ , inserted in eq. (38), for the di-photon decay width we get Γ(η c → γγ) = 6.0(4) keV. 7 We conclude that the usual expression for Γ(η c → γγ) based on factorization approximation [δ = 0 in eq. (39)] leads to the result larger than the experimental value: Γ fact. (η c → γγ) = (6.64 ± 0.27) keV, vs. Γ exp. (η c → γγ) = (5.0 ± 0.4) keV. That discrepancy can be studied in a systematic way by means of non-relativistic QCD expansion, along the lines of ref. [35]. Research in this direction, to elucidate the origin of this discrepancy, would be welcome.

Non-leptonic B decays to charmonia
By using the factorization approximation, the decay rate of the Class-II non-leptonic Bdecays in the Standard Model can be written as where the coefficient a 2 is a combination of Wilson coefficients computed in perturbation theory, encoding the information about the short distance physics. That quantity is considered as a parameter in the generalized factorization [36], that is to be obtained from the experimentally measured one decay mode and then used to describe the other modes of the given Class. By taking the ratios of the above rates, we get where With our result f ηc /f J/ψ = 0.926 (6), one can then compare the measured charged and neutral B-decay modes (3) with eq. (41) and deduce, 6 Even if one assumes this value to be correct, then one could not fit the BaBar data at large Q 2 's by using the expression F ηcγ (Q 2 )/F ηcγ (0) = 1/(Q 2 + m 2 ηc + 2 k 2 ⊥ ) [6], where the perturbative approach is expected to work better. 7 This number remains essentially unchanged if we used the lattice QCD results obtained in ref. [32], V J/ψ→ηc (0) = 1.90 (7)(1) and f J/ψ = 0.405(6) MeV. We get Γ(η c → γγ) = 5.9(5) keV.
These results are consistent with ≈ 1.44, as obtained from the QCDSR calculation near the light cone in ref. [37,38]. They are also consistent with 1.51(3) obtained in the quenched lattice QCD study of ref. [39], but not as well with 1.37 (2) recently obtained in the unquenched lattice study with non-relativistic QCD employed to treat the heavy quark [40]. So this information can be used either to get an idea on the above ratio of the form factors, or as a measure of the deviation with respect to the factorization approximation if the form factors are taken from elsewhere.

Summary
In this paper we presented results of our analysis of four decay constants of the charmonium states. By adopting the strategy of "one resonance + continuum" in the moment QCDSR analysis, we found that the values of the decay constants f J/ψ and f T J/ψ agree quite well with those obtained through the simulations of QCD on the lattice, in the continuum limit. On the other side the QCDSR results for the pseudoscalar meson decay constant f ηc are lower than those obtained on the lattice. Similar holds true for f hc , decay constant of the recently observed J P C = 1 +− charmonium state. Adding more states to the hadronic side of the sum rules helps improving the stability of the sum rules, while the value of the decay constant remains practically unchanged. One reason for disagreement of the QCDSR estimate of f ηc with that obtained on the lattice might be related to the fact that the non-perturbative contribution to the sum rules, proportional to the gluon condensate, has been fixed from the detailed analysis of the vector-vector correlation function. A possible explanation of that discrepancy is that the series of power corrections is truncated and that the higher order terms affect different correlation function differently, which is why f ηc and f hc are not as well reproduced by the QCSR as it is the case with f J/ψ and f T J/ψ . We plan to come back to that issue in the near future. In the one resonance plus continuum setup, the expected precision of the moment QCDSR estimates is of the order 10 ÷ 15%, which is what we observe with our results. Note also that the results presented in this paper for the spectral function and for the gluon condensate contributions to the correlation function of tensor densities are new. We should stress that in view of the approximations made in the method of QCDSR, the agreement of f J/ψ , f T J/ψ and even f ηc with the results obtained from lattice QCD is quite remarkable. The case of f hc is an exception, however. We did not attempt to remedy that discrepancy by adding an extra term to the series of power corrections but we plan to come back to that issue in the future.
Our lattice computation of the same set of decay constants is made in the Wilson regularization of QCD by including the maximally twisted mass term, with N f = 2 dynamical light quarks included in the gauge field configurations. From the simulations made at four different lattice spacings we were able to take the continuum limit. We find that the charmonium decay constants are insensitive to the variation of the mass of the light dynamical quarks. Non-perturbatively computed renormalization constants were implemented in our computation, and our final results are: where we combined the statistical and systematic errors in the quadrature.
With the above results in hands we were able to address two issues of phenomenological interest. First, and by using our f ηc , we get that the standard formula for the decay width of η c → γγ, does not reproduce the experimentally measured width, which might be an indication of the presence of non-factorizable terms. With our values for f ηc /f J/ψ we were able to check on the factorization approximation in the Class-II non-leptonic decays of B-mesons. We found that the most recent lattice results for the B → K indicate the violation of the factorization approximation, whereas those obtained by the QCDSR near the light cone as well as the older lattice results are quite consistent with what we extracted for f B→K from the ratios of non-leptonic decay channels together with our f ηc /f J/ψ . Another lattice QCD estimate of this ratio of form factors would be highly welcome.
Finally, we would like to emphasize that the results for R T J/ψ = f T J/ψ (2 GeV)/f J/ψ as obtained in our QCDSR analysis agree very well with those computed on the lattice:

Appendix: Spectral Functions and gluon condensate contributions
While the leading contributions to the spectral functions, ρ pert i (s), are easy to calculate, the O(α s ) corrections are quite demanding as they require evaluating the two-loop diagrams. To derive perturbative spectral functions ρ pert i (s) = ImΠ pert i /π one needs to calculate the imaginary part of the diagrams shown in fig. 4, with both external currents being either V µ =cγ µ c, or P = 2m c icγ 5 c, or T µν =cσ µν c. By using the standard approach, i.e. multiplying by appropriate projectors and expressing the scalar products in numerators in terms of those in denominators, one performs the tensor decomposition to the basic scalar Feynman integrals. The calculation of the relevant two-loop scalar integrals could be challenging, but since we are only interested in their imaginary part the task becomes much simpler. We computed the scalar integrals in two ways: (i) by the 'cut'-technique using the Cutkosky rules and (ii) by a directly extracting the imaginary part of the integrals from their Feynman parameter representation. Both ways lead to the same results. Since we used dimensional regularization, the above mentioned calculations were performed in ddimensions. Finally, besides the renormalization of α s , m c and the quark field, we accounted for the renormalization of the operators in the MS scheme. The above choice of the pseudoscalar density P is particularly convenient because the anomalous dimension of the icγ 5 c cancels against that of the quark mass, so that P is renormalization group invariant. Therefore, the only correlator in which one should take care of the anomalous dimension is that involving the tensor densities. A standard procedure consists in connection the bare and renormalized current via where Z 2 is the quark field renormalization constant and the anomalous dimension is derived as In the expansion it follows that γ (0) T = 2/3. Therefore after renormalization, the constants f hc and f T J/ψ are µ-dependent quantities.
In the following we give the full expressions for perturbative spectral functions, ρ pert i (s) ≡ ImΠ pert i (s)/π (i = P, V, +, −), written separately for the leading and the next-to-leading term in α s , namely, ρ pert i (s) = ρ We checked that our results for ρ pert P,V agree with those given in the literature, cf. eg. [3]. Our expression for ρ pert + agrees with a similar expression in ref. [3] derived by usingc∂ µ γ 5 c instead of the tensor densitycσ µν c. Expressions for ρ pert − are new. The nonperturbative contributions proportional to the gluon condensate are obtained by computing the diagrams shown in fig. 5. In the notation our results read as follows: 9 The above expressions are obtained by the direct calculation and agree with refs. [14,15].
The result for C G − (Q 2 ) is new.
B(v) is the same as A(u) function in [3]. 9 In the literature there is some discrepancy among results for C G i (Q 2 ) related to a different number of subtractions. For example, in ref. [3], C G P (Q 2 ) and C G S (Q 2 ) (connected with C G + (Q 2 ) here) are obtained by using two-times subtracted spectral function, namely Π i (Q 2 ) − Π i (0) − Q 2 Π i (0) . If we do the same here, our results would agree with their expressions. Similarly, in that way, C G A from ref. [3] would coincide with our C G + .
For calculation of the moments, the integral representation of the above formulas are particularly useful [2,41]. With the help of where ξ = Q 2 /(4m 2 c ), we can express all C G (Q 2 ) as