NNLO QCD corrections to γ + η c ( η 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 η c ( η 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 coeﬃcients 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, integrals over polylogarithms and complete elliptic integrals, 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 − → γ + η c ( η 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.


JHEP01(2018)091
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, 1 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][4][5] and experimental data [6,7], for example in double charmonium production at B-factory, can be indeed amended by the next-to-leading order(NLO) QCD corrections [8,9]. 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 electronposition 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 [10][11][12], 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. [13] and [14] respectively, and a few years ago even the NNNLO result for Υ → e + e − process was obtained [15]. In recent years, at the NNLO level many calculations are performed for quarkonium production and decays [16][17][18]. 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 [19,20]. 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 JHEP01(2018)091 corresponding numerical calculation package, to our knowledge a full numerical calculation of those integrals is still unrealistic and unreliable.
In our calculation, FeynArts [21] 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. 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 summing over the spinor helicities and trace over the dirac chains. Here, s denotes the center of mass energy squared and m Q the heavy quark mass. 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 [22,23], and by virtue of the recent development in choosing basis properly [24], it was found that 86 of the 133 master integrals can be expressed as functions of Harmonic and Goncharov polylogarithms [25]. 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 [26][27][28][29][30], 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 two two-loop sunrise integrals with full massive propagators and those integrals containing them as subtopologies, the other sector comprises the two non-planar two-loop three-point integrals. The typical Feynman diagrams of those two elliptic sectors are shown in figure 2.
In the first sector, there are 2 massive two-loop sunrise integrals and 37 integrals containing massive sunrise integrals as subtopologies. The massive sunrise integrals (E 1 , E 2 ) are shown in figure 2. Note, the problem of massive two-loop sunrise integrals was ever solved in ref. [27], where the results are expressed in terms of one-fold integrals over complete elliptic integrals and polylogarithms. With the basis chosen in ref. [27] and the method proposed in ref. [24], we find the differential equations for all integrals containing massive sunrise integrals can be reduced to compact forms, which are different from [24] and will be presented in detail elsewhere [31].
To be more specific, in our calculation the first sector with 39 integrals E i (i = 1 . . . 39), including 2 massive two-loop sunrise integrals (E 1 , E 2 ) which are already known [27], is treated as subtopologies. With proper basis selection for A i (i = 3 . . . 39), in linear combination of master integrals E i (i = 3 . . . 39), differential equations for A i (i = 3 . . . 39) can be expressed as Here, A is a 37-dimensional basis vector that contains integrals of E i (i = 3 . . . 39); F is a 86-dimensional basis vector that had been calculated in our previous work [25]; W and Y are 37 × 37 and 37 × 86 matrices, respectively; Q i (i = 1, 2) are also 37-dimensional basis vectors; and A 1 is a linear function of the massive sunrise integrals (E 1 , E 2 ), write as

JHEP01(2018)091
with F 1 being the first integral of vector F [25]. Note that (W, Y, Q 1 , Q 2 ) are rational functions of kinematics and hence are free of ǫ. The ǫ 0 term in Q 2 appears only in the inhomogeneous part, and is irremovable. The homogenous parts of the differential equation (3) have been simplified, and the dimensional regularization parameter ǫ in homogeneous parts is factored out from kinematics. With the known inhomogeneous terms containing A 1 [27], the differential equation (3) can be solved recursively in terms of ǫ.
For the second sector, which comprises two non-planar two-loop three-point integrals (B 1 , B 2 ) shown in figure 2, the second order homogenous differential equation for g = (s − 4m 2 Q ) 2 B 1 can be formulated as To obtain the homogeneous solution for above equation, we make a variable transformation . The homogeneous equation for g then turns to Two homogeneous solutions (g 1 (v), g 2 (v)) for (6) are readily obtained and can be expressed as the complete elliptic integrals of the first kind with K(x) being defined as By means of variation of constants, the inhomogeneous solution may be obtained. For details of basis chosen, boundary condition determination, analytical continuation, and solutions for two elliptic sectors, readers are recommended to refer to ref. [31]. After taking the above procedures, all 133 master integrals can be analytically expressed in terms of the Goncharov Polylogarithms, integral over polylogarithms and elliptic functions. Note, in elliptic sectors two-fold iterated integrals are employed in expressing the integrals. Goncharov Polylogarithms can be numerically evaluated with the GINAC implementation [32,33], and can also be transformed into functions of ln, Li n and Li 22 up to weight four, by using the method described in ref. [34]. The numerical evaluation of integrals over polylogarithms and integrals over elliptic functions and polylogarithms can be simply performed by software MATHEMATICA. In order to ensure our analytical results being correct, all of them have been checked against the numerical package yields [35,36], in both physical and non-physical regions. For each integral encountered, we have used at leat two different numerical methods to check it up.
About renormalization, for quark wave function and mass, Z 2 and Z m respectively, we adopt the on-shell scheme [37,38], while for the strong coupling constant the MS scheme

JHEP01(2018)091
is taken up to one-loop order. All ultraviolet divergences are then removed after taking the renormalization procedure. The remaining infra divergences are process independent, and can be canceled by the corrections of quarkonium NRQCD matrix elements under MS scheme, with the results depending on the NRQCD factorization scale ln Λ µ .
Up to NNLO, the cross sections 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. [11,12], and our calculation agree with them. It writes 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 [16].
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 between 1.4 and 1.5 GeV, which agrees with what adopted in [11,12] for η c exclusive production. The bottom quark mass first takes the value in MS scheme, i.e. m b (m b ) = 4.18 GeV in PDG [39], and then converts to the 2-loop pole mass m b = 4.78 GeV [19,20]. In order to evaluate the bottom-quark-mass dependence of the final result, in our numerical calculation the bottom quark mass varies from 4.7 to 4.8 GeV. Of the strong coupling, we first choose the well-determined value at Z-pole, that is α s (M Z ) = 0.1184 [39], and then evolve it to the heavy quark mass scales by employing the program RunDec [40] in four-loop accuracy. The NRQCD matrix elements for heavy quarkonium are taken from fittings in refs. [41,42], 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. [13,19] 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 adopted 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.
In figure 3 we show respectively the LO, NLO, and NNLO cross sections versus the CMS energy, which varies in between 4.0 and 5.5 GeV. From the figure we can read that NLO and NNLO corrections are both negative. With the NNLO correction, the production rate in this region is greatly suppressed, about 3 times reduced from the LO result, and the NNLO cross section varies more smoothly than LO and NLO results with the change of CMS energy, as expected. Recently, BESIII collaboration performed a measurement on e + e − → γ + η c process with CMS energy between 4.01 and 4.60 GeV in data set corresponding to a total integrated luminosity of 4.6 fb −1 [43], 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 1     be a surprise. The NLO and NNLO corrections at B-factory are still negative, and the NNLO cross sections are some half of the LO predictions. Note that though NNLO corrections for e + e − → γ + η c is highly suppressed, it still surpasses the NLO cross section of e + e − → J/ψ+η c process. Therefore, the measurement of exclusive processes e + e − → γ +η c and e + e − → η c + l + + l − , l denoting leptons, are highly desired to further examine the quarkonium production mechanism and high-order QCD calculation reliability.
In figures 4 and 5 we show the renormalization scale dependence of NLO and NNLO Kfactors in γ +η c and γ +η b production at B-factory respectively. It is noticeable that though the NLO and NNLO K-factors are all negative, for γ +η c production the NNLO corrections are large and give no help to diminish the final result renormalization scale 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, integrals over polylogarithms and complete elliptic integrals. The calculation details of all the master integrals can be found in refs. [25,31]. 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.