Inclusive decays of $\eta_c$ and $\eta_b$ at NNLO with large $n_f$ resummation

Based on the nonrelativistic QCD factorization theorem, we resum QCD corrections to the inclusive decay rate of $\eta_c$ and $\eta_b$ in the large-$n_f$ limit using bubble chain resummation. By employing dimensional regularization, we show explicitly the cancellation of the infrared renormalon ambiguity in the factorization formula at leading order in $v$ in the large-$n_f$ limit, where $v$ is the typical heavy quark velocity inside the meson. We also make predictions of the ratio of the inclusive decay rate to the decay rate into two photons. By comparing our results with a fixed-order calculation we conclude that resummation of QCD corrections is crucial in making an unambiguous prediction. We also find significant corrections beyond the large-$n_f$ limit for the decay of $\eta_c$, which may imply that QCD corrections need to be resummed beyond the large-$n_f$ limit to make an accurate prediction of the decay rate.


I. INTRODUCTION
Nowadays there is much experimental effort devoted to investigating the nature of heavy quarkonium states. Precision measurements of the properties of heavy quarkonia such as their masses, decays or transitions can be done at future and ongoing experiments like Belle II, BESIII, and LHCb at CERN. A good understanding of the nature of heavy quarkonium states is essential in exploring other processes that involve heavy quarkonia such as their production in high energy collisions. Following the substantial success of the B factories [1], the Belle II experiment at KEK in Japan is going to collect 50 times more data of what Belle obtained and will be able to investigate the properties of pseudoscalar bottomonium η b state and the corresponding charmonium η c state.
One of the basic observables regarding heavy quarkonium is its inclusive decay rate. Theoretical predictions of the inclusive decay rate of η c and η b have long been pursued using nonrelativistic QCD (NRQCD), which is an effective theory providing a factorization formalism that separates the perturbative short-distance contributions from the nonperturbative long-distance ones [2]. The NRQCD factorization formula for the decay rate of a heavy quarknoium is a sum over products of NRQCD long-distance matrix elements (LDMEs) and the corresponding short-distance coefficients (SDCs). The NRQCD power counting attributes to the LDMEs a specific scaling with v, where v is the relative velocity between the heavy quark and the heavy antiquark in the quarkonium. The SDCs may be computed in perturbation theory. Therefore, the sum in the NRQCD factorization formula is an expansion in powers of v and α s .
There are two major obstacles in achieving theoretical predictions of the decay rates of η c and η b with high accuracy. First, even though the expressions for the SDCs that contribute to the decay rates are available through order v 7 (relative order v 4 ) [3,4], 1 the NRQCD LDMEs are generally not known very well beyond the one at leading order in v. Second, the perturbative corrections to the SDCs, which are currently known to next-to-next-to-leading order (NNLO) in the strong coupling constant α s , are uncomfortably large, hinting at a possible failure of the convergence of the perturbation series. Especially, the nonconvergence of the perturbation series may correspond to renormalon ambiguities that arise when 1 Here, we adopt the power counting rules in Ref. [2]. computing diagrams in dimensional regularization and resumming the perturbation series by using the Borel transform.
These difficulties can be partially overcome by considering R, the ratio of the inclusive decay rate to the electromagnetic decay rate to two photons. A considerable simplification occurs in R: not only it is independent of the leading-order LDME, but also the correction of relative order α 0 s v 2 cancels in R. The NRQCD factorization scale dependence, which arises in the corrections of relative order α 2 s and α s v 2 , also cancel in R. Finally, the renormalon ambiguities associated with the loop corrections to the initial heavy quark-antiquark states also cancel in the ratio up to relative order v 2 . Nevertheless, the renormalon ambiguities associated with the final-state gluons in the inclusive decay rate survive in R, so the perturbation series for R still suffers from nonconvergence.
In this paper, we consider the resummation of perturbative corrections to the ratio R that are associated with the chain of vacuum-polarization bubbles in the final-state gluons in the inclusive decay rate of η Q , where Q = c or b. In Ref. [5], the resummation has been performed by imposing an infrared cutoff in the calculation of perturbative QCD corrections, which allows one to avoid renormalon ambiguities that appear in the resummation if dimensional regularization is used instead. In this work, we employ dimensional regularization to compute corrections in perturbative QCD in the limit of large number of active quark flavors n f and show explicitly in this limit the appearance of renormalon ambiguities in the perturbative QCD amplitude. We also show, by computing the perturbative corrections in NRQCD, that perturbative NRQCD reproduces exactly the large n f leading renormalon ambiguity in perturbative QCD, and therefore, the NRQCD factorization formula is free of this kind of renormalon ambiguities. We argue that, due to the limited knowledge on NRQCD LDMEs of higher orders in v, using a hard cutoff instead of dimensional regularization to regularize the ultraviolet divergences in NRQCD leads to an expression for the decay rate of η Q that is more useful for phenomenological applications. We also combine our resummed calculation with the perturbative calculation of R, which is currently known to NNLO in α s and provide updated numerical results for both η c and η b .
The paper is organized as follows. In Sec. II we consider the resummation of vacuumpolarization bubble chains in the η Q decay rate and obtain a resummed expression for R.
We also combine the resummed result with the perturbative calculation of R. In Sec. III we present our numerical results for the ratio R and compare them to experimental data [6].
We conclude in Sec. IV.

II. RESUMMATION OF VACUUM-POLARIZATION BUBBLE CHAINS IN R
In this section, we present the SDCs that contribute to the inclusive decay rate of η Q .
To this end, we first shortly discuss the NRQCD factorization formula for the decay rate of η Q . Then we use the factorization formula to compute the decay rate of a perturbative QQ state in perturbative QCD and perturbative NRQCD. Finally, the SDCs are obtained by comparing the expressions for the decay rate computed in QCD and NRQCD.
A. NRQCD factorization for the decay rate of η Q The NRQCD factorization formula for the decay rate of η Q , valid through relative order v 3 , reads [3] Γ η Q = 2 Im where m is the pole mass of the heavy quark Q. The four-quark operators O 1 ( 1 S 0 ), P 1 ( 1 S 0 ) and O 8 ( 3 S 1 ) are given by Here, ψ and χ are the Pauli spinor field operators that annihilates a heavy quark and creates a heavy antiquark, respectively. The operator ← → D is the difference between the covariant derivative acting on the spinor to the right and on the spinor to the left, so that The LDMEs η Q |O 1 ( 1 S 0 )|η Q , η Q |P 1 ( 1 S 0 )|η Q , and η Q |O 8 ( 3 S 1 )|η Q are nonperturbative quantities that correspond to the probabilities to find QQ pairs in specific color and angular-momentum states in the η Q state. According to the power counting of Ref. [2], η Q |O 1 ( 1 S 0 )|η Q , which scales like v 3 , is the LDME at leading order in v; η Q |P 1 ( 1 S 0 )|η Q is suppressed by v 2 compared to the leading-order LDME [the suppression comes from the two powers of derivatives in the operator P 1 ( 1 S 0 )] and scales like v 5 ; η Q |O 8 ( 3 S 1 )|η Q , which scales like v 6 , is suppressed by v 3 compared to η Q |O 1 ( 1 S 0 )|η Q . The suppression of the LDME η Q |O 8 ( 3 S 1 )|η Q occurs because the operator O 8 ( 3 S 1 ) annihilates and creates QQ in a color-octet state through a spin-flip process [2].
Power counting rules that are more conservative than those of Ref. [2] have been suggested in Refs. [7,8]. In that power counting, which assumes Λ QCD to be larger than Moreover, there are two additional color-octet LDMEs that scale like v 3 Λ 2 QCD /m 2 and Λ 2 QCD v 3 respectively, which are given by Therefore, if we adopt the power counting rules in Refs. [7,8], Eq. (1) should include the above matrix elements to be valid through relative order Λ 2 QCD /m 2 . Considering, however, that in our numerical results, we will ignore the contribution from the color-octet LDME η Q |O 8 ( 3 S 1 )|η Q and account for its effect only in the uncertainties, and that the additional color-octet LDMEs ignored in Eq. (1) do not affect the calculation of the renormalon ambiguities that we consider in this paper, we conclude that we may consistently neglect also the additional color-octet LDMEs η Q |O 8 ( 1 S 0 )|η Q and η Q |O 8 ( 1 P 1 )|η Q , whose effect is included in the uncertainties. This is equivalent to assuming mv ≫ Λ QCD and restricting our calculation to a precision of relative order v 2 .
The imaginary parts of the SDCs 2 Im[f 1 ( 1 S 0 )/m 2 ], 2 Im[g 1 ( 1 S 0 )/m 4 ] and 2 Im[f 8 ( 3 S 1 )/m 2 ] can be computed in perturbation theory. A general method to compute the SDCs is to consider Eq. (1) with the nonperturbative meson state replaced by the perturbative QQ state with definite color and angular momentum, Here, n denotes the color and angular momentum state of the QQ. We can compute the Γ QQ on the left-hand side in perturbative QCD and compute the LDMEs on the right-hand side in perturbative NRQCD for the QQ states in various color, spin and orbital angular momentum states. Then, the SDCs can be determined by comparing the expressions on the left-and right-hand sides. In fixed-order perturbation theory, all SDCs in Eq. (1) appear from order at leading order (LO) and next-to-leading order (NLO) in α s has been computed in Refs. [9,10], and the corrections at next-to-next-to-leading order (NNLO) in α s have been calculated recently in Ref. [11]. The SDC 2 Im[g 1 ( 1 S 0 )/m 4 ] at LO in α s has been computed in Ref. [12], and the corrections at NLO in α s have been calculated in Ref. [13]. The SDC 2 Im[f 8 ( 3 S 1 )/m 2 ] has been computed up to NLO accuracy in α s in Ref. [14].
The analogous NRQCD factorization formula for the decay of η Q into two photons, valid through relative order v 3 , reads 2 where the electromagnetic operators O EM ( 1 S 0 ) and P EM ( 1 S 0 ) are given by Here, |0 is the QCD vacuum. The SDC 2 Im[f EM ( 1 S 0 )/m 2 ] have been computed up to NNLO in α s in fixed-order perturbation theory [9,15,16], and 2 Im[g EM ( 1 S 0 )/m 4 ] is available up to NLO in α s [12,13,17]. The electromagnetic LDMEs η Q |O EM ( 1 S 0 )|η Q and η Q |P EM ( 1 S 0 )|η Q can be related to the color-singlet LDMEs η Q |O 1 ( 1 S 0 )|η Q and η Q |P 1 ( 1 S 0 )|η Q by using the vacuum-saturation approximation, which holds up to corrections of relative order v 4 [2], Putting Eqs. (1) and (4) together, the NRQCD expression for the ratio R, valid up to relative where the second term in the square brackets corresponds to the correction at relative order v 2 , and the last term on the right-hand side gives the order-v 3 contribution. The order-v 2 correction to R vanishes at LO in α s ; this is because the tree-level Feynman diagrams for QQ → gg and QQ → γγ are same. The order-α s v 2 correction to R can be obtained from the order-α s v 2 corrections to Γ η Q and Γ η Q →γγ . The correction at order α s v 2 is numerically small for both η c and η b , and it is comparable to the nominal size of the order-v 3 correction, which is often neglected (and included in the uncertainties) because the color-octet matrix element η Q |O 8 ( 1 S 0 )|η Q is not known. We will follow this approach also here when providing numerical results (see Sec. III), but we will keep the color-octet matrix element tions that are associated with the running of α s , where a factor of α s is accompanied by a factor of the QCD beta function. One way to resum (partially) such corrections is to consider chains of vacuum-polarization bubbles, which reproduce fixed-order perturbation theory in the limit where the number of active quark flavors n f is large [18][19][20]. In the , the perturbative corrections at large n f that arise from initial-state virtual gluons cancel [5]. Therefore, in R, it suffices to consider only the perturbative corrections to Im[f 1 ( 1 S 0 )/m 2 ] at large n f that arise from the final-state gluons.
In the next part, we resum the QCD corrections to 2 Im[f 1 ( 1 S 0 )/m 2 ] that are associated with the final-state gluons in the large n f limit.
The series for QCD corrections corresponding to bubble-chain diagrams in general do not converge. If one attempts to make use of the Borel transform to carry out the resummation of such series, the nonconvergence manifests itself through singularities in the Borel plane. The Following Ref. [5], we employ two methods to carry out the bubble-chain resummation.
One method is naïve non-Abelianization (NNA), where we consider corrections to the gluon propagator from n f light quark loops, and we promote the light-quark part of the one-loop QCD beta function to the full one-loop QCD beta function β 0 = (33/2−n f )/(6π) [22]. That is, we make the following replacement in the gluon propagator: and Here, d is given by where, in the MS renormalization scheme, C = −5/3 and µ is the renormalization scale.
The strong coupling constant α s is also computed in the MS scheme at the scale µ. Another method is the background-field gauge (BFG) method, where the corrections to the gluon propagator from the gluon and the ghost loops that are gauge dependent are also taken into account [23]. In the BFG method in the R ξ gauge, d is given by where ξ is the gauge-fixing parameter. The choice ξ = 1 corresponds to the Feynman gauge.
If we set ξ 2 = 7/3, we reproduce the NNA method. Hence, the NNA expression for the gluon propagator may be interpreted as a special case of the BFG expression for ξ 2 = 7/3.
In order to examine the dependence on the gauge-fixing parameter ξ, we employ both the NNA method, which is equivalent to the BFG method for ξ 2 = 7/3, and the BFG method in the Feynman gauge (ξ = 1).
In the bubble-chain resummation, the left-hand side of Eq. (3) occurs from order α s through the decay into a single bubble-chain gluon. In order to decay into a virtual gluon, the QQ pair must be in a color-octet state. If we take the QQ pair to be in the color-octet spin-triplet state and take the relative momentum between the Q and theQ to vanish, the occurs from order α 0 s , and the matrix elements In order to compute the SDC 2 Im[f 1 ( 1 S 0 )/m 2 ], we consider the left-hand side of Eq. (3) at order α 2 s . At this order, Γ QQ(n) occurs through the decay into two bubble-chain gluons. If we take the QQ pair to be in the color-singlet spin-singlet S-wave state, the ma- occurs from order α s . If we take the relative momentum q between the Q and theQ to be zero, the matrix element  The appearance of the renormalon ambiguity in Γ QQ(n) and the cancellation of the ambiguity in the SDC 2 Im[f 1 ( 1 S 0 )/m 2 ] can be seen explicitly by computing the SDC in dimensional regularization. In this section, we compute the SDC 2 Im[f 1 ( 1 S 0 )/m 2 ] using bubble-chain resummation, by considering Eq. (3) where the the QQ is in the color-singlet spin-singlet S-wave state with vanishing relative momentum between the Q and theQ. We regulate the infrared divergence using dimensional regularization.

B. Computation in perturbative QCD
We compute the decay rate of a QQ pair in the color-singlet 1 S 0 state into two bubblechain gluons. We use nonrelativistic normalization for the QQ states. We set the momentum of the Q and theQ to be p. To project the QQ pair onto the color-singlet spin-singlet state, we replace the spinors by where Π 1 (p, p) and Λ 1 are the spin-singlet and color-singlet projectors, respectively [24,25], Here, 1 is the SU(N c ) unit matrix. A straightforward calculation of the diagrams in Fig. 1 gives where k and ℓ are the momenta of the final-state gluons, g = , and the trace is over the color and gamma matrices. Even though we employ dimensional regularization, it suffices to work in four dimensions because in the current calculation, we encounter no divergences that require regularization. Then, We change the integration variables k 0 , ℓ 0 to x and y, so that from x = (k 2 0 − k 2 )/(4m 2 ) and y = (ℓ 2 0 − ℓ 2 )/(4m 2 ), we obtain dk 0 = 2m 2 dx/k 0 and dℓ 0 = 2m 2 dy/ℓ 0 . Using the threemomentum components of the delta function to eliminate the integral over ℓ, we replace ℓ with k, and we obtain Here, k 0 = √ 4m 2 x + k 2 and ℓ 0 = 4m 2 y + k 2 . The lower limits of the integrals over x and y are set by the fact that the imaginary parts of K(x) and K(y) vanish for negative values of x and y, respectively. The upper limits of the integrals over x and y are set by the fact that the maximum invariant mass of a final-state particle is equal to the invariant mass of the QQ in the initial state. Since we have chosen a root of the square root function such that k 0 > 0 and ℓ 0 > 0, we can drop θ(k 0 ) and θ(ℓ 0 ). The remaining delta function in Eq. (17) constrains where The sum over n 1 and n 2 corresponds to insertions of n 1 and n 2 vacuum polarization bubbles to the two final-state gluon lines. If we use the relation which is valid for a generic function F (x, y), we obtain where We compute T (t, τ ) in the Appendix. The summation in Eq. (21) can be rewritten in integral form by using the Borel summation formula: the Borel sum of ∞ n=0 a n x n is given by the integral where φ(x) = ∞ n=0 a n x n /n!. Using this formula, we rewrite Eq. (21) as where Note that w depends on the scale µ through dependence on α s (µ) and d. This dependence, however, cancels at the level of one-loop running of α s (µ). The prefactor 2πC F α 2 s /m 2 corresponds to Γ QQ 1 ( 1 S 0 ) at leading order in α s in fixed-order perturbation theory.
The function T (t, τ ) is regular for 0 ≤ t < 1 and 0 ≤ τ < 1. For t ≥ 1 and τ ≥ 1, there are singularities in T (t, τ ) that make the value of the integral in Eq. (24) ambiguous. The singularities in T (t, τ ) that occur for the smallest values of t or τ are at t = 1 or τ = 1, see Eqs. (A20) and (A21). These singularities give the leading renormalon ambiguities in One way to estimate the size of the leading renormalon ambiguity is to inspect the difference between the results for the integral over t and τ when the integration contour is above the renormalon singularity and below the renormalon singularity [21]. From the residue theorem, the estimated ambiguity in Γ QQ 1 ( 1 S 0 ) that arises from the leading renormalon sin- For the case of η c , the numerically estimated size of the leading renormalon ambiguity is of relative order one compared to Γ QQ 1 ( 1 S 0 ) at leading order in α s in fixed-order perturbation theory. This implies that for η c , the value of the perturbation series Γ QQ 1 ( 1 S 0 ) has an ambiguity of order one. Even for the case of η b , the estimated ambiguity can be of relative order 10 −1 , which is comparable to the nominal size of the order-v 2 corrections to the decay rate. Therefore, in order to make an accurate theoretical prediction of the η Q decay rate, it is crucial to have a factorization formula where such ambiguities are absent.
The renormalon ambiguities that arise from the singularities in T (t, τ ) are located at t = t 0 or τ = t 0 with t 0 = 1 being the smallest and involves a factor e −t 0 w . If we consider only the one-loop running of α s , this factor can be written as Therefore, renormalon ambiguities that arise from the singularities in T (t, τ ) that are located at larger values of t or τ are suppressed by powers of Λ QCD /µ. We can estimate the renormalon ambiguity from the first subleading singularities in T (t, τ ) which are located at The estimated renormalon ambiguity in Γ QQ 1 ( 1 S 0 ) from the first subleading singularities is This ambiguity is of relative order 10 −1 for η c and is of relative order 10 −3 for η b . For η c , the ambiguity is comparable to the nominal size of the order-v 4 correction to the decay rate, Because we set the relative momentum between the Q and theQ to be zero, on the righthand side of Eq. (3), only the matrix element Since we only consider the NRQCD operators of the lowest mass dimensions, the contributions from the right-hand side of Eq. (3) will only reproduce the leading renormalon ambiguity in Eq. (24).
At leading order in α s , the color-singlet LDME is given by The color-octet matrix element vanishes at order α 0 s , but it receives contributions at order α s from the insertion of the σ · B vertices to the quark and antiquark lines. The corresponding Feynman diagrams are shown in Fig. 2. The sum of the diagrams gives where Here, we use the following shorthand notation we can deform the the contour for the integration over k 0 so that We first integrate over k 0 . The result is where Here, F (a, b; c; z) is the hypergeometric function. Because we are matching QCD with NRQCD, we expand in 1/m and keep only the contribution at leading power in 1/m [2], In dimensional regularization, the integral over k is scaleless, and hence vanishes, where in the second equality we split the integral over |k| so that the first (second) integral corresponds to the region where |k| is small (large). The first (second) integral is finite only when t < 1 (t > 1). After integrating over |k|, we use analytical continuation to extend the region of t to the whole complex plane. Since the first term in the parenthesis comes from the region where |k| is small, and the second term originates from the region where |k| is large, the first and second terms in the parenthesis correspond to the IR and UV renormalon singularities of the LDME O 8 ( 3 S 1 ) QQ 1 ( 1 S 0 ) , respectively. Then, we can write the right-hand side of Eq. (3) as The subscripts IR and UV denote the origins of the IR and UV renormalon singularities, respectively. Note that to derive Eq. (40) we rewrote Eq. (12) as 2 Im and symmetrized in t and τ . By comparing Eq. (40) with Eq. (24), we can see that the infrared renormalon singularities in T DR 8 (t, τ ) [terms proportional to 1/(1 − t) IR and 1/(1 − τ ) IR ] reproduce the leading renormalon singularities in T (t, τ ) at t = 1 or τ = 1, and, therefore, Eq. (40) reproduces the leading renormalon ambiguity in Eq. (24). Then, the is free of the leading infrared renormalon ambiguity. On the other hand, the UV renormalon singularities in T DR 8 (t, τ ) [terms proportional to 1/(1 − t) UV and 1/(1 − τ ) UV ] has no counterpart in perturbative QCD [Eq. (24)]; therefore, the SDC 2 Im[f 1 ( 1 S 0 )/m 2 ] has a UV renormalon ambiguity. Since the UV renormalon ambiguity is absent in Eq. (24), the UV renormalon ambiguities in the SDC 2 Im[f 1 ( 1 S 0 )/m 2 ] and the LDME cancel in the factorization formula Eq. (3). Since the UV renormalon ambiguities are of ultraviolet origin, the nonperturbative LDME η Q |O 8 ( 3 S 1 )|η Q has the same UV renormalon ambiguities as the perturbative LDME O 8 ( 3 S 1 ) QQ 1 ( 1 S 0 ) , with the perturbative QQ states replaced by the nonperturbative meson state. Therefore, the ambiguity is absent in the factorization formula for the inclusive decay rate of η Q [Eq. (1)].
Even though the UV renormalon ambiguities cancel in the factorization formula, it is still necessary to define the SDC 2 Im[f 1 ( 1 S 0 )/m 2 ] and the LDME η Q |O 8 ( 3 S 1 )|η Q unambiguously in order to compute the inclusive decay rate. An unambiguously defined LDME will lead to an expression for 2 Im[f 1 ( 1 S 0 )/m 2 ] that is free of UV renormalon ambiguities; however, different definitions will lead to different expressions for the SDC 2 Im and the LDME η Q |O 8 ( 3 S 1 )|η Q . Especially, the differences between different definitions of the LDME can be of the size of the UV renormalon ambiguity, which can be estimated from We note that the renormalon ambiguity scales like (Λ QCD /µ) 2 , which is different from velocity-scaling rules of the LDMEs [21]. Hence, there is a possibility that the renormalon ambiguity in the LDMEs can spoil the expansion in powers of v in case the renormalon ambiguity of an LDME exceeds its nominal size.
We can define NRQCD LDMEs so that the LDMEs are free of UV renormalon ambiguities and also respect the velocity-scaling rules by regulating the UV divergences using a cutoff regulator. In perturbative calculations, it is most convenient to apply a hard cutoff Λ on the size of the spatial momentum of the gluon. The cutoff Λ should be large enough so that it encompasses the relevant momentum regions in NRQCD, while Λ < m so that the expansion in powers of 1/m is valid. Hence, it is customary to choose Λ ∼ mv. While the NRQCD LDMEs in hard-cutoff regularization are free of renormalon ambiguities, they depend on the cutoff Λ.
If we regularize the UV divergences in NRQCD with a hard cutoff Λ for the spatial momentum of the gluon, the integral over k in Eq. (39) now becomes which yields Here, the superscript (Λ) denotes that a UV cutoff was used. Then, the right-hand side of As we have discussed in Sec. II A, the SDC 2 Im[g 1 ( 1 S 0 )/m 4 ] does not appear in Eq. (47) because the contribution from the LDME QQ 1 ( 1 S 0 )|P 1 ( 1 S 0 )|QQ 1 ( 1 S 0 ) to the right-hand side of Eq. (3) vanishes at order α 2 s if the QQ is in the color-singlet 1 S 0 state and the relative momentum between the Q and theQ is zero. The singularities of T which reproduce the leading renormalon singularities in T (t, τ ) at t = 1 or τ = 1. It is clear that, from the expression for T Since the functions T (t, τ ) and T (Λ) 8 (t, τ ) have same singularities at t = 1 or τ = 1 [Eqs. (26,49)], those singularities cancel in T (t, τ ) − T (Λ) 8 (t, τ ). Therefore, the leading renormalon ambiguities in Γ QQ 1 ( 1 S 0 ) that originate from the singularities at t = 1 or τ = 1 are absent in Together with our results for 2 Im[f 1 ( 1 S 0 )/m 2 ] and the perturbative expression for 2 Im[f EM ( 1 S 0 )/m 2 ] at LO in α s [2], we obtain the resummed expression for R at leading order in v including resummed QCD corrections in the large n f limit, where Here e Q is the fractional electric charge of the heavy quark Q. As previously discussed, the order-v 2 correction to R vanishes at LO in α s . We neglect the correction at order α s v 2 that was computed in fixed-order perturbation theory, because it was found to be small numerically [13,17] and is of comparable size to the contribution of order v 3 . The order-v 3 contribution to R can be written as Since it is not known how to compute the color-octet LDME η Q |O 8 ( 3 S 1 )|η Q reliably, we ignore R 8 and instead consider its effects in the uncertainties.
We can combine our results for R Resum with fixed-order calculations of R, so that the corrections at NLO and NNLO in α s are valid beyond the large n f limit. By using the expressions for Γ η Q and Γ η Q →γγ valid to NNLO in α s , we obtain [11,15,16] T R n f , and n H is the number of heavy quark flavors.R (2) is known as a function of n f for the case n H = 1 only.R (2) = 117.144 for n f = 3 andR (2) = 86.421 for n f = 4. The large n f limit of R (2) is given in Ref. [11] as lim n f →∞R (2) /n 2 f = 0.37581 (3). The full n f dependence ofR (2) can be obtained from Refs. [11,16] aŝ whereR (2) lbl = 0.7313 × C F × f (e f /e Q ) 2 + 0.6470 × C F n H is the "light by light" contribution to the two-photon decay rate that occurs through QQ → gg → γγ via a light quark loop.
Here, the sum is over n f light quark flavors, and e f is the fractional charge of a light quark of flavor f .
When n H = 1, the heavy quark Q contributes to the renormalization scale dependence of α s , which cancels the explicit renormalization scale dependence of R Pert from the logarithms of µ/m. It is possible to decouple the heavy quark Q from the running of α s by using the decoupling relations between α s for n f and n f + 1 active quark flavors [26], so that the heavy quark Q does not affect the renormalization scale dependence of R Pert for µ < m. By using Eq. (25) of Ref. [26], we decouple the heavy quark, and then we obtain whereR ′(1) = 199 18 − 13π 2
The expression for R Pert in Eq. (54) is obtained from the electromagnetic and inclusive decay rates of η Q that were calculated in the MS renormalization scheme. In order to make Because Eq. (54) is computed by using dimensional regularization, all power divergences are absent in Eq. (54). In the fixed-order perturbation theory calculation using the cutoffregularization scheme, the color-singlet contribution to R receives power-divergent contribution from the one-loop correction to the color-octet LDME, which is of relative order α s Λ 2 /m 2 . If we set Λ ∼ mv, this contribution is of relative order α s v 2 . For α s ∼ v, this is the same size as the color-octet contribution to R, whose relative size is of order where g 2 is defined by the integral We evaluate this integral numerically to find g 2 = −3.22467022 (9). It can be seen that the bubble-chain resummation reproduces the fixed-order perturbation series in the large n f limit by comparing the coefficients of (α s n f ) n for n = 1 and 2 in Eqs.
This expression is valid for an arbitrary number of light quark flavors n f , except that we only consider three light quark flavors for the contribution to "light by light" contribution to the two-photon decay rate at NNLO in α s that occurs through QQ → gg → γγ via a light quark loop, which is proportional to the sum of squares of light quark fractional charges [16]. For the decay of η b , where we consider four light quark flavors in the light by light contribution, we obtain On the other hand, δR Resum gives, for both η c and η b , Here, we chose the gauge-fixing parameter in the BFG method to be ξ = 1, which corresponds to the Feynman gauge. As expected, δR Resum reproduces R Pert only in the large n f limit. While it is not at all surprising that δR Resum does not reproduce R Pert beyond the large n f limit, the size of the coefficients of order α 2 s n f and α 2 s n 0 f are quite large in R Pert . If the large n f limit does not provide a good approximation to the fixed-order calculation, the perturbative convergence of R Pert − δR Resum can be spoiled. To inspect the agreement between R Pert and δR Resum explicitly, we consider the perturbation series of R Pert and δR Resum with n f = 3 and 4 light quark flavors for the case of η c and η b , respectively. For the decay of η c with n f = 3, we obtain δR Resum ηc,NNA = R 0 1 + 5.76 δR Resum ηc,BFG = R 0 1 + 7.76 For the decay of η b with n f = 4, we obtain δR Resum η b ,BFG = R 0 1 + 7.33 In all cases, agreement between R Pert and δR Resum is poor at NNLO in α s , even though the difference vanishes in the large n f limit. Hence, perturbative corrections may not still be in control because of the large radiative corrections in R Resum+∆NNLO beyond the large n f limit.

III. NUMERICAL RESULTS
We now discuss our numerical results, based on our expression of R in Eq. (59). We first describe our numerical inputs. We take the heavy-quark mass m to be 1.5 GeV for the charm quark, and 4.6 GeV for the bottom quark. These values are numerically close to the one-loop pole mass. We have a freedom in choosing the values of the heavy-quark mass due to the ambiguity in the pole mass that is of the order of Λ QCD . The quantity R Resum depends on the heavy-quark mass only through log µ 2 4m 2 . Hence, the dependence on the choice of the value of the heavy-quark mass m is beyond our accuracy. We take the number of light quark flavors to be n f = 3 for η c and n f = 4 for η b . In evaluating R Resum , we consider both the NNA and the BFG method. The NRQCD cutoff Λ must be chosen between mv and m, where v 2 ∼ 0.3 for η c and v 2 ∼ 0.1 for η b . Accordingly, we set the central value for the NRQCD cutoff to be 1 GeV for η c , and 2 GeV for η b . We take the central value for µ to be the heavy-quark mass. We compute α s in the MS renormalization scheme using the Mathematica package RunDec [26]. We use α = 1/137.036. In the BFG method, we set We list the sources of uncertainties that we consider. We vary Λ by ±25% of its central value. We vary µ between 1 and 2 GeV for η c , and between 2 and 6 GeV for η b . Because our expression for R in Eq. (59) only depends on the heavy-quark mass through log µ 2 4m 2 , where µ is the renormalization scale in the MS scheme, the dependence of R on m is very mild. Also, the change in R from varying m is equivalent to the change in R from varying µ. Since we already take into account the uncertainty from the dependence on µ, we ignore the uncertainty from the dependence on m. Finally, we estimate the color-octet LDME by using the perturbative estimate given in Ref. [5], where we choose v 2 = 0.3 for η c , and v 2 = 0.1 for η b . The uncertainty from ignoring the color-octet contribution is then estimated by ±|R 8 |, where R 8 is given by Eq. (53). Note that the color-octet matrix element η Q |O 8 ( 3 S 1 )|η Q depends on Λ; hence, there is a correlation between errors from varying Λ and ignoring the color-octet matrix element. We, however, ignore the correlation and add the uncertainties in quadrature.

A. Decay of η c
We first present our numerical results for the ratio R for η c . In Fig. 3 For the case of NNA, R Resum+∆NLO develops a stronger dependence on µ from the fixed-order corrections of relative order α s in R Pert − δR Resum . In R Resum+∆NNLO , the renormalization scale dependence is slightly worse than R Resum+∆NLO , because of the subleading logarithm of the form α 2 s log µ 2 4m 2 in R Pert − δR Resum . For µ > m, there are also uncanceled leading logarithms in R Pert − δR Resum that are proportional to n H . For the case of the BFG method, R Resum , R Resum+∆NLO and R Resum+∆NNLO all depend on µ very mildly. This is because in the BFG method in the Feynman gauge, there is almost exact cancellation in the fixed-order corrections R Pert − δR Resum , which contain most of the dependence on µ.
In Fig. 4 we show the dependence on the NRQCD cutoff Λ of R Resum , R Resum+∆NLO and R Resum+∆NNLO at µ = m. In all cases, the dependence on Λ is mild, and the numerical values of R rise slowly with increasing Λ. As we will see later, the uncertainty estimated from varying Λ is smaller than the estimated uncertainty from neglecting the color-octet contribution.
For µ = m, the fixed-order corrections in R Pert − δR Resum are positive for both the contributions of relative order α s and order α 2 s . In NNA, R Resum+∆NLO is larger than R Resum by about 16% of the central value of R Resum+∆NNLO , and R Resum+∆NNLO is larger than R Resum+∆NLO by about 20% of the central value of R Resum+∆NNLO . In the BFG method, R Resum+∆NLO is larger than R Resum by about 4% of the central value of R Resum+∆NNLO , and R Resum+∆NNLO is larger than R Resum+∆NLO by about 6% of the central value of R Resum+∆NNLO .
While the effects of the fixed-order corrections appear less dramatic than the effects of the radiative corrections in the fixed-order calculation in Ref. [11], the fact that the corrections are larger at NNLO than at NLO in α s implies that the perturbative corrections may still not be under control. As discussed in the previous section, this is related to the large perturbative corrections in the fixed-order corrections that go beyond the large n f limit; because the treatment of the renormalon ambiguities in this work is only valid in the large n f limit, we have little or no control over the convergence of the perturbation series beyond the large n f limit.
where the first uncertainty is from µ, the second from Λ, and the third uncertainty is from the neglected color-octet contribution. For BFG, we obtain where the uncertainties are as in NNA. Our numerical results for the NNA and the BFG methods are compatible within uncertainties. We note that R Resum+∆NNLO in the NNA method has a large uncertainty from its strong renormalization-scale dependence for small µ.
In estimating the uncertainties in our numerical results we have neglected the possibility that the convergence of the fixed-order corrections in R Pert − δR Resum may not be in control.
We roughly estimate the uncertainty from this nonconvergence by comparing our numerical We can compare our numerical results with measurements. PDG reports two values for the η c branching ratio to two photons [6]. The first PDG value Br(η c → γγ) = (1.59 ± 0.13) × 10 −4 is from a constrained fit of partial widths. If we take the inverse, we obtain R exp fit = (6.29 +0.56 −0.48 ) × 10 3 . The second PDG value is from an average of measurements, which gives Br(η c → γγ) = (1.9 +0.7 −0.6 )×10 −4 . Taking the inverse gives R exp average = (5.3 +2.4 −1.4 )×10 3 . The two PDG values are compatible with each other due to the large uncertainties in R exp average . The uncertainty in R exp fit is smaller than the uncertainty in R exp average or the uncertainties in our numerical results for R. Note that R exp average is compatible with our results for R in Eqs. (66) and (67). There is, however, a tension between R exp fit and our numerical results. We also note that our calculation of R also applies for the η c (2S) state as well. The PDG value for the η c (2S) branching ratio to two photons is Br(η c (2S) → γγ) = (1.9 ± 1.3) × 10 −4 , which is compatible with the PDG values for the Br(η c → γγ) and our results for R in Eqs. (66) and (67).

B. Decay of η b
We now present our results for η b . In Fig. 5, we show the dependence on the renormalization scale µ of R Resum , R Resum+∆NLO and R Resum+∆NNLO at Λ = 2 GeV. Just like for the case of η c , R Resum has some dependence on µ because the renormalization scale dependence in R Resum only cancels at the one-loop level. For the case of NNA, R Resum+∆NLO develops a stronger dependence on µ from the fixed-order corrections of relative order α s in R Pert −δR Resum . This dependence is partially canceled by the corrections of relative order α 2 s , which contain logarithms of µ, and R Resum+∆NNLO depends on µ mildly. For the case of the BFG method, R Resum , R Resum+∆NLO and R Resum+∆NNLO all depend on µ mildly. This is again because the choice ξ = 1 minimizes the size of the fixed-order corrections R Pert − δR Resum , which contain most of the dependence on µ. In Fig. 6 we show the dependence on the NRQCD cutoff Λ of R Resum , R Resum+∆NLO and R Resum+∆NNLO at µ = m, and 1.5 GeV ≤ Λ ≤ 2.5 GeV. In all cases, the dependence on Λ is mild, and the numerical values of R rise very slowly with increasing Λ. Just like for the case of η c , the uncertainty estimated from varying Λ is smaller than the estimated uncertainty from neglecting the color-octet contribution.
For NNA at µ = m and Λ = 2 GeV, the fixed-order corrections in R Pert − δR Resum are positive for both the contributions of relative order-α s and order-α 2 s . In NNA, R Resum+∆NLO is larger than R Resum by about 11% of the central value of R Resum+∆NNLO , and R Resum+∆NNLO is larger than R Resum+∆NLO by about 3% of the central value of R Resum+∆NNLO . In the BFG method, at µ = m and Λ = 2 GeV, R Resum+∆NLO is larger than R Resum by about 2% of the central value of R Resum+∆NNLO , and R Resum+∆NNLO is smaller than R Resum+∆NLO by about 3% of the central value of R Resum+∆NNLO . The effects of the fixed-order corrections are much less dramatic than the corrections to the η c decay rate, thanks to the smaller size of α s and larger n f . We estimate the uncertainties by varying µ between 3 and 6 GeV, and by varying Λ between 1.5 and 2.5 GeV. We also include the uncertainty for ignoring the color-octet contribution. For NNA, we obtain where the first uncertainty is from µ, the second from Λ, and the third uncertainty is from the neglected color-octet contribution. For BFG, we obtain It is not yet possible to compare our results in Eqs. (68) and (69) with measurements because the partial width Γ η b →γγ has not been observed yet. In Ref. [27], the authors made use of the heavy-quark spin symmetry to extract η b LDMEs from the Υ LDMEs and made the prediction Γ η b →γγ = 0.512 +0.096 −0.094 keV. If we multiply this result to our results for R, we obtain Γ η b = 11.9 +2.3 −2.2 MeV for NNA, and Γ η b = 12.4 +2.4 −2.3 MeV for the BFG method. Reference [28] makes use of the ratio of the leptonic decay rate of the Υ to the decay rate of η b into two photons in the potential NRQCD effective field theory to predict Γ η b →γγ from the measured value for Γ Υ→e + e − . The prediction in Ref. [28] is given by Γ η b →γγ = 0.54±0.15 keV, which is compatible with the prediction in Ref. [27]. If we use this prediction, we obtain Γ η b = 12.5 ± 3.5 MeV for NNA and Γ η b = 13.0 ± 3.7 MeV for the BFG method. These predictions for the η b decay rate are compatible with the PDG value for the η b decay width, which is given by Γ η b = 10 +5 −4 MeV.

C. Comparison with previous results
We now compare our numerical results with previous results for R. In Ref. [5], the authors also considered resummation of bubble-chain contributions to R for the decay of η c . The results of Ref. [5] are equivalent to R Resum+∆NLO , except that in Ref. [5], the SDC was computed by imposing a hard IR cutoff which affects both the virtuality and the spacial momentum of the gluon in the perturbative QCD calculation [Eq. (18)]. The authors of Ref. [5] identified the contribution from the momentum region that was neglected in the perturbative QCD calculation as the contribution from perturbative NRQCD which is regulated by a hard UV cutoff on the gluon momentum. The hard cutoff that was used in Ref. [5] is given by k 2 ≤ 4m 2 δ and k 2 ≤ m 2 (2 √ δ − δ) 2 , where k is a gluon momentum, and δ = 0.1. If we take m = 1.5 GeV, we obtain k 2 ≤ 0.9 GeV 2 and |k| ≤ 0.8 GeV. Numerically, the hard cutoff imposed on |k| is similar to the hard UV cutoff Λ that we have employed in this paper, although in this work, there is no cutoff on the virtuality of the gluon. The main advantage of this work compared to Ref. [5] is that in this work, the appearance and The authors of Ref. [11] presented their numerical results for Br(η c → γγ), which is equal to R −1 , based on their fixed-order calculation of the inclusive decay rate of η c and the decay rate of η c into two photons in Ref. [16] to NNLO accuracy in α s . The result in Ref. [11] is based on the perturbation expansion of Br(η c → γγ) = R −1 to NNLO in α 2 s . By varying the renormalization scale µ from 1 GeV to 3 times the charm quark mass, the authors of Ref. [11] obtained Br(η c → γγ) = (3.1-3.3)×10 −4 , which gives R ηc = (3.0-3.2)×10 3 . This result is compatible with our result in Eq. (67) in the BFG method, but is smaller than our result in Eq. (66) in NNA. Also, the uncertainty in the result in Ref. [11] is smaller than the uncertainties in our results, due to the cancellation of the renormalization-scale dependence at the two-loop level in the fixed-order calculation. Moreover, the uncertainty from the coloroctet contribution at relative order v 3 has been neglected in Ref. [11]. One can obtain a different numerical result if one considers the perturbation series of R ηc = [Br(η c → γγ)] −1 , which is given by Eq. (54). If we use Eq. (54) we obtain R ηc = 4.9 × 10 3 at µ = m.
This disagrees with the numerical results in Ref. [11], and the discrepancy is much larger than the uncertainties estimated in Ref. [11] by varying the renormalization scale µ. The difference between the numerical results based on the perturbation series of R ηc and the one based on the perturbation series of Br(η c → γγ) shows that the nonconvergence of the perturbation series generates a sizable ambiguity. This is consistent with our estimate of the leading renormalon uncertainty in Eq. (27). Our results in Eqs. (66) and (67) also suffer from nonconvergence of the fixed order corrections in R Pert − δR Resum , because we have no control over the convergence of the perturbation series beyond the large n f limit.
We have roughly estimated the uncertainty from this nonconvergence by comparing our numerical results with the series expansion of Br(η c → γγ) = 1/R Resum+∆NNLO ηc in powers of α s through NNLO accuracy. We have found that our rough estimate of the uncertainty from the nonconvergence does not exceed the uncertainties in our numerical results.
In Ref. [11], the authors also presented their numerical results based on the perturbative expression of Br(η b → γγ) to NNLO accuracy. They obtained Br(η b → γγ) = (4.8 ± 0.7) × 10 −5 . If we take the inverse we obtain R η b = (2.1 +0.4 −0.3 ) × 10 4 , which is compatible with our numerical results in Eqs. (68) and (69) within uncertainties. The uncertainties in Eqs. (68) and (69) are much smaller than the result in Ref. [11] because the bubblechain resummation reduces the dependence on the renormalization scale µ compared to the fixed-order calculation. If we use the perturbative expression of R in Eq. (54), which is valid to NNLO accuracy in α s , we obtain R η b = 2.39 × 10 4 at µ = m, which agrees with the numerical result in Ref. [11] within uncertainties. The relative discrepancy between the numerical result from the perturbative expression of the branching ratio into two photons and the numerical result from the perturbative expression of R is smaller in the case of η b compared to the case of η c . This can be understood from our estimate of the leading renormalon ambiguity [Eq. (27)] : since the decay of η b occurs at a higher energy scale than the decay of η c , the renormalon ambiguity is suppressed compared to the case of η c .
Nevertheless, the ambiguity is still sizable compared to the estimated uncertainties in our numerical results in Eqs. (68) and (69). Therefore, even for the case of η b , resumming large perturbative corrections is crucial in obtaining a reliable theoretical prediction.
In Ref. [29], the authors applied the principle of maximal conformality (PMC), which is a method for choosing the renormalization scale µ for a given perturbation series, to the perturbative expression for R to NLO accuracy. The authors of Ref. [29] claim that, by applying the PMC, the β function appearing in the perturbation series, which are associated with the running of α s , is absorbed into the running coupling, and the convergence of the perturbation series is improved. When they include the relative order-α s and order-α s v 2 corrections to R, they obtain, after applying the PMC, R = (6.09 +0.21 4m 2 ) n to all orders in α s . It is worth noting that, unlike the expressions for R at LO and NLO accuracies, the perturbative expression for R at NNLO accuracy no longer suffers from the severe dependence on the renormalization scale [11]. Even if we consider a wide range of the renormalization scale, as the authors of Ref. [11] have done, it is not possible to obtain a value of R that is close to the result of Ref. [29] if one uses the expression for R at NNLO accuracy. Shortcomings of the PMC approach have been discussed in Ref. [30].

IV. SUMMARY AND DISCUSSION
In this paper we have presented an analysis of the ratio R of the inclusive decay rate of the η Q meson to the partial decay rate into two photons, where Q = c or b. In the calculation of the short-distance coefficients, we resum large perturbative corrections in the form (α s β 0 ) n to all orders in α s by including contributions from bubble chain insertions in the gluon propagator. This bubble-chain resummation reproduces fixed-order perturbative calculations in the large n f limit. This resummation has been done in Ref. [5] by imposing an infrared cutoff in the perturbative calculations. In this work, we regulate the infrared divergences using dimensional regularization, so that the appearance of the renormalon ambiguity in the perturbative QCD calculation and the cancellation of the ambiguity in the factorization formula can be seen explicitly. We use naïve non-Abelianization and the background-field gauge method to carry out the resummation, which are unambiguous procedures for resumming bubble chains.
We confirmed that, by using the factorization formula valid to relative order v 3 , the leading renormalon ambiguity of infrared origin that arises from the perturbative QCD calculation is reproduced in the perturbative NRQCD calculation, and therefore, the shortdistance coefficients are free of infrared renormalon ambiguities. We also showed that, if we use dimensional regularization to regulate the ultraviolet divergences in NRQCD, the coloroctet LDME suffers from renormalon ambiguities of ultraviolet origin, but the ambiguity cancels in the factorization formula. Since it is not known how to compute the color-octet LDME reliably, and it is known that the color-octet LDME is suppressed by v 3 compared to the leading-order color-singlet LDME, the color-octet contribution is often neglected in the factorization formula. However, in a resummed calculation, the neglect of the coloroctet contribution results in a sizable ambiguity in the factorization formula. We argued that, for phenomenological applications, we obtain a more useful factorization formula if we use hard-cutoff regularization to regulate the ultraviolet divergences in NRQCD where such ambiguity no longer appears.
Our result for the resummed calculation of R is shown in Eq. (51). We combine our result with the calculation in fixed-order perturbation theory to next-to-next-to-leading order accuracy in α s [Eq. (54)] [9-11, 15, 16]. The expression for the combined result is shown in Eq. (59). We use the expression in Eq. (59) in our numerical analysis.
In our numerical analysis, we estimated uncertainties by varying the renormalization scale and the NRQCD ultraviolet cutoff. We also included the effect of the color-octet contribution by estimating the size of the uncalculated color-octet LDME. Our numerical results for the ratio R for the decay of η c are given in Eqs. Again, the numerical results in Eqs. (68) and (69) agree within uncertainties. Since the decay of η b into two photons is yet to be measured, we cannot compare our results for R directly with measurements for the case of η b . By using predictions of Γ η b →γγ in Refs. [27,28], we have obtained predictions of Γ η b that is compatible with the current measurement of the η b decay rate.
We have compared our numerical results with previous calculations of R in Ref. [5], where the authors also considered bubble-chain resummation, and the results based on fixed-order perturbation theory in Ref. [11]. In Ref. [5], the authors made predictions of R for the decay of η c by combining the resummed result, which was computed by using a fixed infrared cutoff, with the fixed-order calculation valid to next-to-leading order in α s . Our numerical results agree with the results in Ref. [5] for the background-field gauge method [Eq. (67)], but there is tension in the result in naïve non-Abelianization [Eq. (66)]. This discrepancy is mostly from the inclusion of the fixed-order corrections R Pert − δR Resum at next-to-next-toleading order in α s [Eq. (59)]. While the authors of Ref. [5] included the effect of the coloroctet contribution in the uncertainties, the uncertainty from the dependence on the infrared cutoff was neglected. The numerical results for R for η c in Ref. [11] is also compatible with our results in the background-field gauge method [Eq. (67)], but disagrees with our results in naïve non-Abelianization [Eq. (66)]. The uncertainties in the result of Ref. [11] is smaller than the uncertainties in our results because in the fixed-order calculation, the dependence on the renormalization scale cancels at two-loop accuracy, and the uncertainty from the uncalculated color-octet contribution is neglected. Also, the fixed-order calculation in Ref. [11] suffers from a sizable uncertainty from the nonconverging perturbation series, which, for the case of η c , can be of relative order one. The authors of Ref. [11] also made a prediction of R for the decay of η b , which agrees with our results in Eqs. (68) and (69) within uncertainties. The uncertainties in our results are smaller than the uncertainty in the prediction from fixed-order perturbation theory in Ref. [11]. Although in the case of η b , the estimated renormalon ambiguity in the perturbation series for R is smaller than the case of η c , our estimate of the ambiguity is larger than the uncertainties in our numerical results in Eqs. (68) and (69). Therefore, we conclude that resummation is necessary in order to obtain an accurate prediction of R for the decay of η b .
In Ref. [29], the authors applied the principle of maximal conformality to the perturbative expression of R for the decay of η c valid to next-to-leading order in α s . The authors of Ref. [29] claimed that a resummed perturbative expression can be obtained, where the β function that is associated with the running of α s are absorbed into the coupling, by using the principle of maximal conformality. However, we find that our resummed result disagrees with the result in Ref. [29]. The result in Ref. [29] also disagrees with the result from fixed-order perturbation theory valid to next-to-next-to-leading order in α s in Ref. [11].
It is noticeable that the uncertainty estimated from the neglected color-octet contribution is quite significant for both η c and η b . This suggests that in order to have a more precise prediction of R, it is necessary to include color-octet contributions in R. Including the coloroctet contribution may also reduce the uncertainty from the NRQCD cutoff dependence, because the dependence on the NRQCD cutoff cancels in the factorization formula between the color-singlet and color-octet contributions. Since currently it is not known how to calculate the color-octet matrix element reliably, it would be important to develop new ideas to investigate the nature of the color-octet matrix element in NRQCD and other effective field theories such as potential NRQCD, which may help constrain the color-octet contribution.
In our numerical results we included corrections from fixed-order calculations to next-tonext-to-leading order in α s . While the bubble-chain resummation reproduces the fixed-order corrections in the large n f limit, the fixed-order corrections are still found to be significant beyond the large n f limit; even after the bubble-chain resummation, the numerical results for R for the decay of η c suggest nonconvergence of perturbative corrections to persist beyond the large n f limit. Therefore, in order to gain control over the perturbation series of R for the decay of η c , it may be necessary to consider renormalon ambiguities beyond the large n f limit. By inspecting the n f -dependence of the fixed-order corrections to the electromagnetic decay rates Γ ηc→γγ and Γ J/ψ→e + e − , which are available up to two [9,15,16] and three loops [31][32][33][34][35], respectively, we find that the fixed-order corrections are also significant beyond the large n f limit in those electromagnetic decay rates. Hence, the bubblechain resummation calculations of those decay rates in Ref. [36] seem to fail to reproduce the fixed-order calculations.
We have examined a method that is often employed for computing renormalon singu-larities in the heavy-quark pole mass described in Refs. [37,38], where the renormalon ambiguities that scale like powers of Λ QCD are subtracted from the divergent perturbation series. This method has an advantage that it does not rely on the large-n f limit. We have found that a naïve application of the method in Refs. [37,38] to the electromagnetic decay rates Γ ηc→γγ , Γ J/ψ→e + e − and the inclusive decay rate of η c lead to estimates of the perturbative series that are in poor agreement with the fixed-order corrections.
By combining the resummed result for R and the fixed-order calculations valid up to next-to-next-to-leading order in α s , we have obtained precise predictions of R for the decay of η b with uncertainties that could be as small as 5%. Therefore, the measurement of Γ η b →γγ in ongoing and future experiments is highly anticipated. We also look forward to improved experimental measurements for the decay rate of η b , as well as the total and partial decay rates of η c . with f (x, y) given in Eq. (19). We need to calculate the derivatives of T (t, τ ) at t = τ = 0.
Plugging Eq. (A3) in Eq. (A2) we obtain T (t, τ ) = 1 π 2 sin(πt) sin(πτ ) where f 0 and f 1 can be any constants or any functions of t and τ that are analytic in the vicinity of the origin of the complex-t and complex-τ planes. Then, we write Eq. (A13) as where . (A19) Setting f 0 and f 1 to unity, the above expression for T (t, τ ) is convergent for all t < 1. In particular, one can verify that T (t, τ ) does not have any singularity at t = 1 2 although it contains a factor of Γ(t − 1/2). One can also show that lim τ →1 (1 − τ )T (t, τ ) = − 3 π sin(πt) , and by symmetry argument