The $P$-wave charmonium annihilation into two photons $\chi_{c0, c2}\rightarrow \gamma\gamma$ with high-order QCD corrections

In this paper, we present a new analysis on the $P$-wave charmonium annihilation into two photons up to next-to-next-to-leading order (NNLO) QCD corrections by using the principle of maximum conformality (PMC). The conventional perturbative QCD prediction shows strong scale dependence and deviates largely from the BESIII measurements. After applying the PMC, we obtain a more precise scale-invariant pQCD prediction, which also agrees with the BESIII measurements within errors, i.e. $R={\Gamma_{\gamma\gamma}(\chi_{c2})} /{\Gamma_{\gamma\gamma}(\chi_{c0})}=0.246\pm0.013$, where the error is for $\Delta\alpha_s(M_\tau)=\pm0.016$. By further considering the color-octet contributions, even the central value can be in agreement with the data. This shows the importance of a correct scale-setting approach. We also give a prediction for the ratio involving $\chi_{b0, b2} \to\gamma\gamma$, which could be tested in future Belle II experiment.

Charmonium decays have been widely used to explore the interplay between the perturbative and nonperturbative dynamics due to its relatively clean platform, and they also play important roles in establishing the asymptotic freedom of quantum chromodynamics (QCD) [1,2]. Among them, many attentions have been paid for the electromagnetic decays χ c0,c2 → γγ. They have been measured by the CLEO and BESIII collaborations [3,4]; especially, in year 2017, the BESIII collaboration issued their measured value for the R-ratio, R exp = Γ χc2→γγ Γ χc0→γγ = 0.295 ± 0.014 ± 0.007 ± 0.027, (1) where the errors are statistical, systematic, and the associated errors of the branching fraction B(ψ(3686) → γχ c0,c2 ) and the total decay width Γ χc0,c2 , respectively. On the other hand, they have been calculated by using various approaches, such as the nonrelativistic potential model, the nonrelativistic QCD theory (NRQCD), the relativistic quark model, and the lattice QCD theory, respectively, c.f. Refs. [5][6][7][8][9][10][11][12][13][14][15][16] and references therein. Within the framework of NRQCD factorization theory, one has observed that the leading-order (LO) prediction is close to the experimental measurements, but this process is extremely sensitive to high-order QCD corrections and relativistic corrections, due to the fact that the typical magnitude of strong coupling constant and the squared relative velocity of the charm quark in charmonium, α s (m c ) ∼ v 2 c ∼ 0.3, are comparatively large. It is important to finish as more perturbative terms as possible so as to achieve a more accurate pQCD prediction. And in order to obtain a convincing fixed-order prediction, the influence of high-order correction on the χ c0,c2 → γγ must be carefully analyzed.
At the present, the QCD corrections to the S-wave heavy quarkonium electromagnetic/leptonic electromagnetic decays have been calculated up to next-to-nextto-leading-order (NNLO) level. The spin-singlet heavy quarkonium decays η c → γγ and η b → γγ have been calculated up to NNLO level in Refs. [17,18]; and the spin-triplet heavy quarkonium decays J/ψ → e + e − and Υ → e + e − decay have been calculated up to NNLO level by Refs. [19,20]. In year 2016, the NNLO QCD corrections to the P -wave charmonium χ c0,c2 → γγ have been done by Ref. [16], which however show large renormalization scale dependence, and the predicted R-ratio cannot explain the above BESIII value. It is important to show what's the reason for such discrepancy.
Within the framework of NRQCD, one can factorize the decay width into non-perturbative matrix elements and perturbatively calculable short-distance coefficients, and the R-ratio becomes where the helicity amplitude of this process is expressed by A χcJ λ1,λ2 , λ = |λ 1 − λ 2 |, λ 1,2 = ±1, J = 0, 2. The amplitude of the P -wave quarkonium electromagnetic decays χ c0,c2 → γγ can be expressed as [16] A χc0,c2 where the color-singlet P -wave long-distance matrix element can be related to the first derivative of the radial wave function at the origin, where N c = 3 is the SU c (3) color number, R ′ χc0,c2 (0) are first derivatives of χ c0,c2 radial wavefunctions at the origin. The spin-splitting effect for the radiation wavefunctions of χ c0,c2 are small, and by defining the Rratio (2), the uncertainties caused by the matrix elements can be greatly suppressed. The perturbative part C χc0,c2 λ (m c , µ r , µ Λ ) up to NNLO level can be read from Ref. [16], where µ r and µ Λ are renormalization and factorization scales, respectively.
In dealing with the perturbative series of R-ratio, one usually sets µ r = m c so as to eliminate large logarithmic terms in powers of ln µ 2 r /m 2 c , and then varies it within certain range to ascertain the uncertainty. This simple method causes the mismatching of α s with its perturbative coefficients at each order, breaks the renormalization group invariance [21], and leads to conventional renormalization scheme-and-scale ambiguities. Such ambiguities could be softened to a certain degree by including higher-order terms. However due to its complexity, the exact NNNLO corrections of χ c0,c2 → γγ shall not be available in near future, thus it is important to find a correct way to achieve a reliable and accurate prediction by using the known NNLO series.
The renormalization scale-setting problem is one of the most important issues for pQCD theory, which has a long history, cf. the review [22]. To solve it, we adopt the single-scale approach [23] of the principle of maximum conformality (PMC) to analyze the decay width of χ c0,2 → γγ up to NNLO QCD corrections. By using the renormalization group equation recursively, the PMC determines the precise α s value of the process by using the non-conformal β-terms in pQCD series [24][25][26][27][28]. After applying the PMC, the resultant pQCD series becomes conformal, the magnitude of α s and the perturbative coefficients become well matched, and then we obtain exact values for each order. The PMC predictions are renormalization scheme-and-scale independent [29], thus the conventional scale-setting ambiguities is eliminated. Due to the perturbative nature of the pQCD theory, there is residual scale dependence because of unknown higherorder terms [30]; For the PMC perturbative series, such residual scale dependence shall generally be highly suppressed even for lower-order predictions [31].
One can rewrite the R-ratio (2) of χ c0,2 → γγ as the following perturbative form, where a s = α s /4π and Ω = |R The perturbative coefficients r i can be derived into conformal terms r i,0 and non-conformal terms r i,j =0 by using the degeneracy relations among different orders [32], i.e., where β 0 = 11 − 2 3 n f , representing the one-loop βfunction, in which n f is the active flavor numbers. Using the NNLO results given in Ref. [16], we obtain (10) where , and e c represents the charm-quark electric charge.
After applying the standard procedures of the PMC single-scale approach [23] to the pQCD series (5), we obtain the following conformal series, where Q * is the PMC scale, which is obtained by requiring all non-conformal items vanish. It is the effective scale which replaces the individual PMC scales at each order in PMC multi-scale approach [27,28] in the sense of a mean value theorem. At present, by using the known NNLO perturbation series, the PMC scale can be fixed at the leading-log accuracy: wherer i,j = r i,j | µr =mc . It is noted that Q * is independent to any choice of renormalization scale µ r , together with the scale-invariant conformal coefficients, the conventional renormalization scale ambiguity is eliminated. Using Eq. (12), we obtain Q * = 0.250 GeV, which is close to the QCD asymptotic scale. Because the effective momentum flow Q * of the process is close to the QCD asymptotic scale Λ, we need to choose a low-energy model for α s so as to achieve a reliable prediction. In the literature, a variety of low-energy α s models have been suggested [33][34][35][36][37][38][39][40][41]. A comparison of various low-energy α s models has been given in Ref. [42]. In the present paper, for clarity, we adopt the CON model to do our discussion. It is derived from continuum theory [41] and uses the exchanged gluons with an effective dynamical mass m g , and determines the non-perturbative dynamics of gluons by using the Dyson-Schwinger equation. More explicitly, the CON low-energy α s model is expressed as follows: where , and m g = 500 ± 200 MeV [35]. The asymptotic scale Λ can be fixed by using the α s measured at a typical energy scale. More definitely, by using α MS s (M τ ) = 0.325 ± 0.016 [43], which leads to Λ| n f =3 = 0.383 +0.029 −0.031 GeV, Λ| n f =4 = 0.324 +0.029 −0.029 GeV, and Λ| n f =5 = 0.223 +0.022 −0.023 GeV. Fig.1 shows the α s -running behavior at different scales, the α s CON-model with m g = 500 +200 −200 MeV is adopted in low-energy region which is shown by shaded band. The smooth connection between the low-energy region and the large-energy region is obtained by using the matching scheme proposed in Ref. [44]. More explicitly, by requiring the first derivatives of α s to be the same at the crossing point of the two energy regions, the α s transition scale is determined to be 0.933 +0.183 −0.191 GeV.  To do the numerical calculation, we take the c-quark pole mass m c = 1.68 GeV [45], and set the factorization µ Λ = 1 GeV. In Table I, we give the contributions of each loop terms of R c under conventional and PMC scale-setting approaches, respectively. The conventional predictions are highly µ r -dependent for each loop terms and the total contributions of R c , e.g. as shown by the following Eqs. (14)(15)(16), the renormalization scale uncertainty for R Conv. c,total within the range of µ r ∈ [1GeV, 2m c ] is about ( −99% +57% ): R Conv. c,total | µr=mc = 0.067 GeV, R Conv. c,total | µr =2mc = 0.105 GeV.
Those values deviate from the BESIII measurement [4] by at least ∼ 3.9σ. Moreover, the separate scale un-certainties for NLO-terms and NNLO-terms are ( +59% −27% ) and ( −60% +12% ), respectively. After applying the PMC, the conventional scale uncertainty is removed, i.e. we obtain R PMC c,total ≡ 0.246 for any choice of µ r , which agrees with the BESIII measurement within errors. By further taking the α s shift, ∆α s (M τ ) = ±0.016, into consideration, we obtain ∆R Conv. c,total | µr =1 GeV = +0.026 −0.029 , ∆R Conv. c,total | µr =mc = ±0.013, ∆R Conv. c,total | µr =2mc = ±0.007, and ∆R PMC c,total = ±0.013.  Fig. 2 shows explicitly how the R c -ratio changes with different choices of µ r . It shows that after applying the PMC, the perturbative series is independent to any choice of µ r , and the conventional renormalization scale ambiguity is removed. The PMC conformal series ensures the scheme independence of the pQCD prediction, and together with the scale invariance, its behavior indicates the intrinsic perturbative behavior of the series. The NNLO conformal coefficient r 2,0 = 126.74 is larger than the conventional coefficient r 2 = −59.93 +50.84 −67.92 for µ r ∈ [1GeV, 2m c ], this explains why the PMC magnitude of the NNLO-terms is larger than the conventional one. It is noted that the smaller conventional NNLO coefficient r 2 is due to accidental cancelation of conformal and non-conformal terms, which is however highly scale dependent, leading to scale dependent series.
A physical observable should not dependent on the choices of manly introduced parameters such as the renormalization scale and the factorization scale. In the above, we have shown that by applying the PMC, the NNLO R-ratio removes the conventional renormalization dependence. Due to the heavy quark spin symmetry, the matrix elements for χ c0 and χ c2 are the same, and R-ratio avoids the uncertainties caused by different choices of matrix elements. However, because the NNLO-coefficient r 2,0 is explicitly µ Λ -dependent, the R-ratio shows explicit µ Λ dependence, whose magnitude is large [16]. We observe that such large factorization scale dependence could be removed by taking the evolution effects of the matrix element into consideration; That is, because the anoma- lous dimensions for χ c0 and χ c2 are different [48], those two matrix elements are different at different µ Λ , which could compensate the µ Λ -dependence from the hard part.
Using the evolution equation [9], we obtain ln . (18) As required, the difference is an O(α 2 s )-order effect. If taking O 1 ( 3 P J )| µΛ0=1GeV = 0.107 GeV 5 [49] as the initial value, we obtain O 1 ( 3 P 0 )| µΛ=3GeV = 0.139 GeV 5 and O 1 ( 3 P 2 )| µΛ=3GeV = 0.124 GeV 5 . Though the differences are small, we observe that the net factorization scale dependence of R c can be removed. This can be explicitly shown by Fig.3, in which the flat lines for R c versus µ Λ indicates that the R c -ratio is independent to any choice of µ Λ . Fig.3 shows that after applying the PMC, the predicted R c becomes more closer to the BESIII value, which still has a slight gap from the data.
According to Refs. [13,50,51], contributions from the color-octet (CO) components in charmonium should not be ignored. For example, it has been argued that the color-octet components shall shift the decay width by ∆Γ CO χc0 ≃ −0.3 GeV and ∆Γ CO χc2 ≃ −0.227 GeV [13]. If taking those extra color-octet contributions into consideration, we obtain R PMC c,total = 0.246 to 0.299. Fig. 4 shows the R c -ratio with or without the color-octet contributions, where the error bars are for ∆α s (M τ ) = ±0.016 and µ r ∈ [1GeV, 2m c ]. These results show a better match with the experimental data. Thus the CO contributions should be taken into consideration for a sound prediction.   Three typical µr are adopted, and the PMC prediction is independent to those choices of µr.
As a final remark, the above analysis can be conveniently applied for the P -wave bottomonium decays to two photons, which could be measured by the future high precision Belle-II experiment. We present the NNLO R b -ratio under conventional and PMC scale-setting ap-proaches in Table II. To do the numerical calculation, we take the b-quark pole mass m b = 4.78 GeV [45]. Due to α s (m b ) ∼ 0.1, a better pQCD convergence is observed for the bottomonium case, and the scale uncertainty is smaller, i.e. the net NNLO scale error is about 29% for µ r ∈ [1 GeV, 2m b ].
As a summary, in the present paper, we have studied the R c -ratio up to NNLO accuracy. Under conventional scale-setting approach, the renormalization scale uncertainty is large, which is about ( −99% +57% ) for µ r ∈ [1 GeV, 2m c ]. By applying the PMC, we obtain a more accurate pQCD prediction without renormalization scale uncertainty, R c | PMC = 0.246 ± 0.013, whose error is caused by ∆α s (M τ ) = ±0.016. This prediction agrees with the latest BESIII data within errors.