NNLO QCD Corrections to $\gamma + \eta_c(\eta_b)$ Exclusive Production in Electron-Positron Collision

Based on the NRQCD factorization formalism, we calculate the next-to-next-to-leading order QCD corrections to the heavy quarkonium $\eta_c (\eta_b)$ production associated with a photon at electron-positron colliders. By matching the amplitudes calculated in full QCD theory to a series of operators in NRQCD, the short-distance coefficients up to NNLO QCD radiative corrections are determined. It turns out that the full set of master integrals that we obtained could be analytically expressed in terms of Goncharov Polylogarithms, Chen's iterated integrals, and elliptic functions, which mostly do not exist in the literature and could be employed in the analyses of other physical processes. In phenomenology, numerical calculations of NNLO K-factors and cross sections of $e^+e^-\rightarrow \gamma + \eta_c(\eta_b)$ processes in BESIII and B-factory experiments are performed, which may stand as a test of the NRQCD higher order calculation while confronting to the data.

The study of heavy quarkonium production and decay is very important in the understanding of Quantum Chromodynamics (QCD), and heavy quarkonium as well. In the analysis of hadron, to overcome the obstacle of nonperturbative QCD is inevitable.
Fortunately, the advent of Nonrelativistic Quantum Chromodynamics (NRQCD) factorization formalism enables the calculation of quarkonium production and decay in a solid footing [1], rather than models. NRQCD can systematically factorize the nonperturbative sectors out and leave others perturbative QCD(pQCD) calculable. Up to now NRQCD achieves a great success in Born level and first order pQCD calculations [2], though some unsatisfactory still remains in phenomenological study.
Since the energy scales set by heavy quarks, the charm and bottom quarks, are moderate, the strong interaction in this regime is perturbative expandable, whereas does not converge quickly. The large discrepancy between leading order theoretical calculation [3] and experimental data [4,5], for example in double charmonium production at B-factory, can be indeed amended by the next-to-leading order(NLO) QCD corrections [6,7]. Higher order calculation itself is critical, challenging and meaningful theoretically, in addition to its phenomenological impact.
The processes of η c (η b ) exclusive production in accompany with a photon at electron-position colliders is right now an interesting topic, because it can be measured accurately and calculated precisely. Furthermore, naive estimation indicates that e + e − → γ + η c cross section overshoots that of e + e − → J/ψ + η c process by an order of magnitude, and hence e + e − → γ + η c (η b ) processes may be well measured in Belle II experiment soon. On theory side, the leading-order calculation and next-to-leading order corrections are ready [8][9][10], where the radiative corrections are negative.
In calculating high order corrections, the central issue is to calculate the multi-loop Feynman integrals, namely master integrals (MI). The first complete NNLO analytical calculations for quarkonium production and decay, the J/ψ(Υ) → e + e − and e + e − → J/ψ(Υ) processes, were achieved in Refs. [11] and [12] respectively, and a few years ago even the NNNLO result for Υ → e + e − process was obtained [13]. In recent years, at the NNLO level many calculations are performed for quarkonium production and decays [14][15][16]. It is interesting to note that in recently, the NNLO corrections to γγ * → η c (η b ) transition form factor and η c (η b ) → light hardrons were calculated numerically [17,18]. Numerical calculation is an unique and promising way for higher order radiative corrections, nevertheless right now it experiences the shortage of proper numerical packages, especially for the kinematics in physical region.
In this work, we are going to calculate analytically the next-to-next-to-leading order(NNLO) QCD corrections to exclusive η c (η b ) production associated with a photon in electron-positron collision, the multi-scale e + e − → γ + η c (η b ) processes. Numerical evaluation for BESIII and Belle II experiments will be presented. Note, since the kinematics of the master integrals in our calculation lie obviously in the physical region, due to the lack of corresponding numerical calculation package, to our knowledge a full numerical calculation of those integrals is still unrealistic and unreliable.
In our calculation, FeynArts [19] is employed to generate the corresponding Feynman diagrams and amplitudes up to NNLO, among them both massive and massless "light by light" diagrams are taken into account. Typical LO, NLO, and NNLO Feynman diagrams of the e + e − → γ + η c (η b ) processes are shown in Figure 1. We calculate the amplitude and reduce the loop integrals in full QCD, before applying the NRQCD projector to quarkonium stte. Up to any loop, the fermion chains in the amplitude of concerned process can be generally expressed as Here, F 1 is a function of kinematics, coupling, renormalization and factorization scales.
In above expression, the leptonic current is implied being factorized out, and the ellipses represent terms with the number of γ matrix less than 3. Of the e + e − → γ * → γ+η c (η b ) processes, only the first term in (1) contributes. In order to get F 1 , we multiply (1) with the projector and sum over the spinor helicities and trace over the dirac chains. Note that by means of this procedure, one can avoid the ambiguity in the definition of heavy quarkonium projector and γ 5 in the dimensional regularization scheme.
With the reduction of the obtained two-loop amplitudes, we find all related scalar integrals may be attributed to a set of 133 master integrals. Using the method of differential equation [21,22], and by virtue of the recent development in choosing basis properly [23], it was found that 86 of the 133 master integrals can be expressed as functions of Harmonic and Goncharov polylogarithms [20]. The calculation of massive multi-loop Feynman integrals is a challenging task due to the appearance of new mathematical structures, which cannot be reduced to the well-known multiple polylogarithms yet. Thanks to the very recent developments in differential equation technique in dealing with multiscale Feynman integrals [24][25][26][27], we are now able to calculate those massive multi-loop Feynman integrals of our concern beyond multiple polylogarithms, which are found can be classified into two categories, the elliptic sectors. One elliptic sector contains the two-loop sunrise integrals with full massive propagators and those integrals containing them as subtopologies, the other sector comprises the non-planar two-loop three-point integrals. The typical Feynman diagrams of those two elliptic sectors are shown in Figure 2.
The problems of massive two-loop sunrise integrals have been solved in Ref. [25], where the results are expressed in terms of one-fold integrals over complete elliptic integrals and polylogarithms. We find that with the basis chosen for sunrise integrals in Ref. [25] and the method for properly choosing basis proposed in Ref. [23], differential equations for all integrals containing massive sunrise integrals as sub-topologies can be reduced to a compact form, different from [23], that are ready to be solved recursively.
In our calculation, we solve them order by order up to weight four, the maximum for two-loop integrals. For the non-planar two-loop three-point integrals, we find that the corresponding homogenous second order differential equations can be transformed linearly to the exact form of the second order differential equation for complete elliptic integrals. Then we can readily obtain the homogeneous solutions and the full solution by solving the conventional second order differential equations. Details about the solutions for elliptic sectors will be presented elsewhere.
By taking the above procedures, all the 133 master integrals can be analytically expressed as Goncharov Polylogarithms, Chen's iterated integrals, and integrals over elliptic functions and polylogarithms. Note, in elliptic sectors two-fold iterated integrals are employed to express the integrals. In order to ensure our analytical results being correct, all of them have been checked against the numerical package yields [28,29] in both physical and non-physical regions. For each integral encountered, we have used at leat two different numerical methods to check. Up to NNLO, the cross section of e + e − → γ * → γ + η c (η b ) processes can be formulated as Here, the LO cross section is well-known and can be reproduced easily, i.e., with O 1 ( 1 S 0 ) being the NRQCD matrix element. The NLO correction coefficient c 1 was obtained in Refs. [9,10], and our calculation agree with them. It reads, The NNLO correction coefficient c 2 obtained in this work is a function of (s, m Q , µ r , µ Λ ), which can be formulated as in which the color factor C F = 4 3 , s represents the center-of-mass-system(CMS) energy squared, m Q denotes the heavy quark mass, µ r is the renormalization scale, and µ Λ for the factorization scale. The function f ( s m 2 Q ) is quite lengthy, containing multiple polylogarithms, iterative integrals and elliptic integrals, which is not suitable to show in a journal paper, but may be provided in electronic form upon request. To guarantee our result being reliable, we have also reproduced the NNLO correction result for η c → γγ process [14].
Before performing numerical calculation, we need first to fix the input parameters.
The number of active quarks n f is taken to be 4 for η c production, and 5 for η b process.
The charm quark mass lies in the region of 1.4 to 1.5 GeV, and the bottom quark mass in 4.7 to 4.8 GeV, to account for the heavy quark mass dependence of final result.
Of the strong coupling, we first choose the well-determined value at Z-pole, that is α s (M Z ) = 0.1184 [32], and then evolve it to the heavy quark mass scales by employing the program RunDec in four-loop accuracy [33]. The NRQCD matrix elements for heavy quarkonium are taken from fittings in Refs. [34,35], i.e., For γ + η c (1S) exclusive production in BESIII experiment, the renormalization scale µ r is set to be at 2m Q . It was found in Refs. [11,17] that the choice of factorization scale µ Λ = 1 GeV results in a better convergence in perturbative expansion rather than the conventional choice of µ Λ = m c for charmonium. We find this conclusion still holds for γ + η c exclusive production, and hence µ Λ = 1 GeV is taken in our numerical evaluation. While for bottomonium, the NRQCD factorization appears to work reasonably well, and the conventional choice of µ Λ = m Q is adequate to get a good convergence in pertubative expansion. and 4.60 GeV in data set corresponding to a total integrated luminosity of 4.6 fb −1 [36], and found no evidence of direct γ + η c (1S) production, which is in accord with the higher order QCD predictions. That is to say, the NNLO result also prefers the enhancement in e + e − → γ + η c process between 4.23 and 4.36 GeV being from exotic particle Y (4260) decay.
Quarkonium study in B-factory has an unique importance, where some critical measurement had ever been achieved. With the SuperKEKB/BelleII run in near future, it is reasonable to expect more from it on quarkonium physics, including the measurement on γ +η c (η b ) exclusive production. In Table I    dependence, even with the choice of factorization scale at µ Λ = 1 GeV. Whereas for γ + η b production, the NNLO corrections are moderate, the NNLO corrections reduce the renormalization scale dependence evidently. These observations imply that the perturbative expansion converges well for γ + η b production, but not for its counterpart in charm sector, which casts a shadow over the higher order corrections and further suggests that NRQCD factorization works more soundly for bottomonium sytem. In summary, we calculate the NNLO QCD corrections for γ + η c (η b ) exclusive production in electron-positron collision within the NRQCD formalism. The result tells that the NRQCD factorization still works well at the two-loop level. Thanks to the recently development in differential equation method, we find that all the master integrals can be analytically expressed in terms of Goncharov Polylogarithms, Chen's iterated integrals, and elliptic functions, which were mostly given in Ref. [20]. To the best of our knowledge, this is the first time that a NNLO calculation for multiscale process about heavy quarkonium has been achieved analytically. Numerical results indicate that the NLO and NNLO cross sections are both negative in charmonium and bottomonium production. In BESIII regime, the cross section is greatly suppressed with NNLO QCD correction, which agrees with the experimental measurement and implies the Y (4260) → γ + η c (1S) decay mode within data. Hopefully, in the forthcoming SuperKEKB run and Belle II measurement on γ + η c and γ + η b production, quarkonium production mechanism might be further elucidated.