${\mathcal O}(\alpha_s)$ corrections to $e^+e^-\rightarrow J/\psi+\eta_{c2}(\chi_{c1}^\prime)$ at $B$ factories

We investigate the ${\cal O}(\alpha_s)$ correction to $e^+e^-\to J/\psi+\eta_{c2}$ in the NRQCD factorization approach. A detailed comparative study between $e^+e^-\to J/\psi+\eta_{c2}$ and $e^+e^-\to J/\psi+\chi^\prime_{c1}$ at $B$ factory energy is also carried out. After incorporating the ${\cal O}(\alpha_s)$ correction, we predict the cross section for the former process to be around 0.3 fb, while that of the latter about 6 times greater. The outgoing $J/\psi$ is found to be dominantly transversely-polarized in the former process, while longitudinally-polarized in the latter. These features may provide valuable guidance for the future experiment to examine the ${}^3P_1$ or ${}^1D_2$ charmonium option of the X(3872) meson through the exclusive double-charmonium production processes. The observation potential of $e^+e^-\to J/\psi+\chi^\prime_{c1}$ looks bright for the current data sample of the \textsc{Belle} experiment, provided that the $\chi^\prime_{c1}$ is indeed the narrow X(3872) state. In the appendix, we also identify the coefficients of the double logarithms of form $\ln^2(s/m_c^2)$ associated with all the relevant next-to-leading order Feynman diagrams, for the helicity-suppressed double-charmonium production channels $e^+ e^- \to J/\psi+ \eta_{c2}(\eta_c,\chi_{c0,1,2})$.


I. INTRODUCTION
Since a number of double charmonium production processes were first discovered at B factory a decade ago [1][2][3], this topic has spurred a widespread interest [4]. Firstly, this production environment provides a unique and powerful means to search for the new Ceven charmonium states, especially those X, Y , Z states, by fitting the recoil mass spectrum against the J/ψ (ψ ′ ). The most famous examples are the discovery of the X(3940) [5] and the X(4160) [6] by this way.
Another appealing reason to study double charmonium production is that it provides a new stage to sharpen our understanding toward perturbative QCD, especially toward the application of the light-cone approach [7,8] and the nonrelativistic QCD (NRQCD) factorization approach [9] to hard exclusive reactions involving heavy quarkonium.
The most famous double-charmonium production process is perhaps e + e − → J/ψ + η c . The original lowest order (LO) NRQCD predictions to this process [10,11] is about one order of magnitude smaller than the measurement [1]. This acute discrepancy has triggered a great amount of theoretical investigations in both NRQCD and light-cone approaches [12][13][14][15][16][17][18][19][20][21][22]. One crucial element in alleviating the discrepancy between the NRQCD prediction and the data is the significant and positive next-to-leading order (NLO) perturbative corrections [16,17]. By contrast, owing to some long-standing theoretical obstacles, the NLO correction to this helicity-suppressed process in the light-cone approach has never been successfully worked out. As a consequence, despite some shortcomings, the NRQCD approach seems to be the only viable method which is based on the first principle and also systematically improvable.
Recently, the O(α s ) corrections to the double-charmonium production processes e + e − → J/ψ + χ c0,1,2 have also been investigated in NRQCD approach [23,24]. For some production channels, the effect of the NLO perturbative corrections can be important.
In this work, we plan to carry out a comprehensive study of the NLO perturbative corrections to the processes e + e − → J/ψ + η c2 and e + e − → J/ψ + χ ′ c1 , again in the NRQCD factorization framework. The calculation for the former process is new, while that for the latter can be readily adapted from Ref. [24]. This work should be considered as a sequel of Refs. [21,24].
Our interest in conducting such a comparative study is strongly motivated by the controversy recently raised at the quantum number of the X(3872) meson, whether it being 1 ++ or 2 −+ [25]. Consequently, the canonical charmonium options for the X(3872) would be the χ ′ c1 and η c2 , respectively. In literature, there have already existed a number of studies about critically distinguishing the properties of χ ′ c1 and η c2 , such as the radiative transitions η c2 (χ ′ c1 ) → J/ψ(ψ ′ ) [26][27][28], the η c2 (χ ′ c1 ) hadroproduction rates [29], or η c2 (χ ′ c1 ) production rates in B decay [30]. By confronting the established properties of X(3872), all these studies tend to disfavor the 1 D 2 assignment of the X(3872) meson. We hope that doublecharmonium production can also be added to the above list as a valuable means to help clarify the situation.
It turns out that the production rate of e + e − → J/ψ + χ ′ c1 is about 6-7 times greater than that of e + e − → J/ψ + η c2 . Based on the 1 ab −1 data currently accumulated at the Belle experiment, it looks promising to observe the former process, if the χ ′ c1 can indeed be regarded as the X(3872) particle with very narrow width. This seems to constitute the strong enough incentive for experimentalists to make an updated analysis of the double charmonium production at B factory.
The rest of the paper is organized as follows. In Section II, we specify the helicity selec-tion rule suited for the hard exclusive reaction e + e − → J/ψ + η c2 , and give the definition for the dimensionless, reduced helicity amplitudes. In Section III, we present the leading order expressions for all the independent helicity amplitudes in the NRQCD factorization framework. In Section IV, we first review some key technical issues about the O(α s ) calculation, then present the asymptotic expressions for the NLO perturbative corrections to all the encountered helicity amplitudes. The pattern of the double-logarithmic scaling violation is confirmed once again. In Section V, a comparative study is performed for both unpolarized and polarized cross sections between the processes e − e + → J/ψ +η c2 and e − e + → J/ψ +χ ′ c1 . This study may shed some light on unveiling the quantum number of the X(3872) meson in the future double-charmonium production experiments. Finally we summarize in Section VI. In Appendix A, we tabulate the coefficients of the double logarithm ln 2 (s/m 2 c ) associated with all the relevant NLO Feynman diagrams, for all the double-charmonium production channels we have studied so far, i.e., e + e − → J/ψ + η c2 (η c , χ c0,1,2 ).

II. HELICITY SELECTION RULE AND REDUCED HELICITY AMPLITUDES
It is often desirable to glean more information than simply present the unpolarized cross section for a hard exclusive reaction, especially for the double-charmonium production process considered in this work. It is of particular advantage to making explicit predictions to various J/ψ + η c2 production rates for different helicity configurations. The underlying reasons of carrying out such detailed studies are two-fold. On the experimental ground, once the sufficient statistics is achieved, the helicity amplitudes themselves in principle can be measured by studying the angular distributions of these charmonia and their decay products; from the theoretical perspective, it is also instructive to stay with the helicity amplitudes. The reason is that, for a hard exclusive reaction, the relative importance of the polarized cross section in a given helicity channel is dictated by the celebrated helicity selection rule (HSR) [31].
We are interested in the hard-scattering limit √ s ≫ m c ≫ Λ QCD , where √ s stands for the center-of-mass energy of the e + e − collider, m c for the charm quark mass, and Λ QCD for the intrinsic QCD scale. In this limit, the asymptotic behavior of the production rate for J/ψ + η c2 in a definite helicity configuration follows from the HSR [10]: where λ 1 , λ 2 represent the helicities carried by the J/ψ, η c2 , respectively. v denotes the characteristic velocity of charm quark inside a charmonium. Equation (1) implies that the helicity state which exhibits the slowest asymptotic decrease, thus constitutes the "leadingtwist" contribution, i.e., σ ∼ 1/s 3 , is (λ 1 , λ 2 ) = (0, 0). We have chosen by default to work in the e + e − center-of-mass frame. Let |P| signify the magnitude of the momentum carried by the J/ψ (η c2 ), and θ denote the angle between the moving directions of the J/ψ and the e − beam. Following the steps elaborated in [24], the differential rate for polarized J/ψ + η c2 production in e + e − annihilation can be expressed as where A λ 1 ,λ 2 is the helicity amplitude associated with the virtual photon decay into J/ψ+η c2 carrying the helicity component (λ 1 , λ 2 ). Parity invariance can be invoked to reduce the number of independent helicity amplitudes: hence the two helicity amplitudes related by flipping the helicities of two chamonia bear the equal magnitude. An immediate consequence of (3) is that, the virtual photon decay into the longitudinally-polarized J/ψ and η c2 is strictly forbidden by parity invariance. Integrating (2) over the polar angle θ and including all the allowed helicity states, it is then straightforward to obtain the unpolarized cross section: In conformity with the constraint |λ 1 − λ 2 | ≤ 1, as demanded by angular momentum conservation, there are totally 4 independent helicity amplitudes for γ * → J/ψ + η c2 (Recall that A 0,0 = 0 owing to parity invariance). In equation (4), we have also retained a factor of 2 explicitly to account for the contributions from those helicity-flipped states. We have also adopted the approximation 2|P| √ s ≈ √ 1 − 4r by assuming M J/ψ ≈ M η c2 ≈ 2m c . In the NRQCD factorization framework, the product of two nonperturbative factors, i.e., the (second derivative of) wave functions at the origin for the charmonia J/ψ, and η c2 : R J/ψ (0), and R ′′ η c2 (0), ubiquitously enters every helicity amplitude, thereby it appears convenient to define a reduced dimensionless helicity amplitude, of which these nonperturbative factors are pulled out. We introduce the reduced helicity amplitude, a λ 1 ,λ 2 , which is related to the standard helicity amplitude A λ 1 ,λ 2 as follows: where e c e = 2 3 e is the electric charge of the charm quark, r ≡ 4m 2 c /s signifies a dimensionless mass ratio. To make the HSR manifest, we have explicitly stripped off a factor r 1 2 (|λ 1 +λ 2 |−1) in (5), so that the dimensionless helicity amplitude a λ 1 ,λ 2 is expected to scale as O(r 0 ). Plugging (5) back into (4), we obtain the NRQCD predictions to the polarized production rate for the J/ψ(λ 1 ) + η c2 (λ 2 ) state:

III. LO PREDICTIONS FOR THE HELICITY AMPLITUDES
The reduced helicity amplitude a λ 1 ,λ 2 can be viewed as the NRQCD short-distance coefficient, which encodes the contribution solely stemming from the momentum region between m c and √ s, so can be computed reliably in perturbation theory. It is convenient to parameterize it as FIG. 1: One sample LO diagram and five sample NLO diagrams that contribute to γ * → J/ψ +η c2 .
Our central task in this work is to decipher the O(α s ) correction to the reduced amplitude, a First we recapitulate the LO calculation. To proceed, it it most convenient to consider the quark amplitude γ * → cc( 3 S 2 ) using the covariant projection technique [10,32]. At LO in α s , there are only four Feynman diagrams that contribute, one of which is depicted in Fig. 1a) (For simplicity, we have neglected the QED fragmentation diagrams, whose effect appears to be modest).
It is straightforward to project out 4 corresponding helicity amplitudes from the decay amplitude for γ * → J/ψ + η c2 . One then follows equation (5) to read off each of the LO reduced helicity amplitudes: For some accidental reason, only the J/ψ(±1) + η c2 (0) channel has the non-vanishing cross section at LO in α s . Substituting (8) into (6), we find agreement with the QCD part of the LO prediction for the unpolarized cross section first given in Ref. [10].

IV. NLO PERTURBATIVE CORRECTIONS TO THE HELICITY AMPLITUDES
We start this section by first sketching some technical issues about the NLO perturbative calculations, followed by presenting the asymptotic expressions of the O(α s ) corrections to all the reduced helicity amplitudes.

A. Outline of the calculation
At NLO in α s , there are 20 two-point, 20 three-point, 18 four-point, and 6 five-point oneloop diagrams for the process γ * → cc( 3 S (1) 2 ), some of which have been illustrated in Fig. 1. The calculation is quite similar to our preceding works on double charmonium exclusive production, i.e., e + e − → J/ψ + χ c0,1,2 [24] and the e + e − → J/ψ + η c [21], so here we will only present a very brief description.
We adopt dimensional regularization to regularize both UV and IR singularities. We follow the 't Hooft-Veltman prescription for γ 5 , as detailed in Ref. [33]. In projecting out the D-wave orbital angular momentum state, one needs expand the relative quark momentum q to the quadratic order. Our strategy is to follow the method of region [34] via directly deducing the NRQCD short-distance coefficients, rather than resorting to the much more expensive matching calculation, through making the expansion in q before carrying out the loop integration.
The only nontrivial technical problem is that one may encounter some unusual one-loop integrals which generally contain the propagators of cubic power, as a consequence of taking the derivative over q twice. The Mathematica package FIRE [35] and the code Apart [36] are utilized to reduce these unconventional higher-point one-loop tensor integrals into a set of masters integrals. With the aid of the integration-by-part algorithm and partial fractioning technique, it turns out that all the encountered master integrals are nothing but the standard 2-point and 3-point one-loop scalar integrals, whose analytic expressions have already been tabulated in Ref. [17].
When adding the contributions of all the diagrams, and after renormalizing the charm quark mass and the QCD coupling constant, we finally end up with both UV and IR finite NLO expressions for the decay amplitude γ * → cc( 3 S (1) 2 ). Since everything becomes finite, we can safely return to the 4 spacetime dimensions, and readily project out all the required helicity amplitudes.
In Ref. [37], an all-order-in-α s proof for exclusive quarkonium production has been outlined in the NRQCD factorization context. It argues that at lowest order in v and to all orders in α s , NRQCD factorization holds for the exclusive production of a S-wave quarkonium plus any higher-orbital-angular-momentum quarkonium in e + e − annihilation. Our explicit calculation confirms that the NRQCD short-distance coefficients affiliated with S-wave plus D-wave charmonia production are indeed IR finite at NLO in α s , which is compatible with what is asserted in [37].

B. Analytic expressions of NLO helicity amplitudes
The a (1) λ 1 ,λ 2 are complex-valued, whose analytic expressions are in general quite lengthy. Rather than reproducing their cumbersome-looking expressions here, we are content with presenting their numerical values over a wide range of r, as shown in Fig. 2.
As a matter of fact, it seems much more illuminating to know the asymptotic behaviors of the helicity amplitudes in the limit √ s ≫ m c . At NLO in α s , one anticipates to see the logarithmic scaling violation to the naive power-law HSR as indicated in (1). Furthermore, conducting the asymptotic expansion in NRQCD short-distance coefficients is theoretically appealing, since it is equivalent to disentangling the contributions occurring at the "hard" scale (virtuality ∼ s) from the "lower-energy" collinear/soft sectors (virtuality ∼ m 2 c ), by which one can intimately link the NRQCD factorization approach and the light-cone approach [38][39][40][41][42]. Such an asymptotic expansion on the NRQCD hard coefficients has been carried out for a number of exclusive double-quarkonium production processes, either in e + e − annihilation [21,24,42], or from bottomonium decay [42][43][44][45], and some general pattern about logarithmic scaling violation has been recognized [24,42]. According to the definition in (6), we find the asymptotic expressions of four reduced NLO helicity amplitudes for γ * → J/ψ + η c2 to be where µ is the renormalization scale, and β 0 = 11 3 C A − 2 3 n f is the one-loop coefficient of the QCD β function, and n f = 4 denotes the number of active quark flavors. In contrast with (8), all the four reduced helicity amplitudes receive non-vanishing O(α s ) corrections.
Note that the β 0 ln(4µ 2 /s) term only resides in the (±1, 0) channel, since this is the only channel that has a non-vanishing tree-level amplitude. For the sake of comparison, all these asymptotic results of the reduced helicity amplitudes are also shown in Fig. 2, in juxtapose with the corresponding exact NLO results. These asymptotic results appear to converge with the exact ones decently well even at relatively lower √ s, say, at B factory energy.
From Eq. (9), one sees that the leading scaling violation is due to the double logarithm ln 2 r. This constitutes another example to further corroborate the earlier conjecture: the occurrence of ln 2 r in the one-loop NRQCD short-distance coefficients is always affiliated with the helicity-suppressed channels in exclusive double charmonium production processes [24,42].

V. PHENOMENOLOGY
Aside from a comprehensive analysis on the process e + e − → J/ψ + η c2 at the B factory, in this section we also target at a detailed investigation on the process e + e − → J/ψ + χ ′ c1 , where the essential elements have already been set up in [24]. This study is largely motivated by the recent concern about the nature of the X(3872) meson, whether its quantum number is 1 ++ or 2 −+ . The canonical charmonium options for the X(3872) are the χ ′ c1 and η c2 , respectively. With this consideration in mind, our study may provide some useful guidance on unveiling the quantum number of the X(3872) through the dedicated double-charmonium experiment at future B factory. Before making concrete predictions to the production rates for e + e − → J/ψ + η c2 , we first address a subtlety about utilizing (6). For the J/ψ(±1) + η c2 (0) helicity state, the canonic way of identifying the O(α s ) correction is by replacing |a ±1,0 | 2 with 2 αs π ℜ[a (0) ±1,0 ]. Because of the null tree-level amplitudes for the remaining helicity configurations, the O(α s ) correction to the total cross section is solely from the (±1, 0) channels. On the other hand, the polarized cross section for each helicity state is also a physical observable by itself. Since the non-vanishing amplitudes are first generated at O(α s ) for the remaining helicity states, it is thereby also sensible to interpret |a λ 1 ,λ 2 | 2 = αs π 2 a (1) λ 1 ,λ 2 2 for (λ 1 , λ 2 ) = (±1, 0). Though formally of order α 2 s , these new pieces do constitute the leading contributions to the respective polarized cross sections, thereby still being consistent 1 . We will follow this viewpoint in presenting our NLO predictions.
In the numerical analysis, we take √ s = 10.58 GeV, and the charm quark pole mass m c = 1.4 ± 0.2 GeV. The fine structure constant is chosen as α( √ s) = 1/130.9 [19]. The running QCD strong coupling constant is evaluated by using the two-loop formula with Λ (4) MS = 0.338 GeV [16,17]. The nonpertubative input parameters, i.e., (the derivatives of ) the wave function at the origin for J/ψ, χ ′ c1 , and η c2 , also suffer from a large amount of uncertainty. Their values have been compiled in Ref. [46], which were deduced from several different potential models. We choose to use those given by the Buchmüller-Tye potential model [47]: (0)| 2 = 0.102 GeV 5 , and |R ′′ η c2 (0)| 2 = 0.015 GeV 7 . Another important source of uncertainty for the NLO predictions stems from the scale affiliated with the strong coupling constant. As is well known, the scale ambiguity is a typical nuisance of NRQCD factorization approach, reflecting the fact that two disparate hard scales,  In the rightmost column, we also give the K factor for the unpolarized cross sections. √ s and m c , are entangled together in NRQCD short-distance coefficients. In fact, the lesson gained from the light-cone approach strongly suggests that it is rather implausible to set the scales entering all α s in NLO short-distance coefficient to be unanimously around √ s [42].
Although we are unable to circumvent the scale ambiguity problem in the confine of NRQCD approach, we may allow the µ to float in the range between 2m c and √ s, hoping that the most trustworthy prediction may interpolate in between. In Fig. 3, we show the total cross sections both for e + e − → J/ψ+η c2 and e + e − → J/ψ+χ ′ c1 as a function of µ, which also take into account the error due to the uncertainty of m c . One sees clearly that, after incorporating the NLO perturbative correction, the scale dependence of the cross section has been reduced. In Fig. 4, we also plot the LO and NLO total cross sections for e + e − → J/ψ + η c2 as a function of √ s. Since the leading-twist (0, 0) channel is absent in both of the processes, the cross sections drop very rapidly (∝ 1/s 4 ) as √ s increases.
In Tables I and II, we also list the predictions for the polarized cross sections from each individual helicity channel as well as the total cross sections for e + e − → J/ψ + η c2 and e + e − → J/ψ + χ ′ c1 , with m c fixed at 1.4 GeV but with the renormalization scale chosen . We have taken m c = 1.4 GeV. In the rightmost column, we also list the K factor for the unpolarized cross sections.  fb, respectively. The production rate of the latter process is about six times more copious than that of the former.
Inspecting Tables I and II, one observes an interesting hierarchy of the polarized cross sections for the above processes. The overwhelming contribution for the J/ψ+η c2 production comes solely from the (±1, 0) state, with J/ψ transversely polarized, while the dominant helicity channel for producing J/ψ + χ ′ c1 is instead the (0, ±1) state, with J/ψ longitudinally polarized. This characteristic of J/ψ polarization may serve as a benchmark for the future double-charmonium production experiment to unveil the nature of the X(3872), provided that the prospective Super B factory can observe sufficient number of X(3872) events recoiling against J/ψ.
We can ascertain the observation potential of these two exclusive double charmonium production processes. Thus far, Belle experiment has accumulated about 1000 fb −1 data. Taking σ[e + e − → J/ψ + η c2 ] ≈ 0.23 − 0.33 and σ[e + e − → J/ψ + χ ′ c1 ] ≈ 1.53 − 2.16 fb from Tables I and II, we estimate that roughly 230 − 330 J/ψ + η c2 events, and 1500 − 2200 J/ψ + χ ′ c1 events have been produced at Belle at around √ s = 10.58 GeV. Since neither of the masses and decay patterns of η c2 and χ ′ c1 is known, it seems experimentally impossible to simultaneously reconstruct the J/ψ and the η c2 (χ ′ c1 ) signals. The viable experimental way is to only reconstruct the J/ψ → l + l − event, and fit the recoil mass spectrum against the J/ψ to estimate the number of η c2 (χ ′ c1 ) peak events. This method does not depend on the concrete decay modes of η c2 (χ ′ c1 ), and is particularly suitable for the limited statistics of signal events like in our case. In fact, this method has already been used routinely by the Belle collaboration to search for the double charmonium production processes e + e − → J/ψ + C-even charmonium.
In fitting the recoil mass spectrum, the net detection efficiency for e + e − → J/ψ + η c2 (χ ′ c1 ) may reach around 4% (Similar for η c2 and χ ′ c1 , with the reconstruction efficiency for J/ψ → l + l − included [48]). As a very crude estimate, the number of observed e + e − → J/ψ+η c2 (χ ′ c1 ) events are expected to reach (230 − 330) × 4% = 9 − 13, and (1500 − 2200) × 4% = 60 − 90, respectively. Since only the J/ψ is reconstructed, the background level can be a little higher. With only around 9−13 observed J/ψ+η c2 events, it appears difficult to observe a significant signal with current 1 ab −1 Belle data. If χ ′ c1 is indeed the X(3872) meson, thanks to the very narrow width of the X(3872), it seems possible to observe the 60-90 J/ψ + χ ′ c1 signal events with the current Belle full data sample, which may provide a strong incentive for updating their earlier e + e − → J/ψ+ charmonium analysis with only 673 fb −1 data [5,49]. This study may hopefully help one to better understand the mechanism about the X(3872) production. On the other hand, in case that χ ′ c1 is not X(3872) and its width is no longer narrow, it may be difficult to observe a significant signal at the moment.
In any event, a much larger data samples seem to be called for to arrive at some definite conclusion. In the prospective Super B factory, which may reach an integrated luminosity of 50 ab −1 by 2022, it seems feasible that the processes e + e − → J/ψ + η c2 (χ ′ c1 ) will eventually be observed, and the quantum number of X(3872) also hopefully be nailed down.

VI. SUMMARY
In this work we have calculated the complete NLO perturbative corrections to e + e − → J/ψ + η c2 within the NRQCD factorization framework, and found that the O(α s ) correction can be sizable for a larger value of the renormalization scale. Our calculation indicates that the dominate contribution to the total cross section comes from the helicity states (±1, 0).
We have also carried out a comparative study for e + e − → J/ψ+η c2 and e + e − → J/ψ+χ ′ c1 at B factory energy, with the specific motivation that these double charmonium production processes may provide useful means to unveil the quantum number of the X(3872) meson in the future B factories. It turns out that the production rate of the latter process is about 6-7 times greater than the former, therefore it is possible to observe the J/ψ + χ ′ c1 signals based on the current 1 ab −1 Belle data sample, if the χ ′ c1 is indeed the very narrow X(3872) particle. The dominantly transverse polarization of the J/ψ is also a useful indicator for identifying the X(3872) with χ ′ c1 . A necessary extension of our current work is to further incorporate the leading relativistic correction to the above double-charmonium production processes, whose effects might be as important as the radiative corrections. The study along this direction is of some theoretical interest, especially regarding that, until today the O(v 2 ) calculation hardly exists for the production processes involving the P , D-wave quarkonium.
It has been recently realized that the O(α s ) NRQCD short-distance coefficients for the helicity-suppressed double charmonium-production processes are often plagued with the double logarithm of form ln 2 (s/m 2 c ). This symptom is likely intertwined with some long-standing failure of applying the light-cone approach to the "higher-twist" hard exclusive reactions beyond tree level. To expedite a better understanding and controlling of these double logarithms, we have scrutinized all the NLO diagrams for the various exclusive doublecharmonium processes that have been studied so far, e.g., e + e − → J/ψ + η c2 (η c , χ c0,1,2 ), and singled out those that contain the double logarithm, and enumerate their coefficients for each individual helicity states. No simple and general pattern has been recognized yet.
However, the majority of phenomenologically relevant double-charmonia production processes are of "higher-twist" nature, exemplified by γ * → J/ψ +η c2 (η c , χ c0,1,2 ), for which some of or all of possible helicity channels are of HSR-suppressed type. Until today, it is still not clear how to consistently compute the NLO radiative correction to this type of processes in the light-cone framework, due to some long-standing problem such as the inevitable emergence of the endpoint singularity. By contrast, NRQCD factorization approach, which is based on the velocity expansion rather than the twist expansion, serves as the only viable tool to investigate the O(α s ) corrections to this kind of helicity-suppressed hard exclusive reactions.
An abnormal feature about the O(α s ) NRQCD short-distance coefficients for the helicitysuppressed double charmonium production processes seems to be the frequent emergence of the double logarithm ln 2 (s/m 2 c ), which seems to seriously deteriorate the convergence of the perturbative expansion in NRQCD factorization 2 . This double logarithm appears to result from the overlap between the "collinear" and "soft" regions in the NRQCD shortdistance coefficient and seems deeply related to the afore-mentioned endpoint singularity problem. Although it is feasible to reproduce the closed form of the double logarithms for each concrete process, the lack of a thorough understanding of their behavior in higher order in α s prevents one from systematically controlling them, i.e., to resum them to all orders in α s .
A useful starting point is to anatomize all the double logarithms of form ln 2 (s/m 2 c ) from the existing NLO calculations for the various exclusive double charmonium processes. For this purpose, we examine all the NLO Feynman diagrams for the processes e + e − → J/ψ + η c2 (η c , χ c0,1,2 ), sort out those that contain the double logarithm, and enumerate their coefficients.
A careful examination reveals that (almost) all the relevant NLO diagrams that contribute the double logarithm ln 2 (s/m 2 c ) (in Feynman gauge) for the above processes are included in Fig. 5, which was first identified in [42].
In Table III, we tabulate the double logarithms that come from each individual diagram in Fig. 5, with the color-structure specified, for the processes e + e − → J/ψ + η c2 (η c ).
The reaction e + e − → J/ψ(λ 1 ) + η c (λ 2 ) only possesses one independent helicity configuration (±1, 0). Each diagram contains a non-vanishing double logarithm. Summing up all their contributions, and multiply by 2 to account for the contribution from the chargeconjugated diagrams, we recover the net double logarithm in the NLO short-distance coefficient: C (1) C (0) asym = 13 24 ln 2 s m 2 c , as given in equation (37) of Ref. [42]. All the four polarized production channels e + e − → J/ψ(λ 1 ) + η c2 (λ 2 ) are helicitysuppressed. Adding up the contributions from all the diagrams in Fig. 5, and multiply by 2, we recover the net double logarithm in the NLO reduced helicity amplitudes a (1) λ 1 ,λ 2 asym , as given in equations (9). Note that the double logarithms from the (±1, 0) channel has the similar diagram-by-diagram structure as its counterpart for e + e − → J/ψ + η c . Neverthe-2 In fact, one will encounter the double logarithm even in the single-quarkonium exclusive production process such as γ * → η c + γ, once the corresponding O(α s ) NRQCD short-distance coefficient is expanded to NLO in m 2 c /s. This can be readily checked by using the analytic O(α s ) short-distance coefficient given in [38,50].  Fig. 5 from the various helicity states (λ 1 , λ 2 ) in the processes γ * → J/ψ + η c2 (η c ). The two types of color factors are denoted by C 1 = C 2 F and C 2 = C F (C F − 1 2 C A ), where C F = 4 3 and C A = 3.
less, for the remaining three helicity channels, a great simplification arises that many of the diagrams in Fig. 5 do not contain double logarithm.