Revisiting the heavy vector quarkonium leptonic widths

We revisit the heavy quarkonium leptonic decays and using the Bethe-Salpeter method. The emphasis is on the relativistic corrections. For the decays, the relativistic effects are %, %, %, % and %, respectively. For the decays, the relativistic effects are %, %, %, % and %, respectively. Thus, the relativistic corrections are large and important in heavy quarkonium leptonic decays, especially for the highly excited charmonium. Our results for are consistent with the experimental data.


Introduction
R(D ( * ) ) R(K ( * ) ) As it gives a clean experimental signal, the dilepton annihilation decay of the heavy vector quarkonium plays an important role in determining the fundamental parameters such as the strong coupling constant [1,2], heavy quark masses [1,[3][4][5], heavy quarkonium decay constants [2,[6][7][8], etc. Its decay amplitude is a function of the quarkonium wave function, and this process can be used to test various theories such as the quark potential model, non-relativistic Quantum Chromodynamics (NR-QCD), QCD sums rules, lattice QCD, etc. The Standard Model prediction of the universality of lepton flavor is questioned by the measured ratios and [9][10][11][12][13][14], and the quarkonium leptonic decay is another way to test the lepton flavor universality.

Υ ′
The vector quarkonium leptonic decays have been studied since a long time [15][16][17][18][19][20][21]. With the progress in computer science and experimental technology, many advances have been reported in literature. For example, one can find the lattice QCD predictions of the leptonic decays of the ground-state Υ and its first radial excitation in [22]; Ref. [23] reported the next-to-leading non-perturbative prediction and Ref. [24] the next-to-leading-log perturbative QCD (pQCD) prediction; in Ref. [25], the Υ(1S ) two-loop QCD correction was computed; Ref. [26] studied the inclusive leptonic decay of Υ up to the next-tonext-to-leading order (NNLO) by including the re-summation of the logarithms (partly) up to the next-to-nextto-leading logarithmic (NNLL) accuracy; the NNNLO corrections have been discussed by various groups [27][28][29][30][31]. A pQCD analysis of the leptonic decay up to NNNLO using the principle of maximum conformality (PMC) [32][33][34][35] was presented in Refs. [36,37], where the renormalization scale ambiguity of the decay width is eliminated with the help of the renormalization group equation.
Even though considerable improvements have been made, there are still deviations between the theoretical predictions and the experimental data for the heavy vector quarkonium leptonic decays. There are two sources which may cause such deviations. The first are the unknown higher order perturbative QCD corrections. By using PMC, the conventional pQCD convergence of the series can be greatly improved by elimination of the divergent renormalon terms, and a more accurate decay width can be obtained. However, there are still large errors due to unknown high-order terms [36,37]. The second source is the relativistic correction, which could be large. However, almost all pQCD predictions are calculated using NRQCD, in which the decay constant of quarkonium, or its wave function at the origin, is treated simply in the non-relativistic approximation. One may argue that the relativistic correction is small for a heavy quarkonium, since the relative velocity among the heavy constituent quarks is small, e.g. [38]. However, many analyses in literature have found that the relativistic effect could be large. For example, Bodwin et al. computed the coefficients of the decay operators for the heavy quarkonium decay into a leptonic pair and found large relativistic correction [39]; Gonzalez et al. pointed out that large relativistic and QCD corrections of the quarkonium leptonic decays are necessary to fit the experimental data [7]; Geng et al. studied the meson semileptonic decays into charmonium and also found that the relativistic corrections are large [40], especially for highly excited charmonium states. Moreover, from the experimental standpoint, Ref. [41] showed that a careful study of leptonic decays is still needed for highly excited charmonium states.
In this paper, we focus on the leptonic decays of charmonium and bottomonium, including their excited states, using the relativistic method. In a previous short letter [42], we presented a relativistic calculation of the quarkonium decays into , where the results disagree with the experimental data. As a step forward, we revisit this topic in more detail, and include the decays into and as well as the ratios . We present the relativistic effects in these quarkonium decays, and discuss the universality of lepton flavor.
The paper is organized as follows. The general equation of the quarkonium leptonic decay width is given in Sec. 2. In Sec. 3, we give a brief review of the Bethe-Salpeter equation, and its instantaneous version, the Salpeter equation. We then show in Sec. 4 in detail how to solve the full Salpeter equation and obtain the relativistic wave function for a vector meson. The calculation of the decay constant in the relativistic method is given in Sec. 5. Finally, in Sec. 6, we give the numerical results and a discussion. A summary is presented in Sec. 7.

The quarkonium leptonic decay width nS
The leptonic partial decay rate of a vector charmonium or bottomonium state V is given by α em e Q e Q = +2/3 is the fine structure constant, is the electric charge of the heavy quark Q in units of the electron charge, for the charm quark and for the bottom quark, is the mass of the state quarkonium, is the lepton mass, is the decay constant of the vector meson that is defined by the following matrix element of the electromagnetic current where P is the quarkonium momentum, and ϵ is the polarization vector.
In the non-relativistic method, the well-known formula for the decay constant is where means the non-relativistic (NR), and is the non-relativistic wave function evaluated at the origin. In the NR method, there is only one radial wave function, and the vector meson and its corresponding pseudoscalar have the same radial wave function and the same decay constant. However, in the relativistic method, they have different wave functions and different decay constants, and more than one radial wave function gives a contribution to the vector meson decay constant.
In the relativistic method, the decay constant is not related to the wave function at the origin, but in the full region. In the following, we focus on the calculation of in the relativistic method.

The Bethe-Salpeter equation and the Salpeter equation
In this section, we briefly review the Bethe-Salpeter (BS) equation [43], which is a relativistic dynamic equation describing the two-body bound state, and its instantaneous version, the Salpeter equation [44]. The BS equation for a meson, which is a bound state of a quark, labelled as 1, and anti-quark, labelled as 2, can be written as [43] where is the relativistic wave function of the meson, is the interaction kernel between the quark and anti-quark, are the momenta and masses of the quark and anti-quark, P is the momentum of the meson, q is the relative momentum between quark and anti-quark. The momenta and satisfy the relations, and , where and . In the case of quarkonium, where , we have .
In the general case, the BS equation is hard to solve due to the complex interaction kernel between the constituent quarks. For the doubly heavy quarkonium considered here, the interaction kernel between the two heavy constituent quarks can be treated as instantaneous, leading to a simpler version of the BS equation. In this case, it is convenient to divide the relative momentum q into two parts, , where and Chinese Physics C Vol. 44, No. 6 (2020) 063104 , M is the mass of the bound state, and we have . Then, we have two Lorentz invariant variables, and . When , that is in the meson center-of-mass frame, they reduce to the usual components and , and .
With this notation, the volume element of the relativistic momentum k can be written in an invariant form , where and is the azimuthal angle. Taking the instantaneous approximation in the center-of-mass frame of the bound state, the kernel changes to . We introduce the three-dimensional wave function and the notation The BS eqution Eq. (4) is then rewritten as where and are propagators of quark 1 and anti-quark 2, respectively, which can be decomposed as Here, we have defined the constituent quark energy and the projection operators , where and 2 for quark and anti-quark, respectively.
Using the projection operators, we can divide the wave function into four parts with the definition . Here, and are called the positive and negative energy wave functions of the quarkonium. q P After integrating over on both sides of Eq. (7) using contour integration, we obtain the famous Salpeter equation [44]: .
Equivalently, the Salpeter equation can be written as four independent equations using the projection operators: The normalization condition for the BS wave function reads Note that usually in literature, it is not the full Salpeter equation Eq. (10) that is solved (or equivalently, the four Eqs. (11)- (14)), but only Eq. (11), which involves only the positive wave function. There is a good reason why such an approximation is made: it is its effective range. The numerical value of in Eq. (11) is much smaller than of in Eq. (12), which means that the positive wave function is dominant, and that the contribution of the negative wave function can be safely neglected. However, we point out that if only Eq. (11) for is considered, then not only is the contribution of the negative wave function neglected, but so are the relativistic effects of these wave functions. The reason is that the number of eigenvalue equations limits the number of radial wave functions, and as is shown below, only the four coupled equations Eqs. (11)- (14) can provide sufficient information to derive a relativistic wave function.

Relativistic wave function and the kernel
Although BS or the Salpeter equation is the relativistic dynamic equation describing the two-body bound state, the equation cannot by itself provide the information about the wave function. This means that we need to provide an explicit form of the relativistic kinematic wave function as input, which can be constructed using all allowable Lorentz and γ structures.

−
From literature, we have the familiar form of the nonrelativistic wave function for the vector meson, e.g.
where M, P and ϵ are the mass, momentum and polarization of the vector meson, is the relative momentum between the quark and anti-quark. There is only one unknown wave function in Eq. (16), which can be obtained numerically by solving Eq. (11) or the non-relativistic Schrödinger equation. The relative momentum is related to the relative velocity between the quark and anti-quark in the meson, . A relativistic wave function should depend on the relative velocity or momentum separately, not merely on the radial part , because the radial part is in fact or equally .

J pc J p J pc
To obtain the form of the relativistic wave function, we start from of a meson, because or are in any case good quantum numbers, where J is the total angular momentum, and p and c are the parity and the charge conjugate parity of the meson. The parity transform Chinese Physics C Vol. 44, No. 6 (2020) 063104 changes the momentum into , so for a meson, after applying the parity transform, the four-dimensional wave function changes to , where p is the eigenvalue of parity. The charge conjugate transform changes to , where c is the eigenvalue of charge conjugate parity, is the charge conjugate transform operator, and T is the transpose transform. Since the Salpeter equation is instantaneous, the input wave function is also instantaneous, and the general form of the wave function for the vector meson can be written as [45, 46] There are in total 8 radial wave functions with , which obviously cannot be obtained by solving only one equation, e.g. Eq. (11), but can be obtained by solving the full Salpeter Eqs. (11)- (14). The above expression does not include the terms with , since the condition of instantaneous interaction is . There are also no higher order terms like , , , etc., because the even powers of can be absorbed into the radial part of , while the odd powers of can be changed to lower power, for example, . By the way, if we delete all terms except those inside the radial wave functions, then the wave function Eq. (17) reduces to . If we further set , the wave function reduces to the non-relativistic case, e.g. Eq. (16). Thus, the terms with , , , , and in Eq. (17) are all relativistic corrections.
When the charge conjugate parity is taken into account, the terms with and vanish because of the positive charge conjugate parity , and the general instantaneous wave function for the quarkonium becomes Before moving on, we would like to discuss the interaction kernel . We know from Quantum Chromodynamics that the strong interaction between a quark and antiquark is given by the exchange of gluon(s), and that the basic kernel contains a short-range vector interaction plus a long-range linear confining λr scalar interaction suggested by the lattice QCD calculations [47]. In the Coulomb gauge and in the leading order, the kernel is the famous Cornell potential where λ is the string tension, is a free constant appearing in the potential to fit the data, and is the running coupling constant. In order to avoid infrared divergence and incorporate the screening effects, an exponential factor is added to the potential [48], i.e.
αr ≪ 1 It is easy to check that when , Eq. (20) reduces to Eq. (19). In the momentum space and in the rest frame of the bound state, the potential takes the form: where Here, is the running coupling of the one loop QCD correction, and . The constants λ, α, and are the parameters which characterize the potential, and for the system, for the system. The reader may wonder why we have chosen a simple basic kernel, and not a relativistic one [47,49] which includes details of the spin-independent potential and the spin-dependent potential, like the spin-spin interaction, spin-orbital interaction, tensor interaction, etc. The reason is that in our relativistic method, with a relativistic wave function for the bound state, we only need the basic potential and not a relativistic one, otherwise we would have double counting. To explain this, let us show how the relativistic potential is obtained: the potential between a quark and anti-quark is constructed from the on-shell scattering amplitude in the center-of-mass frame motivated by single gluon exchange, where the gluon propagator is given in the Coulomb gauge. The basic non-relativistic vector potential is obtained at leading-order from the amplitude (usually in the momentum space). To obtain the relativistic corrections of the potential, the on-shell Dirac spinors of the quark and anti-quark are expanded in quantities like the mass, momentum, etc. The relativistic potential is then obtained, and the relativistic corrections from the free spinors (wave functions for a bound state) are moved to the potential. The corresponding wave function becomes non-Chinese Physics C Vol. 44, No. 6 (2020) 063104 063104-4 relativistic.
In our case, we have a relativistic wave function and the potential is non-relativistic. If both of them are relativistic, then there is double counting. In general, a relativistic method should have a relativistic wave function with a non-relativistic potential, or a non-relativistic wave function with a relativistic potential. In principle, a halfrelativistic wave function with a half-relativistic potential is also permitted, but one has to be careful to avoid double counting. The method with a non-relativistic wave function and a relativistic potential is usually good for calculating the mass spectrum of bound states, while the method with a relativistic wave function and a non-relativistic potential is not only good for calculating the mass spectrum as an eigenvalue problem, but is also good for calculating the transition amplitude.
With the kernel Eq. (21) and the relativistic wave function Eq. (17) or Eq. (18), we are ready to solve the coupled Salpeter equation Eqs. (11)- (14). Substituting the wave function Eq. (18) into Eq. (13) and Eq. (14), taking the trace on both sides, multiplying with the polarization vector on both sides, e.g. or , and then using the completeness of the polarization vector, we obtain the relations where we have used for a quarkonium state. We now have only four independent unknown radial wave functions, , , , , whose numerical values can be obtained by solving Eq. (11) and Eq. (12). Substituting the wave function Eq. (18) into Eq. (11) and Eq. (12) and taking the trace again, we finally obtain four coupled equations where we have used the relation for a quarkonium, and , . Since we have four coupled equations, the four independent radial wave functions can be obtained numerically, and the mass spectrum obtained simultaneously as the eigenvalue problem.

The decay constant in the Salpeter method
The relativistic decay constant in Eq. (2) for a vector quarkonium can be calculated in the BS method as where is the color number, and Tr is the trace operator. We note that when calculating the decay constant, the Salpeter wave function , and not merely the positive wave function gives a contribution. For a vector quarkonium with the relativistic wave function Eq. (18), we obtain the relativistic decay constant The input parameters can be fixed by fitting the mass spectra of charmonium and bottomonium. We choose , , GeV, and GeV [42] 1) . We also choose and GeV for the charmonium system, and and GeV for the bottomonium system.
The mass spectra of vector charmonium and bottomonium are shown in Table 1. The theoretical predictions are consistent with the experimental data given by the Particle Data Group (PDG). As an example of the wave functions, we present four radial wave functions in Fig. 1: the dominant radial wave functions and and the two minor ones and 2) . From now on, we use the symbols and for simplicity.
As described in Sec. 4, the terms with radial wave functions and in the total wave function Eq. (18), are non-relativistic, while all the others are relativistic corrections. Figure 1 shows that the relativistic wave functions and are small and could be safely neglected, but in fact this is not the case. Figure 1 only shows the relative importance of the wave functions in the region of small q, and to see the relative importance of the wave functions in the whole q region, we plot the ratio in Fig. 2. It can be seen that in the large q region, the value of is only a few times larger than of . Thus, the terms which are proportional to ( and others) may have a sizable contribution in the large q region, leading to possibly important relativistic corrections.

Charmonium leptonic decay widths
Our results for are shown in Table 2, where in the second column, 'NR ', the non-relativistic decay rates are shown, meaning that in Eq. (28) the term is ignored, so that the only contribution is from the term. The third column, 'Re ' , show the relativistic results including the contributions of and . One can see that for charmonium the relativistic results are different from the non-relativistic ones. To see this clearly, we add the fourth column in Table 2 with the ratio (NR-Re)/Re, whose value can be called the 'relativistic effect'. Table 2 indicates that the relativistic effect is about 22% for the decay, which is consistent with the usual power relation for the relativistic terms, e.g.
. For the excited states, the relativistic effects are much larger than for the ground state. For the 2S, 3S, 4S and 5S states, the relativistic effects are about 34%, 41%, 52% and 62%, respectively. These results are consistent with our previous study of the semi-leptonic decays , where higher excited charmonium states were shown to have larger relativistic effects 1) In this previous Letter, we have chosen different for charmonium and bottomonium. Since this parameter appears only in , which depends on , and for more convenience of fitting the data, we choose the same for the two systems. 2) Here we show the curves of and other than and , because they always appear in such a combined form in the applications.
Chinese Physics C Vol. 44, No. 6 (2020) 063104 [40]. This conclusion can also be obtained qualitatively from the plots of the radial wave functions. We mentioned that the relative momentum q concerns the relative velocity between the quark and antiquark in quarkonium, . As shown in Fig. 1, two nonrelativistic radial wave functions always dominate over the relativistic wave functions in the whole q region, leading to a small relativistic correction. For the excited states, see Fig. 3 as an example of the radial wave functions of , the non-relativistic wave functions still dominate in the small q region, but there is a node structure in each curve where the wave function changes sign. The contributions in the low q region may cancel each other, and the wave functions for large q ( ) may give sizable contributions, resulting in large relativistic corrections.
There are other methods for considering the relativistic effects in heavy quarkonium decays. For example, Bodwin et al. [39] and Brambilla et al. [51] and the corrections of the decay rate of quarkonium in the framework of NRQCD. In the case of [39], the predicted relativistic effect is % for , and is 23.0% for . These values are consistent with our prediction of 22.3%.
In Table 2, we also show the theoretical uncertainties caused by the choice of input parameters. We vary all parameters simultaneously within % of their central values, and take the largest variation as the uncertainty. With the errors, most predictions are much larger than the experimental data. The only exception is the channel , which has large uncertainties 1) . We note that a calculation of the leptonic decay in lattice QCD with fully relativistic charm quarks was reported in Ref.
[52], and gave 5.48(16) keV , consistent with the experimental data. This indicates that the disagreement of our results with the experimental data may be due to the lack of QCD corrections.
We note that in a recent paper, Soni et al. [53] calcu-   1) The reason is that the mass is only a slightly heavier than that of two , so the phase space of this channel is very sensitive to the variation of parameters.
Chinese Physics C Vol. 44, No. 6 (2020) 063104 lated the quarkonium leptonic decay using the Cornell potential in a non-relativistic version and with the pQCD correction up to NLO. Their results for charmonium are neither consistent with the experimental data nor with our results, while for bottomonium their results are comparable with ours (see below). Also, Badalian et al. [54] calculated the decay rates of the leptinic decays with the QCD correction at NLO using the Cornell potential and the semi-Salpeter equation, and obtained 5.47 keV, 2.68 keV, 1.97 keV, and 1.58 keV , respectively, which are smaller than our charmonium results. These studies indicate that the relativistic corrections and QCD corrections are large for the charmonium system.

Bottomonium leptonic decay widths
We present the non-relativistic and relativistic results of the bottomonium leptonic decay widths in Table 3. Similarly to charmonium, the relativistic corrections are sizable. For the ground state , the relativistic effect is about 14%, and for the excited states , they vary from 20% to 28%. These predictions agree with those in literature. For example, Bodwin et al. [39] predicted the relativistic effect of 13.2% for using NRQCD up to the accuracy, and a lattice QCD prediction indicated that the relativistic effects are about 15%-25% [22] for and up to the accuracy.
We should point out that the above large relativistic effects are specific for bottomonium leptonic decays , and are not universal for processes involving a bottomonium. In the di-lepton decays, the amplitude is proportional to the wave function as , i.e. the wave function is to the power of one. For other processes, such as meson A to meson B semileptonic decays, the amplitude is proportional to the overlapping integral of the wave functions for the initial and final states . Because the wave functions are large in the small q region, the product of two wave functions is suppressed in the large q region compared to the case with one wave function, and the contributions from the relativistic terms are greatly suppressed. In Table 3, we also give the theoretical uncertainties, which are obtained by varying all parameters simultaneously within % of the central values, and the largest variations are taken as the errors. Our relativistic results agree well with the experimental data. We also note that for , PDG gives two different results: directly listed is keV, but a branching ratio % is also given, leading to keV using the full width keV [50]. The second value is the same as our relativistic result. Similarly, in the case of , PDG directly lists keV [50], but from the branching ratio also given in PDG, we get keV. We hope PDG will update the data in the near future. Table 3 shows that all relativistic results for are consistent with the experimental data. Our predictions also agree with the lattice QCD predic-  tion [22], keV and keV, with the NRQCD prediction [26], keV, and with the NR-QCD prediction with NNNLO pQCD corrections [30], keV.

Lepton flavor university
To test the lepton flavor university, we give the ratios and in Table 4. Their definitions are similar, for example,   Table 4 shows the ratios calculated with the 'Re' values. The uncertainties of the ratios are from the variation of the input parameters. In the case of charmonium, the ratios are quite different from each other, since the charmonium mass is a bit higher than of the two τ . For the same reason, we get a large uncertainty. For bottomonium, since the τ mass is much smaller than the bottomonium mass, we get almost the same values for all ratios Their uncertainty is also very small due to the cancellation between the numerator and denominator.
Even though all central values of the ratios are smaller than 1, they are consistent with the existing experimental data within errors.
To cancel the model dependence of the theoretical predictions, we give in Table 5 and Table 6 the ratios and . For the decay, we obtain the same central values for the e and τ final states, so we only present the ratio in Table 5 and Table 6, which are calculated using the final states listed in Table 3. Table 5 shows that the ratio is larger but close to the experimental data, while the ratios for highly excited states are much larger than the experimental data. Table 6 shows the bottomonium leptonic decay ratios. In the row 'Exp1', the value of keV is used, which is directly listed in PDG. In the row 'Exp2', keV is used, which is calculated using the branching ratio of given in PDG. For , the results outside the brackets were obtained using keV from PDG, while the results inside the brackets used keV obtained from the PDG branching ratio. It can be seen that all our theoretical predictions are consistent with the experimental data.     In this paper, we studied the leptonic decays of heavy vector quarkonia. For the charmonium decays, not all states are consistent with the experimental data, while for the bottomonium decays, almost all S wave states are in good agreement with the data.
Theoretical results of the ratios and were given in Ref. [55], where the potential model was used including the relativistic corrections and pQCD corrections at NLO. These results are comparable with ours, i.e. the charmonium leptonic decay widths are not consistent with the experimental data and the bottomonium leptonic widths are in good agreement with the data. This situation was also observed in Ref. [56]. It seems that the same theoretical tool cannot provide satisfactory results for both the charmonium and bottomonium systems [57]. There are several possible reasons for this difference in our study. It may be that the instantaneous approximation works well for bottomonium, but is not good enough for charmonium. An improvement of the Cornell potential may be needed, and more importantly, the perturbative QCD corrections may have larger effect in charmonium decays than in bottomonium decays. Since the BS equation is an integral equation, the QCD corrections from the gluon ladder diagrams are already included, but other QCD corrections in the kernel or in the quark propagators may need to be improved in future calculations.