Next-to-Next-to-Leading-Order Radiative Corrections to $e^+e^-\to\chi_{cJ}+\gamma$ at $B$ factory

Within the nonrelativistic QCD (NRQCD) factorization framework, we have computed the ${\mathcal O}(\alpha_s^2)$ corrections to the exclusive production of $P$-wave spin-triplet charmonia $\chi_{cJ}(J=0,1,2)$ accompanied by a hard photon at $B$ factories. For the first time, we have explicitly verified the validity of NRQCD factorization for exclusive $P$-wave quarkonium production at two-loop order. Unlike the $\chi_{cJ}$ electromagnetic decay processes, the $\mathcal{O}(\alpha_s^2)$ corrections are found to be smaller than the $\mathcal{O}(\alpha_s)$ corrections in all three channels $e^+e^-\to \chi_{c0,1,2}+\gamma$. In particular, the ${\mathcal O}(\alpha_s^2)$ corrections appear moderate for $\chi_{c1}$ case, and insignificant for $\chi_{c0}$. In addition, the next-to-next-to-leading order (NNLO) predictions for the production rates of $\chi_{c0,1}+\gamma$ are insensitive to the renormalization and factorization scales. All of these features may indicate that perturbative expansion in these two channels exhibits a decent convergence behavior. By contrast, both the ${\mathcal O}(\alpha_s)$ and ${\mathcal O}(\alpha_s^2)$ corrections to the $\chi_{c2}+\gamma$ production cross section are sizable, which even reduces the Born order cross section by one order of magnitude after including the NNLO perturbative correction. Taking the values of the long-distance NRQCD matrix elements from nonrelativistic potential model, our prediction to $\chi_{c1}+\gamma$ production rate is consistent with the recent {\tt Belle} measurement. The NNLO predictions for the $\chi_{c0,2}+\gamma$ production rates are much smaller than that for $\chi_{c1}+\gamma$, which seems to naturally explain why the $e^+e^-\to \chi_{c0,2}+\gamma$ channels have escaped experimental detection to date.


Introduction
Heavy quarkonia, the tightly-bound systems composed of a heavy quark and a heavy antiquark, are generally viewed as the simplest hadrons in Quantum Chromodynamics (QCD). A peculiar trait of quarkonia is the coexistence of several distinct mass scales, which makes it an interesting and unique laboratory to sharpen our understanding about the interplay between perturbative and nonperturbative aspects of QCD. Since the heavy (anti)quark inside a quarkonium is essentially nonrelativistic, the mainstream theoretical method is firmly based on the modern effective field theory (EFT) doctrine, the so-called Nonrelativistic QCD (NRQCD) factorization approach [1]. This approach allows to systematically disentangle the short-distance and long-distance effects, formalized by a double expansion in QCD strong coupling α s and heavy quark velocity v. In the past two decades, NRQCD factorization has been widely employed to tackle innumerable quarkonium production and decay processes.
Thanks to its enormous luminosity and simplicity of the initial state, B factories have acted as an fertile and clean playground to investigate charmonium production. For instance, in the past two decades, there have emerged a handful of experimental measurements about exclusive charmonium production processes [2][3][4][5][6], together with intensive theoretical investigations using NRQCD approach [7][8][9][10][11] (for a comprehensive list of references, we refer the interested readers to the recent review article [12]).
Among various charmonium production processes, the exclusive production of positive-C-parity charmonium associated with a hard photon, e.g., e + e − → η c (χ cJ )+γ (J = 0, 1, 2), is of special interest. Due to their simplicity, these processes can be regarded as the golden channels to test our understanding of charmonium production mechanism and the utility of NRQCD approach. The leading-order (LO) cross section for e + e − → η c (χ c0,1,2 ) + γ at B factory was predicted in [13]. The next-to-leading-order (NLO) perturbative corrections were subsequently computed in [14,15], where the O(α s ) corrections in the χ c2 + γ channel turn out to be sizable, even exceeding −60% 1 . On the other hand, since charm quark is not decently heavy, one may expect that relativistic corrections might also have important impact. The leading relativistic corrections to e + e − → η c + γ were first considered in [14]. The relativistic corrections to e + e − → χ c0,1,2 + γ were first explored in [18], yet missing the contribution due to the NRQCD operators that explicitly contains the chromoelectric field. Very recently, the complete O(v 2 ) corrections to these P -wave charmonium exclusive production processes have been given in [19]. Unfortunately, the values of various O(v 2 ) NRQCD long-distance matrix elements (LDMEs) are poorly constrained, therefore it is difficult to present accurate predictions for the χ cJ + γ production rates.
Leaving relativistic corrections aside, one has witnessed remarkable progress in deducing the higher-order perturbative corrections for various quarkonium decay and production processes. More than two decades ago, the next-to-next-to-leading order (NNLO) perturbative corrections for the simplest quarkonium electromagnetic decay, Υ(J/ψ) → e + e − , were analytically determined [20,21]. It is worth noting that, the O(α 3 s ) corrections for these channels have also been available recently [22,23]). In recent years, with the advance of numerical and analytic multi-loop technology, a number of NNLO perturbative corrections to more difficult quarkonium decay processes have been accomplished, e.g., η b,c → γγ [24,25], B c → ℓν [26,27], χ c0,2 → γγ [28], and η c,b → light hadrons [29]. For most of the aforementioned processes, the NNLO radiative corrections turn out to be sizable and significantly modify the lower-order NRQCD predictions, especially for charmonia.
The NNLO perturbative corrections to the simplest channel of exclusive quarkonium production, the γγ * → η c,b transition form factor, have also been reported recently [25,30]. Very recently, the NNLO radiative corrections to the famous double-charmonium production process at B factory, e + e − → J/ψ + η c , have also been inferred [31]. In this case, the O(α 2 s ) corrections are observed to have moderate effect. It is encouraging that the state-of-the-art NRQCD prediction is consistent with the BaBAR measurement [3]. Moreover, the NNLO radiative corrections to e + e − → η c + γ has also recently been computed analytically at lowest order in v [32], again with moderate impact 2 . Very recently, the NNLO corrections to this process were also reinvestigated with the renormalization scales chosen by the principle of maximum conformality [33].
In this work, we proceed to further evaluate the NNLO perturbative corrections for e + e − → χ cJ +γ at B factory. This work constitutes the first NNLO perturbative correction calculation for the P -wave quarkonium production process. It is of theoretical curiosity to examine the validity of NRQCD factorization framework for this case. We also wish to examine the convergence of perturbative expansion in this channel, as well as confront our predictions with the latest measurement on e + e − → χ c1 + γ by Belle Collaboration [34].
The remainder of this paper is structured as follows. In section 2, after outlining the 1 When √ s ≫ mc, the collinear logarithm ln s/m 2 c in NRQCD short-distance coefficients can get large, which may potentially ruin the convergence of fixed-order perturbative expansion. For e + e − → ηc + γ, there have been attempts to resum the leading logarithms [16] and next-to-leading logarithms [17] to all orders in αs. 2 Note that the NRQCD short-distance coefficient has a rather lengthy expressions in term of Goncharov polylogarithms, the integrals over polylogarithms and complete elliptic integrals. NRQCD factorization formula for e + e − → χ cJ + γ, we describe the theoretical strategy to deduce the NRQCD SDCs associated with the J = 0, 1, 2 channels. In section 3, we first briefly describe some technicalities encountered in two-loop calculation, then present the numerical results for various NRQCD short-distance coefficients. Section 4 is devoted to the phenomenological analysis and discussion. We summarize in section 5.
2 Theoretical background for P -wave onium exclusive production In accordance with the NRQCD factorization ansatz, the production rates of e + e − → χ cJ + γ can be expressed in the following factorized form: where F 1 ( 3 P J ) (J = 0, 1, 2) represent the corresponding NRQCD short-distance coefficients (SDCs). Owing to asymptotic freedom, these coefficients can be computed in perturbation theory, order by order in powers of the strong coupling constant α s . O( 3 P J ) represent the process-independent NRQCD long-distance matrix elements (LDMEs), which bear a genuinely nonperturbative origin, and are defined as where ψ and χ † denote the Pauli spinor fields annihilating a heavy quark and antiquark, respectively, and Invoking the approximate heavy quark spin symmetry, we have the following simplifying relations: Since NRQCD SDCs are insensitive to the nonperturbative hadronization effects, they can be deduced with the aid of the standard perturbative matching technique. That is, by replacing the physical χ cJ meson with a fictitious onium composed of free cc pair, carrying the quantum number 3 P J , we compute both sides of (2.5) order by order in α s . After this replacement, (2.4) becomes where the subscript cc( 3 P J ) in the NRQCD LDMEs indicates that the hadronic states χ cJ have been replaced by a cc( 3 P J ) pair, which can be accessed in perturbation theory. After computing both sides of (2.5) in perturbative QCD and NRQCD, we are able to solve for the desired NRQCD SDCs order by order in α s . It is worth emphasizing that, for a quarkonium hard exclusive reaction like in our case, the factorization (2.5) actually also holds at the amplitude level.
To facilitate the perturbative calculation, we assign the momenta of the c andc quarks to be where P and q denote the total momentum of the cc pair and the relative momentum, respectively. The on-shell condition p 2 =p 2 = m 2 (with m signifying the charm quark mass) enforces that Since we are only concerned with the SDCs at the lowest order in v, it is legitimate to approximate the square of the invariance mass of the cc pair by 4m 2 .
It is convenient to employ the covariant spin-projector to enforce the cc pair in the spintriplet state. The relativistically normalized color-singlet/spin-triplet projector reads [35]: J ) amplitude can be projected out by differentiating the colour-singlet/spintriplet quark amplitude A with respect to the relative momentum q, followed by setting q to zero: µν denoting the polarization vectors affiliated with J = 0, 1, 2. In order to obtain the unpolarized cross section, we also need sum over all possible polarizations for each J. It is convenient to employ the polarization sum identities given in [36]: where d = 4 − 2ǫ signifies the space-time dimensions, and the polarization tensor Π µν (P ) is defined through Now we have collected all the necessary ingredients to evaluate the quark-level cross sections σ(cc( 3 P J ) + γ) in perturbative QCD. In the meanwhile, the perturbative NRQCD matrix elements O( 3 P J ) cc( 3 P J ) can also be carried out. It is then straightforward to ascertain the SDCs following the matching procedure.
3 The cross sections through O(α 2 s ) In this section, we first describe the computational techniques utilized to determine the various two-loop SDCs F 1 ( 3 P J ), then present their numerical expressions.

NLO NNLO
regular light by light LO Figure 1. Some representative Feynman diagrams for e + e − → cc( 3 P J ) + γ, which are drawn by JaxoDraw [37] We use FeynArts [38] to generate the Feynman diagrams for e + e − → cc + γ and corresponding Feynman amplitude through two-loop order in α s . Some typical Feynman graphs are displayed in Fig. 1. Employing the color-singlet/spin-triplet projector (2.8), following the recipe as specified in (2.9) to single out the P -wave component of the amplitude, we are able to project out the intended e + e − → cc( 3 P J ) + γ amplitude, order by order in α s . We then employ the packages FeynCalc [39,40] and FormLink [41] to carry out the polarization sum according to (2.10).
The leading-order (LO) NRQCD SDCs have long been known [14]: where α represents the electromagnetic coupling constant, e c = 2/3 signifies the charge of charm quark, s corresponds to the squared center-of-mass energy. The dimensionless ratio r ≡ 4m 2 /s is constrained to be less than 1.
Once beyond the LO, we adopt the standard shortcut to directly extract the shortdistance coefficient, i.e., to take the derivative of the amplitude with respect to q prior to conducting the loop integral, which amounts to directly extracting the contribution from the hard region in the context of strategy of region [42]. Throughout the work we employ the dimensional regularization to regularize the UV and IR divergences. For the loop integrals, we utilize the packages Apart [43] and FIRE [44] to conduct partial fraction and the corresponding integration-by-part (IBP) reduction. We end up with 6 one-loop master integrals (MIs) and 174 two loop MIs, most of which are complex-valued integrals. For the one-loop MIs, one can readily work out the analytical expression for all the MIs. We have confirmed the analytic expression of the O(α s ) corrections to the cross section, first accomplished in [14].
It becomes much more challenging to deduce the analytical expressions for all the encountered two-loop MIs. In this work, we are content with high-precision numerical results 3 . We use the modified FIESTA [45] to perform sector decomposition (SD) for the two-loop MIs. For the real-valued MIs, we directly use CubPack [46] to carry out the numerical integration. In contrast to the application of SD to the Euclidean region, the singularities encountered in the physical region lie inside, rather than sit on, the integration boundary, which render the integrals hard to be numerically evaluated. The difficulty can be overcome to a certain extent by deforming the integration contour via the following variable transformation prior to decomposing the sectors [47]: where F denotes the F -term in the α parametrization, λ k is some positive number. Actually, the integration efficiency may vary drastically with λ k . In our calculation, we first choose a set of λ k and utilize CubPack to conduct the first-round rough numerical integration. For those integrals with large estimated errors, we adjust the values of λ k and perform the integration with a fixed number of sample points. The operation will be repeated until we find a optimized values of λ k , which render the integration error endurable.
With the new determined λ k , the integration will be performed once again with the aid of a parallelized integrator HCubature [48] to reach the desired precision. For more detail, we refer the interested readers to Ref. [31].
To eliminate UV divergences, we carry out the renormalization procedure by implementing the O(α 2 s ) expressions of the renormalization constants Z 2 and Z m from [49][50][51]. Nevertheless, the renormalized NNLO squared amplitudes are found to still contain an uncancelled single IR pole. This symptom is a common feature specific to NRQCD factorization, which have been encountered many times in NNLO perturbative calculations involving quarkonium. This IR pole can be factored into the NRQCD LDME, so that the NRQCD SDCs become IR finite. As a consequence, both of the LDMEs and the corresponding two-loop SDCs develop a log µ Λ dependence (µ Λ refers to NRQCD factorization scale), nevertheless their product must be independent of µ Λ . In fact, from the the coefficient of the single IR pole, one can read off the anomalous dimensions associated with the NRQCD bilinear currents carrying the quantum numbers 3 P J in (2.4): Reassuringly, these values exactly agree with those predicted from the renormalization group analysis in velocity NRQCD [52].
As mentioned before, we will be content with only providing the numerical expressions for various NRQCD SDCs. We then substitute these results into (2.5) to predict the unpolarized production rates for e + e − → χ cJ + γ, to the prescribed order in α s . For numerical calculation, we take the B factory center-of-mass energy to be √ s = 10.58 GeV. We choose two typical values for charm mass, m = 1.4 GeV and 1.68 GeV, which correspond to the one-loop and two-loop charm quark pole mass. The production rates for e + e − → χ cJ + γ through NNLO accuracy are then predicted to be  (2) ,  4) and (3.5), β 0 = (11/3)C A − (4/3)T F n f is the one-loop coefficient of the QCD β function, with T F = 1 2 and n f signifying the number of active quark flavors. In this work, we take n f = n L + n H , where n L = 3 labels the number of light quark flavors and n H = 1 indicates the number of heavy quark flavors. In addition, the symbol 'lbl' labels the contributions from the Feynman diagrams with the "light-by-light" topology, which are illustrated by the last diagram in Fig. 1. Note that the occurrence of γ χ cJ ln µ 2 Λ terms is reminiscent of the remnant of the uncancelled single IR pole, while the occurrence of the β 0 ln µ 2 R simply reflects the renormalization group invariance.

Phenomenology
In this section, we apply the formulas obtained in section 3 to make a concrete phenomenological analysis. First we need to fix the various input parameters. We take the running QED coupling constant evaluate at the B factory energy scale, α( √ s) = 1/130.9. The QCD running coupling constant is evaluated to two-loop accuracy with the aid of the package RunDec [53]. The NRQCD LDME for χ cJ is approximated by the first derivative of the Schroödinger radial wave function at origin through The 1P radial wave function at origin for χ c varies with different quark potential models. For instance, |R ′ 1P (0)| 2 = 0.075 GeV 5 in Buchmüller-Tye (BT) potential model, and |R ′ 1P (0)| 2 = 0.1296 GeV 5 in Cornell potential model [54,55]. Substituting these values into (4.1), one immediately obtains the corresponding NRQCD LDME O( 3 P J )) = 0.107 GeV 5 from BT potential model, and O( 3 P J ) = 0.186 GeV 5 from Cornell model 4 .
By setting the renormalization scale µ R = √ s/2 and the NRQCD factorization scale µ Λ = m, we can express the cross sections of e + e − → χ cJ + γ as the following power series: for m = 1.40 GeV, and In Ref. [13], the LDME O( 3 PJ ) is fitted via equating the NRQCD factorization predictions accurate through O(αs) with the measured values for χc0,2 → 2γ compiled by the particle data group [56]. The LDME is determined to be O( 3 PJ ) = 0.060 +0.043 −0.029 GeV 5 [13], which seems to be considerably smaller than the values given by potential models. Nevertheless, the O(α 2 s ) corrections χc0,2 → γγ turns out to be substantial [28]. For consistency, we will not use the fitted value of LDME in [13] in current work. Very recently, the authors in Ref. [57] have also explored the LDME O( 3 PJ ) based on the strong coupled pNRQCD, the values of which are found to be consistent with the potential models within the errors. , we clearly observe that the O(α 2 s ) corrections are less important than the O(α s ) corrections. Therefore, to some extent, the perturbative expansion in α s exhibits a decent convergence behavior, especially for the χ c0 + γ channel, in which the NNLO corrections only plays a minor role.
Assuming the LDME O( 3 P J ) = 0.107 GeV 5 as given by the BT potential model, and adopting two benchmark values of charm quark mass, in Table 1 we tabulate the NRQCD predictions to production rates for e + e − → χ cJ +γ at various level of perturbative accuracy. The uncertainties affiliated with the cross sections are estimated by varying µ R from 2m to √ s, with the central values evaluated at µ R = √ s/2. We notice the cross sections with m = 1.68 GeV are considerably smaller than those with m = 1.4 GeV, which can be attributed to the factor 1/m 3 that arise in the NRQCD SDCs for σ(e + e − → χ cJ + γ), as is evident in (3.1). It is interesting to note that, the perturbative corrections to σ(e + e − → χ c2 + γ), i.e., both the NLO and NNLO corrections, are sizable and negative. Incorporating the NNLO corrections reduces the LO prediction by almost one order of magnitude. In contrast, both the O(α s ) and O(α 2 s ) corrections are moderate for χ c1 + γ production, and have minor effect for the χ c0 + γ. In addition, from Table 1 one may also observe that the production rates for χ c0,1 + γ are insensitive to the factorization scale µ Λ .
Were the value of LDME chosen from Cornell model instead of BT model, our predicted cross sections would be enhanced roughly by a factor of 1.   . NRQCD predictions for σ(e + e − → χ cJ + γ) as a function of µ R at various levels of accuracy in α s . The value of the LDME is taken from the BT potential model. We take m = 1.4 GeV for the figures in the left panel, while m = 1.68 GeV in the right panel. The green band represents the uncertainty band by varying µ Λ from 1 GeV to m. In addition, we also display the Belle measurement by the yellow band for the e + e − → χ c1 + γ channel.
In Figure 2, we plot the cross sections as a function of the renormalization scale µ R at various level of perturbative accuracy, with the value of LDME taken from the BT potential model. We take m = 1.4 GeV on the left panel, and m = 1.68 GeV on the right panel. The green band labeled with "NNLO" is obtained by varying the factorization scale from 1 GeV to m. In order to facilitate comparison, we also demonstrate the Belle data by the yellow band in the plot of χ c1 + γ, where the red-dotted curve corresponds to the central value of the experimental measurement. We observe that the NNLO prediction seems to have a slightly reduced µ R -dependence relative to the NLO prediction. More importantly, our NRQCD prediction to σ(χ c1 + γ) with m = 1.4 GeV agrees well with the Belle measurement within the reasonable range of µ R .

Summary
In summary, in this work we have computed the O(α 2 s ) corrections to σ(e + e − → χ cJ + γ) at B factory. We choose two benchmark values for charm quark mass m = 1.4 GeV and m = 1.68 GeV, which correspond to the one-loop and two-loop charm quark pole mass, respectively. For the first time, we have verified that NRQCD factorization remains valid for exclusive P -wave quarkonium production to two-loop order. The impact of NNLO perturbative corrections is found to be significant for χ c2 + γ, moderate for χ c1 + γ, and rather marginal for χ c0 + γ. Unlike the electromagnetic decays χ c0,2 → γγ, the NNLO perturbative corrections to σ(e + e − → χ cJ + γ) are found to be smaller than the NLO pertubative corrections for all χ cJ + γ channels, which may indicate a satisfactory convergence in perturbative expansion. The NRQCD LDMEs are approximated by the first derivative of the χ cJ wave function at the origin deduced from nonrelativistic potential models. When taking the BT potential model as input, we find the NNLO predictions to the cross section of χ c1 + γ with m = 1.4 GeV is consistent with the Belle measurement within errors. On the other hand, when taking the LDME from Cornell model, the NRQCD prediction with m = 1.68 GeV turns out to be also consistent with the experimental measurement. After including the NNLO perturbative corrections, the production rates for χ c0,2 + γ are still found to be much smaller than that of χ c1 + γ, which may naturally explain why the e + e − → χ c0,2 + γ channels have escaped the experimental detection until today.