Lattice QCD study of the radiative decays J/ψ → ηcγ and hc → ηcγ

We present the results of our lattice QCD study of the hadronic matrix elements relevant to the physical radiative J/ψ → ηcγ and hc → ηcγ decays. We used the twisted mass QCD action with Nf = 2 light dynamical quarks and from the computations made at four lattice spacings we were able to take the continuum limit. Besides the form factors parameterizing the above decays we also computed: (i) the hyperfine splitting and obtained Δ = 112 ± 4 MeV, (ii) the annihilation constant fJ/ψ which agrees with the one inferred from the measured Γ(J/ψ → e+e−).


Introduction
After the observation of η b at BaBar [1,2] the radiative decays of heavy quarkonia received a significant attention in the literature. As Υ(1S) → η b γ is not yet experimentally accessible due to the smallness of phase space, the experimenters turned to studying Υ(2S) → η b γ and Υ(3S) → η b γ modes from which they then extracted m η b . A potential problem in that value is the insufficient control of theoretical uncertainties in the transition matrix elements to guarantee an accuracy at a percent level, especially when the radial excitations are involved. The corresponding transition matrix elements have been computed by using quark models [3][4][5][6]. 1 Indeed the resulting value m η b = 9390.9±2.1 MeV extracted from the BaBar experiment [1,2] was consistent with the value obtained from the similar measurements made at CLEO [8,9], leading to the following value of the hyperfine splitting [10], that turned out to be much larger than the values predicted by methods based on perturbative QCD, namely ∆ b = 44 ± 11 MeV [11], 39 ± 14 MeV [12]. The lattice QCD results are inconclusive on this issue so far, although they seem to point towards the values larger than those obtained in refs. [11,12]. For example, by using the simulations with N f = 2 + 1 staggered quarks and the Fermilab treatment of the heavy quarks on the lattice, the value ∆ latt b = 54 ± 12 MeV was obtained in ref. [13], while the non-relativistic QCD (NRQCD) treatment of the heavy quarks lead to ∆ latt b = 70 ± 9 MeV [14]. Simulations of QCD with 1 For a complete reviews with extensive lists of references please see ref. [3][4][5][6][7].
Proponents of the extensions of the Standard Model involving more than one Higgs doublet speculated that the experimentally established pseudoscalar state η b might be actually a mixture of the true η b and the light parity-odd Higgs boson A 0 [18][19][20][21]. 2 This would solve the puzzle of too large a hyperfine splitting and would give more support to a plausible solution that m A 0 ∼ 9 GeV. However, to give these speculations more support it is essential to check whether or not the hadronic matrix element used to extract m η b from the mentioned experiments coincides with the results obtained by using the methods based on QCD from first principles. For example, by using NRQCD on the lattice, the authors of ref. [23] obtained much larger values for the transition matrix elements than those inferred from the measured Υ(nS) → η b γ (n = 2, 3). Since the direct QCD simulations of the bbsystems are difficult because the lattice spacings are still too large to resolve the b-quark mass, we decided to explore the similar physics processes in the charmed systems and study, J/ψ → η c γ and h c (1P ) → η c γ. The established methodology of this paper will then be used for our future attempt to compute the amplitude for h b (1P ) → η b γ decay on the lattice. Besides methodological issues these decays are physically interesting on their own. One is the so called magnetic dipole (M1) and the other electric dipole decay (E1). Quark models fail to reproduce the measured Γ(J/ψ → η c γ) and, instead, obtain a significantly larger value [3][4][5][6]. On the other hand, h c (1P ) has been discerned from the experimental background only recently and its dominant decay is indeed h c (1P ) → η c γ, the branching fraction of which has been measured accurately.
The first extensive study of the radiative decays of charmonia on the lattice has been reported in ref. [24] where the authors computed relevant matrix elements for a number of decay channels in the quenched approximation of QCD and with one lattice spacing. That computation has been extended to the case of N f = 2 dynamical light quark flavors at single lattice spacing [25]. In this paper we will focus on J/ψ → η c γ and h c (1P ) → η c γ, for which we compute the desired form factors at four lattice spacings and then extrapolate them to the continuum limit. Those results may be used for a cleaner exclusion of the possibility of having very light m A 0 < 2m τ (see e.g. [18][19][20][21]).

Hadronic matrix elements
The transition matrix element responsible for the J/ψ → η c γ * decay reads,

JHEP01(2013)028
where J em µ = Q cc γ µ c is the relevant piece of the electromagnetic current, with Q c = 2/3 in units of e = √ 4πα em . Information regarding the non-perturbative QCD dynamics is encoded in the form factor V (q 2 ) and represents the most challenging part on the theory side. For the physical process, i.e. with the photon on-shell q 2 = 0, the decay rate is given by [24] where ∆ stands for the hyperfine splitting ∆ = m J/ψ −m ηc . When both the initial and final hadrons are at rest the matrix element (2.1) is zero by definition. The smallest momentum that can be given to a hadron on the lattice with periodic boundary conditions is 2π/L, which is very large for the lattices that we work with today and would make q 2 < 0, far from q 2 = 0. As a result we would have to work at several negative q 2 's, then model the q 2 shape of the form factor as to extrapolate to the physical point, q 2 = 0. That methodology has been adopted in refs. [24,25]. In this work, instead, we will use the so called twisted boundary conditions [26][27][28] which allow us to work directly at q 2 = 0. This is achieved by tuning the twisting angle θ 0 via the three momentum given to the pseudoscalar meson that fulfills the condition, where we use q = (1, 1, 1) × θ 0 /L. For that purpose, and for each of our lattices, we first computed the masses of m J/ψ and of m ηc , and then by using eq. (2.3) we determined θ 0 that is then used in the computation of one of the charm quark propagators. In practice this last step is made by "twisting" the gauge links according to on which the quark propagator is computed according to, S θ c (x, 0; U ) = e i θ· x/L S c (x, 0; U θ ) . (2.5) In our notation the quark propagator S c (x, 0) ≡ S c ( x, t; 0, 0) = c(x)c(0) U , and we only in eq. (2.5) we write explicitly the gauge field configuration in the argument, to distinguish U from U θ . In what follows U will be implicit. Similarly, in the computation of the physical h c → η c γ decay, we compute the transition matrix element that is parameterized in terms of two form factors as, 3  [30,31]). Data obtained at different β's are rescaled by using the Sommer parameter r 0 /a, and the overall lattice spacing is fixed by matching f π obtained on the lattice with its physical value, leading to r 0 = 0.440 (12) fm (c.f. ref. [35]). All quark masses are given in lattice units.
The decay rate for the on-shell photon is [24] Γ To reach the physical form factor at q 2 = 0, with h c at rest, the twisted boundary condition applied on one of the charm quark propagators is made with (2.8)

Two-point correlation functions
Similarly to our recent publication [29], we use the gauge field configurations produced by ETM Collaboration [30,31] employing the maximally twisted mass QCD [32], the details of which are summarized in table 1. The masses of charmonia are extracted from the following correlation functions:

JHEP01(2013)028
in which the Dirac structures are chosen to provide the coupling to the charmonium states with quantum numbers J P C = 0 −+ , 1 −− , and 1 +− , for η c , J/ψ and h c , respectively. S c (0, 0; x, t) and S ′ c (0, 0; x, t) refer to the propagators of the charm quark in the doublet ψ(x) = [c(x) c ′ (x)] T entering the maximally twisted mass QCD action on the lattice [32] and therefore the propagator S c (0, 0; x, t) is obtained by inverting the above lattice Dirac operator with the Wilson parameter r, while S ′ c (0, 0; x, t) is obtained by using −r. In practice r = 1 and m cr is the same as the one used in the production of the gauge field configurations [30,31]. Finally, ∇ µ and ∇ * µ are the usual forward and backward derivatives on the lattice. In this study we also implement the Gaussian smearing on one of the sources [36]. In other words, one replaces c(x) by where the smearing operator H is defined via [37] H where U na i,µ is the n a times APE smeared link [38], defined in terms of (n a −1) times smeared link U We chose the parameters κ = 4, n g = 30, α = 0.5, n a = 20 , (3.6) which are kept fixed for all of our lattices. The value of the bare charm quark mass, µ c , at each of our lattices is given in table 1. It has been fixed according to the result of ref. [35] where it was shown that the charm quark computed from the comparison of the lattice results with the physical m ηc fully agrees with the value obtained by using the physical m Ds or m D . Therefore, we can say that m ηc , obtained by numerically solving and then by fitting m eff ηc (t) at large time separations to a constant, is merely a verification that, after a smooth continuum extrapolation, we indeed reproduce m exp. To extract the values of m J/ψ and m hc we proceed along the same line and compute . (3.8) In figure 1 we show all three effective mass plots as obtained by using all four lattice spacings explored in this work and for one value of the sea quark mass which we choose to be the least light ones. We see that the effective masses for the pseudoscalar ad vector charmonia are excellent while the signal for the orbitally excited state, h c , is much more noisy. The effective masses are then combined to . (3.9) The advantage of these ratios is that they have smaller statistical errors than any of the effective meson masses separately. In figure 2 we show one such a ratio. On the plateaus, that we carefully examined for each of our lattices, we then fit R J/ψ,hc (t) to a constant R J/ψ,hc . In table 2 we collect our results for R J/ψ,hc as obtained from all of the lattice ensembles at our disposal. As indicated in the plot in figure 1 the fitting intervals involving the state h c are shorter. The shift of the plateau region to the right is made to account for the different lattice spacings, so that the fit is made at approximately the same physical separation between the interpolating field operators. After a smooth linear extrapolation to the continuum limit, (3.10) we obtain For the reader's convenience we also quoted the values obtained from experiments [10]. In eq. m q ≡ m MS q (2 GeV), while the parameter c J/ψ,hc ≈ 3 % measures the leading discretization effects. Division by a β=3.9 = 0.086 fm is made for convenience. The linear fit (3.10) describes our data very well except for the results obtained at β = 3.8. The results obtained at β = 3.8 can be either excluded from the continuum extrapolation, which is how we got the above results, or a term proportional to a 4 can be added in (3.10) which leads to a result that is fully consistent with the one quoted above. Finally, we should stress that the disconnected, OZI-suppressed, contributions to the correlation functions discussed in this work have been neglected. The fact that our lattice results agree with the experimental values (3.11) can be viewed as a proof that the OZI suppressed contributions to the two-point functions are indeed very small.

Dispersion relation and the sea quark mass dependence
We already mentioned that for the determination of the desired radiative decay form factors one of the charm quark propagators is to be computed with twisted boundary conditions (2.5), with the twisting angle tuned to ensure q 2 = 0, c.f. eqs. (2.3), (2.8). In this section we check on the energy-momentum relation in the case of the pseudoscalar charmonium η c by exploring five different values of the twisting angle, covering the range of the meson's three-momenta, 0 ≤ | p| ≤ √ 2 × 2π/L. We compute with p = θ/L. For | θ| = 0 we obviously get m ηc . For | θ| = 0, we proceed like in eq. (3.7) and fit E eff p (t) to a constant E p . As expected, due to the discretization effects the continuum relativistic formula, E 2 p = m 2 ηc + p 2 , is not accurately verified on the lattice. Instead, our (3.14) The results at β = 4.05 are shown in figure 3. At this stage it is not clear whether or not with higher statistics all of the above C β values would get closer to 1. The finite volume effects on the sea quark mass are unlikely to modify the expression (3.13). That point we could check from our simulations at β = 3.9 and µ sea = 0.0040 where the results obtained at two lattices 24 3 × 48 and 32 3 × 64 are perfectly consistent.
Since we are working with heavy quarks it is tempting to use the non-relativistic energy-momentum relation as well, but accounting for the lattice artifacts that according to ref. [39] can be included by the distinction between the "rest" (m (0) ηc ) and the "kinetic" mass (m (1) ηc ) in, We fit our data to this expression and find that at each lattice spacing the kinetic mass is indeed larger than the rest one but they ultimately converge to the same value in the continuum limit, as they should on the basis of the restored Lorentz invariance, as shown in figure 4. Interestingly, however, we see that the discretization error in both the rest and kinetic masses are of the same size, but of the opposite sign. The above discussion on the dispersion relations is made on the data obtained at fixed value of the sea quark mass. We checked that indeed The same is, however, not true with m J/ψ , although this is hard to see from the mass ratio in eq. (3.10) alone. Instead we examined the hyperfine splitting, and from the fit of our data to,    in good agreement with the experimental result written in brackets [10], and in excellent agreement with the result of lattice QCD computation presented in ref. [40], ∆ = (111 ± 5) MeV. Once more, this implicitly suggests that the contribution of the OZI breaking JHEP01(2013)028 diagrams in the charmonia, neglected in our computations, are very small (c.f. discussion in ref. [41][42][43]). Note also that from the fit of our data to eq. (3.18) we find, b ∆ = 1.0(3) GeV −1 , c ∆ = 0.47 (6) , (3.20) in qualitative agreement with ref. [13] where a tiny decrease of ∆ is found while lowering the sea quark mass. Note, however, that this observation (b ∆ 0) disagrees with earlier findings of ref. [44]. To better appreciate that disagreement we fit our data at each lattice spacing separately as a linear function of the light sea quark mass and plot the resulting curves in figure 5. We see that at each of our β's the hyperfine splitting mildly decrease as the sea quark mass approaches the chiral limit. The above value for ∆ cont. is obtained without including the results at β = 3.8 in the continuum extrapolation. If they are included our final result does not change but the error becomes smaller by 1 MeV. To be consistent with what we quote as a result of the continuum extrapolation in all other quantities computed in this work, we will quote ∆ = (112 ± 4) MeV.

Radiative transition form factors
To extract the desired hadronic matrix element (2.1) we computed the three point correlation functions where P =cγ 5 c ′ , V i =cγ i c ′ are the interpolating operators fixed at t = 0 and t = t y = T /2 (T being the time extension of our lattices). Using the fact that our three-momentum is isotropic, q = −(1, 1, 1) × θ 0 /L, we can average where the last line is valid for the sufficiently separated operators in the correlation function (4.1), which is ensured by the Gaussian smearing procedure. Coupling to the smeared pseudoscalar (vector) interpolating field operator P (V ), is denoted by Z S P (Z S V ). 5 Electromagnetic current, J em j = Z V (g 2 0 )cγ j c, is local and renormalized by using Z V (g 2 0 ) listed in table 1. Notice also that the pseudoscalar source operator has been fixed at t y = T /2, which simplifies the averaging of the signals propagating in both halves of the lattice, 5 The values of Z S P and Z S V are easily computed from the correlators of smeared source operators as to which we include the signal propagating from the opposite end of our lattice. i.e. t ∈ (0, ±T /2). In computing the propagators S c (x; y) we used the stochastic source techniques as explained in ref. [30,31]. Since the signals for the pseudoscalar and vector charmonia are very good, we can proceed to eliminate the sources in two ways. We can either divide the correlator (4.1) by the corresponding two point functions, which we refer to as the numerical ratio, or by the analytic expression of the coupling of the source operators to the lowest states, which we call the semi-analytic ratio. Obviously, the values for Z S P,V , m J/ψ and E ηc are obtained from the study of the two-point correlation functions. The resulting plateaus for all the lattice spacings used in this work are shown in figure 6. We see that thanks to efficiency of the smearing procedure the plateaus of the two ratios indeed coincide over a JHEP01(2013)028  broad range of time-slices. In the plots shown in figure 6 we also indicate the intervals (shaded areas) in which we fitted the ratios R sa,num (t) to a constant. This is then used to obtain the form factor V (0), as indicated in eq. (4.2). In table 2 we present our results for the form factor V (0), as obtained from all of the lattice data-sets and by using the ratio R sa (t). These results should be now extrapolated to the physical limit (m sea ≡ m q → 0, a → 0) by (4.5) from which we then obtain the two following values: 1.941(35) (without β = 3.8) , depending on whether or not we include the results obtained on our coarsest lattices in the continuum extrapolation. A few more comments are in order: • The values of the J/ψ → η c γ form factor computed on our lattices do not depend on the light sea quark mass, b V ≃ 0.
• Although the smallness of discretization errors in tmQCD is guaranteed by construction [they are O(a 2 )], they are quite large for the form factor V (0). More specifically, from the fit of our data to eq. • The above results are obtained by using the semi-analytic ratio R sa . By repeating the same analysis with R num we verify the same features concerning the O(a 2 ) effects and after extrapolating the results computed on all our lattices we obtain V (0) = 1.903 (48), entirely consistent with the results quoted in eq. (4.6).
Our final result is: where the second error is a difference between the central values obtained by using the semi-analytic and numerical ratios discussed above.
We now turn to the discussion of the form factor F 1 (0), relevant to the h c → η c γ decay, as defined in eq. (2.6). To that end we compute the three point correlation functions Keeping in mind that we consider h c at rest and for q = −(1, 1, 1) ×θ 0 /L, withθ 0 tuned as JHEP01(2013)028 indicated in eq. (2.8), we can easily isolate the term proportional to F 1 (0) by combining: +C 232 ( q; t) + C 313 ( q; t) + C 121 ( q; t)] , where, as before, in the last line we show the result of the spectral decomposition when all operators are sufficiently separated. By fitting to a constant we obtain m hc F 1 (0). The results for the form factor F 1 (0), as obtained from all of our lattice data sets, are presented in table 2. These results too need to be extrapolated to the continuum and chiral limits, in a way similar to eq. (4.5).
• The error on F 1 (0) computed on each of our lattices is under 10%, but is nevertheless twice larger than those we have when computing V (0). This is expected as the signal for h c is much harder to tame. For the same reason we could not compute the form factor by employing the numerical method, i.e. by dividing by the two-point correlation functions.
• Contrary to the case of V (0), the discretization effects on F 1 (0) are small (c.f. figure 7). From the fit of our data to the form similar to eq. (4.5) we find c F 1 ≈ 2%. We get: • Like in the case of V (0), within the accuracy of our data, the form factor F 1 (0) is insensitive to the variation of the light sea quark mass.
• Our final result is where the second error is our estimate of the uncertainty due to the method for extracting the form factor from the correlation functions. In the case of V (0) we were able to estimate that error from the difference between the central values obtained by using the semi-analytical (4.4) and numerical ratios (4.3). That difference turned out to be less than 2%. We add the same error to F 1 (0) leaving its sign free.

J/ψ annihilation constant
In addition to the above quantities we also computed the annihilation constant f J/ψ defined as, which enters decisively in the expression for the well measured electronic width Γ(J/ψ → e + e − ). Above, e λ µ stands for the polarization vector of J/ψ. f J/ψ is computed along the same lines discussed in our previous paper [29], namely from the fit of the correlation function C J/ψ ii (t) from eq. (3.1) to the form 14) where the non-perturbatively determined Z A (g 2 0 ) are those listed in table 1. 6 In practice we of course combine the correlation functions in which both interpolating operators are smeared and those in which one operator is local and the other one is smeared. The results are given in table 2, which after extrapolating to the continuum by using a form analogous to eq. (4.5) lead to The continuum extrapolation is smooth (c f = 5(2)%) and the result is obtained by using all our lattices. If we leave out the results obtained at β = 3.8, the resulting f J/ψ is larger by 9 MeV, which is the second error quoted above. Note also that f J/ψ depends very mildly on the light sea quark mass (b f = 0.5(2) GeV −1 ).

JHEP01(2013)028
where for the physical result we used the measured Br(J/ψ → η c γ) = (1.7 ± 0.4)% and the full width Γ J/ψ = 92.9 ± 2.8 keV [10], as well as the physical values of m J/ψ = 3096.92(1) MeV, and ∆ = 116.6 ± 1.2 MeV. Had we used our ∆ = 112 ± 3 MeV, the resulting Γ(J/ψ → η c γ) would have been 10% smaller. Although somewhat lower than the quark model result, Γ(J/ψ → η c γ) = 2.85 keV [3][4][5][6], our lattice result obviously gives a larger J/ψ → η c γ decay rate, and the agreement with experiment is only at 2σ. The effective theory approach of ref. [47] and the QCD sum rule analyses [48,49] succeeded at getting lower value for the decay rate of this decay, but with large uncertainties. Note that the dispersive (model independent) approach of ref. [50] predicted, Γ(J/ψ → η c γ) = 2.2 ÷ 3.2 keV, many years ago. All these results agree with ours too, except that we have smaller and controlled uncertainties. We hope more effort on the experimental side will be devoted to clarify the disagreement among various experiments, including the recent ones. For example, a study of this decay at BESIII would give us a very valuable information. Recent result at KEDR suggested a larger value for the branching fraction Br(J/ψ → η c γ) = (2.34 ± 0.15 ± 0.40) % [51], which would result in Γ(J/ψ → η c γ) = (2.2 ± 0.6) keV, in very good agreement with our result (5.2).
As for the other lattice calculations of this form factor, we note that the quenched result of ref. [24], V (0) = 1.85(4), is only slightly lower than ours, while the one obtained at single lattice spacing with N f = 2 light flavors in ref. [25], V (0) = 2.01(2), is larger than our values at β = 4.05 listed in table 2. Apart from different methodology, a notable difference is that the authors of ref. [25] used the point-split electromagnetic current that does not require renormalization, whereas we use the local current on the lattice that is properly renormalized by non-perturbatively determined Z V (g 2 0 ). Keep in mind that our final results is obtained after the extrapolation to the physical limit of the results obtained at several lattice spacings and for several different values of the light sea quark mass.

h c → η c γ
h c escaped the experimental detection for a long time and only recently CLEO succeeded to isolate this state [52] and observed that its prominent mode is precisely h c → η c γ, the branching fraction of which was later accurately measured at the BESIII experiment, with a result: Br(h c → η c γ) = (53 ± 7)% [53]. We obviously cannot compute the branching ratio on the lattice, but with our form factor result (4.12) we can compute the decay width using eq. (2.7). We get This result can be combined with the measured Br(h c → η c γ) to estimate the width of the h c state. We obtain: where the first error comes from our determination of the form factor F 1 (0), and the second one reflects the experimental uncertainty in the branching ratio. Notice also that JHEP01(2013)028 we symmetrized the error bars. This constitutes a prediction that would be interesting to check against the actual experimental measurement once the latter becomes available. 7 To compare our estimate F 1 (0) = −0.57(2) with other lattice results we convert the value of the value reported in ref. [24] to our dimensionless form factor and obtain F 1 (0) = −0.53 (3), which agrees very well with our result. 8 Similar conversion of the result of ref. [25] would result in F 1 (0) = −0.33(1), much smaller value than ours, whether we compare it with the values we obtain at β = 4.05 or the one in the continuum limit.

Summary
In this paper we presented results of our analysis of the radiative decays of charmonia by means of QCD simulations on the lattice. Using the twisted mass QCD with N f = 2 dynamical flavors at several small lattice spacings we were able to smoothly extrapolate the relevant form factors to the continuum limit.
• We used the twisted boundary conditions to make sure that we extract the physical form factors, i.e. at q 2 = 0. We checked that at every lattice spacing explored in this paper our data indeed reproduce the latticized energy-momentum relation given in eq. (3.13). We also showed that the so-called kinetic and rest masses converge to the same value in the continuum limit, but that at fixed lattice spacing they both have discretization errors that are equal in size but different in sign.
• We computed the hyperfine splitting and obtained ∆ = m J/ψ − m ηc = 112 ± 4 MeV , (6.1) and showed that it depends very mildly on the sea quark mass, with a slope being positive. From our computation of R hc in eq. (3.11) we obtain m hc = 3.537(32) GeV, in good agreement with m exp. hc = 3.525 GeV.
• We computed the hadronic form factor relevant to the J/ψ → η c γ (M1) transition and found V (0) = 1.92(3)(2) , (6.2) which is larger than the one we would infer from the measured Γ(J/ψ → η c γ), although compatible at the 2σ level. We found that the discretization effects are large and negative, but that they are adequately described by the linear function in a 2 .
• Our result for the h c → η c γ (E1) transition form factor is

JHEP01(2013)028
which mildly depends on the lattice spacing. Within the uncertainties quoted above, both our form factors are insensitive to the change of the sea quark mass. After combining our F 1 (0) with the measured Br(h c → η c γ) we deduced the value of the width, Γ hc = 1.37 (22) MeV.
In the above results we do not make any estimate of the size of systematic uncertainty due to the omitted s and c quarks in the sea. In ref. [55] it was claimed that the contributions from dynamical charm might be important. That point will be numerically assessed from the analysis similar to the one presented in this paper but on the set of gauge field configurations that include N f = 2 + 1 + 1 dynamical quark flavors. We also emphasize that our results are obtained without inclusion of the OZI suppressed contributions. Their impact appears to be small in J/ψ → e + e − decay, but their size in the radiative decays is unknown. They were neglected and the associated uncertainty cannot be estimated without actually attempting to compute the corresponding disconnected diagrams on the lattice. Such a computation would be very welcome. The strategy employed in this paper can be applied to compute the much needed F h b →η b 1 (0) which we plan to do in the near future.