Scale-fixed predictions for γ + ηc production in electron-positron collisions at NNLO in perturbative QCD

In the paper, we present QCD predictions for γ + ηc production at an electron-positron collider up to next-to-next-to-leading order (NNLO) accuracy without renormalization scale ambiguities. The NNLO total cross-section for e+ + e− → γ + ηc using the conventional scale-setting approach has large renormalization scale ambiguities, usually estimated by choosing the renormalization scale to be the e+e− center-of-mass collision energy s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document}. The Principle of Maximum Conformality (PMC) provides a systematic way to eliminate such renormalization scale ambiguities by summing the nonconformal β contributions into the QCD coupling αs(Q2). The renormalization group equation then sets the value of αs for the process. The PMC renormalization scale reflects the virtuality of the underlying process, and the resulting predictions satisfy all of the requirements of renormalization group invariance, including renormalization scheme invariance. After applying the PMC, we obtain a renormalization scale-and-scheme independent prediction, σ|NNLO,PMC ≃ 41.18 fb for s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document}=10.6 GeV. The resulting pQCD series matches the series for conformal theory and thus has no divergent renormalon contributions. The large K factor which contributes to this process reinforces the importance of uncalculated NNNLO and higher-order terms. Using the PMC scale-and-scheme independent conformal series and the Padé approximation approach, we predict σ|NNNLO,PMC+Pade ≃ 18.99 fb, which is consistent with the recent BELLE measurement σobs=16.58−9.93+10.51\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\sigma}^{\mathrm{obs}}={16.58}_{-9.93}^{+10.51} $$\end{document} fb at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} ≃ 10.6 GeV. This procedure also provides a first estimate of the NNNLO contribution.


Introduction
Processes involving the production of heavy quarkonium are important for testing Quantum Chromodynamics (QCD) as well as the effective theory of Nonrelativistic QCD (NRQCD) [1]. The framework of NRQCD factorization theory allows the non-perturbative dynamics involving the binding of the heavy quark-antiquark pair in quarkonium to be factored into universal NRQCD matrix elements which can be extracted and fixed via a global fitting of experiments involving heavy quarkonium production. The remaining 'hard' contribution involving higher momentum transfers is then perturbatively calculable. Thus reliable calculations of quarkonium production and its decay now appear viable.
The NRQCD approach has been successfully applied to a number of quarkonium processes, but many challenges and puzzles have remained. At present, most NRQCD results have been done at the next-to-leading order (NLO) level. The next-to-next-to-leading order (NNLO) and higher-order calculations are much more difficult. Thus it is important to test a variety of NNLO predictions before drawing any definite conclusion on the general applicability of NRQCD, especially since the K factors which appear in heavy quarkonium production processes can be very large.
The η c production via the process e + + e − → γ + η c is an important charmonium production process which can be precisely measured at high-energy, high-luminosity electronpositron colliders. For example, the Belle II experiment is expected to produce a sizable number of γ + η c events in the near future, which can be used to make precise comparisons with theoretical predictions. This heavy quarkonium production process has been calculated in pQCD up to NNLO level [3][4][5][6][7]. However, the NNLO calculation performed by Chen, Liang and Qiao [7] displays both a large K factor and large renormalization scale uncertainties. As a cross check, and also as important step forward, we shall recalculate the cross section of e + + e − → γ + η c at the NNLO level, make a detailed discussion on how to eliminate the unnecessary renormalization scale ambiguities, and present a first prediction of the NNNLO contribution.
A physical observable, corresponding to an infinite-order pQCD approximation, doesn't depend on the choice of renormalization scale. If the perturbative coefficients and the strong JHEP01(2021)131 coupling constant α s are not well matched at a fixed order, as is the case of conventional scale-setting approach in which the renormalization scale is simply guessed, one finds significant renormalization scale-and-scheme ambiguities; cf. the reviews [8][9][10]. Any dependence of pQCD prediction on the choice of renormalization scheme violates a fundamental principle of the renormalization group. In fact, predictions based on conventional scale-setting approach are even incorrect for Abelian theory -Quantum Electrodynamics (QED); the renormalization scale of the QED coupling constant α can be set unambiguously by using the Gell-Mann-Low method [11]. It is thus essential to use a rigorous scale-setting approach in order to achieve reliable and precise scale-and-scheme independent fixed-order pQCD predictions.
The Principle of Maximum Conformality (PMC) [12][13][14][15][16] provides a rigorous approach to renormalization scale setting, extending the BLM method [17] to all orders in pQCD. The purpose of the PMC is not to find an 'optimal' renormalization scale, but to achieve renormalization scale independent pQCD prediction. This can be achieved by systematically and rigorously determining the effective value of α s of the process based on the renormalization group equation (RGE), which is free of renormalization scale dependence. Following this way, the scheme-dependent non-conformal {β i }-terms 1 are eliminated in the pQCD series, which matches the corresponding conformal series. Thus after applying the RGE, the resulting pQCD prediction is also independent of the choice of the renormalization scheme. The commensurate scale relations [18] which relate PMC predictions for different observables among each other also ensure the scheme independence. Moreover, the PMC procedure reduces in the N C → 0 Abelian limit to the Gell-Mann-Low method [19]. Thus the PMC eliminates renormalization scale-and-scheme ambiguities simultaneously, satisfying the principles of renormalization group invariance [20,21]. In addition, since the {β i }-terms have been removed, the divergent renormalon terms like α n s β n 0 n! disappear and a convergent perturbative series can be naturally achieved.
In the paper, we shall adopt the PMC single-scale approach [22] for our analysis. The PMC scale can be interpreted as the effective momentum flow within the production process; its value displays stability and convergence with increasing order in pQCD, and any residual scale dependence due to unknown higher-order terms is highly suppressed.
The remaining parts of the paper are organized as follows: in section 2, we give the calculation technology for the total cross section of e + + e − → γ + η c up to NNLO level. We present the numerical results in section 3, and section 4 is reserved for the summary.

Calculation technology
In the following, we shall give a brief description on the calculation techniques, and then present the numerical expressions up to the NNLO level.
As the first step, the package FeynArts [23] is used to generate the Feynman diagrams for the γ +η c production via the process e + +e − → γ +η c and the corresponding amplitudes  up to NNLO level. There are totally 2 LO Feynman diagrams, 8 NLO Feynman diagrams and 120 NNLO Feynman diagrams for the considered process, whose typical ones are presented in figure 1. The evaluation involves the regular two-loop diagrams and the "light-by-light" (lbl) scattering diagrams. The lbl graphs are both gauge-invariant and free of ultra-violet (UV) and infrared (IR) divergences. Second, we employ the packages FeynCalc/FormLink [24,25] to deal with the trace over the Dirac and color matrices. And then with the help of the packages Apart [26] and FIRE [27], we conduct the partial fraction and the integration-by-parts reduction. Eventually, we end up with 6 independent one-loop Master Integrals (MIs) and 174 independent two-loop MIs. Those MIs can be evaluated through performing sector decomposition and subsequent numerical integrations with quadruple precision via the packages FIESTA/CubPack/HCubature [28][29][30]. This explain our main procedures, and some more descriptions can be found in ref. [31].
To implement the renormalization of the quark wave function and quark mass, we adopt the O(α 2 s )-order on-shell renormalization constants Z 2 and Z m , which expressions can be read from refs. [32,33]. For self-consistency, the strong coupling constant under the MS scheme is taken up to order-α 2 s . After the renormalization procedure, the entire UV divergences are eliminated. To deal with infrared (IR) divergence, before conducting the loop integration, we expand the integrand of the quark amplitude in power series of q 2 , where q is the relative momentum between the two heavy constituent quarks. For this performance, the Coulomb singularities are automatically removed, thus we only need to consider the soft divergence. It is noted that the NNLO amplitude contains a single IR pole with the very coefficient as anticipated from the following equations (2.10), (2.11), (2.12), as shown in the following illustrations, those IR terms can be explicitly factorized into the non-perturbative η c NRQCD matrix element according to MS prescription, with log µ Λ manifested as the remnant (µ Λ refers to the NRQCD factorization scale).
According to the NRQCD factorization formalism, the cross section for the process e + + e − → γ + η c up to NNLO in α s can be expressed as where a s ≡ α s /π, f (0) represents the LO short-distance coefficient (SDC), r 1 and r 2 correspond to the NLO and NNLO radiative corrections to the LO SDC, O 1 ( 1 S 0 ) denotes the NRQCD matrix element for the color-singlet and spin-singlet charmonium state. Since SDCs are insensitive to the nonperturbative hadronization effects, they can be deduced with the aid of the perturbative matching technique. That is, by replacing the physical η c meson with a fictitious onium composed of free cc pair, carrying the same quantum number as η c . After this replacement, we have Both the perturbative cross section σ(cc( 1 S 0 )) and the NRQCD matrix element are calculable, and we can solve all the three SDCs f (0) , r 1 and r 2 . Concretely, we can employ the aforementioned techniques to evaluate perturbative cross section σ(cc( 1 S 0 )) up to O(α 2 s )-order. Meanwhile, we are required to carry out the nonperturbative NRQCD matrix element at the two-loop level. Since the contributions from the potential region in the QCD calculation have been removed, to be consistent, we must remove this counterpart in our NRQCD evaluation. There still remain UV as well as IR divergences in NRQCD matrix element at the two loop level. The UV divergence can be renormalized in MS scheme, thus, at the lowest order in velocity expansion, we have which can be directly translated from refs. [34,35]. After hard work, we find that the remained divergence in σ(cc( 1 S 0 )) can be exactly cancelled by the IR divergence in eq. (2.3), which render the SDCs be free of any divergences, however develop a log µ Λ dependence.
In the following, we will directly present the SDCs and the cross section. The LO SDC reads Immediately, we obtain the cross section at O(α 0 s ) For future convenience, we reexpress eq. (2.1) as where the tree-level cross section σ (0) is

JHEP01(2021)131
where e c = +2/3 denotes the electric charge of charm quark, m c denotes the mass of charm quark, α is the fine structure constant, √ s is the center of mass energy, and O 1 ( 1 S 0 ) is the matrix element for η c . The NLO coefficient r 1 takes the form [4] where q = 4m 2 c /s and b = 1 − 4m 2 c /s. The NNLO calculation technology has been described in detail in ref. [36]. The explicit expressions for the NNLO coefficient r 2 is extremely lengthy. For convenience, we define a K factor, Numerical results for the K factor at the e + e − collision energy where µ Λ is the factorization scale, n l is the number of light flavors (u, d and s), and "lbl" denotes the contribution from the light-by-light Feynman diagrams. Though the lbl-terms are implicitly proportional to n l , they are free of ultra-violet (UV) divergence and are irrelevant to the running of α s , so they should be kept as conformal terms when applying the PMC (Thus the n l in lbl-terms is fixed to be 3). Using the above equations, the perturbative coefficients r i can be fixed. As a useful reference, we rewrite the K factors

JHEP01(2021)131
for those typical m c values with explicit color structures in the following (2.14) where C F = 4/3 and C A = 3 are SU c (3) color factors. In order to apply the PMC, the perturbative coefficients r i need to be divided into conformal terms and non-conformal terms [15,16], and for the present NNLO analysis, the coefficients for the perturbative series (2.6) can be written as where β 0 = 11−2/3n f with the active flavor numbers n f = n l +1 (n f = 4 for η c production). The coefficients r i,0 are the conformal coefficients and the r i,j =0 are non-conformal ones. By using the standard PMC scale-setting procedures, the NNLO total cross section (2.6) can be written as a conformal series, where Q stands for the PMC scale which is determined by requiring all non-conformal terms to vanish. The PMC scale can be fixed up to leading-logarithm (LL) accuracy by using the known NNLO pQCD series; i.e., wherer 2,1 = r 2,1 | µr=mc . The PMC scale Q can be regarded as the effective momentum flow of the process, since it is determined by using the RGE and the effective value of α s (Q ) has been fixed.

JHEP01(2021)131
Eq. (2.19) shows that the scale Q is independent of the choice of µ r . Since the conformal coefficients r i,0 are scale-independent, the PMC prediction is exactly free of µ r -dependence. Thus the conventional renormalization scale ambiguity is solved.
The unknown higher-order terms in the perturbative series of ln Q 2 /m 2 c can have some residual scale dependence [37]; however, this dependence is distinct from the conventional renormalization scale ambiguities, and it has both α s -power suppression and exponential suppression. In practice, most applications of the PMC published in the literature show that such residual scale dependence is rather small 2 [8,10].

Numerical results
To do the numerical calculation, we take α=1/130.7, α s (M Z )=0.1181 [39], √ s = 10.6 GeV, and the η c matrix element [40,41]. This value is fixed via a matching at the NLO level without factorization scale dependence. At present, the NNLO evolution equation of the matrix element is missing and we adopt the one-loop evolution equation [1] to fix the matrix element at different scales, there is thus factorization scale dependence for the NNLO total cross section (2.6). For definiteness, we have explicitly set the factorization scale of the matrix element as the usually adopted one, i.e. µ Λ = 1GeV. If setting its scale as another usual choice, e.g. µ Λ = m c , we find that the total cross-section shall only be slightly changed, e.g. less than 5% of the total cross-section. The required two-loop α s running is calculated by using the RunDec program [42].
Assuming m c = 1.5 GeV and the factorization scale µ Λ = 1 GeV, the total cross sections under conventional scale-setting approach (Conv.) for three typical choices of renormalization scale are  Table 1 shows how the magnitude of each loop term to the NNLO K factor changes under different choices of µ r . After applying the PMC, the scale dependence of each loop term is eliminated; however, the pQCD convergence is still poor. This is due to the fact that the conformal coefficient is large and dominates over the total NNLO coefficient r 2 , e.g. r 2 = −79.9 ± 7.4 where r 2,0 ≡ −83. 5   numerically that if one sets µ r = √ 2s, the conventional prediction gives total cross section and perturbative behavior which are close in comparison to the PMC prediction; thus √ 2s can be treated as the optimal renormalization scale of conventional prediction. Figure 2 shows how the net scale uncertainty varies when more loop terms have been included. Contrary to usual expectations, the net NNLO scale uncertainty of the conventional series is larger than that of the NLO series, since cancellations of terms at different orders are absent. In contrast, after applying the PMC, both the NLO and NNLO K factors are free of the renormalization scale ambiguities. In this sense, the scale-invariant PMC conformal series is extremely important for precise pQCD predictions.
After eliminating the renormalization scale uncertainty, there are other error sources such as the factorization scale, the charm quark mass m c , the value of α s (M Z ), and the η c matrix element O 1 ( 1 S 0 ) . The matrix element is an overall parameter, whose error can be determined separately. Thus, in the following, we shall only discuss uncertainties which come from µ Λ , m c and ∆α s (M Z ). When discussing the uncertainty of the PMC prediction from one parameter, the other parameters are set as their central values.
First, to discuss the factorization scale uncertainty, we adopt µ Λ = 1 GeV, m c and 2m c with m c = 1.5 GeV, respectively. The   and PMC scale-setting approaches, respectively. It shows that the K factor decreases with increasing factorization scale. Second, to discuss the m c uncertainty, we adopt m c = 1.4, 1.5, and 1.6 GeV, respectively. By setting µ Λ = 1 GeV, we obtain K PMC = 0.65, 0.64, 0.63, (3.6) accordingly. Figure 4 shows that the m c dependence under conventional and PMC scalesetting approaches, respectively. The K factor decreases with the increment of m c . Third, to discuss the ∆α s (M Z ) uncertainty, we adopt ∆α s (M Z ) = 0.1181±0.0011 [39], and we obtain K PMC = 0.65, 0.64, 0.63, (3.7) accordingly. Figure 5 shows that the ∆α s (M Z ) dependence under conventional and PMC scale-setting approaches, respectively. Recently, the BELLE Collaboration published their measured total cross-section for e + +e − → η c +γ at √ s = 10.58 GeV. There measured Born cross section σ B =11.3 +7.0+1.5 −6.6−1.5 fb [2], and by employing the following formula [43] σ B (s) = (1 + δ(s)) σ obs (s) we inversely obtain σ obs =16.58 +10.51 −9.93 fb. This value is smaller than the theoretical prediction (3.8). The central theoretical value shows a 2.3σ deviation from the data, which increases to 3.1σ when taking the usual choice of µ Λ = 1 GeV. If one takes the uncertainty of the η c matrix element into consideration, e.g.
−0.105 GeV 3 [40,41], we then have an additional uncertainty, +10.94 −10.34 fb, to the total cross section. This gives a slight overlap of the theoretical prediction with the measurement.
As a final remark, as shown by table 1, the poor pQCD convergence is due to the intrinsic nature of the present process, which cannot be improved even after applying the PMC. Thus it is necessary to know the contributions from unknown terms before we draw any definite conclusions. The conventional error estimate obtained by varying the guessed scale over a certain range cannot give a reliable prediction of the unknown terms, since it only partly estimates the non-conformal contribution but not the conformal one.
If one has a renormalization-scale independent conformal series, one can often obtain a reliable prediction of unknown higher-order contributions [44] with the help of the Padé JHEP01(2021)131 resummation [45][46][47]. The diagonal [1/1]-type Padé series is generally preferable for estimating the unknown contributions from a poor pQCD convergent series [48][49][50]; and for the present process, the poor convergence of PMC series is due to large conformal coefficients. Detailed procedures for obtaining a combined Padé +PMC prediction can be found in ref. [50]. By using the diagonal [1/1]-type PAA approximant, the predicted coefficient of the N 3 LO term is -2713.77, which results in an extra 38% suppression from the LO cross-section, leading to a NNNLO prediction in better agreement with the data; i.e., where µ Λ = 1 GeV and the other parameters are set to their central values. This indicates the importance of a strict NNNLO calculation, even though it would be much more difficult than the present NNLO calculation.

Summary
In this paper, we have presented a detailed study of the cross section for γ + η c production in electron-position collisions up to NNLO level. If one used the conventional procedure, the renormalization scale uncertainty for the NNLO total cross-section is estimated as 28% by varying the scale within the range of µ r ∈ [ √ s/2, 2 √ s]. However, after applying PMC scale setting, the conventional scale uncertainty is eliminated. Our NNLO prediction is, where the errors are squared average of the errors caused by varying µ Λ ∈ [1, 3] GeV, m c ∈ [1.4, 1.6] GeV, and ∆α s (m Z ) = ±0.0011. Among the uncertainties from the other input parameters, the factorization scale error is the largest. The central value of the NNLO cross section deviates substantially from the measured data. The poor pQCD convergence of the series indicates the importance of uncalculated NNNLO terms for this process. An initial estimate of the NNNLO terms with the help of the PMC and Padé resummation has been given in section 3; the magnitude of the NNNLO contribution is sufficiently large to explain the data. Even though we need more accurate data to confirm the theoretical results, this application of the PMC shows the importance of a correct renormalization scale setting for a reliable pQCD prediction.

JHEP01(2021)131
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.