${\cal O}(\alpha_sv^2)$ corrections to hadronic decay of vector quarkonia

Within the nonrelativistic QCD (NRQCD) factorization framework, we compute the ${\mathcal O}(\alpha_s v^2)$ corrections to the hadronic decay rate of vector quarkonia, exemplified by $J/\psi$ and $\Upsilon$. Setting both the renormalization and NRQCD factorization scales to be $m_Q$, we obtain $\Gamma(J/\psi\to {\rm LH})= 0.0716\frac{\alpha_s^3}{m_c^2} \langle \mathcal{O}_1({}^3S_1)\rangle_{J/\psi} [1-1.19\alpha_s+(-5.32+3.03\alpha_s)\langle v^2\rangle_{J/\psi}]$ and $\Gamma(\Upsilon\to {\rm LH})= 0.0716\frac{\alpha_s^3}{m_b^2}\langle\mathcal{O}_1({}^3S_1)\rangle_{\Upsilon}[1-1.56\alpha_s+(-5.32+4.61\alpha_s)\langle v^2\rangle_{\Upsilon}]$. We confirm the previous calculation of $\mathcal{O}(\alpha_s)$ corrections on a diagram-by-diagram basis, with the accuracy significantly improved. For $J/\psi$ hadronic decay, we find that the ${\mathcal O}(\alpha_sv^2)$ corrections are moderate and positive, nevertheless unable to counterbalance the huge negative corrections. On the other hand, the effect of ${\mathcal O}(\alpha_sv^2)$ corrections for $\Upsilon(nS)$ is sensitive to the $\mathcal{O}(v^2)$ NRQCD matrix elements. With the appropriate choice of the NRQCD matrix elements, our theoretical predictions for the decay rates may be consistent with the experimental data for $\Upsilon(1S,2S)\to {\rm LH}$. As a byproduct, we also present the theoretical predictions for the branching ratio of $J/\psi(\Upsilon)\to 3\gamma$ accurate up to $\mathcal{O}(\alpha_s v^2)$.


I. INTRODUCTION
The successful predictions of charmonium annihilation decay rates are among the earliest triumph of perturbative QCD [1]. Among the quarkonia families, the spin-triplet S-wave sates with J P C = 1 −− , exemplified by J/ψ and Υ, undoubtedly occupy the central stage, whose various decay channels have been extensively studied both experimentally and theoretically. Among a variety of decay channels of the vector quarkonia, the inclusive hadronic decays are particularly interesting and important. Obviously these are not only the most dominant decay channels, but an ideal place to test our understanding about the interplay between perturbative and nonperturbative aspects of QCD. Historically, Υ → light hadrons has been employed to calibrate the strong coupling constant at the scale of bottom mass [2,3].
At lowest order, the inclusive hadronic decays of vector quarkonia proceed via J/ψ(Υ) → 3g, as demanded by conservation of C parity. This is very similar to the ortho-positronium (o-Ps) annihilation decay into three photons [4], so the corresponding expression can be directly transplanted supplemented with the proper color factor.
Since quarkonia are essentially the bound states composed of slowly-moving heavy quark and heavy antiquark, the consensus is nowadays that quarkonium annihilation decay processes can be reliably tackled in nonrelativistic QCD (NRQCD) factorization framework [5], which effectively organizes the theoretical predictions as double expansion in α s and v.
The O(α s ) perturbative corrections, which turn out to be negative, was first computed by Mackenzie and Lepage in 1981 [2]. In sharp contrast to the significantly negative O(α s ) corrections in J/ψ → 3γ, the O(α s ) correction in J/ψ(Υ) hadronic decay appears to be moderate in magnitude. On the other hand, the leading relativistic corrections to J/ψ(Υ) → 3g, 3γ were first calculated by Keung in 1982 [6]. With reasonable assumption of the relativistic NRQCD matrix elements for J/ψ, the O(v 2 ) corrections appear to be significantly negative for both J/ψ(Υ) inclusive hadronic decay channel and three-photon channel. Bodwin et al. proceeded further to compute the O(v 4 ) corrections to the J/ψ(Υ) hadronic decay [7]. It was observed that the relativistic corrections in the color-singlet channel exhibit decent convergence pattern, but for consistency one should also include contributions from various color-octet operator matrix elements at O(v 4 ). Unfortunately, the actual values of the higher-order NRQCD color-octet matrix elements are rather poorly known. The proliferation of many poorly constrained nonpertubative matrix elements severely hinders the predictive power of NRQCD approach 1 .
It is interesting to recall how a longstanding puzzle for J/ψ → 3γ is resolved. After incorporating both substantially negative O(α s ) and O(v 2 ) corrections, one simply ends up with a negative, hence unphysical, prediction to the partial width! Fortunately, the dilemma is greatly reconciled after including the joint perturtative and relativistic correction, e.g., the O(α s v 2 ) correction, which turns out to be positive and unexpectedly sizable [11]. Incorporating this new piece of correction, with the renormalization scale in the range 1.2 GeV < µ < 1.4 GeV, one can reach satisfactory agreement between the stateof-the-art NRQCD prediction [11] and the latest measurement at BESIII [12] for the rare decay channel J/ψ → 3γ. 1 Although the NRQCD matrix elements of the color-octet production operators have been fitted in various places [8][9][10], there does not exist rigorous relation between the color-octet production matrix elements and the decay matrix elements [5].
Inspired by the important role played by the O(α s v 2 ) correction to J/ψ → 3γ, the aim of this work is to investigate the O(α s v 2 ) correction to the J/ψ(Υ) hadronic widths. Obviously, the corresponding calculation is much more challenging than that for J/ψ → 3γ. We note that this joint radiative and relativistic correction is formally of the comparable magnitude with the O(α 2 s v 0 ) and O(α 0 s v 4 ) corrections. Nevertheless, unlike these two types of corrections, which are either beyond current calculational capability or lacking predictive power, we are equipped with sufficient technicality to tackle the O(α s v 2 ) correction, and can also make some unambiguous and concrete predictions. We hope our calculation can serve to further test the validity of NRQCD factorization approach in these basic vector quarkonia decay channels.
The remainder of the paper is organized as follows. In Section II, we recapitulate the NRQCD factorization formula for the hadronic decay of vector quarkonia and set up the notations. In Sec. III, we sketch the calculational strategy for the intended NRQCD SDCs. In Sec. IV we present the main numerical results. In particular, a detailed comparison with the existing O(α s ) corrections in a diagram-by-diagram basis has also been given. We devote Section V to a comprehensive phenomenological analysis. Finally we summarize in Section VI. We also dedicate Appendix A to describe the technicality in treating multi-body phase space integration in dimensional regularization.

II. NRQCD FACTORIZATION FORMULA FOR HADRONIC WIDTH OF VEC-TOR QUARKONIA
In line with the NRQCD factorization [5], the annihilation hadronic decay rate of a vector quarkonium V (V may refer to Υ or J/ψ) can be expressed as where LH denotes the abbreviation for light hadrons, and m Q is the heavy quark mass with quark species Q = c, b. The four-fermion NRQCD operators in (1) are defined by whose matrix elements are expected to obey the velocity counting rule and scale as v 3 and v 5 , respectively. Hence, the second term constitutes the O(v 2 ) relativistic correction to the first term in (1). The coefficients F and G in (1) are referred to as the NRQCD short-distance coefficients (SDCs), which capture the relativistic (k ∼ 1/m Q ) short-distance effects of QCD. Owing to asymptotic freedom, they can be reliably computed in perturbation theory: The various SDCs can be deduced via the standard perturbative matching technique. The key idea is because the SDCs like F , G are insensitive to the long-distance physics, we can replace the physical vector quarkonium state V in (1) with a fictitious quarkonium state carrying the same quantum number 3 S 1 , yet composed of a free heavy quark-antiquark pair. Concretely speaking, our matching equation reads We can calculate both sides of (4) using perturbative QCD and NRQCD, then iteratively solve for a variety of SDCs to a prescribed order in α s . Our main task in this work is to determine the unknown coefficient G 1 ( 3 S 1 ).

III. TECHNICAL STRATEGY OF CALCULATING NRQCD SDCS
In this section we sketch some important technicalities encountered in the perturbative matching calculation. For the fictitious vector quarkonium appearing in (4), we assign the momenta carried by the heavy quark Q and antiquark Q to be where P and q denote the total momentum and half of the relative momentum of the QQ pair, respectively. The on-shell condition indicates with E = m 2 Q − q 2 > m Q . Moreover, we assign the momentum of each massless parton in the decay products of V to be k i (i = 1, 2, 3 for three particles in the final state, and i = 1, 2, 3, 4 for four-body decay), which is subject to the on-shell condition k 2 i = 0. To facilitate the perturbative QCD calculation on the left-hand side of (4), it is convenient to employ the covariant projector to enforce the QQ pair to be in the spin-triplet and colorsinglet state. The nonrelativistically-normalized spin-triplet/color-singlet projector reads [7] where ǫ µ represents the spin-1 polarization vector. With the aid of (7), the amplitude for the spin-triplet/color-singlet QQ pair to annihilating into massless partons can be obtained through where A denotes the quark amplitude for QQ → LH with the external quark spinors amputated, and it is understood that the trace acts on both spinor and color indices. We proceed to project out the S-wave piece from (8). To our purpose, we need expand the amplitude A through the second order in q. The desired amplitude for QQ( 3 S 1 ) → LH can be expanded into with Here q 2 = −q 2 is defined in the rest frame of QQ pair, and the symmetric tensor is given by The above A S0 and A S2 represent the O(v 0 ) and O(v 2 ) amplitudes respectively. Squaring the S-wave amplitude in (9) and integrating over multi-body phase space, truncating through order-q 2 , we then accomplish the calculation of the quark-level decay rate on left side of (4). Nevertheless, a technical subtlety may be worth mentioning. Firstly, some portion of relativistic correction is hidden in the phase space integral, since the invariant mass of QQ is 4E 2 instead of 4m 2 Q . Secondly, the squared S-wave amplitude explicitly involves factors P · k i , which also depend on the variable E other than m Q . As pointed out in Refs. [11,13], in order to extract the relativistic correction, it may be more convenient to expand the quark amplitude (8) in power series of q 2 E 2 . This strategy guarantees that the relativistic effect is automatically taken into account in the phase space integral. Finally, we return to more conventional expressions by further expanding Another point in calculating the virtual correction is also worth mentioning. Prior to conducting loop integration, we have already expanded the integrand of the quark amplitude in power series of q 2 . This operation amounts to directly extracting the hard matching coefficient in the context of strategy of region [14]. Therefore, we are no longer concerned with Coulomb singularity and other soft contributions ubiquitously arising in virtual correction calculation. As a consequence, calculating the radiative correction to the NRQCD side of (4) is no longer required. It suffices to know the following perturbative NRQCD matrix elements at leading order (LO) in α s : In a sense, rather than literally follow perturbative matching calculation, we take an efficient shortcut to obtain the NRQCD SDCs.

B. Strategy of calculating virtual and real corrections
To coherently unify the treatment of the virtual and real corrections to QQ( 3 S 1 ) inclusive decay, we utilize the optical theorem to organize the calculation for O(α s ) correction. Concretely, we follow [2] to consider the imaginary part of the QQ( 3 S 1 ) forward-scattering amplitude, and place all possible cuts in accordance with Cutkosky's rule.
Throughout the work, we choose Feynman gauge to compute the quark amplitude, and adopt the dimensional regularization to regularize both UV and IR divergences. We use the It is tacitly understood that a Custkosky cut should be placed on all possible intermediate state for each diagram, since we are only interested in the imaginary part of the amplitude. We follow the same classification convention as Ref. [2]. The solid blob in class n) signifies the one-loop vacuum polarization diagrams that are composed of gluon and massless quarks.
package FeynArts [15] to generate Feynman diagrams and the corresponding QQ amplitudes. As illustrated in Fig. 1, all Feynman diagrams are categorized into a dozen of groups with distinct topology. The Cutkosky's cut is implicitly assumed and acts on all possible places. For each side of the cut, we apply the spin/color projector (7) together with (8), (9) and (10) to project out the amplitude for QQ( 3 S 1 ) transitioning into massless partons through O(v 2 ). The packages FeynCalc [16] and Form [17]/FormLink [18] are utilized for tensor contraction to obtain the squared S-wave amplitudes.
We utilize the packages Apart [19] to conduct partial fraction and FIRE [20] for the corresponding IBP reduction. We end up with obtaining 17 MIs for the real corrections and 58 MIs for the virtual corrections. For all the MIs, we use the packages FIESTA [21]/SecDec [22] to perform sector decomposition (SD) 2 . For each decomposed sector, we first use CubPack [26]/Cuba [27] to conduct the first-round rough numerical integration. For those integrals with large estimated errors, a parallelized integrator HCubature [28] is utilized to repeat the numerical integration until the prescribed accuracy is achieved.
Both UV and IR divergences could arise in the NLO correction calculation. The UV divergence is affiliated only to the virtual correction calculation. We implement the on-shell renormalization for the heavy quark wave function and mass, and renormalize the strong 2 To improve efficiency, we have interchanged the order of the operations for contour deformation [23], SD and series expansion in the FIESTA. For more technical details, we refer the interested readers to Refs. [24,25].
coupling constant under MS prescription. The virtual correction then becomes UV finite after standard renormalization procedure. Nevertheless, both virtual and real corrections are still plagued with IR divergences. At O(α s v 0 ), upon summing both contributions, the IR divergences are exactly cancelled [2]. As we shall see in Section IV, a peculiar pattern actually emerges at O(α s v 2 ). Upon summing virtual and real corrections, the IR singularities do not exactly cancel away, but leave a logarithmic IR divergence as the remnant. This single IR pole should in fact be factored into the NRQCD matrix element, so that one finally arrives at IR-finite yet scale-dependent SDC at O(α s v 2 ). Some remarks are in order for treating the multi-body phase space integration. Rather than directly integrating the squared amplitudes over the three(four)-body phase space, which appears to be a daunting task, we employ some modern technique, e.g., the reverse unitarity method [29,30] to simplify the calculation. The trick is to convert a phase-space integral into a loop integral, which is facilitated by the following identity involving the i-th cut propagator: Since the differentiation operation is insensitive to the ±iε in the propagator, one can fruitfully apply the integration-by-parts (IBP) identities, which are widely used in reducing the multi-loop integrals into a set of simpler Master Integrals (MIs), also to reduce the various phase-space integrals into a set of simple MIs. We then perform the corresponding multibody phase space integration numerically with these much simpler MIs. In Appendix A, we provide some details on how to parameterize the multi-body phase space integration in dimensional regularization.

IV. NUMERICAL RESULTS FOR NRQCD SDCS
In this section, we collect the known results of various SDCs in the NRQCD factorization formula for J/ψ(Υ) hadronic decay, and report the major new result of this work, the O(α s v 2 ) correction to SDC.

A. Leading-order NRQCD SDC
At LO in α s , the inclusive hadronic decay of the QQ( 3 S 1 ) pair is characterized by its annihilation into three gluons. It is a straightforward exercise to deduce the LO decay rate in d = 4 − 2ǫ spacetime dimensions: We have introduced the symbol Γ (i,j) to denote the perturbative decay rate computed at a specific order in α s and v expansion, e.g., with the superscript (i, j) signifying the joint O(α i s v 2j ) correction. For future usage, we have deliberately kept the O(ǫ) piece in (14).
Substituting (14) and (12a) into (4), and setting ǫ → 0, we then solve the SDC at O(α 0 s v 0 ): In deriving this, we have replaced E by m Q , which is valid at LO in v 2 . As expected, we have just reproduced the classical result [2,5,7].
Following the shortcut to extract the relativistic correction, as outlined in Section III A, we find the O(α 0 s v 2 ) correction to the perturbative decay rate to be To determine the SDC G 1 , we need also expand E in (16) to first order in q 2 m 2 Q . Substituting (16) and (12b) into (4), with the aid of the known LO SDC F (0) 1 in (15), we then deduce the SDC at O(α 0 s v 2 ): which agrees with Refs. [5][6][7].
The QCD radiative correction to the vector quarkonia hadronic decay, yet at LO in velocity expansion, was first computed by Mackenzie and Lepage nearly forty years ago [2]. Dividing all the cut diagrams into several groups of distinct topology, the authors of [2] tabulate the numerical result of each individual class of diagrams, some of which have rather limited accuracy. In order to make a close comparison with their result, we also adopt the same classification convention of Feynman diagrams as [2].
In Table I, we tabulate the numerical value of Γ (1,0) /(α s Γ (0,0) /π) for each class of cut diagrams. Following [2], we also include in the class c) the diagrams with insertion of heavy quark mass counterterm. Moreover, the entry "C. T. (δα s )" in Table I signifies those tree diagrams with insertion of the counterterm for the strong coupling constant. The results from Ref. [2] are also juxtaposed with our results, which employed the gluon mass λ to regularize IR divergences. To expedite the comparison, we replace log λ 2 /m 2 Q in [2] with 1/ǫ IR + log µ 2 /m 2 Q . As can be clearly observed from the Table I, we find perfect agreement between our and their results for each class of cut diagrams. Since we are equipped with the modern IBP method tailored for multi-loop (multi-body phase space) integrals and the more accurate non-Monte-Carlo integrator, it is conceivable that a much better numerical accuracy can be achieved with respect to that in [2].
Summing up the contributions from each class of cut diagrams, we end up with the IR-finite O(α s v 0 ) correction to the perturbative decay rate:  v 0 ) contribution to the Γ (1,0) /(α s Γ (0,0) /π) from each class of cut diagrams, with our results juxtaposed with the counterparts given in [2]. We have shifted the 't Hooft unit mass according to µ 2 → µ 2 e γ E /(4π), which is equivalent to renormalizing α s in the MS scheme. For simplicity, we temporarily set the renormalization scale µ R = 2E.
Cut diagrams Virtual corr. Real corr. Real+Virtual Corr.  where β 0 = 11/3C A − 2/3n f denotes the one-loop coefficient of the QCD β function, with n f signifying the number of the active light flavors. In this work, we take n f = 3 for J/ψ decay and n f = 4 for Υ decay. µ R in (18) represents the renormalization scale, whose explicit occurrence is demanded by the renormalization group invariance. Substituting (18) and (12a) into (4), setting E = m Q , we then deduce the SDC at O(α 0 s v 0 ): which is compatible with the result in literature [2,5,7], yet bears a much smaller error.
We proceed to determine the numerical value of G In Table II, we tabulate the individual contribution from each class of cut diagrams to the perturbative hadronic decay rate. The numbering convention for the Feynman diagrams is the same as Table I. Upon summing all the individual contribution in Table II, and implementing standard renormalization procedure, we obtain the following inclusive decay Clearly, the occurrence of G (0) 1 ( 3 S 1 )β 0 ln µ R reflects the µ R -independence of the decay rate, which is demanded by the standard renormalization group equation.
A salient trait in (20) is the occurrence of an uncancelled single IR pole. We assign the symbol µ Λ to denote the scale accompanied with this single IR pole, which has very different origin from the QCD renormalization scale µ R .
The uncanceled single IR pole in (20), which arises at the hard region at O(α s v 2 ), is actually not surprising at all. In fact, this symptom was first discovered in the O(α s v 2 ) corrections for J/ψ → e + e − [31] (see also [32]), η c → γγ [33], J/ψ → 3γ [11], and η c total hadronic width [34]. Physically, this unremoved IR pole indicates the breakdown of the color transparency once beyond the LO in v. This IR divergence can be factored in the renormalized NRQCD matrix element, so that the O(α s v 2 ) SDC will explicitly depend on the NRQCD factorization scale µ Λ , which naturally ranges from m Q v to m Q .
Expanding E in (20) around m Q in power series of q 2 /m 2 Q , then comparing (20) with (4), we finally obtain the desired SDC in MS factorization scheme: Equation (21) constitutes the major new finding of this work. It is interesting to note that the SDC at O(α s v 2 ) contains the same ln µ Λ term as the case for J/ψ → 3γ [11]. However, (21) also contains a negative non-logarithmic constant, in sharp contrast with what is found for J/ψ → 3γ [11].

V. PHENOMENOLOGICAL ANALYSIS
Having all the SDCs through O(α s ) at hand, we are ready to conduct a detailed phenomenological analysis based on the NRQCD factorization formula (1).

A. Phenomenology for J/ψ(Υ) → LH
Before launching into concrete numerics, we would like to first develop some intuition on the relative importance of the various corrections. Setting µ R = µ Λ = m Q , we then obtain To condense the notation, we have introduced the dimensionless ratio v 2 V to characterize the relativistic correction: To make concrete predictions, we need further specify various input parameters: the value of heavy quark masses, various NRQCD matrix elements, as well as the running strong coupling constant evaluated around the quarkonium mass. We choose the quark pole mass to be m c = 1.4 GeV and m b = 4.6 GeV [11]. We take the same value of NRQCD matrix elements for J/ψ as in Ref. [35]: We also take the values of the NRQCD matrix elements for Υ(nS) from Ref. [36] 3 : We use the package RunDec [38] to compute the QCD running coupling constant to twoloop accuracy: Our predictions for the hadronic widths of J/ψ and Υ(1S, 2S) are tabulated in Table III As one easily recognizes from Table III, a salient feature for J/ψ → LH is the disquietingly negative O(α 0 s v 2 ) correction, which is even greater than LO prediction in magnitude. Even worse, the O(α s v 0 ) radiative correction further decreases the predicted hadronic width. Although the new O(α s v 2 ) correction tends to increase the prediction, unfortunately its effect is not significant enough to bring the predicted hadronic width to a positive value, so that we are unable to make a meaningful prediction to confront the measurement. same as in (25a), while on the right panel we take v 2 Υ(1S) = 0.056 in (27), which is determined through G-K relation. To closely examine this apparent discrepancy between theory and experiment, we further display in Fig. 2 the dependence of J/ψ hadronic decay rate on the renormalization scale µ R at various level of accuracy in α s and v expansion, with the NRQCD factorization scale µ Λ = m c held fixed. In the reasonable range of 1 GeV < µ R < 2m c , we again observe that the total decay rate is always negative. As mentioned before, the root of this dilemma can be readily traced from (22a), that is, mainly due to the substantially negative O(α 0 s v 2 ) and O(α s v 0 ) corrections and moderately positive O(α s v 2 ) correction. Thus, we conclude that the state-of-the-art NRQCD prediction ceases to yield a physically meaningful result for J/ψ → LH. How to account for the experimental data from the perspective of NRQCD factorization remains an open challenge. As a potentially appealing solution, one may conjecture that the uncalculated O(α 2 s v 0 ) correction might be significantly positive, which would bring the NRQCD prediction closer to measured value. Unfortunately, the bottleneck of the current multi-loop(leg) calculational capability impedes one to tackle this type of NNLO perturbative correction in the foreseeable future. From Table III, one observes that the O(α s v 2 ) correction has rather minor effect for Υ(1S, 2S) → LH, much less significant than the O(α s v 0 ) and O(α 0 s v 2 ) corrections. In Fig. 3 and Fig. 4, the hadronic decay rates of Υ(1S, 2S) are displayed as a function of the renormalization scale µ R . As can be easily visualized, though the LO NRQCD prediction exhibits strong µ R dependence, after including the corrections through O(α s v 2 ), the decay rate becomes much less sensitive to the variation of µ R , especially when µ R > m b .
The relativistic corrections in Υ(1S) hadronic decay appear to be negligible, which can be attributed to the highly suppressed relativistic NRQCD matrix elements, as clearly seen in (25a). We choose two different values for the relativistic NRQCD matrix element to make predictions for Υ(1S) → LH. When choosing the values specified in (25a), the finest NRQCD prediction for Υ(1S) hadronic width is considerably larger than the measured value, which can be readily seen in Table III and the left panel of Fig. 3.
Alternatively, one can estimate the v 2 Υ(1S) according to the Gremm-Kapustin (G-K) relation [40]: once m b = 4.6 GeV is chosen. With such choice, the O(α s v 2 ) correction may reach 10% of the LO decay rate. Consequently, for m b ≤ µ R ≤ 2m b , the hadronic width of Υ(1S) is predicted to be in the range 44 − 51 keV. Interestingly, as one can tell from the right panel of Fig. 3, the finest NRQCD prediction is in satisfactory agreement with the measured value: 44.1 ± 1.4 keV [39]. The knowledge of the O(α s v 2 ) correction enables us not only to present the most complete NRQCD predictions for hadronic widths of vector quarkonia, but present the finest NRQCD predictions for the rare electromagnetic decay J/ψ(Υ) → 3γ.
The O(α s v 2 ) corrections to the partial width of J/ψ(Υ) → 3γ were addressed by the authors some time ago [11]: where e Q represents the electric charge of the heavy quark, and α denotes the fine structure constant. Curiously, the ln µ Λ term in the O(α s v 2 ) SDC is identical to that in G 1 ( 3 S 1 ) in (21). has not yet been observed experimentally. We take v 2 Υ(1S) = 0.056 in (27), which is inferred from G-K relation.
Dividing (28) by (1), and plugging the result for various SDCs through O(α s v 2 ) that are tabulated in Section IV, we then obtain the state-of-the-art NRQCD predictions for the branching fractions of J/ψ(Υ) → 3γ. In the spirit of NRQCD expansion, we further expand the ratio in power series of α s and v 2 and obtain To reduce the theoretical error, we take Br(J/ψ → LH) = 64.1% and Br(Υ(1S) → LH) = 81.7% as experimental input [39]. As is evident in (29), The LO NRQCD matrix element has canceled in the ratio. More interestingly, the O(α 0 s v 2 ) correction has also completely disappear in the ratio. Furthermore, the explicit logarithmic dependence on NRQCD factorization scale µ Λ has also disappeared.
In Fig. 5 we display our predictions for the branching fractions of J/ψ → 3γ (left panel) and Υ(1S) → 3γ (right panel) at various levels of accuracy in α s and v expansion, with µ Λ = m Q . The sensitivity of the branching fractions to µ R seems not to improve much after incorporating the higher order corrections. This may be attributed to the sizeable radiative corrections and the α 3 s factor in the denominator. Setting µ R = m c , the finest NRQCD prediction is Br(J/ψ → 3γ) = 1.45 × 10 −6 , about one order of magnitude smaller than the measured value 1.16 × 10 −5 . As can be seen in Fig. 5, varying the renormalization scale in the range 1 GeV ≤ µ R ≤ 3 GeV, we always observe that Br(J/ψ → 3γ) is far smaller than the experimental measurement. The resolution to this dilemma may need call for including further higher order relativistic corrections [35], or incorporating the O(α 2 s v 0 ) correction. Such a study appears to be extremely challenging on technical ground, which is beyond the scope of this work.
In a similar fashion, with the aid of (29b), we can predict Br(Υ(1S) → 3γ) = 2.06 +0.20 −1.55 × 10 −7 and Br(Υ(2S) → 3γ) = 1.78 +0.11 −0.99 × 10 −7 , where the uncertainties are estimated by varying µ R from m b /2 to 2m b . The uncertainties looks significant due to strong µ R dependence. Alternatively, if we choose v 2 Υ(1S) = 0.056 as indicated by the G-K relation, the branching fraction shifts to Br(Υ(1S) → 3γ) = 2.33 +0.18 −1.44 × 10 −7 . It is a clear sign that the relativistic effect is less important for bottomonium case, and the branching fraction appears to be insensitive to the relativistic corrections. We hope future experiments for Υ rare electromagnetic decay can test our predictions.

VI. SUMMARY
In summary, we have computed the O(α s v 2 ) corrections to J/ψ(Υ) hadronic decays, precisely deducing the corresponding SDCs. Unlike the case in J/ψ(Υ) → 3γ, the corrections are moderate, therefore the expansion in v 2 exhibits a decent convergence behavior. We find that the theoretical predictions for the decay width of J/ψ → LH is negative through O(α s v 2 ), which is certainly unphysical and can be attributed to the sizable and negative relativistic corrections. On the other hand, we find the theoretical predictions for Υ(1S, 2S) hadronic decays are consistent with the experimental data, if the NRQCD matrix elements v 2 Υ(nS) are appropriately chosen. Incorporating the formulas in the Ref. [11], we derive the branching fraction for J/ψ(Υ) → 3γ accurate up to O(α s v 2 ). Although the O(v 2 ) relativistic corrections get exactly cancelled, we observe a strong µ R dependence in the branching ratio. Once again, we find the theoretical prediction for Br(J/ψ → 3γ) can not explain the experimental measurement. In our opinion, the theoretical incompetence to precisely predicting J/ψ decay in NRQCD factorization approach may be attributed to the not-yet-known higher-order relativistic and perturbative corrections.
We start with the Born-order process for J/ψ(Υ) → LH, where the decay product is comprised of three gluons. We define the following Lorentz invariants: where the momenta P and k i (i = 1, 2, 3) have been defined in Section III A. The dimensionless variables x 1 and x 2 are constrained to be 1 ≤ x 1 + x 2 ≤ 2 by energy conservations. It is useful to change variables from (x 1 , x 2 ) to (x, y) so that the integration range of the new variables lie in a square [41]: with 0 ≤ x, y ≤ 1. It is then straightforward to parameterize the phase space integral of three massless particle in d = 4 − 2ǫ space-time dimensions as (A3) Calculation of the real corrections demands considering the J/ψ(Υ) decay into four massless partons. The four-particle phase space integration in dimensional regularization is somewhat intriguing. We employ the parametrization scheme introduced in Ref. [42]. The Mandelstam variables are defined as s ij = k i · k j , which obey s 12 + s 13 + s 14 + s 23 + s 24 + s 34 = 4E 2 constrained by energy conservation. It is convenient to introduce a set of dimensionless variables x i through rescaling the Mandelstam invariants by In terms of the dimensionless variables in (A4), the massless four-body phase space can be parameterized as [42] where λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + xz + yz) is the Källen function, Θ represents the Heaviside step function, and the volume factor C Γ is given by with V (d) = 2π d/2 Γ(d/2) designating the area of the unit sphere imbedded in d-dimensional space.
The presence of Heaviside step function in (A5) generally renders the integration boundary highly irregular, so that one has to resort to the Monte Carlo recipe for numerical integration, with some intrinsic limitation on integration accuracy. In fact, we can manage to eliminate the Heaviside function through a clever trick, so that we can employ more accurate numerical integrator other than the Monte Carlo algorithm.
First let us explicitly write down the Källen function in (A5), abbreviated by λ for simplicity: With the aid of Cheng-Wu theorem [43][44][45], we have the freedom to pull the variables x 1 , x 2 and x 3 outside the δ-function when carrying out phase space integration (A5). Consequently, the integration range of x 1 , x 2 and x 3 then become [0, ∞), essentially unbounded. We can freely rename the variables: Now the Källen function reduces to We can divide the integration range for x 1 and x 2 into two sectors: x 1 ≥ x 2 and x 1 < x 2 . For the first sector x 1 ≥ x 2 , it is convenient to further change the variables: For the second sector x 1 < x 2 , we can also make variable change: The advantage of changing variables is to ensure that the support of the new variables are within a square, e.g., both integration ranges of new x 1 and x 2 lie in [0, ∞). Without loss of generality, we choose the first sector to illustrate our recipe. After changing variables in line with (A10), the Källen function now simplifies to Similarly, we can further break the integration into two sectors: x 1 ≥ x 3 and x 1 < x 3 . For the x 1 ≥ x 3 sector, we are free to continue to rename the variables: And for the x 1 < x 3 sector, we change the variables as Now the integration range of the new x 1 and x 3 variables lie in [0, ∞). We can repeat the game of dissecting the integral into two sub-sectors. For concreteness, let us concentrate on the sub-sector x 1 ≥ x 3 . After changing variables as advised in (A13), the Källen function now reduces to We further rename the variable x 2 as so that λ in (A15) reduces to One more again, we can dissect the integration range into two sectors: x 1 ≥ x 2 and x 1 < x 2 . As indicated by the Heaviside function Θ(−λ) in (A5), the first sub-sector (x 1 ≥ x 2 ) apparently make vanishing contribution. Therefore, only the second sub-sector with x 1 < x 2 survives. We further change the variables as follows: So the integration range of the new x 1 and x 2 variables also lies in [0, ∞). Meanwhile, the Källen function reduces to λ = − x 1 x 2 2 . Since the Heaviside function Θ(−λ) remains 1 in this prescribed integration support, it can be safely dropped. This sequences of manipulation can be recursively applied to any integration sector.
At the last step, we can invoke Cheng-Wu theorem again to put x 1 , x 2 and x 3 back inside the δ-function.
Through the aforementioned manipulations, we are able to eliminate the Heaviside Θ function from the four-particle phase space integral. Since the support of the multidimensional integration variables is inside a regular hypercube, we can resort to any efficient numerical integrator, e.g., the parallelized integrator HCubature [28] for high-precision numerical integration.
As a test, we apply sector decomposition method together with our parametrization to numerically calculate five MIs of the massless four-body phase-space type (labelled by R 4 , R 6 , R 8,a , R 8,b and R 8,r in [42]). These five MIs have been analytically worked out in 2003 [42], with which we find exquisite numerical agreement.