Next-to-leading-order study of $J/\psi$ angular distributions in $e^{+}e^{-} \to J/\psi+\eta_c,\chi_{cJ}$ at $\sqrt{s} \approx 10.6$ GeV

In this paper, we present a detailed next-to-leading-order (NLO) study of $J/\psi$ angular distributions in $e^{+}e^{-} \to J/\psi+\eta_c,\chi_{cJ}$ ($J=0,1,2$) within the nonrelativistic QCD factorization (NRQCD). The numerical NLO expressions for total and differential cross sections, i.e., $\frac{d\sigma}{d\cos\theta}=A+B\cos^2\theta$, are both derived. With the inclusion of the newly-calculated QCD corrections to $A$ and $B$, the $\alpha_{\theta}(= B/A)$ parameters in $J/\psi+\chi_{c0}$ and $J/\psi+\chi_{c1}$ are moderately enhanced, while the magnitude of ${\alpha_\theta}_{J/\psi+\chi_{c2}}$ is significantly reduced; regarding the production of $J/\psi+\eta_c$, the $\alpha_\theta$ value remains unchanged. By comparing with experiment, we find the predicted ${\alpha_\theta}_{J/\psi+\eta_c}$ is in good agreement with the $\textrm{B}\scriptsize{\textrm{ELLE}}$ measurement; however, ${\alpha_\theta}_{J/\psi+\chi_{c0}}$ is still totally incompatible with the experimental result, and this discrepancy seems to hardly be cured by proper choices of the charm-quark mass, the renormalization scale, and the NRQCD matrix elements.


Introduction
In the past twenty years, the BELLE and BABAR Collaborations have independently measured the total cross sections of e + e − → J/ψ + η c , χ c0 at √ s ≈ 10.6 GeV [1][2][3], which significantly overshoot the results [4][5][6][7] calculated at leading order (LO) in α s using the nonrelativistic QCD (NRQCD) factorization [8]. As a breakthrough of the theoretical attempts  to explain this inconsistency, the next-to-leading-order (NLO) QCD corrections [13,16,17] can provide considerable and positive contributions, largely alleviating the discrepancies between theory and experiment.
Besides the total cross sections, based on a data sample of 140 fb −1 , BELLE has also measured the J/ψ angular distributions in e + e − → J/ψ + η c , χ c0 [2], i.e., the α θ parameter in dσ d cos θ = A + B cos 2 θ = A(1 + α θ cos 2 θ). By simultaneously fitting the production-and helicity-angle distributions, the measured α θ reads α θ (e + e − → J/ψ + η c ) = 0.93 +0.57 −0.47 , α θ (e + e − → J/ψ + χ c0 ) = −1.01 +0. 38 −0.33 . (1.1) On the theoretical side, the existing studies of the differential cross section ( dσ d cos θ ) are only accurate to the LO level in α s [4,6]. For J/ψ +η c , the LO calculations using NRQCD give a prediction α θ = 1, which is consistent with the measured value in equation (1.1); however, the NRQCD-based LO predictions of α θ(J/ψ+χ c0 ) 0.25 are fundamentally different from the above BELLE measurement. In order to fill the huge gap between theory and experimental result, spurred by the significant impacts of the QCD corrections on total cross section, it is urgent to carry out a NLO analysis to the differential cross section. Moreover, whether the inclusion of the uncalculated QCD corrections would spoil the existing coincidence of α θJ/ψ+η c with experiment need also to be verified.
At this stage, BELLE and BABAR have not observed any evident event of e + e − → J/ψ + χ c1,2 , owing to the relatively small production rates of the two channels. Fortunately, the recently-commissioning Super-B factories with a high luminosity designed to reach up to about 50 ab −1 by 2022 would bring great opportunities to fulfill the observations, which can aid in further understanding the double-charmonium productions in e + e − annihilation. Taken together, in this paper we will for the first time perform a comprehensive NLO study of the differential cross sections in e + e − → J/ψ + η c , χ cJ with J = 0, 1, 2.
Note that, in the context of NRQCD factorization, due to the inadequate knowledge about the heavy-quarkonium production mechanism, the calculated σ e + e − →J/ψ+ηc,χ cJ suffer severely from the indeterminacy inherent to the nonperturbative NRQCD long distance matrix elements (LDMEs), which would then significantly weaken the predictive power of NRQCD. On the contrary, as a result of the cancellation of the dependences of A and B on LDME, the theoretical result of α θ (= B/A) dose NOT involve this kind of ambiguity. In this sense, the α θ parameter is expected to be a more ideal laboratory than total cross section for the study of exclusive double-charmonium productions in e + e − annihilation. Considering the large uncertainty of α θJ/ψ+η c,χc0 in equation (1.1) and the lack of measured α θJ/ψ+χ c1,2 , we strongly suggest the Super-B factories to reperform with better precision the measurements of α θ , and our state-of-the-art calculation results would pave the way for the future comparisons.
The rest of the paper is organized as follows: Section 2 is an outline of the calculation formalism. Then, the phenomenological results and discussions are presented in Section 3. Section 4 is reserved as a summary.

Leptonic current
L µν can directly be calculated and obtained as where q = p 1 + p 2 and s = (p 1 + p 2 ) 2 . Note that, in calculating total cross section, the standard l µν in equation (2.3) can be replaced equivalently with − 8 3 g µν [34], which, however, is no longer applicable for the case of differential cross section. Therefore, in order to evaluate dσ d cos θ , we will employ a new l µν obtained in our recent paper [35]; the following is a brief description of its derivation.
By integrating H µν over all the final states other than J/ψ, we obtain the hadronic tensor W µν h (p 3 , q), which is dependent only on p 3 and q. Subsequently we decompose W µν h (p 3 , q) as a linear combination of the tensors constituted of −g µν , p 3 , and q as which satisfies the following relation where p 3 is the three momenta of J/ψ, and θ is the angle between p 3 and the spatial momentum of e − (or e + ) in the e + e − center-of-mass frame.
It is easy to verify that, with where A 1 = 1 + cos 2 θ, one can reproduce the results in equation (2.6). Further employing the current conservation will finally yield With the integration over cos θ, the term involving A 2 in equation (2.9) vanishes, and the first term on the right-hand side will reduce to − 8 3 g µν . Comparing the two leptonic tensors in equations (2.3) and (2.9), one can easily find the newly-derived l µν is related only to the hadronic-process momentum (p 3 ), which, by the absence of p 1 and p 2 , would greatly reduce the complication in computing the differential cross sections, especially at the NLO accuracy.

Cross sections
Using the leptonic tensor in equation (2.9), the differential cross sections of e + e − → J/ψ + η c , χ cJ can be written in the following form and then we have The universal factor κ reads (2.12) |R 1S (0)| and |R 1P (0)| are the wave functions at the origin, which can be related to the NRQCD LDMEs by the following formulae: According to the LO Feynman diagrams for e + e − → J/ψ + η c , χ cJ , which are free of divergence and which are representatively shown in figure 1(a), we figure straightforwardly out the coefficients C 1 and C 2 through LO order in α s for various processes, 14) where the dimensionless variable r is defined as Following the relations in equation (2.11), one can directly obtain the analytical expressions of the LO-level A, B and σ concerning e + e − → J/ψ + η c , χ cJ , which can be proved to agree with the results in ref. [6].

NLO
Due to the absence of the real-correction processes in NLO, we need only to calculate the virtual corrections, which include 60 one-loop diagrams and 20 counter-term diagrams, as illustrated in figures 1(b)-1(f). We utilize the dimensional regularization with D = 4 − 2 to isolate the ultraviolet (UV) and infrared (IR) divergences. The on-mass-shell (OS) scheme is employed to set the renormalization constants for the c-quark mass (Z m ) and heavy-quark filed (Z 2 ); the minimal-subtraction (M S) scheme is adopted for the QCD-gauge coupling (Z g ) and the gluon filed Z 3 . The renormalization constants are taken as where γ E is the Euler's constant, N = Γ[1 − ]/(4πµ 2 r /(4m 2 c )) is an overall factor in our calculation, and β 0 = 11 3 C A − 4 3 T F n f is the one-loop coefficient of the β function. n f (= n L + n H ) represents the number of the active-quark flavors; n L and n H denote the number of the light-and heavy-quark flavors, respectively. At √ s = 10.6 GeV, n L = 3 and n H = 1. In SU(3), the color factors are given by T F = 1 2 , C F = 4 3 , and C A = 3. By taking into account the QCD corrections, we acquire the NLO-level C 1 and C 2 , and then A and B of NLO in α s using the relations in equation (2.11). Since the fully analytical expressions of the NLO results are lengthy, we in the following just provide the numerical expressions, among which the choice of charm-quark mass m c = 1.5(±0.1) GeV that is often used in calculating the charmonium-involved processes is chosen.
The NLO expressions of A, B, or σ can be generally written as where the dimensionless coefficients c L , c H , and c for various processes are summarized in tables 1, 2, 3, and 4.   Table 3: The coefficients c L , c H , and c for e + e − → J/ψ + χ c1 at √ s = 10.6 GeV. The unit of m c is GeV.  Table 4: The coefficients c L , c H , and c for e + e − → J/ψ + χ c2 at √ s = 10.6 GeV. The unit of m c is GeV. In our calculations, we use FeynArts [36] to generate all the involved Feynman diagrams and the corresponding analytical amplitudes. Then the package FeynCalc [37] is applied to tackle the traces of the γ and color matrixes such that the hard scattering amplitudes are transformed into expressions with loop integrals. In the next step, we utilize our self-written Mathematica codes with the implementations of Apart [38] and FIRE [39] to reduce these loop integrals to a set of irreducible Master Integrals, which could be numerically evaluated by using the package LoopTools [40].
To check the correctness of our calculations, we simultaneously apply another independent package, Feynman Diagram Calculation (FDC) [41], to perform the QCD corrections and acquire the same numerical results. As another crucial cross check, with the employment of our numerical expressions, one can immediately reproduce the same values of the K(= σ NLO /σ LO ) factors as in refs. [13,16,21,22].

Phenomenological results
In order to do the numerical calculations, we choose m c = 1.5 ± 0.1 GeV, M J/ψ = M ηc = M χ cJ = 2m c , α = 1/137, and employ the two-loop α s running coupling constant. The center-of-mass collision energy is √ s = 10.6 GeV. By using the NRQCD-based results of Γ J/ψ→e + e − and Γ χ c2 →γγ of NLO in α s to respectively match the latest measurements [42], we obtain |R 1S (0)| 2 = 0.941 GeV 3 and |R 1P (0)| 2 = 0.067 GeV 5 . Table 5: Total and differential cross sections at √ s = 10.6 GeV. A and α θ are the coefficients in dσ d cos θ = A(1 + α θ cos 2 θ). m c = 1.5 GeV and µ r = 2m c . The NRQCD predictions for the total and differential cross sections are summarized in tables 5 and 6. Inspecting the two tables, on can observe the NLO corrections are crucial for the total cross sections of e + e − → J/ψ + η c , χ c0 , while moderate in the case of J/ψ  plus χ c1,2 . This can be understood by analyzing the NLO expressions in equation (2.18). Taking m c = 1.5 GeV for instance, the coefficient of c in σ NLO ηc(χ c0 ) , i.e., 12.728(8.2307), would bring forth significant enhancements to the LO results; however, as to J/ψ + χ c1 (χ c2 ), this coefficient is only 0.5627(-0.4205).
From the data in tables 5 and 6, it is apparent that the LO predictions of α θ are independent on the choice of the renormalization scale µ r . With the inclusion of the QCD corrections, α θJ/ψ+η c remains unchanged; however, by varying µ r in [2m c , √ s/2], α θJ/ψ+χ c0 and α θJ/ψ+χ c1 are enhanced by about 9 − 12% and 14 − 23%, respectively, and α θJ/ψ+χ c2 is reduced in magnitude by about 43 − 80%. These enhancing or reducing effects on α θ are also clearly visualized by figures 2 and 3.
It is worth noting that the impacts of the QCD corrections on α θ differ substantially from that in the case of total cross section. From equation (2.18)   BELLE result; however, the increase of α θ in J/ψ + χ c0 stemming from the incorporation of the QCD corrections further intensifies the disturbing discrepancy in existence between theory and experiment. α θJ/ψ+χ c0 and α θJ/ψ+χ c1 exhibit strong stability under the m c variation, while the deviation in m c from 1.5 GeV by ±0.1 GeV would bring about a 30% fluctuation of α θJ/ψ+χ c2 .
In consideration of the noticeable disagreement between the theoretical result of α θJ/ψ+χ c0 and BELLE data, which can hardly be remedied by properly choosing the values of m c , µ r , and NRQCD LDME, it seems to be premature to draw a decisive conclusion concerning the experimental verification of the NRQCD description of e + e − → J/ψ + χ c0 from the coincidence between total cross section and experiment. To shed light on this issue, it would be desirable to extend the existing measurements to include e + e − → J/ψ + χ c1 and e + e − → J/ψ + χ c2 , especially at the Super-B factories which are designed to run with a luminosity up to ∼ 10 36 cm −2 s −1 .

Summary
In order to provide deeper insight into the well-known exclusive double-charmonium productions in e + e − annihilation at B factories, we in this paper carry out the first NLO study of the J/ψ angular distributions in e + e − → J/ψ + η c , χ cJ with J = 0, 1, 2 based on the NRQCD factorization. The numerical results show that the newly-calculated QCD corrections moderately affect the LO predictions of α θJ/ψ+χ c0 and α θJ/ψ+χ c1 , while significantly reduce in magnitude the LO result of α θJ/ψ+χ c2 . Concerning the production of J/ψ + η c , the QCD corrections do not change its α θ value. By comparisons to experiment, we find α θJ/ψ+η c is consistent with the BELLE measurement; however, the existing radical incompatibility between the LO prediction of α θJ/ψ+χ c0 and experiment becomes even worse by including the QCD corrections.