NLO QCD Corrections to $\eta_c+h_c(1P)/\psi_{1,2}(1D)$ Production at B-factories

We calculate the next-to-leading order (NLO) quantum chromodynamics (QCD) corrections to double charmonium production processes $e^+e^-\to\gamma^*\to\eta_c+h_c(1P)/\psi_{1,2}(1D)$ within the non-relativistic QCD (NRQCD) factorization framework. We find that the corrections to $\eta_c+h_c(1P)$ production are positive, while those to $\eta_c+\psi_{1,2}(1D)$ are negative. Unlike the $J/\psi + \eta_c$ case, all the corrections here are not large. Uncertainties in the renormalization scale, quark mass and running energy of center-of-mass are discussed, and the scale dependence of these processes is found to be greatly reduced with the NLO QCD corrections.


Introduction
Once, the double charmonium inclusive(exclusive) processes e + e − → J/ψ + cc(η c ) encouraged many investigations because of large discrepancies between leading order (LO) calculations [1][2][3][4] and experimental results at B-factories [5,6]. It was found that this problem may be somehow remedied by next-to-leading order (NLO) Quantum Chromodynamics (QCD) corrections [7][8][9][10][11] in the framework of non-relativistic QCD (NRQCD) [12], which indicates that the NLO corrections might be significant to LO results. In order to further understand the NLO properties and investigate the convergence of perturbative expansion for heavy quarkonium production and decays in the NRQCD formalism, we need to study more heavy quarkonium production and decay processes.
There have been great advances in recent years in the calculation of radiative corrections to charmonium inclusive and exclusive production and decays. The LO estimations for h c inclusive production at B-factories were given in Refs. [13,14]. Recently, a complete NLO calculation for the P-wave charmonium χ cJ ( 3 P [1] J , 3 S [8] 1 ) inclusive production processes e + e − → χ cJ + cc/gg/qq(q = u, d, s) was carried out [15]. Though so far there are insufficient data on the double excited charmonium exclusive processes e + e − → γ * → H 1 +H 2 * , theoretical studies have already started [3,16], and even the NLO corrections for J/ψ+χ cJ production at B-factories [17,18] have been performed.
The rest of this paper is organized as follows: Section 2 presents our formalism and calculation method, numerical evaluation and some discussion of the results are given in Section 3, and Section 4 gives a summary and conclusions.

Calculation Scheme Description and Formalism
In our calculation the Mathematica package Fey-nArts [25] was applied to generate all the Feynman diagrams and amplitudes of partons. The standard projec-tion operators for charmonia may be expressed as [26]: Here, Π 0 corresponds to the spin-singlet charmonia with ǫ 0 = γ 5 , while Π 1 corresponds to spin-triplet states with ǫ 1 = ǫ * , the spin polarization vector. p = P 2 + q and p = P 2 − q are the momenta of the quark and antiquark within the charmonium, respectively. P denotes the momentum of quarkonium and q is the relative momentum between quarks inside the quarkonium. 1 stands for the unit matrix in color space. Using the method described in Ref. [3], and applying FeynCalc [27] to assist the calculation of amplitudes, one can readily obtain the tree-level results. For oneloop QCD corrections, the representative Feynman diagrams are shown in Fig. 1. We first used FeynCalc to trace the Dirac matrix chains as well as the color matrices, and to perform the derivative on momentum q. Next, the package $Apart [28] was employed to reduce the propagators of each one-loop diagram into linearly independent ones. Then, we applied the package FIRE [29] to reduce all the one-loop integrals into several master-integrals(A 0 , B 0 , C 0 , D 0 ). Finally, package Loop-Tools [30] was employed to evaluate the scalar masterintegrals numerically.
Throughout our calculation, we adopted Feynman gauge, and took the conventional dimensional regularization with D = 4 − 2ǫ to regularize the divergences. The ultraviolet divergences are canceled by the counter terms and the infrared divergences in the short distance coefficients cancel each other, which confirms the NRQCD factorization for e + e − → η c + h c (1P )/ψ 1,2 (1D) processes at NLO level. When handling the counter terms, we found that terms related to Z 3 , the renormalization constant corresponding to the gluon field, vanish in the end. The renormalization constant Z g , corresponding to the strong coupling constant α s , was computed in the modified-minimal-subtraction (M S) scheme, while Z 2 and Z m , corresponding respectively to the quark field and quark mass, were in the on-shell (OS) scheme. Eventually, the expressions for the relevant renormalization constants read:

Numerical results and analysis
Before carrying out numerical calculation, the input parameters need to be fixed. The NLO running coupling constant In numerical evaluation, the effective quark flavor number n f = 4 was adopted and Λ QCD = 0.297GeV [31]. At leading order in relativistic expansion M ηc/hc/ψ 1,2 = 2m c , and hence the charm quark mass m c = 1.7±0.2GeV was taken. The values of NRQCD matrix elements were evaluated from the Bunchmüller-Tye potential model [32], i.e.,  Table 3. Cross sections σ(µ)(fb) of e + e − → ηc + ψ2(1D) at leading order and next-to-leading order at µ = 2mc and √ s/2. (1) For η c +h c (1P ) production, the NLO correction enhances tree-level result, and the K-factor is bigger at large scale µ = √ s/2; (2) For η c + ψ 1,2 (1D) productions, the NLO corrections are negative, and the K-factors are also bigger at large scale µ = √ s/2; (3) At LO, cross sections of these three processes decrease with the scale µ increasing from 2m c to √ s/2. The NLO result of η c + h c (1P ) also decrease with the scale, while the NLO results of η c +ψ 1,2 (1D) have inverse relation with the energy scale in the discussed region; (4) The hierarchy σ ψ 1 < σ hc < σ ψ 2 remains for both LO and NLO results at scales µ = 2m c and µ = √ s/2.
To illustrate the results in Tables 1-3 more clearly, we show the LO and NLO cross sections of η c + h c (1P )/ψ 1,2 (1D) production with different m c versus renormalization scale µ in Figs. 2-4 respectively. One may notice that the scale dependence of the NLO results is obviously depressed, and the convergence of pertubative expansion for these processes works well, whereas the cross sections are quite sensitive to the quark mass.  In Figs. 5-7, we present the LO and NLO cross sections for e + e − → η c + h c (1P )/ψ 1,2 (1D) versus center-ofmass energy (E cm = √ s), respectively, where scale µ and mass m c are fixed. The figures show that NLO correction for η c +h c (1P ) production is positive in the displayed region, while that for η c + ψ 1 (1D) is negative. Interestingly, for the η c + ψ 2 (1D) channel, the NLO correction has an inversion at E cm = 10.36 GeV. Going across this energy point, the NLO correction changes from positive to negative. Generally, all curves go down with the increase of center-of-mass energy. It is worth noting that NRQCD factorization formalism for double charmonium only holds when √ s ≫ m c [33], and it may break down at the double charmonium threshold because of the color transfer effect [34], whereas in drawing Figs. 5-7 the center-of-mass energy has been extended to the threshold region, lower than the typical B-factory energy, just for a schematic display.   Since M ψ 1,2 > M DD , i.e. it is above the open charm threshold, and η c +ψ 1,2 (1D) cross sections are at the order of 1 fb, in experiment it turns out to be harder to detect ψ 1,2 than h c through double charmonium processes at super B-factories. For η c + h c (1P ) production, defining the relative production ratio R as and considering the Belle data of σ[e + e − → η c +J/ψ] = 33 fb [5], we have R = 0.010(0.009) at µ = 2m c ( √ s/2). If the BaBar data σ[e + e − → η c +J/ψ] = 17.6 fb [6] are adopted, we then have R = 0.019(0.017) at µ = 2m c ( √ s/2). If we adopt the NLO results of Ref. [7] for σ[e + e − → η c +J/ψ] = 18.9 fb, then R = 0.018 under the same renormalization scale µ = 2m c . In future super B-factories, this ratio R may stand as a benchmark for the estimation of the possibility of observing the h c (1P ) through double charmonium process.

Summary and Conclusions
In this work, we have studied the double charmonium processes e + e − → η c +h c (1P )/ψ 1,2 (1D) at NLO accuracy under the NRQCD factorization mechanism. Cross sections with varying charm quark mass m c = 1.7±0.2 GeV at typical renormalization scales (µ = 2m c , √ s/2) were analyzed in detail. The magnitudes of cross section versus energy scale µ and center-of-mass energy √ s at LO and NLO were evaluated. We have also estimated the relative production rate of R = σ[e + e − → η c +h c ]/σ[e + e − → η c +J/ψ], which might be helpful for the measurement of double charmonium exclusive production in future super B-factories. Through our study, we find that the NLO corrections for e + e − → η c +h c (1P )/ψ 1,2 (1D) are small, and the convergence of perturbative expansion works well up to NLO. When the scale µ lies in the region [2m c , √ s/2], the NLO correction for η c + h c (1P ) production is positive, yet it is negative for η c +ψ 1,2 (1D). The scale dependence has been clearly reduced when NLO corrections are taken into account. At √ s = 10.6 GeV, the relation σ ψ 1 < σ hc < σ ψ 2 exists at both LO and NLO, no matter whether µ = 2m c or µ = √ s/2. It is worth noting that the cross sections of the three processes considered are sensitive to the charm quark mass, which hence is the main source of uncertainty in our results.
Although h c (1P ) has been measured through the h c → J/ψ + π 0 → (e + e − ) + π 0 process, its branching fraction is not determined, and there has been no signal observed in the other h c primary decay mode h c → J/ψ+2π → (e + e − )+2π [19]. The results in this work may be helpful to h c and NRQCD studies in B-factories in future.