Is the colour-octet mechanism consistent with the double $J/\psi$ production measurement at B-factories?

Double $J/\psi$ production in $e^+e^-$ collisions involving colour-octet channels are evaluated up to order $\alpha^2\alpha_s^3$. Having implemented the variation of the parameters ($m_c$, $\mu_r$ and long-distance matrix elements), we found that the cross sections for producing double $J/\psi$ at B-factories range from $-0.016$fb to $0.245$fb, which are even much smaller than that via the colour-siglet mechanism. Accordingly, this result is consistent with the measurement by the Belle and BABAR Collaborations.


I. INTRODUCTION
The phenomenological study on the nonrelativistic QCD (NRQCD) effective theory [1] is making new progress since the LHC started its running. Copious data not only provides evidences for the colour-octet (CO) mechanism, but also indicates challenges to the theory.
In addition to the fact that the J/ψ hadroproduction data can be well reproduced by the theoretical evaluations within the NRQCD framework [2][3][4], χ c hadroproduction [5,6] gives another strong support. In low transverse momentum (p t ) region, even though the factorization might not hold, the colour-glass-condensate model [7][8][9] associated with NRQCD [10] did a good job in the description of the J/ψ production in proton-proton and proton-nucleus collisions [11,12]. Despite all the successes, we can not overlook the challenges it is facing.
The universality of the NRQCD long-distance matrix elements (LDMEs) has not yet been suggested in all the processes. As an example, the constraint [13] on the CO LDMEs indicated by the QCD next-to-leading order (NLO) study on the J/ψ production at B factories is apparently below the LDME values obtained through the fit of the J/ψ production data at other colliders [3,[14][15][16]. The perspectives of the long-standing J/ψ polarization puzzle still have not converged. Three groups [15,17,18] achieved the calculation of the J/ψ polarization at hadron colliders at QCD NLO; however, with different LDMEs, their results are complete different from one another. Recently, the η c hadroproduction was measured by the LHCb Collaboration [19], which provides another laboratory for the study of NRQCD.
Ref. [20] considers it as a challenge to NRQCD, while Refs. [21,22] found these data are consistent with the J/ψ hadroproduction measurements. Further, with the constraint on the LDMEs obtained in Ref. [22], Ref. [23] discovered some interesting features of the J/ψ polarization, and found a possibility of understanding the J/ψ polarization within the NRQCD framework.
The J/ψ pair production at B factories is another challenge that NRQCD is facing.
Belle [24] and BABAR [25] Collaborations observed the process e + e − → J/ψ+Charmonium, and found no evidence for the J/ψ pair events, while the QCD leading order (LO) calculation based on the colour-siglet (CS) mechanism predicted a significant production rate [26,27].
This was understood by the QCD NLO corrections [28], which contribute a negative value and cancel the large LO cross sections. Ref. [28] only talked about the CS contributions.
However, the Belle and BABAR measurements actually did not exclude the double J/ψ plus light hadron events. Both of the experiments measured the M res spectrum, where M res denotes the invariant mass of all the final states except for the fully reconstructed J/ψ. These distributions exhibited no significant excess in the range of about 300 MeV above the J/ψ mass, which suggested that the cross section for the J/ψ pair plus light hadron (e.g. π 0 , π + π − ) associated production is also too small to observe. To accord with NRQCD, the double J/ψ production cross sections involving the CO channels must not be significant, which, however, is not manifest. Although suppressed by the CO LDMEs, the double J/ψ yield due to the CO mechanism is enhanced by the powers of α s /α, relative to via the CS channels. As is pointed out in Refs. [26,28], at B factories, double cc( 3 S   J . These two processes are suppressed by the CO LDMEs by a factor of v 8 ≈ 0.002, where v is the typical charm-quark velocity in the charmonium rest frame, however, enhanced by the coupling constants by a factor of α 2 s /α 2 ≈ 1000, relative to via the CS channels. The double J/ψ can also be produced through such kind of processes, e + e − → cc(n 1 ) + cc(n 2 ) + g, where g denotes a gluon. When n 1 = 3 S  1 ), this kind of processes are enhanced by the coupling constants by a factor of α 3 s /α 2 ≈ 200 and reduced by the CO LDME by a factor of v 4 ≈ 0.05, relative to the processes involving only the CS channels. In sum, the processes involving CO states are enhanced by a synthetic factor of about 2 ∼ 10, comparing with the processes considered in Refs. [26,28], the LO cross sections of which is large enough to be observed by Belle and BABAR experiments. Accordingly, we need to calculate the cross sections for the J/ψ pair production involving the CO channels to see whether NRQCD can endure this paradox.
In this work, we will present a comprehensive study on the double J/ψ production in e + e − annihilation involving CO channels up to order O(α 2 α 3 s ), and check whether it is consistent with the meausurements by Belle and BABAR Collaborations. The J/ψ plus χ c production at B factories has already been studied in Refs. [29][30][31] J ) production. We also notice that the J/ψ may come from the χ cJ feed down, where the χ cJ can be produced via the 3 S [8] 1 channel. By employing the LDMEs obtained in Refs. [6,23], 1 ) = 1.08 × 10 −2 GeV 3 in association with the branching ratios listed above, we find that this contribution is much smaller than that from the J/ψ directly produced through 1 channel. Similarly, the J/ψ production cross sections via the ψ(2s) feed down is also smaller than that for the directly produced ones. For this reason, we completely omit the discussions on the feed down contributions from both χ c and ψ(2s).
The rest of this paper is organised as follows. In Sec.II, we outline the formalism of the calculation. Sec.III presents the numerical results and related discussions, followed by a concluding remark in Sec.IV.

II. DOUBLE J/ψ PRODUCTION IN NRQCD FRAMEWORK
Following the NRQCD factorization, the total cross sections for the J/ψ pair production can be expressed as J , and the representative Feynman diagrams are illustrated in Fig.(1a, 1b, and 1c). At QCD NLO (O(α 3 s α 2 )), in addition to the virtual corrections (the representative Feynman diagrams for which are shown in Fig.(1d, 1e, and 1f)) to the processes presented in Eq.(1), double cc states in association with a gluon production is also required for consideration, as illustrated in Eq. (2). The real-correction processes to the LO ones are in addition to which, five processes are also at this order, as listed below, and will be calculated in our paper.
e + e − → cc( 1 S The representatvie Feynman diagrams for the final-state-gluon-emission processes are presented in Fig.(1g, 1h, and 1i). However, not all the processes have the three types of diagrams. So, we summarize the possible diagrams for each of the processes in TableI.  The processes are abbreviated to the numbers of the equations, namely Eq. (4,5,6,7,8,9), in which the processes are presented.
Before we present the numerical results, we first address the divergences rising from the processes listed above. First of all, the LO processes are divergence free, and we denote their total cross sections as σ LO = σ LO 3 S corrections to σ LO can be expressed as where c-quark mass and the wave functions of the c-quark and gluon, and modified-mininumsubtraction (MS) scheme for that of the QCD coupling constant, which are coincide with Ref. [33].
where µ r is the renormalization scale, γ E is Euler's constant, β 0 = 11 3 C A − 4 3 T F n f is the one-loop coefficient of the QCD beta function, n f is the number of active quark flavors. The cross sections for the processes listed in Eq.(8 and 9) also have divergences, which, however, can be eliminated through the renormalization of the SDCs for them. We take process 8 as an example. The cancellation of its divergences requires the calculation of the 1 ) . The bare LDME can be expressed as where m c is the c-quark mass, and N c = 3 for SU(3) gauge theory. Here we adopt the µ Λ -cutoff renormalization scheme [6] to subtract the UV divergence. By substituting the relation between the bare and renormalized LDMEs, into Eq. (14), we obtain the renormalized LDME as Then process e + e − → cc( 3 S 0 ) contribute an additional divergent term which cancels the IR sigularities rising from process 8. In this sense, we can redefine the SDC for process 8 aŝ 0 )), (18) whereσ (e + e − → cc( 1 S 2 ) + g) (19) has been implicated in Eq.(18) (the same convention applies to the SDCs with the subscript "renorm").σ renorm is a finite quantity, therefore, we can replace, in Eq.(3), the divergent one by it. The same operation can be done for process 9 as well. Then, we denote all the divergence-free total cross sections for the processes listed in Eq. (5,6,7,8,9) as σ n 1 +n 2 , where n 1 and n 2 are the corresponding cc states.

III. NUMERICAL RESULTS
In our analytic calculation, we use our Mathematica package with the employment of FeynArts [34], FeynCalc [35], FIRE [36] and Apart [37]. As a cross check, we also compute the processes using the FDC package [38], except for process 9. To subtract the IR divergences in the gluon-emission processes, we adopt the two-cutoff slicing strategy [39]. The independence of the cutoff has been checked.  Table II, where the SDCs for the 3 P [8] J channels are defined by multiplying a factor of m 2 c to those defined in Ref. [1], in order to keep the homogeneity of the dimensions (for double 3 P J productions are 60.07 fb/ GeV 3 and 290.9 fb/ GeV 3 , respectively.  Employing the LDMEs obtained in Refs. [22,23], namely we list the cross sections for each channel in Table III. The total cross section for double J/ψ production at B factories up to order O(α 2 α 3 s ) is the sum of those for different channels. Note that the two processes which are symmetric in the sense of switching n 1 and n 2 are only counted once to avoid the double counting. We obtain this value as σ = 0.1046 GeV. Compared with the results obtained in Ref. [28], it is even smaller than that via the CS channels up to QCD NLO.    [22,23].
To investigate the uncertainties brought in by the two scales, we vary m c from 1.2 GeV to 1.7 GeV and µ r from 3.0 GeV to √ s/2 and calculate the corresponding total cross sections.
When one of these scales varies its value, the other is fixed. Note that the LDMEs used in our calculation are obtained with the configuration m c = 1.5 GeV. When investigating the m c dependence, we need to take the scaling rule, O J/ψ (n) ∝ m 3 c , into account. The total cross section as a function of charm quark mass m c is presented in Fig.2 J + g cross sections are negative . The total cross section increases from about 0.08 fb to about 0.12 fb as the m c increases from 1.2 GeV to 1.7 GeV. The µ r dependence of the total cross section is presented in Fig.3, and the σ decreases from 0.104 fb to 0.08 fb as µ r increases from 3.0 GeV to √ s/2. The dependence on the two scales is not severe, which indicates good convergence of the perturbative expansion.
Since there are several parallel extractions of the LDMEs, we need to investigate the uncertainties brought in by the different values of them. As is shown in TABLEIV, the total cross sections obtained by using the LDMEs in Ref. [16] are almost twice of ours, however, still too small to be observed by the experiment.

IV. SUMMARY AND CONCLUSION
We calculated the total cross sections for double J/ψ production in e + e − annihilation at the B-factory energy up to O(α 2 α 3 s ) within the framework of NRQCD. We studied the m c and µ r dependence of the total cross sections, and found that the results ranges from 0.08 fb to 0.12 fb. Also, we investigated the uncertainties by trying different set of the LDMEs.
Even for the largest results, the total cross section is too small for Belle to observe any significant access. This result is consistent with the Belle measurement.  GeV 3 in the calculation.