Revisiting the $P$-wave charmonium radiative decays $h_{c}\rightarrow\gamma\eta^{(\prime)}$ with relativistic corrections

The $P$-wave charmonium decays $h_{c}\rightarrow\gamma\eta^{(\prime)}$ are revisited by taking into account relativistic corrections. The decay amplitudes are derived in the Bethe-Salpeter formalism, in which the involved one-loop integrals are evaluated analytically. Intriguingly, from both the quark-antiquark content and the gluonic content of $\eta^{(\prime)}$, the relativistic corrections make significant contributions to the decay rates of $h_{c}\rightarrow\gamma\eta^{(\prime)}$. By comparison with the leading-order contributions from the quark-antiquark content (one-loop level), the ones from the gluonic content (tree level) are also important, which is compatible with the conclusion obtained without relativistic corrections. Usually, for $\eta$ production processes, the predicted branching ratios are sensitive to the angle of $\eta-\eta^{\prime}$ mixing. As an illustration, using the Feldmann-Kroll-Stech result about the mixing angle $\phi=39.3^{\circ}\pm1.0^{\circ}$ as input, we find that the predicted ratio $R_{h_{c}}=\mathcal{B}(h_{c}\rightarrow\gamma\eta)/\mathcal{B}(h_{c}\rightarrow\gamma\eta^{\prime})$ is much smaller than the experiment measurement. While, with $\phi=33.5^{\circ}\pm0.9^{\circ}$ extracted from the asymptotic limit of the $\gamma^{\ast}\gamma-\eta^{\prime}$ transition form factor, we obtain $R_{h_{c}}=30.3\%$ in consistent with $R_{h_{c}}^{exp}=(30.7\pm11.3\pm8.7)\%$. As a cross-check, the mixing angle $\phi=33.8^{\circ}\pm2.5^{\circ}$ is extracted by employing the ratio $R_{h_{c}}$, and a brief discussion on the difference in the determinations of $\phi$ is given.

Turning to P -wave charmonia decays, the physical picture seems more complex, since the higher Fock-state contributions and the relativistic corrections may become important. It is well known that infra-red (IR) divergences are encountered in the color-singlet state contributions for the inclusive P -wave charmonia decays with the zero-binding approximation [52][53][54][55]. Although these IR divergences can be removed by considering the higher Fock-state contributions from the point of view of nonrelativistic QCD (NRQCD) [56,57], they may imply that the effects beyond those contained in the derivative of the nonrelativistic wave function at the origin R (0) may play a key role. Generally, it should be noted that similar IR divergences do not appear in exclusive P -wave charmonia decays [58][59][60][61]. Nevertheless, as pointed out in Refs. [58][59][60][61][62][63][64], the higher-order contributions, such as the higher Fock-state contributions [58][59][60][61]64] and the relativistic corrections [62][63][64], are still important to exclusive P -wave charmonia decays.
For the exclusive P -wave charmonium decays h c → γη ( ) , there are a few studies [65][66][67] in the theoretical aspect ever since the branching ratios B(h c → γη ) and B(h c → γη) are first measured to be, respectively, (1.52±0.27±0.29)×10 −3 and (4.7±1.5±1.4)×10 −4 by the BESIII Collaboration [28]. In the nonrelativistic limit [65,67], the relativistic corrections related to the internal momentum of the P -wave charmonium h c have been neglected in the calculation of the decay rates, and all the nonperturbative effects are absorbed in R hc (0) with the Taylor expansion of the hard-scattering amplitudes up to the linear terms. Then it is found that the calculations are IR safe and the predicted branching ratios B(h c → γη ( ) ) are much smaller than the experimental measurements. Obviously, this indicates that the relativistic corrections or/and the contributions from the higher Fock-state of h c are highly significant. While, from the point of view of NRQCD, the next-to-leading-order Fock-state contributions are suppressed by a relative factor v 2 cc α s in the decays h c → γη ( ) [65]. So it means that the relativistic corrections are needed in the exclusive P -wave charmonium decays h c → γη ( ) .
One of the major concerns of this paper is to study the relativistic corrections in the exclusive P -wave charmonium decays h c → γη ( ) by performing an explicit calculation. To make clear these relativistic corrections, the Bethe-Salpeter (B-S) framework [68][69][70] is used to calculate the wave function of h c and the decay amplitudes of h c → γη ( ) , where the internal momentum of h c is retained in both the soft bound-state wave function and the hard-scattering amplitude. Here, it is worth noting that there are at least two sources of the relativistic corrections. One is from the kinematical corrections which appear in the annihilation amplitudes, and the other is from the dynamical corrections of bound-state wave function itself. For the final light mesons η ( ) , light-cone distribution amplitudes (DAs) are adopted because of the large momentum transfer.
And the contributions of the quark-antiquark content and those of the gluonic content of η ( ) are both taken into account in our calculations.
In this paper, with the technique of helicity projector, we evaluate analytically the involved one-loop integrals with the internal momentum of h c kept. For the contributions from the quark-antiquark content of η ( ) in the decays h c → γη ( ) , the relativistic effects mainly originate from the kinematic part of the annihilation amplitudes, especially when the internal momentum of h c makes the propagator near on-shell. For the contributions of the gluonic content of η ( ) in the decays h c → γη ( ) , the next-to-leading-order effects related to the internal momentum are not substantially suppressed in the major region of the wave function of h c , and therefore the corresponding relativistic corrections are extremely important. Furthermore, we find out that the gluonic contributions and the quark-antiquark contributions are comparable, unlike the situation in the heavy vector quarkonium decays V → γη ( ) [36,40,71] where the gluonic contributions are strongly suppressed due to the special form of the spin structure of their amplitudes. In addition, it is also unlike the phenomenological fits [72][73][74] where the gluonic content of η can be neglected. This signifies that the decays h c → γη ( ) can be used to test the gluonic content of the η ( ) more efficiently than the decay processes V → γη ( ) . It is worthwhile to point out that the decay rates of h c → γη ( ) are insensitive to the light quark masses and the shapes of the η ( ) DAs [67], so the theoretical uncertainties from the η ( ) DAs are negligible, and the mixing angle of the η − η system could be reliably extracted in our calculations.
The paper is organized as follows. The formalism for the decays h c → γη ( ) is presented in section 2. In section 3, we obtain our numerical results, and the final section is our summary.
The expressions of the numerators involved in section 2 are given in the Appendix.

Bethe-Salpeter equation
It is generally known that the B-S equation [75,76] is an effective relativistic equation describing a bound state and has a solid basis in quantum field theory. So it is a conventional approach to treat various relativistic bound-state problems. In this subsection, we briefly review the formulation of the B-S framework. For charmonia, the B-S equation has the form [77][78][79][80] where K(K, q, q ) represents the interaction kernel between the internal quark and antiquark, and S F (p) = i/( / p −m c + i ) represents the propagator with the effective mass of c quarkm c .
The momenta of the quark and antiquark can be written as where q and K represent the internal momentum and the total momentum of the charmonia respectively.
For convenience, one can divide the internal momentum q into two parts. One part is the transverse componentq withq · K = 0, and the other is the longitudinal component q which is parallel to the total momentum K: Here both q K = q·K M andq 2 = q 2 − q 2 K are Lorentz invariant variables and M is the mass of the charmonia.
Then the B-S wave function can be expressed as where the hadron-quark vertex function reads Using the operators [78][79][80]85] the propagators can be decomposed as . For the axial vector meson h c , the Salpeter wave function can be approximately written as [79,80,[86][87][88] where M and ε(K) are the mass and the polarization vector of h c respectively, and the front factorq · ε(K) indicates that the wave function is of P -wave nature mainly, and f (q 2 ) is a scalar function ofq 2 . In the rest frame of h c , the momenta K andq have the form 11) and the scalar function f (q 2 ) satisfies the harmonic oscillator equation (more details and dis-cussions could be found in Refs. [80,85]). The expression of f (q 2 ) reads where N A is the normalization constant and β A is the harmonic oscillator parameter [80,85].
The normalization equation of f (q 2 ) reads In the rest frame of h c , the amplitude of h c → γg * g * has the form [35, 67, 89-91] where k, k 1 , k 2 and (k), (k 1 ), (k 2 ) stand for the momenta and polarization vectors of the photon and the gluons respectively, the factor √ 3 is included to account for the color properties of the quark-antiquark content, O(f,f ) is the hard-scattering amplitude, and the momenta of the quark and antiquark read with the c quark mass m c .
The light-cone expansion of the matrix elements of η ( ) over quark and antiquark fields reads [3,92,93] where p represents the momentum of η ( ) , the superscript q = u, d, s denotes the flavor of the light quarks, the "· · · " stands for the high twist terms and the decay constants f q η ( ) are defined as Then one can obtain the coupling of g * g * − η ( ) up to twist-3 level [94][95][96]: Here u is the momentum fraction carried by the quark, m q is the mass of the quark (q = u, d, s). The light-cone DA has the form [97] with the asymptotic form of DA φ AS (u) = 6u(1 − u) and the Gegenbauer moments c q n (µ). In Table 1, we list three models of the DAs given in Ref. [97]. Schematically, we also show their shapes at the scale of µ 0 = m c in Fig. 2. As pointed out in Refs. [40,67], the decay rates of h c → γη ( ) barely depend on the shapes of η ( ) DAs (we will estimate them below). It means that the mixing angle of η − η system could be reliably extracted in our calculations due to the negligible uncertainties from η ( ) DAs. Table 1: Gegenbauer moments of three sample models at the scale of µ 0 = 1 GeV. To proceed, the decay amplitude of h c → γη ( ) can be obtained directly by contracting the two couplings A αβµν and M µν , inserting the gluon propagators and integrating over the loop momentum: where the factor 1/2 takes into account that the two gluons have already been interchanged in both A αβµν and M µν . By Lorentz invariance, parity conservation and gauge invariance, one can obtain [35] i.e., there is only one independent helicity amplitude H q QCD : where With the help of the helicity projector [35] one can obtain the helicity amplitude where the dimensionless function H q reads I q (u,q) represents the sum of the loop integrals of all the Feynman diagrams with l = k 1 − k 2 and the denominators of the propagators As shown in Eq. (2.30), the spin structures of the bound-state wave function are absorbed into the loop function I q (u,q), and the expressions of the numerators N 1 , N 2 and N 3 are presented in the Appendix. Since the loop function I q (u,q) has no soft singularities and the dimensionless function H q is very insensitive to the light quark mass m q [40,67], one can take the following simplicity safely: i.e., H q (q = u, d, s) = H 0 . Then the helicity amplitude in Eq. (2.29) can be rewritten as with the effective decay constants By using the algebraic identity (ξ = ±1) with ξ = 1 − 2u and the mass of η ( ) meson m, the loop function I 0 (u,q) can be decomposed into a sum of four-point one-loop integrals When ξ = 1, the denominators of the propagators have the relation D 1 = D 4 , and the loop function I 0 (u,q) becomes And when ξ = −1, the denominators of the propagators have the relation D 1 = D 5 , then the loop function I 0 (u,q) becomes (2.39) With the program P ackage − X [98,99], one can evaluate the above one-loop integrals analytically. Similar to the situations without considering the internal momentum of charmonium [40,67], we find that the loop function I 0 (u,q) is also almost unchanged over the most

The contributions of the gluonic content of η ( )
The gluonic content of η ( ) can directly contribute to the decay processes h c → γη ( ) from tree level. One typical Feynman diagram is exhibited in Fig. 3, and there are other two diagrams from permutations of the photon and the gluon legs. Generally, contributions of gluonic content are supposed to be small [40,72], since gluonic content can be seen as the higher-order effects from the point of view of the QCD evolution of the two-gluon DA, which vanishes in the asymptotic limit. However, as we have pionted out in Ref. [67], these contributions may become with the effective decay constant f 1 One can obtain the corresponding helicity amplitude where the dimensionless function H g has the form Here the denominators of the propagators C 3 and C 4 read (2.44) and the expressions of the numerators N 4 , N 5 and N 6 are given in the Appendix.
In the remaining part of this section, we present a brief discussion about the relativistic effects. To the next-to-leading-order corrections related to the internal momentum from the numerators N 4 , N 5 and N 6 , if we takem c ≈ m c ≈ M/2 and m 2 /M 2 ≈ 0, they exhibit the following behavior: Obviously, the next-to-leading-order contributions are not suppressed enough in the major region † of the integration variableq. Therefore, in the decay processes h c → γη ( ) , the relativistic corrections should be taken into account.

Numerical results
The decay widths of h c → γη ( ) can be expressed as effect of quark-antiquark interaction, are respectively taken asm c = 1490 MeV and β A = 590 MeV, and more discussions could be found in Refs [79,80,85]. As we have already mentioned, the theoretical uncertainties from η ( ) DAs are negligible. So in our calculations, we choose the Model I of the meson DA in Table 1.

Branching ratios
The phenomenological parameters φ, f q and f s have been determined in [17] as  Tables 2, 3 and 4, where the branching ratios B(h c → γη), B(h c → γη ) and their ratio R hc = B(h c → γη)/B(h c → γη ) are presented in the first, second and third lines of these tables, respectively.
The contributions from the quark-antiquark content of η ( ) and those from the gluonic content of η ( ) are presented in Tables 2 and 3, respectively. The total contributions both from the quark-antiquark content and the gluonic content of η ( ) are presented in Table 4. In order to show the contributions from the relativistic effects more clearly, we present the results with the zero-binding approximation † (i.e., without relativistic corrections) [67] and those with the relativistic corrections in the first and second columns of these tables respectively. † We update the previous results [67] with the QCD running coupling constant α s (m c ) = 0.38.    Tables 2 and 3, one can find that the gluonic contributions and the quark-antiquark   contributions are comparable with each other, whatever the relativistic corrections are taken into account or not. It is unlike the situation where the gluonic contributions are strongly suppressed in the heavy vector quarkonium decays V → γη ( ) [36,40,71] because of an additional suppression factor (i.e., m 2 /M 2 ) from the spin structure of their amplitudes. Besides, our results are also unlike the phenomenological fits where the gluonic content can be neglected [72][73][74], especially for the meson η. By comparing with the results of the zero-binding approximation, we find out that the relativistic corrections are significant for both the quark-antiquark contributions and the gluonic contributions. Intriguingly, the importance of the gluonic contributions in the decay processes h c → γη ( ) indicates that these two decay processes can test the gluonic content of η ( ) more efficiently than the decay processes V → γη ( ) .
Comparing the results listed in Tables 2 and 3 with those listed in Table 4, we find that the branching ratios both B(h c → γη) and B(h c → γη ) are greatly enhanced with the constructive interference of the quark-antiquark contributions and the gluonic contributions, whatever the relativistic corrections are taken into account or not. Unfortunately, the branching ratio B(h c → γη) is still much smaller than its experimental value [28], even though it is substantially enhanced by the relativistic corrections. Moreover, it is thought-provoking that the ratio R hc with containing more dynamical corrections from the initial meson h c hardly change † and is still smaller than the experimental value [28]. These might imply the set of parameters of FKS results is not unquestionable.
As we have already mentioned, there are some obvious discrepancies in the determinations of the mixing angle [17, 34, 40, 43-45, 103, 105, 106]. Besides the value of the mixing angle φ ∼ 40 • (see, e.g., Refs [17,43,44,103,105,106]), a smaller value of the mixing angle φ ∼ 34 • is also usually obtained in many methods [34,40,45,103]. Therefore, it would be interesting to show the results of the branching ratios B(h c → γη), B(h c → γη ) and their ratio R hc with the different sets of the phenomenological parameters.
With the set of the parameter values [103] φ = 33.5 • ± 0.9 • , f q = (1.09 ± 0.02)f π , f s = (0.96 ± 0.04)f π , (3.4) extracted from the transition form factor (TFF) F γ * γη (+∞), which is in accord with the BABAR measurement in the timelike region at q 2 = 112 GeV 2 [111], we present the numerical results in Tables 5, 6 and 7 with Γ hc = 0.70 MeV.   certainly sensitive to the mixing angle of η − η mixing, and the main reason is that the decay amplitude is proportional to the factor ( √ 2f q cos φ − f s sin φ), which lead to large cancellations in the matrix elements. More interestingly, after taking into account the relativistic corrections, not only the ratio R hc but also the individual branching ratios B(h c → γη) and B(h c → γη ) are in very nice agreement with their experiment data [28]. Furthermore, we find that relativistic corrections increase the individual branching ratios B(h c → γη) and B(h c → γη ) by about a factor of 2, which are independent on the choice of the phenomenological parameters. This may imply that relativistic corrections are also important in other decay processes of the P -wave h c .
In addition, the gluonic contributions are comparable with the quark-antiquark contributions, and this is also independent on the choice of the phenomenological parameters.

η − η mixing
As a cross-check, we give a prediction of the mixing angle φ in the B-S formalism. By the ratio and the experimental measurements [28,102,112] R exp hc = (30.7 ± 11.3 ± 8.7)%, Γ exp (η → γγ) = 4.34 (14) keV, one can obtain the mixing angle where the uncertainty comes mainly from the R exp hc . It is worth noting that the recent lattice QCD calculations give the following values: φ = 34 • ± 3 • from the UKQCD collaboration [45] and φ = 38.8 • ±3.3 • from the ETM collaboration [107]. Schematically, we show the dependence of the ratio R hc on the mixing angle φ in Fig. 4. Obviously, the prediction of the mixing angle in the B-S framework is consistent with the value φ = 33.5 • ± 0.9 • extracted from η TFF [103].
Moreover, it is also in good agreement with our previous determinations φ = 33.9 • ± 0.6 • [40] and φ = 33.8 • ± 2.5 • [67], which are obtained by the ratios R J/ψ and R hc without considering the relativistic corrections respectively.
It is worth noting that there are discrepancies in the determinations of the mixing angle φ, and it may imply an incomplete understanding of η − η mixing. Last but interestingly, in the same framework of perturbative QCD, we obtain a consistency check about the mixing

Summary
In this work, we have revisited the P -wave charmonium radiative decays h c → γη ( ) in the B-S formalism, where the internal momentum of h c has been retained in both the soft wave function ψ(q) and the hard-scattering amplitude O(q). The B-S wave function is employed to describe the bound-state properties of the P -wave charmonium h c , while the light-cone DAs are adopted for the light mesons η ( ) . And then the involved one-loop integrals are carried out analytically. It is found that the relativistic corrections from both the quark-antiquark content and the gluonic content of η ( ) make significant contributions to the decay rates of h c → γη ( ) .
In addition, the predicted branching ratios B(h c → γη ( ) ) are insensitive to the shapes of the η ( ) DAs, and the gluonic contributions as well as the quark-antiquark contributions are both important in these two decay processes. What is more, by the ratio R hc , the widths Γ(η ( ) → γγ) and their experimental values, we have obtained a consistency check about the mixing angle φ in the framework of perturbative QCD.
As aforementioned discussion, since the decay amplitude of the η channel is proportional to the factor ( √ 2f q cos φ − f s sin φ) which could lead to large cancellations in the matrix elements, the predicted branching ratio B(h c → γη) is sensitive to the angle of η − η mixing. This means that the branching ratios are generally hard to predict precisely, but would be more efficient to determine the mixing angle in η production and decay processes. On the other hand, by the comparison between the results without the relativistic corrections and these obtained in this work, we find that the relativistic corrections are rather significant in the exclusive P -wave charmonium decays h c → γη ( ) . This may imply that the relativistic effects should be taken into account in the production and decay processes of the higher excited charmonia, especially for the radially-excited states with nodes contained in their wave functions. Further investigations about these issues are certainly deserved.