Revisiting the radiative decays $J/\psi \rightarrow \gamma\eta^{(\prime)}$ in perturbative QCD

In the framework of perturbative QCD, the radiative decays $J/\psi\rightarrow\gamma\eta^{(\prime)}$ are revisited in detail, where the involved one-loop integrals are evaluated analytically with the light quark masses kept. We have found that the sum of loop integrals is insensitive to the light quark masses and the branching ratios $\mathcal{B}(J/\psi\rightarrow\gamma\eta^{(\prime)})$ barely depend on the shapes of $\eta^{(\prime)}$ distribution amplitudes. With the parameters of $\eta-\eta^{\prime}$ mixing extracted from low energy processes and $J/\psi\rightarrow\gamma\eta^{(\prime)}$ by means of nonperturbative matrix elements $\langle0|G_{\mu\nu}^a\tilde{G}^{a,\mu\nu}|\eta^{(\prime)}\rangle$ based on $U_{A}(1)$ anomaly dominance argument, we could not give the ratio $R_{J/\psi}$ in agreement with experimental result. However, using the parameters, especially the mixing angle $\phi=33.5^{\circ}\pm0.9^{\circ}$, extracted from $\gamma^{\ast}\gamma-\eta^{\prime}$ transition form factor measured at $q^{2}=112~\mathrm{GeV}^{2}$ by BaBar collaboration, we obtain $R_{J/\psi}=4.70$ in good agreement with $R_{J/\psi}^{exp}=4.65\pm0.21$. As a crossing check, with $\Gamma^{exp}(\eta^{(\prime)}\rightarrow\gamma\gamma)$ and our results for $J/\psi\rightarrow\gamma\eta^{(\prime)}$, we get $\phi=33.9^{\circ}\pm0.6^{\circ}$. The difference between the determinations of $\phi$ is briefly discussed.


Introduction
Since its discovery [1,2], the J/ψ meson has always been an active topic in particle physics [3][4][5]. The heavy quarkonium physics might be crucially important to improve our understanding of quantum chromodynamics (QCD), especially, the interplay of perturbative quantum chromodynamics (pQCD) with nonperturbative QCD. The Okubo-Zweig-Iizuka (OZI)-forbidden radiative decays of J/ψ → γη ( ) are expected to proceed predominantly via two virtual gluons which subsequently convert to η ( ) , with the photon emitted from the initial charm quarks. Such decays have been the subject of many experimental and theoretical studies [3][4][5]. Theoretically, these decays provide a clean environment to study the conversion of gluons into hadrons. In this respect, the radiative decays J/ψ → γη ( ) are of particular interest, since they are also closely related to the issues of η − η mixing, which are important ingredients for understanding many interesting phenomena related to the η and η mesons.
In the literature, the exclusive radiative decays J/ψ → γη ( ) have been studied in different approaches. The decay widths Γ(J/ψ → γη ( ) ) were calculated by Novikov et al. [6], with the assumption that these decays occur as a consequence of the U A (1) anomaly and are, therefore, controlled by the nonperturbative gluonic matrix elements 0|G a µνG a,µν |η ( ) , where G a µν is the gluon field strength tensor andG a,µν = 1 2 µναβ G a αβ its dual tensor. Then the ratio of the decay widths takes the form [6] R where the matrix elements 0|G a µνG a,µν |η ( ) can be calculated with the QCD sum rules and other approaches. For example, Chao [7,8] and Kuang et al. [9] derived the expressions of the matrix elements 0|G a µνG a,µν |η ( ) in the large-N c approach and QCD multipole expansion, respectively.
In the Feldmann-Kroll-Stech (FKS) scheme for η−η mixing system [10], the nonperturbative matrix elements 0|G a µνG a,µν |η ( ) can also be expressed as function of the phenomenological parameters f q , f s , φ, and the ratio R J/ψ reads where φ denotes the mixing angle of η − η system. Compared R J/ψ with its experimental value, the mixing angle is found to be φ = 39.3 • ± 1.0 • [10]. However, as it had been pointed out that the above equation was calculated in the approximation with the assumption of ground state dominance and neglection of continuum contributions to the dispersion relations [6,11].
It is also noticed that the matrix elements 0|G a µνG a,µν |η ( ) induced by the U A (1) anomaly are a higher twist effect(twist-4) [12]. The 0|G a µνG a,µν |η ( ) would give the main contributions to the radiative decays J/ψ → γη ( ) , only in the case of the leading twist contributions were strongly suppressed. Although the leading twist contributions from the gluonic content of η ( ) are suppressed by a factor of m 2 η ( ) /M 2 J/ψ [13], the assumption that the matrix elements 0|G a µνG a,µν |η ( ) dominate the radiative decays J/ψ → γη ( ) could be broken, because the leading twist contributions from the quark-antiquark content of η ( ) are not suppressed much at the energy scale of M J/ψ . The first pQCD investigation of these decays were carried out by Körner et al. [14] about thirty years ago. They took the annihilation of cc quarks to be a short-distance process described by pQCD, and nonperturbative dynamics of the bound states factorized to wave functions. It has been argued that pQCD asymptotic behaviors may be expected at the scale of M J/ψ [12,[15][16][17]. However, in this pioneer work [14], the nonrelativistic quark model with the weak-binding approximation, has been taken for both J/ψ and η ( ) . Whereafter, the nonrelativistic approximation for η ( ) wave functions in [14] was improved by Kühn [18] with light-cone expansion. However, in the calculation of the loop integrals, the approximation of m 2 η ( ) /M 2 J/ψ ≈ 0 was made [18]. In our calculation, we would keep m η ( ) and m u,d,s , which result in our final analytical expression of the loop function is much complicated than the ones in Refs. [14,18]. We notice that our results can reproduce the one in Ref. [18] in the limit of m 2 η ( ) /M 2 J/ψ → 0, and detail comparison is presented in the Appendix. Recently, several groups have revisited J/ψ → γη ( ) in the framework of pQCD [19][20][21][22][23]. In works [20,21], the light-cone distribution amplitudes (DAs) were adopted for η ( ) . However, in Ref. [21], the decay widths of J/ψ → γη ( ) were found to be very sensitive to the light quark masses of η ( ) involved in the loop integrals, which is needed to be made clear, since the sensitivities of the loop integrals to light quark masses usually point to the possible infra-red (IR) divergences.
Besides the aforementioned QCD approaches, the radiative decays J/ψ → γη ( ) have also been studied with phenomenological models, such as the approach with an effective lagrangian [24] and the approach considering the η c − η ( ) mixings [7,8,25]. Generally, predictions compatible with the experimental measurements could be obtained.
In this paper, we present a detail calculation of these decays in the framework of pQCD.
The bound-state property of J/ψ is parameterized by its Bethe-Salpeter (B-S) wave function, while the η and η are described by their light-cone DAs. Then the loop integrals are evaluated analytically with the light quark masses kept.
In our calculation, the loop function is found to be insensitive to the light quark masses, which differs from the results in Ref. [21]. Moreover, our results of the B(J/ψ → γη ( ) ) are also insensitive to the shapes of the light meson DAs. The theoretical uncertainties due to choices of different η ( ) DAs available in the literature are negligible in the prediction for the ratio R J/ψ , so that, the mixing angle of η − η mixing could be reliably extracted. In addition, the corrections from QED processes J/ψ → γ * → γη ( ) are also considered in our calculation.
The paper is organized as follows. In section 2, we present the formalism for calculating the decay amplitudes of J/ψ → γη ( ) . Numerical results are presented in section 3, and the final part is our summary. The analytical expressions for the dimensionless key function H 0 and discussions of its properties are presented in the Appendix. In the rest frame of J/ψ, the amplitude of J/ψ → γg * g * can be decomposed into hardscattering part and the B-S wave function of J/ψ [26] where χ(q c ) is the B-S wave function of J/ψ and √ 3 is the color factor. K and ε(K) stand for the momentum and the polarization vector of the J/ψ, k, k 1 , k 2 and (k), (k 1 ), (k 2 ) stand for momenta and polarization vectors of the photon and the gluons, respectively. The momenta of the c andc quarks are parameterized as where q c is the relative momentum between the c andc quarks. Since the B-S wave function χ(q c ) is sharply peaked when q c ∼ 0 for the heavy J/ψ in the nonrelativistic limit, one may neglect the q c -dependence in the hard-scattering amplitudeÔ(q c ) in the leading order approximation, then the amplitude can be simplified to where the Salpeter function is defined as For the vector meson J/ψ, the Salpeter function with leading order Dirac structures reads [27,28] where φ(q 2 c ) is a scalar function of q 2 c , and M represents the mass of J/ψ. With the help of the definition of the S-wave wave function evaluated at the origin [29] one can obtain [14] A And the hard-scattering amplitudeÔ(0) can be written aŝ where we have made the nonrelativistic approximation M ≈ 2m c .
The light-cone expansion of the matrix elements of the meson η ( ) over quark and antiquark fields reads [12] η where the high twist terms are omitted. With this definition, one can obtain the coupling of [30][31][32]: where u is the momentum fraction carried by the quark andū = 1 − u, m q is the mass of the with the asymptotic form of DA φ AS (u) = 6u(1 − u) and c q n (µ) the Gegenbauer moments. In Table 1, we list three models of the DAs discussed in Ref. [33]. Their shapes are shown in  Table 1: Gegenbauer coefficients of three sample models at the scale of µ 0 = 1 GeV. The decay amplitude of J/ψ → γη ( ) can be obtained by multiplying the above two cou-plings, inserting the gluon propagators and performing the loop integrations where the factor 1/2 takes into account that the two gluons have already been interchanged both in A αβµν and M µν . Using parity conservation, Lorentz invariance and gauge invariance, it can be proved that so there is only one independent helicity amplitude H q QCD [14] T With the help of the helicity projector [14] one can obtain the helicity amplitude The dimensionless function H q reads is the sum of the loop integrals of the six Feynman diagrams for the decays Here the expressions of the denominators are given by 20) and the six numerators N i read The third and the sixth terms in I(u) are five-point loop integrals, and the other terms are fourpoint loop integrals since the denominator C 1 is independent of the loop momentum q. Before going to calculate the integral, we would present a short analysis of its IR properties. When one of the gluons goes on-shell, i.e., q → p, the denominators D 2 , D 3 (with u = 1 − m q /m), D 5 tend to zero. Following to the procedure in Ref. [34], we can find that the one-loop integrals in the individual Feynman diagram for J/ψ → γη ( ) have soft singularities, which depend on the momentum fraction u, and can make the convolution integral over u become sensitive to the shapes of the η ( ) DAs. Moreover, the divergent term due to light quark pole would result in the final numerical results strongly dependent on the light quark masses. However, summing up the six Feynman diagrams, one can obtain When q = p + λ (λ is a small quantity), the numerator, which arises from the summation of all Feynman diagrams, 23) and the on-shell propagators For the ultrasoft gluon(λ → 0), the contributions to the loop integral have the form It means the sum of the loop integrals is IR safe.
With algebraic identities the loop function I(u) can be decomposed into a sum of four-and three-point one-loop integrals which can be analytically calculated with the technique proposed in Ref. [35] or the computer program P ackage − X [36,37]. Performing the convolution integral between the loop function I(u) and the DA φ q (u), we find that the function H q in Eq. (2.18) is very insensitive to the light quark masses. Specifically, the change of the absolute value of the function H q does not exceed 3% when the value of m q goes from 0 to 100 MeV for all the three kinds of η ( ) DAs in Fig. 2. Actually, when the light quark propagator is near its mass shell, i.e., q ∼ ±(u −ū)p, the factor (k · pq 2 − k · qp · q) in the numerator of Eq. (2.22) tends to zero and cancels the light quark propagator pole. In Fig. 3, we show the dependence on m q of the functions H η q = H q | m=mη and H η q = H q | m=m η with a asymptotic DA. It is noticed that the one-loop QCD contribution to the decay J/ψ → γπ 0 vanishes as a consequence of the antisymmetrical flavor wave function of π 0 , even for m u = m d , which disagrees with the result in Ref. [21]. For simplicity, we can take the following limit safely that is H q (q = u, d, s) = H 0 . The expression of I 0 (u) is presented in the Appendix. Numerically, we find that the loop function I 0 (u) is quite steady over the most region of u. In Fig. 4, we show the u dependence of I η 0 (u) = I 0 (u)| m=mη and I η 0 (u) = I 0 (u)| m=m η with the range u ∈ (0, 1). Unlike the result in the limit of m 2 /M 2 → 0 [18], the loop functions I η 0 (u) and I η 0 (u) change slowly over the momentum fraction u near the endpoints. As a consequence, the convolution integral between the loop function I 0 (u) and the DA becomes insensitive to the shapes of η ( ) DAs. For example, the difference between the convolution integral of I 0 (u) with a "narrow" DA (model I in Fig. 2) and the one with a "broad" DA (model II and III in Fig. 2) is less than 2%.
After these analyses, we return to the remaining calculations. The helicity amplitude H q QCD in Eq. (2.17) can be simplified to with the effective decay constants Besides the one-loop QCD contributions, QED processes J/ψ → γ * → γη ( ) can also contribute to the decays J/ψ → γη ( ) . The corresponding Feynman diagrams are depicted in Fig. 5, and the contribution reads where the dimensionless function is and the Q q represents the light quark charge.

(3.4)
For numerical calculations, all the values of meson masses, Γ J/ψ and f π are quoted from the PDG [44]. We use the FKS scheme for the η − η mixing [10], and then the effective decay constants f q η ( ) can be parameterized as For the three phenomenological parameters, i.e., the mixing angle φ and the decay constants f q , f s , they have been determined in different methods [10,[45][46][47][48]. Here we take the up-to-date values from Refs. [45,48] which are recapitulated in Table 2. The values in the first line are Table 2: The values of φ, f q and f s obtained with three phenomenological approaches [45,48]. [45] 40.6 ± 0.9 1.10 ± 0.03 1.66 ± 0.06 ηTFF [48] 40.3 ± 1.8 1.06 ± 0.01 1.56 ± 0.24 η TFF [48] 33.5 ± 0.9 1.09 ± 0.02 0.96 ± 0.04 extracted from the low energy processes(LEPs) V → η ( ) γ, η ( ) → V γ(V = ρ, ω, φ). The second and the third lines are the results extracted with rational approximations for the η ( ) transition form factor(TFF) F γ * γη ( ) (Q 2 ). The results of both the first and the second lines are generally consistent with the known FKS results [10] where the ratio R J/ψ with the approximation of nonperturbative matrix elements 0|G a µνG a,µν |η ( ) was adopted. While in the third line, the parameters are extracted with rational approximations for the η TFF F γ * γη (Q 2 ), which is in accord with the BaBar measurement in the timelike region at q 2 = 112 GeV 2 [49].
The Gegenbauer moments c q 2 (µ), c q 4 (µ) still have large uncertainty as depicted in Table 1. Fortunately, the dimensionless function H q in Eq. (2.18) is insensitive to the shapes of the η ( ) DAs as we have shown in section 2.1. Furthermore, both H g and h q are an order of magnitude smaller than the H q , therefore, the uncertainty of the Gegenbauer moments impacts our numerical calculations lightly (less than 2% among the results with three different models in Table 1). So in the following numerical calculations, we choose the Model I with As known, it is very hard to give precise predictions for individual decay width. With the nonperturbative matrix elements 0|G a µνG a,µν |η ( ) , the decay widths of J/ψ → γη ( ) have been given in Ref. [6] Γ(J/ψ → γη ( ) ) = 2 3 5 2 3 10 where the factor m −8 c will bring very large uncertainty. While, in the pQCD approach employed in this paper, there also exist large uncertainty due to the factor α 4 s (µ). For comparison, in Table 3, we present results of these two methods with m c = 1.5 GeV and α s = 0.34 as benchmarks. For the results in the second column, the matrix elements 0|G a µνG a,µν |η ( ) are evaluated with the updated values of φ, f q and f s in the first line of Table 2. Generally, one may expect the order of magnitude could be correctly predicted by the both approaches. However, we find the pQCD estimation of B(J/ψ → γη) is too small to be comparable to its experimental one. The reason may be due to our choice of the inputs φ, f q and f s extracted from low energy processes, or our understanding of η − η mixing scheme is incomplete which is beyond the scope of this paper. Table 3: The branching ratios B(J/ψ → γη ( ) ) obtained with nonperturbative gluonic matrix elements and pQCD approaches. Recently, the BaBar collaboration [49] has made the measurements of the η and η TFFs at q 2 = 112 GeV 2 , which have challenged the theoretical prediction very much. It is true that precise theoretical estimation of the TFFs is hard, due to uncertainties in the φ, f q , f s , the DAs of η ( ) and even the mixing scheme at high energy. In the literature, there are extensive discussions of this hot topic [33]. In Ref. [48], Escribano et al. have presented a treatment of η and η TFFs with rational approximations and extracted these parameters in Table 2. With these inputs and α s = 0.34, we show our results in Table 4. From the Table, we find R J/ψ is in good agreement with its experimental value R exp J/ψ = 4.65 ± 0.21 [44], only when we adopt the set of parameter values of η TFF. For other two sets of parameter values, B(J/ψ → γη) is estimated to be too small, which results in R J/ψ two orders higher than R exp J/ψ . In Table 5, we present results with the contributions due to the gluonic content of η ( ) and the one from QED processes J/ψ → γ * → γη ( ) . We can find that such contributions enhance B(J/ψ → γη) much. However, B(J/ψ → γη) is still far from its experimental value for both LEPs and ηTFF values of φ, f q and f s . Only η TFF set of parameter values can give B(J/ψ → γη ( ) ) and R J/ψ comparable with their experimental data. In the remaining part of this section, we would present a determination of φ without the input values in Table 2. The ratio R J/ψ in our calculation is Since the scalar functions are insensitive to the shapes of the η ( ) DAs, the ratio R J/ψ mainly depends on the angle φ and the ratio f s /f q . However, with the help of the ratio and the experimental values [44,53] Γ exp (η → γγ) = 4.36 (14) KeV, Γ exp (η → γγ) = 0.516(18) KeV, (3.10) the ratio R J/ψ becomes a function which only depends on the mixing angle φ. In Fig. 7, we show the dependence of the ratio R J/ψ on the mixing angle φ. As displayed by the horizontal dashed lines in Fig. 7 of the experimental measurement R exp J/ψ = 4.65 ± 0.21 [44], one can find which is in good agreement with η TFF result φ = 33.5 • ± 0.9 • [48] † , but in clear disagreement † It is noticed that the predicted η TFF is lim Q 2 →+∞ Q 2 F ηγ * γ (Q 2 ) = (0.160 ± 0.024) GeV, which is not in line with the BaBar measurement q 2 |F ηγ * γ (q 2 )| q 2 =112 GeV 2 = (0.229 ± 0.031) GeV, while the predicted η TFF, lim Q 2 →+∞ Q 2 F η γ * γ (Q 2 ) = (0.255 ± 0.004) GeV, is in accord with the BaBar measurement q 2 |F η γ * γ (q 2 )| q 2 =112 GeV 2 = (0.251 ± 0.021) GeV. More discussions could be found in Ref. [48].
Last but not least, we would compare our result of the loop function with the one obtained by Kühn [18]. Since, besides the approximation M ≈ 2m c , no further approximation is made, so our loop function I 0 (u), as shown in the Appendix, is much more complicated than the one in Ref. [18]. However, in the limit of m 2 /M 2 → 0, our I 0 (u) is reduced to the weight function W(u) in the Eq. (11) of Ref. [18].

Summary
In this paper, we have revisited the radiative decays J/ψ → γη ( ) in detail in the framework of pQCD. Comparing to the pioneer work [14], we do not take the weak-binding approximation for the final mesons and evaluate the involved one-loop integrals analytically with the light quark masses kept. Moreover, we also consider the contributions from the QED processes and the gluonic contents of η ( ) . Different from the results obtained in Ref. [21], our numerical results are insensitive to the light quark masses. In addition, we find that the B(J/ψ → γη ( ) ) are insensitive to the shapes of the η ( ) DAs.
Using three sets of values for φ, f q and f s , namely LEPs [45], ηTFF [48] and η TFF [48], we have presented our numerical results in Table 4 and 5, where B(J/ψ → γη) is too small to be comparable with the experimental one for the parameters extracted from LEPs and η TFF.
Only with the values extracted from η TFF, which is in accord with the BaBar measurement in the timelike region at q 2 = 112 GeV 2 [49], our results of B(J/ψ → γη ( ) ) and R J/ψ are in good agreement with the experimental data.
As a crossing check, we use our calculation of R J/ψ and Γ exp (η ( ) → γγ) as inputs to extract the value of φ, and find φ = 33.9 • ± 0.6 • , which is in good agreement with the η TFF one φ = 33.5 • ±0.9 • remarkably. However, such a small φ differs too much from the φ = 40.6 • ±0.9 • extracted from low energy processes and J/ψ → γη ( ) with nonperturbative matrix elements 0|G a µνG a,µν |η ( ) due to U A (1) anomaly dominance argument. The difference may arise from g * g * − η ( ) TFF used in our calculation, in like manner, γ * γ − η TFF in the extraction of φ from the BaBar measurement of q 2 |F η γ * γ (q 2 )| q 2 =112 GeV 2 . Anyhow, the physics under the difference is interesting and certainly worth further investigations. , In the above expressions, the dilogarithm function Li 2 (x) and the logarithm function ln(x) without "i " prescription are defined as Li 2 (x) = lim →0 + Li 2 (x − i ), ln(x) = lim →0 + ln(x + i ).
While the matrix elements 0|G a µνG a,µν |η ( ) are twist-4 effects which are suppressed by the factor of m 2 /M 2 relative to the leading twist terms. That is the reason why the mass dependence of the decay widths obtained in this work differs from that obtained in Ref. [6] by an additional factor M 4 /m 4 . In additions, when the factor m 2 /M 2 is not very small, the suppression for the absorptive part is no longer operative. For example, the dimensionless function H η 0 = H 0 | m=mη is dominated by its dispersive part, while the dispersive part and the absorptive part of the H η 0 = H 0 | m=m η are comparable.