The heavy quarkonium inclusive decays using the principle of maximum conformality

The next-to-next-to-leading order (NNLO) pQCD correction to the inclusive decays of the heavy quarkonium $\eta_Q$ ($Q$ being $c$ or $b$) has been done in the literature within the framework of nonrelativistic QCD. One may observe that the NNLO decay width still has large conventional renormalization scale dependence due to its weaker pQCD convergence, e.g. about $(^{+4\%}_{-34\%})$ for $\eta_c$ and $(^{+0.0}_{-9\%})$ for $\eta_b$, by varying the scale within the range of $[m_Q, 4m_Q]$. The principle of maximum conformality (PMC) provides a systematic way to fix the $\alpha_s$-running behavior of the process, which satisfies the requirements of renormalization group invariance and eliminates the conventional renormalization scheme and scale ambiguities. Using the PMC single-scale method, we show that the resultant PMC conformal series is renormalization scale independent, and the precision of the $\eta_Q$ inclusive decay width can be greatly improved. Taking the relativistic correction $\mathcal{O}(\alpha_{s}v^2)$ into consideration, the ratios of the $\eta_{Q}$ decays to light hadrons or $\gamma\gamma$ are: $R^{\rm NNLO}_{\eta_c}|_{\rm{PMC}}=(3.93^{+0.26}_{-0.24})\times10^3$ and $R^{\rm NNLO}_{\eta_b}|_{\rm{PMC}}=(22.85^{+0.90}_{-0.87})\times10^3$, respectively. Here the errors are for $\Delta\alpha_s(M_Z) = \pm0.0011$. As a step forward, by applying the Pad$\acute{e}$ approximation approach (PAA) over the PMC conformal series, we obtain approximate NNNLO predictions for those two ratios, e.g. $R^{\rm NNNLO}_{\eta_c}|_{\rm{PAA+PMC}} =(5.66^{+0.65}_{-0.55})\times10^3$ and $R^{\rm NNNLO}_{\eta_b}|_{\rm{PAA+PMC}}=(26.02^{+1.24}_{-1.17})\times10^3$. The $R^{\rm NNNLO}_{\eta_c}|_{\rm{PAA+PMC}}$ ratio agrees with the latest PDG value $R_{\eta_c}^{\rm{exp}}=(5.3_{-1.4}^{+2.4})\times10^3$, indicating the necessity of a strict calculation of NNNLO terms.

The next-to-next-to-leading order (NNLO) pQCD correction to the inclusive decays of the heavy quarkonium ηQ (Q being c or b) has been done in the literature within the framework of nonrelativistic QCD. One may observe that the NNLO decay width still has large conventional renormalization scale dependence due to its weaker pQCD convergence, e.g. about ( +4% −34% ) for ηc and ( +0.0 −9% ) for η b , by varying the scale within the range of [mQ, 4mQ]. The principle of maximum conformality (PMC) provides a systematic way to fix the αs-running behavior of the process, which satisfies the requirements of renormalization group invariance and eliminates the conventional renormalization scheme and scale ambiguities. Using the PMC single-scale method, we show that the resultant PMC conformal series is renormalization scale independent, and the precision of the ηQ inclusive decay width can be greatly improved. Taking the relativistic correction O(αsv 2 ) into consideration, the ratios of the ηQ decays to light hadrons or γγ are: R NNLO ηc |PMC = (3.93 +0. 26 −0.24 ) × 10 3 and R NNLO

I. INTRODUCTION
The heavy quarkonium, being a common bound state of Quantum Chromodynamics (QCD) which consists of a pair of heavy quark and antiquark, has been continuously studied either experimentally or theoretically. For example, the η c decays to light hadrons and γγ have been measured by the BES III detector [1], which also gives the evidence for η c → γγ. In year 2018, the Partical Data Group (PDG) [2] issued the most recent information for heavy quarkonium from various measurements. The heavy quarkonium processes involve both perturbative and non-perturbative effects, and these processes are important tests of QCD factorization theories.
The nonrelativistic QCD (NRQCD) factorization theory provides us an effective framework to deal with heavy quarkonium processes [3], which factorizes the pQCD approximant into the non-perturbative but universal longdistance matrix elements (LDMEs) and the perturbatively calculable short-distance coefficients. More explicitly, the short-distance coefficients can be expressed as a perturbative series over the strong coupling constant (α s ) and the relative velocity of the heavy quarkonium (v); e.g. the decay width of the heavy quarkonium η Q * Electronic address: yuq@cqu.edu.cn † Electronic address: wuxg@cqu.edu.cn ‡ Electronic address: zengj@cqu.edu.cn § Electronic address: hxud@cqu.eud.cn ¶ Electronic address: yuhm@cqu.edu.cn can be factorized as the following form [4] where Q = c or b, respectively, F 1 ( 1 S 0 ) and G 1 ( 1 S 0 ) are short-distance coefficients. The perturbative series is arranged by the velocity scaling rule. The NRQCD matrix elements η Q |O 1 ( 1 S 0 )|η Q and η Q |P 1 ( 1 S 0 )|η Q refer to the possibility of observing the specific color and angular-momentum state of the QQ pair, where the 4fermion operators P 1 and O 1 produce or annihilate a QQ pair in the Fock state |η Q . It is convenient to compare the following ratio such that to avoid/suppress the uncertainties from the bound state parameters, which is defined as where LH stands for light hadrons. The next-to-leadingorder (NLO) QCD corrections to the perturbative part of η c → LH and η c → γγ have been done in Refs. [5,6], and the relativistic O(α s v 2 )-corrections have also been done in Refs. [7][8][9]. Recently, the next-to-next-to-leading order (NNLO) QCD corrections, including the relativistic corrections, have been finished by Feng et al. [10] and Brambilla et al. [11], which also show that factorization scale dependence of the R-ratio cancels at the NNLO level. Thus we are facing the chance of achieving a more accurate pQCD prediction on the R-ratio.
There is still large renormalization scale dependence for the NNLO pQCD approximant of the R-ratio under conventional scale-setting approach. That is, conventionally, people adopts the "guessed" typical momentum flow of the process such as m Q as the renormalization scale with the purpose of eliminating the large logarithmic terms or minimizing the contributions of the higher-order loop diagrams [12,13]. Such a naive treatment, though conventional, directly violates the renormalization group invariance [14] and does not satisfy the self-consistency requirements of the renormalization group [15], leading to a scheme-dependent and scale-dependent less reliable pQCD prediction in lower orders.
To eliminate the unwanted scale and scheme ambiguities caused by using the "guessed" scale, the principle of maximum conformality (PMC) scale-setting approach has been suggested [16][17][18][19][20]. The purpose of PMC is not to find an optimal renormalization scale but to find the effective coupling constant (whose argument is called as the PMC scale) with the help of renormalization group equation (RGE). By using the PMC, the effective coupling constant is fixed by using the β-terms of the pQCD series, which are arranged by the general degeneracy relations in QCD [21] among different perturbative orders. Since the effective coupling is independent to the choice of the renormalization scale, thus solving the conventional scale ambiguity. Furthermore, after using the PMC to fix the α s behavior, the remaining coefficients of the resultant series match the series of the conformal theory, leading to a renormalization scheme independent prediction. At the present, the PMC approach has been successfully applied for various high-energy processes.
For the PMC multi-scale method [16,19], we need to absorb different types of {β i }-terms into the coupling constant via an order-by-order manner. Different types of {β i }-terms as determined from the RGE lead to different running behaviors of the coupling constant at different orders, and hence, determine the distinct PMC scales at each perturbative order. The precision of the PMC scale for higher-order terms decreases at higherand-higher orders due to the less known {β i }-terms in those higher-order terms. The PMC multi-scale method has thus two kinds of residual scale dependence due to the unknown perturbative terms [22], which generally suffer from both the α s -power suppression and the exponential suppression, but could be large due to possibly poor pQCD convergence for the perturbative series of either the PMC scale or the pQCD approximant [13] 1 .
Lately, in year 2017, the PMC single-scale method [24] has been suggested to suppress those residual scale dependence. The PMC single-scale method replaces the individual PMC scales at each order by an overall scale, which effectively replaces those individual PMC scales derived under the PMC multi-scale method in the sense of a mean value theorem. The PMC single scale can be regarded as the overall effective momentum flow of the process; it shows stability and convergence with increasing order in pQCD via the pQCD approximates. It has been demonstrated that the PMC single-scale prediction is scheme independent up to any fixed order [25], thus its value satisfies the renormalization group invariance. Moreover, it has been found that the first kind of residual scale dependence still suffers α s and exponential suppression and the second kind of residual scale dependence can be eliminated by applying the PMC single-scale method; thus the residual scale dependence emerged in PMC multi-scale method does greatly suppressed. In the following, we shall adopt the PMC singlescale method to do our discussions.
The remaining parts of the paper are organized as follows. In Sec.II, we give the calculation technology for the R-ratio under the PMC single-scale method. In Sec.III, we give the NNLO numerical results for η c and η b Rratios and the approximate NNNLO predictions. Sec.V is reserved for a summary.

II. CALCULATION TECHNOLOGY
The pQCD approximant for the R-ratio has been calculated up to NNLO level under the MS-scheme where those short-distance coefficient All the coefficients f 1 , f 2 , g 1 , f γ , and g γ can be read from Refs. [8,10,26], the logarithm L = ln µ 2 r 4m 2 and the non-logarithmic constant and we adopt v 2 c ≃ 0.430/m 2 c and v 2 b = −0.009 [27,28]. The factorization scale µ Λ = m and the active quark flavor n f = n L + n H with n H = 1. Here n L = 3 for η c decay and n L = 4 for η b decay, and nL i=1 e 2 i /e 2 Q sums up the fractional charges of the light flavors involved in those two decays. The fractional charge e Q equals to 2 3 for η c and − 1 3 for η b , respectively. Before applying the PMC to the R-ratio, we first transform the MS-scheme pQCD series (3) into the one under the minimal momentum space subtraction scheme (mMOM) [29,30]. This transformation avoids the confusion of distributing the n f -terms involving the threegluon or four-gluon vertexes into the β-terms [31][32][33][34]. This transformation can be achieved by using the relation of the running coupling under the mMOM-scheme and the MS-scheme, e.g.
where a s (µ r )=α s (µ r )/(4π), and under the Landau gauge (ξ = 0), the first three coefficients D 1 , D 2 and D 3 are [35] where ζ 3 and ζ 5 are usual Riemann zeta functions. Then, Eq. (3) can be rewritten as ) are short notations whose explicit forms can be found in Ref. [13], and r i,0 are conformal coefficients and r i,j =0 are non-conformal coefficients which can be related to the known r i coefficients by using the standard PMC formulae and by using the degenerated relations among different orders [19,20], i.e. r 1 = r 1,0 , r 2 = r 2,0 + pr 2,1 β 0 , The usual β-function is defined as where the β i -functions under the MS-scheme up to five loop level can be found in Refs. [36][37][38][39][40][41][42][43][44]. The β 0 -and β 1 -functions are scheme independent, and the β i>1functions for the mMOM-scheme can be related to the MS-scheme ones via the relation, The β i -functions under the mMOM-scheme up to four loop level can be found in Ref. [35]. The coefficients r i,j are functions of the logarithm ln(µ 2 r /4m 2 Q ). By using the RGE, these coefficients can be expressed as where the coefficientsr i,j = r i,j | µr =2mQ and the combination coefficients C k j = j!/k!(j − k)!. Substituting Eq. (20) into Eq. (18) and by requiring all the RGE-involved non-conformal terms to zero, one can determine the effective running coupling of the process and hence the optimal scale Q ⋆ of the process, e.g. i+j≤3 i≥1,j≥1,0≤k≤j Thus we obtain If we have known the NNNLO pQCD series, the PMC scale Q ⋆ can be fixed up to next-next-to-leading-log (NNLL) accuracy, e.g. where T 1 = (p + 1)(r 2,0r2,1 −r 1,0r3,1 ) pr 2 1,0 T 2 = (p + 1) 2 (r 1,0r2,0r3,1 −r 2 2,0r2,1 ) + p(p + 2)(r 1,0r2,1r3,0 −r 2 1,0r4,1 ) p 2r3 Using the present known NNLO pQCD seriers, we can fix the PMC scale Q ⋆ only up to the NLL accuracy. One may observe that the effective scale Q ⋆ is explicitly independent of the choice of the renormalization scale µ r at any fixed order, thus there is no renormalization scale ambiguity for the PMC prediction R| PMC . Therefore, the precision of the pQCD approximant can be greatly improved by using the PMC.
It is helpful to estimate how the uncalculated higherorder terms contribute to the pQCD series. Many attempts have been tried in the literature, all of which are based on the known perturbative terms. The usual er-ror estimate obtained by varying the scale over a certain range is unreliable, since it only partly estimates the non-conformal contribution but not the conformal one. Moreover, we should point out that predictions using the conventional renormalization scale-dependent pQCD series (3) cannot be reliably used for the purpose; It could have negligible net renormalization scale dependence for the whole pQCD approximant by including enough high-order terms due to the large cancellation among the scale-dependent terms at various orders, but the large scale dependence for each perturbative term cannot be eliminated. On the contrary, the PMC predictions are renormalization scheme-and-scale independent, highly precise values at each order can thus be achieved. As has been pointed out by Ref. [45], by using the renormalization scheme-and-scale independent conformal series (21), one can reliably predict how the uncalculated NNNLO-terms contribute to the R pQCD series by using the Padé approximation approach (PAA) [46][47][48].
The PAA was introduced for estimating the (n + 1) thorder coefficient in a given n th -order perturbative Talyor series and feasible conjectures on the likely high-order behaviour of the series. It has been previously demonstrated their applicability to the QCD problems with the help of some resummation methods [49][50][51][52]. For the Padé approximation of a general pQCD approximant ρ n (Q), its [N/M ]-type form can be defined as where p is the starting α s -order, and the parameter M ≥ 1, N +M +1 = n. For the present case of R-ratio, we have n = 3. There is a one-to-one correspondence between coefficients C i in Eq. (27) and the conformal coefficientŝ r i,0 in Eq. (21). The one-order higher coefficient C n+1 can be expressed by those coefficients b i∈[0,N ] and c j∈ [1,M] which can also be related to the known coefficients C i . More explicitly, if [N/M ] = [n − 2/1], we have In the following, we adopt the PMC series (21) together with the PAA to estimate the NNNLO contribution to the R-ratio. It has been found that for the divergent pQCD series, such as the conventional pQCD series which has renormalon divergence (the dominant factor for the n th -order coefficient is proportional to n!β n 0 ), the diagonal PAA series is preferable [53,54]; And if the pQCD series is much more convergent, such as the PMC conformal series which is free of renormalon divergence, the preferred PAA series should be consistent with the GM-L method [45].
More explicitly, for the present considered R-ratio, the predicted magnitude of the NNNLO-term, is

III. NUMERICAL RESULTS
To do the numerical calculation, we take the quark pole mass from PDG [2]: the c-quark pole mass m c = 1.67 GeV and the b-quark pole mass m b = 4.78 GeV.  As mentioned above, we shall calculate the R-ratio under the mMOM scheme, and the Landau gauge is adopted 2 . The asymptotic scale Λ n f QCD for the mMOM scheme is fixed by using the α s -value at the reference point, α s (M Z ) = 0.1181 ± 0.0011 [2] and by using the three-loop RGE. The results are presented in Table I   We present the R ηc and R η b ratios up to NNLO level as a function of µ r under conventional and PMC scalesetting approaches in Figs. (1, 2). Figs. (1, 2) show that under the conventional scale-setting approach, there are convex behaviors for R ηc and R η b ratios within the scale range of [m Q , 4m Q ], whose peak are ∼ 1.6m c for R ηc and ∼ 1.7m b for R ηc . Contributions from each 2 A detailed discussion on the gauge dependence of the mMOM scheme before and after applying the PMC is in preparation [55]. loop term for R ηc and R η b up to NNLO level are presented separately in Tables II and III, where the errors are for µ r ∈ [m Q , 4m Q ]. Tables II and III show that under conventional scale-setting approach, the renormalization scale dependence for the total NNLO prediction becomes small, whose magnitude is ∼ +4%

−34%
for R ηc , and ∼ +0.0 −9% for R η b . Such a smaller scale dependence for a NNLO level prediction is caused by the cancellation of scale dependence among different terms, and the scale dependence for each loop term is still very large for conventional series. For example, for the case of R ηc , the scale errors are −40%

−1092%
for the NNLO-term, respectively. On the other hand, by fixing the effective coupling α s (Q ⋆ ) of the process with the help of PMC, the renormalization scale dependence for each loop term and hence the total NNLO prediction can be eliminated simultaneously. Thus by applying the PMC, a more accurate pQCD prediction without renormalization scale dependence can be achieved; Basing on the scale-invariant PMC series, it is helpful to predict the unknown NNNLO contribution, as shall be done in next subsection.
Tables II and III indicate that under conventional series, even though the NNLO-terms are highly scale dependent, their magnitudes sound more convergent than the PMC series, which are due to accidentally cancellation between the large conformal terms and the divergent non-conformal terms at the NNLO level. This is different from previous PMC examples whose pQCD series converges much more quickly than conventional pQCD series due to the elimination of divergent renormalon terms (and also due to the reason of the magnitude of the conformal terms are usually moderate). By using the PMC, the RGE-involved non-conformal terms have been elimi-nated and have been adopted to fix the renormalization scale-independent strong coupling constant of the process, α s (Q ⋆ ). The resultant pQCD series is conformal and scheme independent, one can thus conclude that the PMC series shows the intrinsic property of the pQCD approximant and shows the correct convergent behavior of the pQCD series. At present, the PMC scales for R ηc and R η b can be determined up to NLL-accuracy: and ln Q 2 which are at the order of O(m Q ). The second lines of those two equations show that the second terms are about 35% and 6% of the first terms, indicating the perturbative series of ln Q 2 ⋆ /Q 2 has good convergence at the NLL level for both cases. The effective PMC scale Q ⋆ is physical, which is renormalization scale independent and determines the correct value of the strong running coupling and hence the correct momentum flow of the process. The heavy quark mass (m Q ) provides a natural hard scale for the heavy quarkonium decays into light hadrons or photons, and O(m Q ) is usually chosen as the renormalization scale. The PMC scale-setting approach provides a reasonable explanation for this conventional "guessing" choice.
B. The PAA prediction of the contribution from the uncalculated NNNLO-terms Table II shows that the NNLO PMC prediction of R ηc is 3.93 × 10 3 , which is only ∼ 65% of the NLO PMC prediction, R ηc = 6.09 × 10 3 [31], and is also smaller than the PDG central value, . This is due to the large negative NNLO-term 3 , as shown by Table II. Because the pQCD series of R ηc shows a slowly convergent behavior, the facts of the previous NLO prediction agrees with the data and the present NNLO prediction does not agree surely do not indicate the failure of the PMC or the pQCD factorization. We should at least know the magnitude of the NNNLO-term before drawing any definite conclusions.
A strict NNNLO calculation is unavailable in near future due to its complexity. Because the conventional pQCD series cannot be adopted for a reliable prediction, since its known terms are both scale dependent and scheme dependent, and in the following, we shall give an approximation of NNNLO prediction by applying the PAA to the PMC scheme-and-scale independent conformal series.   Table IV shows that if taking the approximate NNNLO-term into consideration, the R NNNLO ηc -ratio agree well with the PDG value within errors. The Rη c -ratio under various approaches. "EC" is the exact prediction by using the known NLO or NNLO PMC series, and "PAA+PMC" is the PAA prediction by using the PMC NNLO series. The error of "PAA+PMC" at the NNNLO level (N=4) is caused by ∆αs(M 2 Z ) = ±0.0011. As comparisons, the PDG value, the NNLO predictions (N=3) of Feng et al. [10], NNA [11] and BFG [11] are also presented.
the η c,b decays up to NNLO level in the large-n f limit by further using the bubble chain resummation: I) Using the naive non-Abelianization resummation [56], they obtained R ηc (NNA) = (4.28 +1. 38 −0.72 ) × 10 3 and R η b (NNA) = (23.2 +0.8 −0.9 ) × 10 3 ; II) Using the background field gauge resummation (BFG) [57], they obtained R ηc (BFG) = (3.39 +0.61 −0.64 ) × 10 3 and R η b (BFG) = (24.1 +0.9 −1.0 ) × 10 3 . We should point out that even though those two predictions are consistent with the PDG value, they only partly resum the large renormalon-terms with the purpose of improving the pQCD convergence along [11], which however by using the guessed scale cannot get the correct magnitude of the running coupling and there are still large renormalization scale errors. Thus those predictions cannot be treated as precise pQCD predictions.
Our present PMC predictions are based on the NNLO fixed-order result of Ref. [10], which includes the important relativistic O(α s v 2 )-contribution and has negligible factorization scale dependence. We present the R ηc -ratio and R η b -ratio under various approaches in Figs. (3, 4). "EC" is the exact prediction by using the known NLO or NNLO PMC series, and "PAA+PMC" is the PAA prediction by using the PMC NNLO series. The PDG value, the NNLO predictions of F.Feng [10], NNA [11] and BFG [11] are also presented as comparisons. Different to the previous theoretical predictions whose uncertainties are mainly caused by the renormalization scale dependence, there is no renormalization scale dependence in PMC prediction, and we give a prediction of the error of "PAA+PMC" approach at the NNNLO level by taking ∆α s (M 2 Z ) = ±0.0011.

IV. SUMMARY
In the paper, we have studied the R ηQ -ratio up to NNLO level. By using the RGE, the scale-independent α s -value can be achieved at any perturbative order, and the correct momentum flow of the process can be determined up to NLL accuracy; and because of the eliminating of the RGE-involved {β i }-terms, the resultant pQCD series is conformal and scheme independent. Thus we achieve an accurate NNLO fixed-order prediction which is free of conventional renormalization scheme and scale dependence. Fig.3 shows that the NNLO PMC prediction is somewhat smaller than the PDG value. Table II shows that the resultant PMC conformal series of R ηc | PMC has a slowly convergent behavior, thus before drawing definite conclusions, it is important to make a prediction on the contribution of the uncalculated NNNLO-terms. By using the Pade approximation approach, we give the approximate NNNLO predictions: R ηc | PMC = 5.66 +0.65 −0.55 × 10 3 and R η b | PMC = 26.02 +1.24 −1.17 × 10 3 . The approximated NNNLO PMC prediction of R ηc -ratio agrees with the PDG value within errors, indicating the necessity of finishing a strict NNNLO calculation.