QCD analysis of the P -wave charmonium electromagnetic Dalitz decays h c → η ( ′ ) ℓ + ℓ −

The P -wave charmonium electromagnetic Dalitz decays h c → η ( ′ ) ℓ + ℓ − ( ℓ = e, µ ) with large recoil momentum are investigated in the framework of perturbative QCD, and the contributions from the small recoil momentum region are described by the overlap of soft wave functions. The transition form factors f h c η ( ′ ) ( q 2 ) and the normalized transition form factors F h c η ( ′ ) ( q 2 ) in full kinematic region are derived for the first time. It is noticed that there are no IR divergences at one-loop level, and the transition form factors with the relativistic corrections from the internal momentum of h c are insensitive to both the shapes of η ( ′ ) distribution amplitudes and the invariant mass of the lepton pair in the large recoil momentum region. Intriguingly, unlike the situation in the S -wave charmonium decays J/ψ → η ( ′ ) ℓ + ℓ − , we find the contributions from the small recoil momentum region are comparable with those from the large recoil momentum region in the P -wave charmonium decays h c → η ( ′ ) ℓ + ℓ − . By employing the obtained F h c η ( ′ ) ( q 2 ), we give the predictions of the branching ratios B ( h c → η ( ′ ) ℓ + ℓ − ), which may come within the range of measurement of present or near-future experiments.

The P -wave charmonium h c cannot be directly produced in e + e − collisions because of the quantum numbers J P C = 1 +− , but it can be produced through ψ(2S) → π 0 h c [31,32].In recent years, many more decay modes of h c have been searched for at BESIII [32][33][34][35][36][37][38].Due to the negative C parity, the h c most likely decays into a photon and a pseudoscalar meson η c or η (η ′ ), in which the radiative decays h c → γη (′) have first been observed by the BESIII Collaboration using about 0.4 billion ψ(2S) events [33].So far, there are around 3 billion ψ(2S) events collected with the BESIII detector [39][40][41], and it represents about an order-of-magnitude increase in statistics.This provides a good opportunity to study the EM Dalitz decays h c → η (′) ℓ + ℓ − , and their branching ratios can be reached by present or near future experiments, especially for the η ′ channels (reaching 10 −5 ).These EM Dalitz decays could not only offer useful information to constrain theoretical models (as mentioned above) in the charmonium region, but also shed light on the transition mechanism of h c → η (′) and the η − η ′ mixing effects [25,42,43] in different kinematic regions.Besides, it is more interesting for the P -wave charmonia decays.Generally, inclusive P -wave charmonia decays suffer from IR divergences in the color-singlet state contributions with the zero-binding approximation [44][45][46]; while the similar IR divergences do not appear in the exclusive P -wave charmonia decays [47][48][49][50].This may imply that the effects beyond those contained in the derivative of the nonrelativistic wave function at the origin play a key role.Recently, it has been pointed out that the relativistic corrections from the internal momentum of h c are extremely important in the decays h c → γη (′)  [51].This indicates that the relativistic corrections may also be important in the EM Dalitz decays h c → η (′) ℓ + ℓ − due to the same EM structure arising at the h c → η (′) transition vertex.
In this paper, one of the major concerns is to clarify the dynamical picture of the EM Dalitz decays h c → η (′) ℓ + ℓ − in different kinematic regions.Phenomenologically, there exist three types of contributions in the decay processes h c → η (′) ℓ + ℓ − : (i) In the large recoil momentum region, i.e., the square of the invariant mass of the lepton pair q 2 = m 2 ℓ + ℓ − ≃ 0, the transition mechanism of h c → η (′) could be described by the perturbative QCD approach, which has been reliably employed to treat the corresponding radiative decays h c → γη (′) [51,52].And we call this transition mechanism the hard mechanism.(ii) In the small recoil momentum region, i.e., the square of the invariant mass of the lepton pair q 2 ≃ q 2 max = (M hc − m η (′) ) 2 , the transition mechanism of h c → η (′) is governed by the overlapping integration of the soft wave functions, and this transition mechanism is the so-called wave function overlap [53][54][55][56][57][58][59][60][61](i.e., the soft mechanism).The TFFs f hcη (′) (q 2 ) account for the size effects from the spatial wave functions of the initial-and final-state hadrons.(iii) In resonance regions, such as q 2 ≃ m 2 ρ , m 2 ω , m 2 ϕ , the transition mechanism of h c → η (′) can be universally described by the vector meson dominance (VMD) model [62], in which the resonance interaction between photons and hadrons is predominant.However, on the one hand the contributions from VMD are negligibly small due to the narrow widths of resonances (see Ref. [15] for more details), and on the other hand there are still some open questions for the VMD model, such as the sign ambiguity in the amplitude from the intermediate vector mesons and the off-mass-shell effects of the coupling constants [63].
To make the dynamical picture of the EM Dalitz decays h c → η (′) ℓ + ℓ − clear, we will mainly present the detailed discussions about the hard mechanism and the soft mechanism in the later parts of this paper.In the large recoil momentum region, by employing the Bethe-Salpeter (B-S) framework [51,[64][65][66][67][68], we work out the B-S wave function of h c , in which the internal momentum is retained.Considering a large momentum transfer, one can adopt the light-cone distribution amplitudes (DAs) to describe the internal dynamics of the final light mesons η (′) , and the involved quark-antiquark and gluonic contents of η (′) are taken into account in our calculations.By an analytic calculation of the involved one-loop integrals, we find that the TFFs f hcη (′) (q 2 ) are UV and IR safe, and they barely depend on the shapes of the light meson DAs.Furthermore, the gluonic contributions and the quark-antiquark contributions are comparable in the TFFs f hcη (′) (q 2 ).It is compatible with the situation in the corresponding radiative decays h c → γη (′) [51,52].In the small recoil momentum region, the TFFs are calculated phenomenologically by the wave function overlap.Through a detailed calculation, we obtain the TFFs f hcη (′) (q 2 ) in the whole kinematic region for the first time.It is worthwhile to point out that the contributions from the soft mechanism and those from the hard mechanism are comparable with each other in the branching ratios B(h c → η (′) ℓ + ℓ − ), unlike the situation in S-wave charmonium EM Dalitz decays J/ψ → η (′) ℓ + ℓ − [15] where the soft contributions are suppressed because of the special form of the spin structure of their amplitudes.In order to remove the main uncertainties arising from the bound-state wave functions, we use the normalized TFFs F hcη (′) (q 2 ) ≡ f hcη (′) (q 2 )/f hcη (′) (0) to obtain the predictions of the branching ratios The paper is organized as follows: The theoretical framework for the EM Dalit decays h c → η (′) ℓ + ℓ − is presented in detail in Sec. 2. In Sec. 3 we show our numerical results and some phenomenological discussions, and Sec. 4 is our summary.In the large recoil momentum region of η (′) , the EM Dalitz decays h c → η (′) ℓ + ℓ − can be described by the perturbative QCD approach.The leading order Feynman diagrams for the quark-antiquark content of η (′) arise from one-loop QCD processes.One of them is illustrated in Fig. 1, and the other five diagrams come from permutations of the photon and the gluon legs.

THEORETICAL FRAMEWORK
Here and in what follows, the involved kinematical variables are labeled in Fig. 1, where u and ū are the momentum fractions carried by the light quark and the light antiquark, respectively.
According to the Feynman diagrams, one can obtain the amplitude of where A αβ represents the amplitude of h c → η (′) γ * ; K and ε(K) stand for the momentum and polarization vector of h c , respectively; q stands for the momentum of the virtual photon, and q 2 = m 2 ℓ + ℓ − is the square of the invariant mass of the lepton pair; l 1 and l 2 stand for the momenta of the leptons ℓ − and ℓ + , respectively.It is convenient to convert the amplitude of h c → η (′) γ * into two parts: the effective coupling of the process h c → g * g * γ * and that of the process g * g * → η (′) .By multiplying the two parts, inserting the gluon propagators and performing the loop integrations, one can obtain the final amplitude of h c → η (′) γ * .
In the rest frame of the P -wave charmonium h c , the amplitude of h c → g * g * γ * can be written as [69-71] where Ψ(K, k) represents the B-S wave function of h c ; O(f, f ) represents the hard-scattering amplitude; √ 3 represents the color factor; ϵ(q) represents the polarization vector of the virtual photon; k 1 , k 2 and ϵ(k 1 ), ϵ(k 2 ) represent the two gluons' momenta and polarization vectors; f and f represent the momenta of the quark c and antiquark c, and they read with k the relative momentum between the quark c and antiquark c, i.e., the internal momentum of the P -wave charmonium h c .Here M is the mass of h c .For convenience in subsequent calculations, we divide the internal momentum of h c into two parts: the transverse component k with k • K = 0 and the longitudinal component k ∥ with k ∥ • k = 0, i.e., where both k K = k•K M and k2 = k 2 − k 2 K are Lorentz invariant variables.Considering the rest frame of h c , one can easily know that k involves 3 degrees of freedom (namely, the component k) orthogonal to the total momentum K, and k K contains the remaining 1 degree of freedom (namely, the component k 0 ).Now the volume element of the internal momentum k can be expressed in the form d 4 k = d 3 kdk K .Furthermore, with a more relevant treatment k 0 ≪ M , we obtain the momenta and the hard-scattering amplitude and this treatment maintains the gauge invariance of the hard-scattering amplitude [72].
By employing the B-S equation [64,65] of the P -wave charmonium h c , one can reduce the B-S equation to the Salpeter equation under the covariant instantaneous ansatz (CIA) [66][67][68].
The Salpeter wave function is defined as Subsequently, we obtain an analytic Salpeter wave function of h c by solving the Salpeter equa-tion (more details can be found in our recent investigation [51]), where mc is the effective mass of c quark, and the front factor k • ε(K) indicates that the wave function is the nature of P -wave, and the scalar function f ( k2 ) reads (2.9) with N A the normalization constant and β A the harmonic oscillator parameter.And the normalization equation of f (q 2 ) reads Using Eqs.(2.6) and (2.7), we can rewrite the amplitude of h c → γg * g * : where the hard-scattering amplitude O( k) reads with the c quark mass m c .
By contracting the above two amplitudes, inserting the gluon propagators and integrating over the loop momentum, we obtain the decay amplitude of Considering parity conservation, Lorentz invariance, gauge invariance, and current conservation, one knows

17)
Then the h c → η (′) γ * TFFs can be defined by With the help of the projection operator the TFFs can be rewritten as Here we show the expression of the TFFs more clearly with u and ū = (1 − u) the momentum fractions arising from the light mesons η (′) .The expressions of the denominators read and the expressions of the numerators N 1 ∼ N 6 are presented in the Appendix.
With the help of the algebraic identity (u ̸ = 0, 1) ) can be decomposed into a sum of four-point and three-point one-loop integrals.When u = 0, 1, the denominators of the propagators have the relation and the TFFs can also be decomposed into four-point or three-point integrals.Then one can analytically evaluate these one-loop integrals with the technique proposed in Refs.[86][87][88] or the computer program PACKAGE-X [89,90].It is found that the TFFs are UV and IR safe.
Similarly to the situation in the radiative decays h c → γη (′) [51,52], the TFFs are insensitive to the DAs of η (′) .Our numerical results show that the change of the modulus of the TFFs does not exceed 1% with the different models of the DAs.Therefore, the theoretical uncertainties from the DAs are ignorable in our calculations of the TFFs, and we choose model I of the meson DA in Table 1 of Ref. [52].

Contributions of the gluonic content of η (′)
Generally speaking, the contributions of the gluonic content of η (′) are expected to be small because the gluonic content cloud be seen as higher-order effects from the point of view of the QCD evolution of the two-gluon DA, which vanishes in the asymptotic limit.For example, the gluonic contributions are strongly suppressed by the factor m 2 η (′) /M 2 J/ψ in the radiative decays J/ψ → γη (′) [76].However, there is no suppression factor in the P -wave charmonium radiative decays h c → γη (′) [51,52], due to the special form of the spin structure in their amplitudes.So the gluonic contributions may become important in the radiative decays of the P -wave charmonium h c .In fact, as pointed out in Refs.[51,52], the gluonic contributions and the quark-antiquark contributions are comparable with each other in the radiative decays h c → γη (′) .Obviously, this situation should be found in the large recoil momentum region of η (′)  of the decays h c → η (′) ℓ + ℓ − , due to the same spin structures in their hadronic matrix elements.
The corresponding Feynman diagram is depicted in Fig. 2, and the other two diagrams arise from permutations of the photon and the gluon legs.
with the gluonic content of η (′) .The kinematical variables are labeled.
At the leading twist level, the matrix elements of the mesons η (′) over two-gluon fields in the light-cone expansion can be written as [77,91,92] where n = 1 √ 2 (1, − p |p| ) is a lightlike vector along the opposite direction of the mesons η (′) [91], ) are the effective decay constant, and the gluonic twist-2 DA is [77,92,93] (2.25) After a series of calculations, we obtain the corresponding TFFs where the expressions of the denominators read and the expressions of the numerators N 7 ∼ N 9 are presented in the Appendix.
Performing the integral calculations of the TFFs f Q hcη (′) (q 2 ) and f G hcη (′) (q 2 ), we find that the modulus square of these TFFs is very insensitive to the dilepton invariant mass m ℓ + ℓ − (or, q 2 ).
In Fig. 3, the m ℓ + ℓ − dependence of the modulus square |f Q,G hcη (′) (q 2 )| 2 is shown.Schematically, we can clearly see that the modulus square |f Q,G hcη (′) (q 2 )| 2 has only negligible changes in the range of (0 − 1000) MeV.Comparing the quark-antiquark contributions from |f Q hcη (′) (q 2 )| 2 with the gluonic contributions from |f G hcη (′) (q 2 )| 2 , we find that the former is about twice that of the latter.In other words, the gluonic contributions and the quark-antiquark contributions are both important in these decay processes.It is compatible with the situation in the corresponding radiative decays h c → γη (′) [51,52].
Based on the foregoing discussions, in the large recoil momentum region of η (′) , the h c → η (′) γ * TFFs can be obtained by which includes the dynamical structure information from the quark-antiquark content and the gluonic content of η (′) .In addition, the relativistic corrections related to the internal momentum of h c are taken into account in the TFFs.Specifically, there exist the kinematical corrections from the annihilation amplitudes and the dynamical corrections from the boundstate wave function of h c .Since the physical picture in the large recoil momentum region of the electromagnetic Dalitz decays h c → η (′) ℓ + ℓ − is the same with the one in the corresponding radiative decays, these electromagnetic Dalitz decays can be calculated by the non-relativistic quark model with zero-binding approximation, as well as the Bethe-Salpeter formalism.Similarly to the situations in the radiative decays h c → η (′) γ [51,52], we find that these relativistic corrections also play an important role in the decay processes there is about a 2-times enhancement of |f H hcη (′) (q 2 )| 2 over the results with the zero-binding approximation.

Soft mechanism
Generally speaking, a perturbative QCD approach will become invalid in the three-body decay processes h c → η (′) ℓ + ℓ − with a small recoil momentum of η (′) , and a special handling is needed in principle.To deal with this issue properly, a picture of the soft wave function overlap is proposed in Ref. [15], where the picture has been proved valid in the decays J/ψ → ηe + e − with a small recoil momentum.Because of the similar physical picture, the TFFs of the h c → η (′) γ * in the small recoil momentum region, which are dependent on the recoil momentum |p η (′) | and reflect the size effects from the spatial wave functions of the initial-and final-state hadrons, can be adopted the empirical form in the rest frame of h c [53,57,60,94]: Here − γ * coupling and can be determined by the continuity condition of the TFFs between the large and the small recoil momentum regions, and the parameter β is an experiment-related quantity.We adopt β = 400 MeV, which is compatible with the fitted value in the J/ψ decays [15].And it is worth noting that the decay amplitude A αβ in soft mechanism has the form ef S hcη (′) (q 2 ) g αβ − q α K β /q • K , in which the spin structure is determined just by the quantum number of the initial-and final-state hadrons.
In the whole recoil momentum region of η (′) , the TFFs can be given by (2.30) Incidentally, the recoil momentum |p η (′) | = λ 1 2 (M 2 , m 2 , q 2 )/(2M ) is a monotonically decreasing function of the square of the invariant mass of the lepton pair q 2 .The recoil momentum is above 1 GeV when q 2 ≤ 1 GeV 2 .It is commonly asserted that perturbative QCD is self-consistent when the recoil momentum is above 1 GeV [95][96][97][98].Namely, the transition to perturbative QCD appears at about q 2 = 1 GeV 2 , and the hard mechanism begins to dominate as the q 2 decreases.
On the contrary, the contributions from the soft mechanism would become important with the q 2 increases.Although we could obtain the hard contributions from the large recoil momentum region with the perturbative QCD approach and the soft ones from the small recoil momentum region with the overlapping integration of the soft wave functions, how to precisely match these two contributions in the intermediate recoil momentum region is still an open question and needs further investigation.Even so, our description of the EM Dalitz decay processes h c → η (′) ℓ + ℓ − may constitute an important step forward toward a satisfactory description.
In the numerical calculations, all the values of the involved meson masses, quark masses, decay widths and decay constant are quoted from the Particle Data Group [99].By employing the two-loop renormalization group equation, we obtain the strong coupling constant α s (m c ) = 0.38.The effective mass of the c quark and the harmonic oscillator parameter appearing in the bound-state wave function are taken as mc = 1490 MeV and β A = 590 MeV respectively, and more discussions can be found in Refs.[100][101][102].For the Gegenbauer moments from η (′) DAs, we just adopt model I in Table 1 of Ref. [52] due to the negligibly small uncertainties from these Gegenbauer moments, which have been mentioned in the foregoing discussions.For the phenomenological parameters, i.e., the mixing angle ϕ and the decay constants f q(s) , we adopt the set of values [82] ϕ = 33.5 • ± 0.9 • , f q = (1.09± 0.02)f π , f s = (0.96 ± 0.04)f π (3.3) extracted from the TFF F γ * γη ′ (+∞), which is in excellent agreement with the BABAR measurement [103].And more discussions about these phenomenological parameters can be found in Refs.[11,51,52,76,80,104].

Hard mechanism Soft mechanism Total
We now proceed with a full calculation of the branching ratios B(h c → η (′) ℓ + ℓ − ), and our nu-merical results are shown in Table 1.Here we do not present the theoretical uncertainties which come mainly from B exp (h c → ηγ) = (4.7±2.1)×10−4 or B exp (h c → η ′ γ) = (1.5±0.4)×10−3 [33], and they are expected to be 30% ∼ 50%.In the second column, the results from the hard mechanism, which can be described by the perturbative QCD approach, are presented; in the third column, we show the results from the soft mechanism, which is governed by the overlapping integration of the wave functions of the initial-and final-state hadrons.And the total contributions from both the hard mechanism and the soft mechanism are presented in the last column.
First of all, it is noticed that the soft contributions in the decay processes h c → η (′) e + e − are equal to the ones in the decay processes h c → η (′) µ + µ − with an accuracy of more than three (two) significant digits.This is because the differential branching ratios are proportional to (1 − O(m 4 ℓ /q 4 )) in large q 2 region (see Eq. (3.1) or Eq.(3.2)).Therefore, the difference caused by the lepton mass m ℓ is ignorable in the small recoil momentum region.While, the hard contributions of the µ + µ − channel are about five times smaller than the ones of the e + e − channel, since the phase space of the former shrinks in the small q 2 region.Besides, we find that the hard contributions and the soft ones are comparable with each other.It is unlike the situation where the soft contributions are negligibly small in the S-wave charmonium EM Dalitz decays J/ψ → η (′) ℓ + ℓ − because of a suppression of the kinematic factor (i.e., |p η (′) | 3 /q) [15].
It is worthwhile to point out that there are around 3 billion ψ(2S) events collected with the BESIII detector so far [39][40][41], and this may imply that our predictions of the branching ratios B(h c → η (′) ℓ + ℓ − ) may come within the range of measurement of present or near-future experiments, especially for the η ′ channels (reaching 10 −5 ).
Last but not least, we concentrate on the q 2 -dependent TFFs f hcη (′) (q 2 ), which could provide the dynamical information of the EM structure arising at the h c → η (′) transition vertex and offer a powerful probe of the intrinsic structure of the P -wave charmonium h c .Experimentally, one is interested in the normalized TFFs F hcη (′) (q 2 ), because their modulus square |F hcη (′) (q 2 )| 2 can be directly extracted by comparing the measured invariant mass spectrum of the lepton pairs from the Dalitz decays with the point-like QED prediction [10,62].On the other hand, there exists a large uncertainty from sources such as the bound-state wave function and the QCD running coupling constant α s in the TFFs f hcη (′) (q 2 ).However, the dependence of the normalized TFFs F hcη (′) (q 2 ) (i.e., the ratios f hcη (′) (q 2 )/f hcη (′) (0)) on the bound-state wave function and the QCD running coupling constant is cut down to a large extent.So, one can expect that the predictions of the normalized TFFs are more reliable.Besides, the q 2 dependence of the TFFs f hcη (′) (q 2 ) is still retained in the normalized TFFs F hcη (′) (q 2 ) because of the constants f hcη (′) (0).In Fig. 4, we present the q 2 dependence of the modulus square of the normalized TFFs |F hcη (′) (q 2 )| 2 in their full kinematic region.Here it should be noted that, even though the normalized TFFs |F hcη (′) (q 2 )| 2 are independent of the lepton mass, the kinematic region of the e + e − channel is different from the one of the µ + µ − channel.Specifically, the value of the dilepton invariant mass m ℓ + ℓ − starts at 2m e in the e + e − channel and it starts at 2m µ in the µ + µ − channel.As shown schematically in Fig. 4, the difference in the modulus squares |F hcη (q 2 )| 2 and |F hcη ′ (q 2 )| 2 mainly arise from their phase space.One can find that the modulus square |F hcη (′) (q 2 )| 2 is quite steady in small q 2 region and increasing rapidly in large q 2 region, which is compatible with the situation in the EM Dalitz decay processes Jψ → η (′) ℓ + ℓ − [1,15].Present or near future experimental measurement is expected to provide tests for these predictions.

SUMMARY
In this paper, we investigate the P -wave charmonium EM Dalitz decays h c → η (′) ℓ + ℓ − with a QCD analysis.In the large recoil momentum region of η (′) , these decay processes are described by the perturbative QCD approach.For the primary heavy meson h c , we work out its B-S wave function in the framework of B-S equation, and its internal momentum is retained in both the wave function and the hard-scattering amplitude; for the final light mesons η (′) , the light-cone DAs are adopted due to a large momentum transfer.By an analytic calculation of the involved one-loop integrals, we find that the TFFs are UV and IR safe, and the gluonic contributions and the quark-antiquark contributions are both important in the TFFs.In the small recoil momentum region of η (′) , the picture of the soft wave function overlap is adopted to describe the transition mechanism of h c → η (′) .By relating to their radiative decay processes, the branching ratios B(h c → η (′) ℓ + ℓ − ) are obtained.Intriguingly, the contributions from the soft mechanism and those from the hard mechanism are comparable with each other, unlike the situation in S-wave charmonium decays J/ψ → η (′) ℓ + ℓ − [15] where the soft contributions are suppressed because of the special form of the spin structure of their amplitudes.Furthermore, the q 2 -dependent TFFs are analyzed briefly, and we obtain the q 2 dependence of the modulus square of the normalized TFFs |F hcη (′) (q 2 )| 2 in their full kinematic region.Lastly, it should be pointed out that there are around 3 billion ψ(2S) events collected with the BESIII detector so far [39][40][41], and this may imply that our predictions of the branching ratios B(h c → η (′) ℓ + ℓ − ) may come within the range of measurement of present or near-future experiments, especially for the η ′ channels.