${\mathcal O}(\alpha_s v^2)$ correction to $J/\psi$ plus $\eta_c$ production in $e^+e^-$ annihilation at $\sqrt{s}=10.6GeV

Based on the nonrelativistic QCD factorization approach, ${\mathcal O}(\alpha_s v^2)$ corrections to $\jpsi$ plus $\eta_c$ production in $e^+e^-$ annihilation at $\sqrt{s}=10.6 \gev$ is calculated in this work. The numerical results show that the correction at $\alpha_s v^2$ order is only about a few percent for the total theoretical result. It indicates that the perturbative expansions for the theoretical prediction become convergence and higher order correction will be smaller. The uncertainties from the long-distance matrix elements, renormalization scale and the measurement in experiment are also discussed. Our result is in agreement with previous result in ref [1].


I. INTRODUCTION
Study on heavy quarkonium decay and production is a very important and interesting issue to understand quantum chromodynamics (QCD), the fundamental theory of strong interactions.Many experimental and theoretical reserches have been performed since the discovery of the J/ψ charmonium meson in 1974 followed by the Υ bottomonium meson in 1977, for reviews see Ref. [2].In experimental side, it is easy to detect J/ψ and Υ signal.In theoretical side, quarkonium bound states offer a solid ground to probe QCD, due to the high scale provided by the large mass of the heavy quarks, which make the QCD factorization possible in the related calculation.To explain the large discrepancy on the transverse momentum distribution of chromonium hadroproduction between the experimental measurement and theoretical prediction as well as to arrange the infrared divergence cancellation in p-wave quarkonium related calculation, the nonrelativistic QCD (NRQCD) factorization approach [3] has been introduced.It allows consistent theoretical prediction to be made and to be improved perturbatively in the QCD coupling constant α s and the heavy-quark relative velocity v in heavy quarkonium rest frame.
In last five years, most of the important theoretical studies on heavy quarkonium based on NRQCD are calculated at next-to-leading order (NLO) of QCD and many of them are also calculated at next-to-leading order of v.Among them, the J/ψ polarization puzzle at hadron colliders is still unclear after the important progresses at QCD NLO [4], It seems that the inclusive J/ψ production at B-factories can be explained by just color singlet contribution at QCD NLO [5,6], but it causes the problem for the color-octet long distance matrix elements [7].The theoretical calculation with NLO QCD and relativistic correction can cover the experimental measurements on exclusive double chromonium production at B-factories although the corrections are very large.For theoretical prediction based on perturbative expansion, the convergence of the expansion is a very important issue.Therefore, it is important to test the calculation at higher order when the NLO correction is large.Usually, higher order calculation is much more complicate, so far there are only a few simple processes whose O(α s v 2 ) corrections are calculated [8][9][10].
For the exclusive double chromonium production at B-factories, its higher order calculation is studied, so we give a detailed review on it.The exclusive production cross section of double charmonium in e + e − → J/ψη c at √ s = 10.6 GeV measured by Belle [11,12] is ) fb, where B ηc [≥ 2] denotes the branching fraction for the η c decaying into at least two charged tracks.Meanwhile, the NRQCD LO theoretical predictions in the QCD coupling constant α s and the charm-quark relative velocity v, given by Braaten and Lee [14], Liu, He and Chao [15], and Hagiwara, Kou and Qiao [16] are about 2.3 ∼ 5.5 fb, which is an order of magnitude smaller than the experimental results.Such a large discrepancy between experimental results and theoretical predictions brings a challenge to the current understanding of charmonium production based on NRQCD.Many studies have been performed in order to resolve the problem.From treatments beyond NRQCD, Ma and Si [17] treated the process by using light-cone method, a similar treatment was performed by Bondar and Chernyad [18] and Bodwin, Kang and Lee [19], possible contribution from intermediate meson rescatterings was considered by Zhang, Zhao, and Qiao [20], it was also studied in the Bethe-Salpeter formalism by Guo, Ke, Li, and Wu in Ref [21].Based on NRQCD, Braaten and Lee [14] have shown that the relativistic corrections would increase the cross section by a factor of about 2, and the NLO QCD correction of the process has been studied by Zhang, Gao and Chao [22] and Gong and Wang [23], which can enhance the cross section with a K factor (the ratio of NLO to LO) of about 2, again the relativistic corrections have been studied by Bodwin, Kang, Kim, Lee and Yu [24] and by He, Fan and Chao [5], which is significant.More detailed treatment, such as including the resummation of a class of relativistic correction, has been taken into consideration by Bodwin and Lee and Yu [25].In another way, Bodwin, Lee and Braaten [26] showed that the cross section for the process e + e − → J/ψ + J/ψ may be larger than that for J/ψ + η c by a factor of 1.8, in spite of a suppression factor α 2 /α 2 s that is associated with the QED and QCD coupling constants.They suggested that a significant part of the discrepancy of J/ψ + η c production may be explained by this process.Hagiwara, Kou and Qiao [16] also calculated and discussed this process.In 2004, a new analysis of double charmonium production in e + e − annihilation was performed by Belle [27] based on a 3 times larger data set and no evidence for the process e + e − → J/ψ + J/ψ was found.Both the NLO QCD corrections and relativistic corrections to e + e − → J/ψ + η c give a large K factor of about 2. It is obvious that these two types of corrections to e + e − → J/ψ + J/ψ should be studied to explain the experimental results.In fact, they have been studied by Bodwin, Lee and Braaten for the dominant photon-fragmentation contribution diagrams [28].The results show that the cross section is decreased by K factor of 0.39 and 0.78 for the NLO QCD and relativistic corrections respectively.A more reliable estimate, 1.69 ± 0.35 fb, was given by Bodwin, Lee, Braaten and Yu in ref. [29].And light-cone method is used in ref. [30] by V.V. Braguta.Gong and Wang performed a complete NLO QCD calculation on e + e − → J/ψ + J/ψ [31] and the results show that the cross section would be much smaller than the rough estimate in Ref. [28].Therefore it is easy to understand why there was no evidence for the process e + e − → J/ψ + J/ψ at B-factories.
It is easy to see that both the QCD correction (α s ) and relativistic correction (v 2 ) are very large for e + e − → J/ψη c at B-factory energy, and with these corrections the experimental measurement can be explained.Therefore it is natural to ask the question, how is the situation for the higher order corrections beyond α s and v 2 correction ?α 2 s correction is very difficult to do, but recent progress make it available to do α s v 2 correction already.It is very interesting to see that the α s v 2 correction, given in a recent work [1], is a small contribution.It convinces us in some sense (with α 2 s correction absent) that the double expansions in NRQCD converges quite well on this problem.Since the calculation is quite complicate and plays an important role to convince us the convergence on the theoretical predication which can explain the experimental data, in this paper we performed an independent calculation on it by using the our package Feynman Diagram Calculation (FDC) [32] with the built-in method to calculate relativistic correction.The remainder of this paper is organized as follows.Base on the NRQCD frame, we briefly introduce theoretical formulism for the calculation of heavy quarkonium production and give the corresponding results in perturbative NRQCD in Sec.II.The details in perturbative QCD are summarized in Sec.III.We give the numerical results of α s v 2 corrections and some discussion in Sec.IV.Finally, in Sec.V, we present a brief summary.

II. NRQCD FACTORIZATION FORMULA UP TO v 2 ORDER
According to NRQCD effective theory, the production of the charmonium are factorized into two parts, the shortdistance part and the long-distance part.The long-distance parts are related to the four fermion operators, characterized by the velocity v of the charm quark in the meson rest frame.The long-distance matrix elements can be estimated by lattice calculations or phenomenological models, or determined by fitting experimental data.The production cross section up to v 2 order is expressed as with the long-distance matrix elements being defined by using related operators as for J/ψ and for η c where m c is the charm quark mass.It is the basic point that the NRQCD factorization for hadron related process will also hold when the hadron state are replaced by Q Q states with exactly the same quantum numbers as the corresponding hadron state.In this way, the short-distance coefficients c 00 , c 01 and c 10 can be obtained in perturbative calculation through the matching condition, and they are calculated up to QCD next-to-leading (NLO) order.In order to obtain the short-distant coefficients, the matrix elements of the operators for quantum states need to be calculated perturbatively, and there are where there are N c = 3 for SU (3) group and E q = m 2 c + q 2 .From the NRQCD effective Lagrangian, we could easily get the Feynman rules.Therefore we have calculated order α s v 2 corrections to the leading order O 1 2s+1 Ss in perturbative NRQCD with the dimensional regularization and defined the renormalization constants Z MS O of the operator by using the MS scheme [3,33].
At last we could easily give the perturbative NRQCD results. σ(e

III. DETAILS OF PERTURBATIVE QCD CALCULATION
For a Q(p) Q(p) quantum state, we denote P as the total momentum and q as the relative momentum between Q and Q pair.Therefore, there are where m Q is the mass of the heavy quark Q, and the Q and Q are on their mass shells.
To do the perturbative calculation in related process for the quantum states, we could obtain the projectors for each quantum states.The spin-singlet and spin-triplet components of each Q Q state can be projected out by making use of the spin projectors.After multiplying corresponding Clebsch-Gordan coefficients to the spin component of the outer product of the spinors for each Q Q pair, one can find that Π1 and Π3 are the spin-singlet and spin-triplet projectors of the Q Q production, respectively.The spin projectors that are valid to all orders in the relative momentum can be found in Refs [34].
where Π 1 and Π 3 are projectors for spin 0 and spin 1 s-wave quantum states respectively, and ǫ * (λ) is the polarization vector of the spin-triplet state.
For process e + (p 1 )e − (p ), the production matrix element is expressed as where we have used the following relation As for the expansion of q, we should consider the effect that the external momentum and polarization vector may be the implicit function of q.From the momentum conservation and on-shell conditions, p 2 3 = 4E 2 q1 , p 2 4 = 4E 2 q2 , we could find that p 3 , p 4 are implicit functions of q 1 , q 2 respectively.However,it is obvious that the short-distance coefficients, to be obtained in the perturbative calculation, are functions of the independent variables which are the invariant mass s of the e + and e − system and cos θ. θ is the angle between J/ψ and the electron.Where s and cosθ are independent of the relative momentum q.
Since the final results are Lorentz invariance and irrelevant to the reference frame, we choose to do the calculation in the center-of-mass of this system where p 1 + p 2 = p 3 + p 4 = ( √ s, 0, 0, 0) is the explicit expression of the momentum conservation.Therefore the following results are obtained: We choose two vectors r 1 = (0, − → r 1 ) and r 2 = (0, − → r 2 ) with − → r 1 and − → r 2 being unit vectors, while − → r 1 , − → r 2 and − → p 3 are perpendicular to each other, i.e r 1 • r 2 = 0, p 3 • r 1 = 0, p 3 • r 2 = 0. Then vector dp 3 dq 2 1 can be expressed as linear combination of four independent vectors as From the following conditions together with previous conditions in Eq.( 13), we can easily obtain the solution For the ǫ * (λ), the polarization four-vector of the |Q Q( 3 S 1 ) with helicity λ, there are the relation dǫ * (±1) dq 2 1 = 0 since θ is independent of the relative momentum q.It is easy to obtain dǫ * (0) dq 2 1 Therefore, we obtain the relation between the polarization four-vector and q as dǫ * (λ) The treatment about q 2 is similar to these.Considering the two body phase space, we need to expand it.
need to be expanded since cos θ and s are independent of q 1 and q 2 .Then there is where . We could square the amplitude, integrate over the phase space, and expand in powers of q in order to obtain the desired perturbative result Most of the steps in this section are realized in a small program in FDC package, and the final Fortran source for numerical calculation are prepared by using FDC package together with the small program for q 2 expansion.
Since there is no O(α s ) real process in NLO, we only need to calculate virtual corrections.Dimensional regularization has been adopted for isolating the ultraviolet(UV) and infrared(IR) singularities.UV divergences are cancelled upon the renormalization of the QCD gauge coupling constant, the charm quark mass and field, and the gluon field.A similar renormalization scheme is chosen as in ref. [35] except that both light quarks and charm quark are included in the quark loop to obtain the renormalization constants.The renormalization constants of the charm quark mass Z m and field Z 2 , and the gluon field Z 3 are defined in the on-mass-shell(OS) scheme while that of the QCD gauge coupling Z g is defined in the modified-minimal-subtraction(MS) scheme: where γ E is Euler's constant, OS cancel each other, thus the result is independent of renormalization scheme of the gluon field.

IV. RESULTS
The final results are obtained by using the matching method with the UV and IR divergences being cancelled.
σ LO , σ N LO(αs) , σ N LO(v 2 ) , σ N LO(αv 2 ) are the contributions from the leading order, the next leading order in α s , the next leading in v 2 and the next leading in αv 2 .Then the production rate up to O(α s v 2 ) order is expressed as where there are e c = 4 3 , r = and u r is the renormalization scale.Therefore, the obtained analytic expressions of the v 2 correction is in agreement with that in the paper [1].At the same time, the analytic expression of f 3 (r) in the results of the α s correction is also in agreement with that in the paper [23].Since the analytic expressions of f 4 (r) and f 5 (r) in the O(α s v 2 ) correction are so lengthy that we just give the numerical results for them.In the numerical calculation, there are f 1 = 0.97466, f 2 = 1.3080, f 3 = 12.358, f 4 = 3.8382, f 5 = 3.2537, f or r = 4×1.4 2 10.58 2 ; f 1 = 0.87465, f 2 = 1.2098, f 3 = 11.806,f 4 = 2.0543, f 5 = 2.6668, f or r = 4×1.5 2 10.58 2 .
And we take √ s = 10.58GeV and µ Λ = m c .The running strong coupling constant is evaluated by using the two-loop formula with Λ (4) MS = 0.338 GeV as used in Ref [23].Our results are presented in the table I with parameters given in table caption.The results are in agreement with that in ref [1].The contribution from the O(α s v 2 ) order  is small.There are some differences for the results in the table II if the long-distance matrices and QED coupling constant are chosen as in Ref [5], the correction at O(α s v 2 ) order is also small.we also give the relation of the cross sections and the µ r in the FIG. 1.There are about 10 percents differences in the total cross sections between m c = 1.5 and m c = 1.4.We find the uncertainty of the total cross sections from m c is not small.If we choose µ r = 2m c , we could give the relations of the ratios of different parts to total cross sections with the √ s in the FIG. 2. We will find that the contributions from O(α s ) and O(α s v 2 ) become important and the one from LO becomes small when √ s is large, but there are just about 10 percents contribution in O(α s v 2 ).The contribution from the O(α s v 2 ) order is small once again.

V. SUMMARY
In this work we have calculated the O(α s v 2 ) correction in detail for the processes e + e − → J/ψ + η c within the frame of NRQCD.The result at O(α s v 2 ) order give about 6 percent contribution to the total theoretical prediction while the O(α s ) correction and O(v 2 ) are about 40 percent and 14 percent contribution respectively.It indicates that the convergence in the double perturbative expansions in QCD α s and relativistic v 2 are very well for the theoretical calculation on the production rate of the process e + e − → J/ψ + η c .Up to O(α s v 2 ) order, the theoretical prediction with quite large uncertainty from charm quark mass and renormalization scale can describe the experimental measurement.
the one-loop coefficient of the QCD beta function and n f is the number of active quark flavors.There are three massless light quarks u, d, s, and one heavy quark c, so n f =4.In SU (3) c , color factors are given byT F = 1 2 , C F = 4 3 , C A = 3.And β ′ 0 ≡ β 0 + (4/3)T F = (11/3)C A − (4/3)T F n lf where n lf ≡ n f − 1 = 3is the number of light quarks flavors.Actually in the NLO total amplitude level, the terms proportion to δZ 3

FIG. 1 :
FIG. 1: The cross section as a function of the µr at √ s = 10.58GeV.The black and blue solid curves are the cross sections in the mc = 1.4 and mc = 1.5 respectively.The red and green bands represent the measured cross sections by the Belle and BaBar experiments, with respective systematic and statistical errors.