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

Abstract

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\pm 0.013$, where the error is for $\Delta \alpha _s(M_\tau )=\pm 0.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} \rightarrow \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 non-perturbative 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 $\chi _{c0, c2}\rightarrow \gamma \gamma$. 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,

$$\begin{aligned} R_\mathrm{exp}= & {} \frac{\Gamma _{\chi _{c2\rightarrow \gamma \gamma }}}{\Gamma _{\chi _{c0\rightarrow \gamma \gamma }}}=0.295\pm 0.014\pm 0.007\pm 0.027, \end{aligned}$$

(1)

where the errors are statistical, systematic, and the associated errors of the branching fraction $\mathcal{B}(\psi (3686)\rightarrow \gamma \chi _{c0, c2})$ and the total decay width $\Gamma _{\chi _{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, $\alpha _s(m_c)\sim v^2_c\sim 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 $\chi _{c0, c2}\rightarrow \gamma \gamma$ 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-next-to-leading-order (NNLO) level. The spin-singlet heavy quarkonium decays $\eta _{c}\rightarrow \gamma \gamma$ and $\eta _{b}\rightarrow \gamma \gamma$ have been calculated up to NNLO level in Refs. [17, 18]; and the spin-triplet heavy quarkonium decays $J/\psi \rightarrow e^+ e^-$ and $\Upsilon \rightarrow 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 $\chi _{c0, c2}\rightarrow \gamma \gamma$ 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

$$\begin{aligned} R_c = \frac{\left( |A^{\chi _{c2}}_{1,1}|^2+|A^{\chi _{c2}}_{1,-1}|^2\right) }{5(|A^{\chi _{c0}}_{1,1}|^2)}, \end{aligned}$$

(2)

where the helicity amplitude of this process is expressed by $A^{\chi _{cJ}}_{\lambda _1,\lambda _2}$, $\lambda =|\lambda _1-\lambda _2|$, $\lambda _{1,2}=\pm 1$, $J=0,2$. The amplitude of the P-wave quarkonium electromagnetic decays $\chi _{c0, c2}\rightarrow \gamma \gamma$ can be expressed as [16]

$$\begin{aligned} A^{\chi _{c0, c2}}_{\lambda _1,\lambda _2} = C^{\chi _{c0, c2}}_{\lambda }(m_c,\mu _r,\mu _\Lambda ) \frac{\left\langle 0|\chi ^{^\dag }K_{3_{P_{0,2}}}\Psi | \chi _{c_{0,2}} \right\rangle _{\mu _\Lambda }}{m^{3/2}_c}, \end{aligned}$$

(3)

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,

(4)

where $N_{c}=3$ is the $\mathrm{SU}_{c}(3)$ color number, $\overline{R}^{\prime }_{\chi _{c0, c2}}(0)$ are first derivatives of $\chi _{c0, c2}$ radial wavefunctions at the origin. The spin-splitting effect for the radiation wavefunctions of $\chi _{c0, c2}$ are small, and by defining the R-ratio (2), the uncertainties caused by the matrix elements can be greatly suppressed. The perturbative part $C^{\chi _{c0, c2}}_{\lambda }(m_c, \mu _r, \mu _{\Lambda })$ up to NNLO level can be read from Ref. [16], where $\mu _r$ and $\mu _\Lambda$ are renormalization and factorization scales, respectively.

In dealing with the perturbative series of R-ratio, one usually sets $\mu _r = m_c$ so as to eliminate large logarithmic terms in powers of $\ln \mu _r^2/m_c^2$, and then varies it within certain range to ascertain the uncertainty. This simple method causes the mismatching of $\alpha _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 $\chi _{c0, c2}\rightarrow \gamma \gamma$ 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 $\chi _{c0,2}\rightarrow \gamma \gamma$ up to NNLO QCD corrections. By using the renormalization group equation recursively, the PMC determines the precise $\alpha _s$ value of the process by using the non-conformal $\beta$-terms in pQCD series [24,25,26,27,28]. After applying the PMC, the resultant pQCD series becomes conformal, the magnitude of $\alpha _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 higher-order 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 $\chi _{c0,2}\rightarrow \gamma \gamma$ as the following perturbative form,

$$\begin{aligned} R_c= & {} \Omega [1+r_{1}a_s(\mu _r)+r_{2}a^2_s(\mu _r)+\mathcal{O}(a_s^3)], \end{aligned}$$

(5)

where $a_s={\alpha _s}/{4\pi }$ and $\Omega =|{\overline{R}^{\prime }_{\chi _{c2}}} / {\overline{R}^{\prime }_{\chi _{c0}}}|^2 \simeq 1$. The perturbative coefficients $r_i$ can be derived into conformal terms $r_{i,0}$ and non-conformal terms $r_{i,j\ne 0}$ by using the degeneracy relations among different orders [32], i.e.,

$$\begin{aligned} r_1= & {} r_{1,0} , \end{aligned}$$

(6)

$$\begin{aligned} r_2= & {} r_{2,0}+r_{2,1}\beta _0 , \end{aligned}$$

(7)

where $\beta _0=11-\frac{2}{3}n_f$, representing the one-loop $\beta$-function, in which $n_f$ is the active flavor numbers. Using the NNLO results given in Ref. [16], we obtain

$$\begin{aligned} r_{1,0}= & {} -4.40967C_F, \end{aligned}$$

(8)

$$\begin{aligned} r_{2,0}= & {} 47.0754C_{F}^2-\frac{9.23656C_F T_F}{e^{2}_{c}}\nonumber \\&-(22.7351)C_F n_H T_F\nonumber \\&+C_A C_F\left( -34.8488+123.208T_F\right) \nonumber \\&+C_F(-2.84012\times 10^{-15}C_A\nonumber \\&+37.8993C_F)\ln \frac{\mu _\Lambda }{m_c}, \end{aligned}$$

(9)

$$\begin{aligned} r_{2,1}= & {} -33.6021C_F T_F-4.40967C_F\ln \left[ \frac{\mu ^2_r}{m^2_c}\right] , \end{aligned}$$

(10)

where $n_H=1$, $C_{F}=\frac{N^{2}_{c}-1}{2N_{c}}$, $C_{A}=N_{c}$, $T_F=1/2$, 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,

$$\begin{aligned} R_c= & {} \Omega \left[ 1+r_{1,0}\alpha _{s}(Q_*)+r_{2,0}\alpha ^{2}_{s}(Q_*)\right] , \end{aligned}$$

(11)

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:

$$\begin{aligned} \ln \frac{Q^2_*}{m^2_c}= & {} -\frac{\hat{r}_{2,1}}{\hat{r}_{1,0}}+\mathcal{O}(\alpha _s), \end{aligned}$$

(12)

where $\hat{r}_{i,j}=r_{i,j}|_{\mu _{r}=m_c}$. It is noted that $Q_*$ is independent to any choice of renormalization scale $\mu _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 $\Lambda$, we need to choose a low-energy model for $\alpha _s$ so as to achieve a reliable prediction. In the literature, a variety of low-energy $\alpha _s$ models have been suggested [33,34,35,36,37,38,39,40,41]. A comparison of various low-energy $\alpha _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 $\alpha _s$ model is expressed as follows:

$$\begin{aligned} \alpha ^\mathrm{CON}_{s}(\mu )= & {} \frac{\pi }{\beta _{0}\ln {\frac{4M^{2}_{g}+\mu ^{2}}{\Lambda ^{2}}}}, \end{aligned}$$

(13)

where $M^{2}_{g}=m^{2}_{g}\left( \frac{\ln {(\mu ^{2}/\Lambda ^{2} +4m^{2}_{g}/\Lambda ^{2})}}{\ln (4m^{2}_{g}/\Lambda ^{2})}\right) ^{-12/11}$, and a lattice QCD calculation gives $m_{g}\in [650-700]~\mathrm MeV$ [43], which results in $\alpha ^\mathrm{CON}_s(Q_{*})\in [0.535-0.566]$. The asymptotic scale $\Lambda$ can be fixed by using the $\alpha _s$ measured at a typical energy scale. More definitely, by using $\alpha ^{\overline{\mathrm{MS}}}_{s}(M_{\tau })=0.325\pm 0.016$ [44], which leads to $\Lambda |_{n_{f}=3}=0.383^{+0.029}_{-0.031}$ GeV, $\Lambda |_{n_{f}=4}=0.324^{+0.029}_{-0.029}$ GeV, and $\Lambda |_{n_{f}=5}=0.223^{+0.022}_{-0.023}$ GeV.

Fig. 1

The strong coupling constant $\alpha _s$ within different energy scales. In low-energy region, the shaded band is for CON-model with $m_{g}\in [650$–$700]~\mathrm MeV$. In large-energy region, the conventional $\overline{\mathrm{MS}}$ $\alpha _s$ running behavior is adopted. $\Lambda |_{n_{f}=3}=0.383$ GeV

Figure 1 shows that the $\alpha _s$-running behavior at different scales, the $\alpha _s$ CON-model with $m_{g}\in [650$–$700]~\mathrm MeV$ is adopted in low-energy region which is shown by a 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. [45]. More explicitly, by requiring the first derivatives of $\alpha _s$ to be the same at the crossing point of the two energy regions, the $\alpha _s$ transition scale is determined to be 1.070–1.115 GeV.

Table 1 Contributions from each loop terms for $R_c$ up to NNLO level under conventional (Conv.) and PMC scale-setting approaches, respectively. The PMC predictions are scale invariant, and the Conv. predictions are highly scale dependent, whose central values are for $\mu _r=m_{c}$ and the errors are for $\mu _r\in [1\;\mathrm{GeV}, 2m_c]$. $\mu _{\Lambda }=1$ GeV

To do the numerical calculation, we take the c-quark pole mass $m_c=1.68$ GeV [46], and set the factorization $\mu _{\Lambda }=1$ GeV, the CON model parameter $m_{g}=700$ MeV, corresponding to $\alpha ^\mathrm{CON}_{s}(Q_{*})=0.535$. In Table 1, we give the contributions of each loop terms of $R_c$ under conventional and PMC scale-setting approaches, respectively. The conventional predictions are highly $\mu _r$-dependent for each loop terms and the total contributions of $R_c$, e.g. as shown by the following Eqs. (14–15), the renormalization scale uncertainty for $R^\mathrm{Conv.}_\mathrm{c, total}$ within the range of $\mu _r\in [1\mathrm{GeV}, 2m_c]$ is about $(^{-99\%}_{+57\%})$:

$$\begin{aligned} R^\mathrm{Conv.}_\mathrm{c, total}|_{\mu _{r}=1\;\mathrm{GeV}}= & {} 0.001~\mathrm{GeV}, \end{aligned}$$

(14)

$$\begin{aligned} R^\mathrm{Conv.}_\mathrm{c, total}|_{\mu _{r}=m_c}= & {} 0.067~\mathrm{GeV}, \end{aligned}$$

(15)

$$\begin{aligned} R^\mathrm{Conv.}_\mathrm{c, total}|_{\mu _{r}=2m_c}= & {} 0.105~\mathrm{GeV}. \end{aligned}$$

(16)

Those values deviate from the BESIII measurement [4] by at least $\sim 3.9\sigma$. Moreover, the separate scale uncertainties 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^\mathrm{PMC}_\mathrm{c, total}\equiv 0.246$ for any choice of $\mu _r$, which agrees with the BESIII measurement within errors. By further taking the $\alpha _s$ shift, $\Delta \alpha _s(M_\tau )=\pm 0.016$, into consideration, we obtain $\Delta R^\mathrm{Conv.}_\mathrm{c, total}|_{\mu _{r}=1\;\mathrm{GeV}} =\left( ^{+0.026}_{-0.029}\right)$, $\Delta R^\mathrm{Conv.}_\mathrm{c, total}|_{\mu _{r}=m_c}=\pm 0.013$, $\Delta R^\mathrm{Conv.}_\mathrm{c, total}|_{\mu _{r}=\mathrm{2m_c}}=\pm 0.007$, and $\Delta {R^\mathrm{PMC}_\mathrm{c, total}}=\pm 0.013$.

Fig. 2

The NLO order and NNLO order $R_c$-ratios under conventional and PMC scale-setting approaches as a function of $\mu _r$. $\mu _\Lambda =1$ GeV

Figure 2 shows explicitly how the $R_{c}$-ratio changes with different choices of $\mu _r$. It shows that after applying the PMC, the perturbative series is independent to any choice of $\mu _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 $\mu _{r}\in [1\mathrm{GeV}, 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 cancellation of conformal and non-conformal terms, which is however highly scale dependent, leading to scale dependent series.

Fig. 3

The NNLO $R_c$-ratio versus the factorization scale $\mu _\Lambda$ under conventional and PMC scale-setting approaches. Three typical $\mu _r$ are adopted, and the PMC prediction is independent to those choices of $\mu _r$. The shaded band represents the BESIII data [4] with the dashed line being its central value

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 $\chi _{c0}$ and $\chi _{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 $\mu _\Lambda$-dependent, the R-ratio shows explicit $\mu _\Lambda$ 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 anomalous dimensions for $\chi _{c0}$ and $\chi _{c2}$ are different [47], those two matrix elements are different at different $\mu _\Lambda$, which could compensate the $\mu _\Lambda$-dependence from the hard part. Using the evolution equation [9], we obtain

$$\begin{aligned} \ln {\frac{\langle \mathcal{{O}}_{1}(^{3}P_{0})\rangle |_{\mu _{\Lambda }}}{\langle \mathcal{{O}}_{1}(^{3}P_{0})\rangle |_{\mu _{\Lambda 0}}}}= & {} \frac{C_{F}(4C_{F}+C_{A})}{6} \alpha ^{2}_{s}\ln {\frac{\mu ^{2}_{\Lambda }}{\mu ^{2}_{\Lambda 0}}}, \end{aligned}$$

(17)

$$\begin{aligned} \ln {\frac{\langle \mathcal{{O}}_{1}(^{3}P_{2})\rangle |_{\mu _{\Lambda }}}{\langle \mathcal{{O}}_{1}(^{3}P_{2})\rangle |_{\mu _{\Lambda 0}}}}= & {} \frac{C_{F}(13C_{F}+C_{A})}{60} \alpha ^{2}_{s}\ln {\frac{\mu ^{2}_{\Lambda }}{\mu ^{2}_{\Lambda 0}}}. \end{aligned}$$

(18)

As required, the difference is an $\mathcal{O}(\alpha _s^2)$-order effect. If taking $\mathcal{{O}}_{1}(^{3}P_{J})|_{\mu _{\Lambda 0}=1\mathrm GeV}=0.107~\mathrm GeV^{5}$ [48] as the initial value, we obtain $\mathcal{{O}}_{1}(^{3}P_{0})|_{\mu _{\Lambda }=\mathrm 3 GeV}=0.139~\mathrm GeV^{5}$ and $\mathcal{{O}}_{1}(^{3}P_{2})|_{\mu _{\Lambda }=\mathrm 3 GeV}=0.124~\mathrm 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 $\mu _\Lambda$ indicates that the $R_c$-ratio is independent to any choice of $\mu _\Lambda$. Figure 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.

Fig. 4

The NNLO $R_c$-ratios with or without color-octet (CO) contributions under conventional and PMC scale-setting approaches. The error bars are for $\Delta \alpha _{s}(M_\tau )=\pm 0.016$ and $\mu _r\in [1\mathrm{GeV}, 2m_c]$. The dashed line is the central value of BESIII and the shaded band represents its error [4]

According to Refs. [13, 49, 50], 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 $\Delta \Gamma ^\mathrm{CO}_{\chi _{c0}}\simeq -0.3$ GeV and $\Delta \Gamma ^\mathrm{CO}_{\chi _{c2}} \simeq -0.227$ GeV [13]. If taking those extra color-octet contributions into consideration, we obtain $R^\mathrm{PMC}_\mathrm{c, total}=0.246$ to 0.299. Figure 4 shows the $R_c$-ratio with or without the color-octet contributions, where the error bars are for $\Delta \alpha _{s}(M_\tau )=\pm 0.016$ and $\mu _r\in [1\mathrm{GeV}, 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.

Table 2 Contributions from each loop terms for $R_b$ up to NNLO level under conventional (Conv.) and PMC scale-setting approaches, respectively. The PMC predictions are scale invariant, and the Conv. predictions are highly scale dependent, whose center values are for $\mu _r=m_{b}$, $\mu _{\Lambda }=m_{b}$. and the errors are for $\mu _r\in [1\;\mathrm{GeV}, 2m_b]$

Fig. 5

The NNLO $R_b$-ratios versus the factorization scale $\mu _\Lambda$ under conventional and PMC scale-setting approaches. Three typical $\mu _r$ are adopted, and the PMC prediction is independent to those choices of $\mu _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 approaches in Table 2 and Fig. 5. To do the numerical calculation, we take the b-quark pole mass $m_b=4.78$ GeV [46]. Due to $\alpha _s(m_b)\sim 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 $\mu _r\in [1\;\mathrm{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 $\mu _r\in [1\;\mathrm{GeV}, 2m_c]$. By applying the PMC, we obtain a more accurate pQCD prediction without renormalization scale uncertainty, $R_c|_\mathrm{PMC}=0.246\pm 0.013$, whose error is caused by $\Delta \alpha _s(M_\tau )=\pm 0.016$. This prediction agrees with the latest BESIII data within errors.