Lattice QCD estimate of the $\eta_{c}(2S)\to J/\psi\gamma$ decay rate

We compute the hadronic matrix element relevant to the physical radiative decay $\eta_{c}(2S)\to J/\psi\gamma$ by means of lattice QCD. We use the (maximally) 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. The value of the mass ratio $m_{\eta_c(2S)}/m_{\eta_c(1S)}$ we obtain is consistent with the experimental value, and our prediction for the form factor is $V^{\eta_{c}(2S)\to J/\psi\gamma}(0)\equiv V_{12}(0)=0.32(6)(2)$, leading to $\Gamma(\eta_c (2S) \to J/\psi\gamma) = (15.7\pm 5.7)$ keV, which is much larger than $\Gamma(\psi (2S) \to \eta_c\gamma)$ and within reach of modern experiments.


Introduction
Recent progress in simulations of QCD on the lattice allowed to solve several long standing problems in hadronic physics. One such a problem was a conflict between theoretical predictions of the radiative decay width Γ(J/ψ → η c γ) and its value experimentally measured in 1986 [1]. That early measurement turned out to be too small when confronted with theoretical predictions based on various quark models [2,3,4], dispersion relations [5] and QCD sum rules [6,7]. Only in 2009 the CLEO Collaboration [8] was able to provide a new measurement of this decay width and found it to be over 1σ larger than the one measured using the Crystal Ball detector in 1986 [1], but still somewhat smaller than predicted. A very recent measurement at the KEDR experiment confirmed the CLEO result in that Γ(J/ψ → η c γ) is large [9], and reported a value 1.4σ larger than that by CLEO. The charm factory at BESIII is expected to provide a new experimental determination of this decay width and close this issue.
On the theory side, as we just mentioned, the numerical simulations of QCD on the lattice helped solving this problem since the corresponding form factor was computed at several lattice spacings, thus allowing to take the continuum limit and obtain a viable physical result. The effects of light dynamical quarks were included in simulations. In particular, by using the maximally twisted mass QCD on the lattice with N f = 2 dynamical light quarks we confirmed in ref. [10] the discrepancy between theory and the old measurement [1]. Our finding was soon corroborated by a completely independent simulations made by the HPQCD Collaboration, in which the effects of N f = 2 + 1 staggered light quark flavors were included [11]. Both lattice results implemented the non-perturbative renormalization procedure and, in the continuum limit, exhibited quite a remarkable agreement for Γ(J/ψ → η c γ), also in agreement with the most recent experimental findings [9]. Furthermore, recent improvement of the effective theory approach, based on potential non-relativistic QCD (pNRQCD) [12], lead to a very good agreement with lattice QCD results [13]. Therefore, as of today, the theoretical estimate of Γ(J/ψ → η c γ) is very solid. Another important motivation for a more dedicated experimental study of this decay mode lies in the fact that this decay rate could be sensitive to the CP-odd light Higgs boson if its mass were very light, i.e. close to that of the η c -meson [14].
In this paper we discuss another class of so-called magnetic dipole (M1) transitions, namely a decay of a radially excited charmonium to a ground state. In particular, we will focus on η c (2S) → J/ψγ, the width of which is yet to be measured and the lattice QCD computation of its hadronic matrix element should be regarded as a clear theoretical prediction. Quark model predictions of the hadronic matrix element governing this decay are difficult to control because the leading term vanish due to orthogonality of the meson wave functions and the relativistic corrections are highly sensitive to the form of the used potential. Phenomenology of this decay mode has not been considered in the effective theory approach [12]. It is often assumed to be similar in size to Γ(ψ(2S) → η c γ) which has been measured and found to be small. Even if the physical form factors V ηc(2S)→J/ψγ (0) and V ψ(2S)→ηcγ (0) were equal, the width of η c (2S) → J/ψγ would still be about three times larger due to different spin of the initial state. However, there is a dynamical reason why Γ(ψ(2S) → η c γ) is suppressed with respect to Γ(η c (2S) → J/ψγ). It was first observed in ref. [2] that a substantial part of relativistic corrections cancel in ψ(2S) → η c γ, whereas they add up in the case of η c (2S) → J/ψγ to further enhance its decay rate. In the effective field theory approach [12], that point has been recently emphasized in ref. [13], and can be compactly written as: where , m c is the pole charm quark mass, and the matrix elements of O(1/m 2 c ) corrections are for shortness written in the form A|O|B ≡ A O B . The two decays differ in sign of the terms involving the spin dependent potential, which is a peculiarity of these, so-called hindered, M1 transitions, in contrast to the allowed ones (e.g. ψ(nS) → η c (nS)γ) for which the spin dependent corrections do not occur in a first few terms of the 1/m c -expansion. Assuming that J/ψ V S 2 ( r) ηc(2S) ηc V S 2 ( r) ψ(2S) and is positive, the second decay width will be enhanced with respect to the first one. To compute the hadronic matrix elements of the spin dependent part of the potential one should either rely on quark models or attempt computing them by means of NRQCD on the lattice. In pNRQCD, in situations in which the charmonium states with principal quantum number larger than 1 are involved, the computation of the matrix elements in eq. (1) cannot be handled analytically and should be computed by a non-perturbative method.
In this paper we will compute the hadronic matrix element relevant to η c (2S) → J/ψγ by using QCD on the lattice, without relying on NRQCD, and show that the value of the corresponding form factor is indeed significantly larger than the one governing the ψ(2S) → η c γ decay. This is the first time that such a computation is conducted and on the basis of our result we obtain The remainder of this paper is organized as follows: in sec. 2 we define the matrix element and the corresponding form factor, and discuss the strategy to ensure q 2 = 0; in sec. 3 we discuss the two-point correlation functions and the method used to isolate the radially excited state the efficiency of which we test on the mass splitting between the radially excited and the lowest lying states; in sec. 4 we describe the computation of the threepoint correlation functions, extract a desired form factor and discuss the phenomenological consequences of our result; we conclude in sec. 5

Hadronic Matrix Element
The hadronic matrix element governing the radiative decay η c (2S) → J/ψγ * decay can be parameterized in terms of the form factor V 12 (q 2 ) as, 2 where the relevant part of the electromagnetic current is J em . The form factor V 12 (q 2 ) can be computed at various values of q 2 ≡ q 2 γ = (p − k) 2 but the one relevant to the physical η c (2S) → J/ψγ decay rate (on-shell photon) should be obtained at q 2 = 0, viz.
In our previous paper we showed that the charmonium decays J/ψ → + − and J/ψ → η c γ do not depend on the sea quark mass. In this paper we again work with the charmonium states that are all bellow the D ( * )D( * ) production threshold and therefore the dependence on the light (sea) quark should remain negligible. 3 For that reason, in this study, we focus on a subset of gauge field configurations considered in ref. [10] and study one value of the light sea quark mass per lattice spacing but we increase the statistics in order to be able to isolate the radially excited state from our correlation functions. We rely on the gauge field configurations produced by the ETM Collaboration [18] in which the effect of N f = 2 mass-degenerate dynamical light quarks has been included by using the maximally twisted QCD on the lattice [17]. Using the same action we then compute the quark propagators and correlation functions needed for the physical problem discussed in this paper. Details concerning the lattice ensembles and the main results of this paper are listed in tab. 1.   Table 1: Summary of the lattice ensembles used in this work (more information can be found in ref. [18]).
Lattice spacings and bare charm quark masses have been determined in ref. [20]. µ sea and µ c are the bare quark masses and are given in lattice units. n g is the smearing parameter, cf. eq. (15). Values of the mass ratios and the form factors obtained at each lattice spacing are also given. q 2 0 is specified in eq. (29). t sep is the separation between the source operators chosen in computation of the three-point correlation functions. # meas. is written in terms of a number of independent gauge field configurations × a number of time sources used to compute propagators.

Two-point correlation functions
A crucial step in extraction of the matrix element between the lowest lying vector charmonium and a radially excited pseudoscalar one is to reliably project out the radially excited state. That is made through a careful study of two-point correlation functions which will be discussed in this section.
To compute the mass of η c (2S) we use a set of interpolating field operators, P 1 , . . . , P N , each coupling to a tower ofcc-states with J P C = 0 −+ , and build a N × N matrix of twopoint correlation functions: Spectral decomposition of each correlation function can be written as where the sum runs over η c (nS) states. In the above decomposition we neglected the multi-particle states which is legitimate since we consider a few lowest lying states, below the D ( * ) D ( * ) production threshold. Z j (nS) in eq. (8) denotes the hadronic matrix element, Z j (nS) ≡ 0|P j |η c (nS) . We wish to find a linear combination of operators P i that couples optimally to a state n, viz.
can be obtained by solving the Generalized Eigenvalue Problem (GEVP) [21], and from the resulting eigenvectors v our desired solution. The parameter t 0 in eq. (9) should be chosen large enough so that the correlation functions C ij (t 0 ) are dominated by the lightest n states, η c (nS). The role of C (t 0 ) is to optimize the problem and help us to better isolate the lowest n-states. The above GEVP can be solved for each time-slice t, and the corresponding eigenvectors v (n) (t) are expected to be independent of t when focusing onto the lowest n states. In practice v (n) (t) are independent on t, up to the effects of statistical noise which can be reduced.
In this paper we define c i P i that couples optimally to the state η c (nS). The corresponding correlation function, at larger values of t, is dominated by η c (nS), the mass of which is then extracted from In the present study we use a basis of three operators, 4 where the symmetric covariant derivative on the lattice is defined as: Note that in eq. (13) we use c = Hc, to distinguish the smeared quark field from the local one, with the smearing operator H given by [22]: The links U na n;i entering the smearing operator and the covariant derivative (14) are n a -times APE smeared [23], i.e. they are obtained from the (n a − 1)-times smeared link U (na−1) n;i and its surrounding staples, denoted by V In the present study we use α = 0.5 and n a = 20 for all our lattices. The other smearing parameters that appear in H are n g and κ. We keep κ = 4 for all our lattices, but the number of steps n g as given in tab. 1.
With the above set of operators (13) we focus on the first three states, our main target being the radial excitation η (2S). We checked that the inclusion of more operators, in a way it has been done in e.g. ref. [24], does not lead to any improvement of the signal for the state η (2S). With the statistical quality of our data more operators included in GEVP would not help improving the extraction of higher excited states either. The values of coefficients c (n) 1,2,3 we find for the first two states, as well as the values we take for t 0 and t opt , are given in tab. 2.
To determine simultaneously the mass m ηc(nS) and the matrix element Z (n) (nS) we fit the correlation functions C ηc(nS),ηc(nS) (t), defined in eq. (11), over an appropriate time interval using eq. (12). In fig. 1 we show the effective mass m eff ηc(nS) (t) defined as: 5 m eff ηc(nS) (t) = log C ηc(nS),ηc(nS) (t) C ηc(nS),ηc(nS) (t + 1) , together with the results of the fits to a constant for three lowest lying states, and for all four lattice spacings discussed in this work. Note that at β = 4.2 the statistical quality  i determined by solving the GEVP in eq. (9) in the basis of operators listed in eq. (13). n = 1 refers to the optimal coupling to the lowest lying state, and n = 2 to its first radial excitation. We also give the values of t 0 and t opt that we chose while solving the GEVP [cf. text following eqs. (9,10)].
of our data did not allow us to distinguish the second radial excitation. The mass of the lowest lying state is improved with respect to the results presented in our previous papers [10,25], but the overall error bars remain the same since it is entirely dominated by the error in lattice spacing. Instead of looking for absolute values of the meson masses, we prefer to compute the ratio of the radial excitation with respect to the ground state, thus eliminating the error on lattice spacing from the discussion. In fig. 1 we show the plateaux for the first two states that are pronounced and of good quality. After fitting each effective mass to a constant we were able to extract m ηc(2S) /m ηc(1S) , in an obvious notation η c ≡ η c . The results are reported in tab. 1.
Strictly speaking m ηc is not a lattice result. It is just a cross-check because the mass of the charm quark (µ c in tab. 1) has been tuned in ref. [20] in such a way as to reproduce the correct m exp ηc = 2980.3 MeV, and was then checked to result in a correct physical m D (s) in the continuum limit. The results for m ηc(2S) /m ηc , instead, are clean lattice QCD results. To get a physically relevant result we need to make the continuum extrapolation, which we do by using where we account for the dominant O(a 2 ) discretization effects [17]. The above form appears to be adequate to describe our data, and the result of that extrapolation is shown in fig. 2. We obtain R cont 2 = 1.230 (18),  in very good agreement with the experimentally established R exp 2 = 1.220(1) [26]. The parameter X R measures the shift of the continuum value with respect to the one obtained at the lattice with a = 0.085(3) fm, which appears to be in the range of 3 ÷ 5%. The above quoted errors are statistical only. By modifying (enlarging) the plateau region and including 2 more points, we end up with the fully compatible results, which then in the continuum limit give R cont 2 = 1.226 (18). In view of the fact that our lattice QCD result has a much larger error than the corresponding physical result, we will not further dwell on systematics but simply conclude that the lattice results obtained by solving the GEVP are adequately described by eq. (19) and the result obtained in the continuum limit is fully compatible with the physical R exp 2 = 1.220(1). This is not the first lattice determination of R 2 but it is the first in which the maximally twisted mass QCD on the lattice has been used for its computation. In ref. [ (1) MeV, despite the fact that they worked at one lattice spacing only [24]. Another important observation comes from the comparison of our result, R cont 2 = 1.230 (18), with the physical R exp 2 = 1.220 (1). In that respect the inclusion of non-local operators P 1 and P 3 in the set of operators used to solve the GEVP in eq. (13) is crucial. In a preliminary study we used only a set of P 2 operators that differ between each other by a choice of smearing parameters. Despite the fact that we optimized the smearing parameters in a way that the coupling to η c (2S) is larger/smaller, the resulting splitting between η c (2S) and η c (1S), in the continuum limit, was much larger (by about 250 MeV) than the physical one. That observation depends on the physical quantity we consider. For example, the results for the form factor V 12 (q 2 ) are more robust and remain fully compatible in both situations: (i) in which we use the operator basis (13), (ii) when the basis consists of P 2 operators only, differing from each other by the amount of smearing implemented.

Transition Form Factor for η c (2S) → J/ψγ
We now turn to the extraction of the form factor relevant to η c (2S) → J/ψγ. We first compute the three-point correlation function, where P (2) is the operator optimally interpolating the η c (2S) state obtained in the previous section, V i =cγ i c is the smeared operator interpolating the J/ψ, and J em j =cγ j c is the local vector current, renormalized by using Z V (g 2 0 ) given in tab. 1. We benefit from the time reversal symmetry that relates the photon emission and the photon absorption processes which in terms of our correlation functions means, Since the computation of correlation function for the latter process requires less propagator inversions, in the following we will discuss the right hand side of eq. (22). The corresponding Wick contraction reads, where, instead of explicitly injecting the momentum q to the correlation function, we use the twisted boundary condition on one of the charm quark propagators labelled by the superscript 'θ'. In our computation of the quark propagator S c (x; y) we use the stochastic source technique described in ref. [18]. In the present study we neglect the disconnected contractions arising in the computation of all correlation functions, which is equivalent to studying these processes in a theory that contains a doublet of charm quarks and we focussed on the non-singlet states. Such an approximation is expected to have a small impact on physical observables which is what we observed in our previous paper [10] (see also ref. [11]). Following the discussion made in sec. 3, the explicit expressions of the interpolating field operators relevant to η c (2S) and J/ψ states, at a given time-slice t, are: where the coefficients c to η c (2S) |J em µ |J/ψ , i.e.
where η c (2S) is at rest, and the couplings Z are given by: In practice, t sep = T /2 that was suitable for the study of J/ψ → η c γ, may be too large for extraction of the matrix element involving a radially excited state. To increase the region in which we can extract the matrix element η c (2S) |J em j |J/ψ( q, i ) we choose t sep < T /2. The values of t sep are also given in tab. 1. Furthermore, since we take our three-momentum to be isotropic, q = (1, 1, 1) × ϑ 0 π/L, we can average over six non-zero contributions, namely, The matrix element is then obtained after dividing the source operators from the threepoint function (21), namely, where the values of Z V,P (2) , E J/ψ and m ηc are obtained from the study of two-point correlation functions. In fig. 3 we illustrate the plateau of R 3 (t) which is then fit to a constant (shaded area in fig. 3) that corresponds to the matrix element, J/ψ|J em µ |η c (2S) , from which we then get the form factor V 12 (0), c.f. eq. (3). Furthermore, by replacing the coefficients c 1,2,3 in eq. (24) and from the same correlation function (23), we also get the matrix element relevant to J/ψ → η c γ * decay, i.e. the form factor V 11 (q 2 0 ). Since the value of ϑ 0 has been chosen to ensure that the emitted photon in η c (2S) → J/ψγ is on shell (q 2 = 0), after a trivial algebra one gets which is convenient because m ηc is in our study used as input to fix the value of the bare charm quark mass, and therefore m ηc computed on each of our lattices is equal to the physical m exp ηc = 2.981(1) GeV by construction. Instead, R J/ψ = m J/ψ /m ηc , discussed in our previous paper [10], and R 2 = m ηc(2S) /m ηc , computed in this work, vary with lattice spacing, and therefore q 2 0 also changes from one lattice spacing to another. We checked that by using V 11 (0) = V 11 (q 2 0 ) exp[|q 2 0 |/(16b 2 )], with b = 0.54(1) GeV or b = 0.58(2) GeV, as found in refs. [28] and [29] respectively, we reproduce our results for V 11 (0) presented in ref. [10].
Finally, we need to extrapolate our results for V 12 (0) to the continuum limit. To that end we use the expression similar to eq. (19) and fit our data to That extrapolation is shown in fig. 4 and in the continuum limit we obtain V 12 (0) = 0.32(6) , We see that the final error on V 12 (0) is rather large and the small effects due to fixing the charm quark mass and of the overall lattice spacing are completely immaterial at this stage. To account for a more important source of systematic uncertainty, we performed the continuum extrapolation by removing either the finest or the coarsest lattice and obtained V 12 (0) = 0.31 (8), and V 12 (0) = 0.34(8), respectively. We can then take the spread of central values as an estimate of the error due to extrapolation to the continuum limit. To evaluate the impact of of higher excited states to our matrix element extraction, we also shortened the fitting region of R 3 (t), but keeping the points on the left of the plateaus (see fig. 3) which are more likely to be sensitive to the higher excited states, and after the continuum extrapolation we obtain V 12 (0) = 0.30 (7). As our final result we quote

Phenomenological discussion
Let us first remind the reader of the value of the form factor V 21 (0), that parameterizes the hadronic matrix element describing ψ(2S) → η c γ in a way completely analogous to eq. (3). The decay branching fraction is given by, which can be combined with the experimental values [26], and α(m c ) = 1/134. [30], to get much smaller than V 12 (0) = 0.32 (6), that we obtained after extrapolating our lattice QCD results to the continuum limit. We attempted computing the form factor V 21 (0) at single lattice spacing and found it to be very small, consistent with zero within our error bars. It would take a huge statistics of the lattice data sample to be able to compute V 21 (0) comparable with experimental accuracy. The value of V 12 (0), instead, was found to be large at every lattice spacing. Knowing that Γ η c (2S) = (11.4 ± 3.1) MeV, much larger than Γ ψ(2S) , the branching fractions of two decay modes become similar in size. More specifically, with we get very similar to the measured B(ψ(2S) → η c γ), and should be within reach at BESIII, KEDR, LHCb or Belle-2.
As we mentioned in introduction, the fact that V 12 (0) is much larger than V 21 (0) might come as a surprise because they differ by 1/m 2 c -corrections and higher, that are expected to be reasonably small. However, from the quark model picture we know that a dominant contribution to these (hindered) transitions are absent due to orthogonality of the wave functions of initial and final states, and therefore the decay rates are almost entirely determined by the size of power corrections. In the effective field theoretical treatment of this problem, within pNRQCD, the terms O(1/m 2 c ) have been identified [12,13]. Unlike the allowed M1 transitions, such as J/ψ → η c γ and ψ(2S) → η c (2S)γ, the hindered processes depend on the spin-spin interaction of the heavy quark potential which is precisely the one that affects differently η c (2S) → J/ψγ and ψ(2S) → η c γ, cf. eq. (1). Assuming that all three operators (O), electromagnetic radius (r 2 ), typical velocity (p 2 /m 2 c ) and the spin operator (V S 2 ( r)) in eq. (1), satisfy J/ψ O ηc(2S) = ηc O ψ(2S) , from the difference of the two amplitudes we can estimate J/ψ V S 2 ( r) ηc(2S) ≡ V S 2 ( r) . More specifically, which can be useful for phenomenology based on pNRQCD as a first estimate of the corresponding decay rates in the case of bottomia.

Summary
In this paper we made the first lattice computation of the form factor needed for a theoretical estimate of the radiative decay η c (2S) → J/ψγ, and showed that its value is larger than the one entering the similar ψ(2S) → η c γ decay. The explanation of that phenomenon can be understood in the pNRQCD description of these processes because they both involve a spin-spin interaction term which in the former decay enhances the decay rate and in the latter cancels the dominant power correction term.
To be able to extract the matrix element for η c (2S) → J/ψγ on the lattice we needed to solve the GEVP and construct an interpolating field operator that couples mostly to η c (2S). We checked that in the continuum limit our result for m ηc(2S) /m ηc = 1.23 (2), is fully consistent with the measured (m ηc(2S) /m ηc ) exp. = 1.22. Our computations are made by using the (maximally) twisted mass QCD on the lattice by including N f = 2 dynamical light quarks and at four different lattice spacings. After taking the continuum limit and assuming the dependence on the light (sea) quark mass to be negligible, which we showed in our previous work to be the case for similar charmonium decays [10], we obtain that the decay width is Γ(η c (2S) → J/ψγ) = (15.7 ± 5.7) keV .
The available experimental information on B(ψ(2S) → η c γ) and B(ψ(2S) → η c (2S)γ) together with our lattice results for the form factors describing J/ψ → η c γ [10] and η c (2S) → J/ψγ, confirm the expected pattern: the form factor is indeed large in the case of the allowed M1 transitions, while it is small for the hindered ones. Lattice QCD helped to solve the nonperturbative QCD effects where the experimental information was poor or not available, and we now have: The value for the form factor V 12 (0) presented here can be improved by increasing statistics and by verifying that the form factor does not depend on the mass of the light sea quark. Furthermore a computation of V 12 (0) by using a different lattice QCD discretization scheme would be very welcome as well.