Revisit the Heavy Vector Quarkonium Leptonic Widths

In this paper, we revisit the heavy quarkonium leptonic decays $\psi(nS) \to \ell^+\ell^-$ and $\Upsilon(nS) \to \ell^+\ell^-$ using the Bethe-Salpeter method with QCD correction. The emphasis is paid to the relativistic correction. For $\psi(1S-5S)$, we obtained the relativistic effects are $22\%$, $34\%$, $41\%$, $52\%$ and $62\%$, respectively. For $\Upsilon(1S-5S)$, the values are $14\%$, $23\%$, $20\%$, $21\%$ and $28\%$. We conclude that the relativistic corrections of the heavy quarkonium leptonic decays are large and important even for bottomonium, especially for higher excited states. Our results of $\Upsilon(nS) \to \ell^+\ell^-$ consist with experimental data.


I. INTRODUCTION
Heavy vector quarkonium dilepton annihilation decay as a clean experimental signal process played an important role in determining fundamental parameters such as the strong coupling constant [1,2], heavy quark masses [1, [3][4][5], decay constants [2,[6][7][8], etc.Since the corresponding decay amplitude is a function of the quarkonium wave function, this process can be used to test theories such as the quark potential model, the non-relativistic Quantum Chromodynamics (NRQCD), the QCD sums rules, the Lattice QCD, and etc.On the other hand, the Standard Model prediction of the university of lepton flavor is questioned by the ratios of R(D ( * ) ) and R(K ( * ) ) etc, and the processes of leptonic decays of quarkonium provide another place to test the lepton flavor universality [9][10][11][12][13][14].
The vector quarkonium leptonic decays have been studied long time ago [15][16][17][18].Recently, along with the progresses both on the computer science and experimental technology, these topics have been studied extensively.For examples, there are Lattice QCD predictions [19] on the leptonic decays of the ground-state Υ and its first radial excitation Υ ′ ; Ref. [20] gave the next-to-leading nonperturbative calculation and Ref. [21] the next-to-leading-log perturbative QCD calculation; Ref. [22] computes the two-loop QCD correction to the quarkonium leptonic decays; Ref. [23] studied the inclusive leptonic decay of bottomonium ground state at next-to-next-to-next-to-leading order (NNNLO) and includes the resummation of logarithms at next-to-next-to-leading logarithmic (NNLL) accuracy; the NNNLO corrections [24] and logarithm NNNLO corrections [25,26] are also studied by different groups; The latest progresses are made by Ref. [27], in which the third order corrections for bottomonium have been presented; Ref. [28] also gave the three loop corrections, Ref. [29] applied the principle of maximum conformality (PMC) scale-setting approach [30][31][32][33] to set the optimal renormalization scale of the decay width and obtained more convergence perturbative QCD series, Ref. [19] provided a nonperturbative calculation using full lattice QCD method with perturbative corrections.
We can see that large improvements have been made on the heavy vector quarkonium leptonic decays.However, there are still deviations between theoretical predictions and experimental data.There are two sources which may cause such deviations: 1) The first one is the unknown higher order perturbative QCD corrections, since the perturbative QCD convergence is questionable under conventional scale setting approach; By using PMC, the pQCD convergence is greatly improved and a more reasonable decay width can be achieved, which however still has large error due to unknown high-order terms [29]; 2) The second one is the relativistic correction, which could be large.However almost all the existing pQCD predictions are calculated by using the NRQCD, in which the decay constant of quarkonium, or its wave function at origin is treated simply within the non-relativistic approximation.
One may argue that the relativistic corrections are small for heavy quarkonium.While Ref. [34] computed the coefficients of the decay operators for the decay of 3 S 1 heavy quarkonium into a leptonic pair, found large relative velocity correction in order v 2 ∼ 0.1 − 0.3, means large relativistic corrections.Ref. [35] studied the B c meson semileptonic decays to charmonium and found that the relativistic corrections are also large, especially for highly excited charmonium states.So it is reasonably to assume that the relativistic corrections for highly excited bottomonium states are also large.From the experimental standpoint, Ref. [36] showed that careful study of leptonic decays is still needed for excited charmonium states.
So in this paper, we focus on the leptonic decays of charmonium and bottomonium, including excited states, using a relativistic method with perturbative QCD corrections.
Because the nonperturbative effects are rather moderate, and the perturbative contributions are dominated compared with nonperturbative corrections [37].In a previous short Letter [38], we have provided a relativistic calculation of the quarkonium leptonic decays to e + e − .In this work, we did not include the QCD corrections, and the results disagree with the experimental data.In this paper, we shall revisit this topic, supplement the decays to µ + µ − and τ + τ − as well as applying the existing high order perturbative QCD corrections and relativistic corrections.Then, we shall discuss the universality of lepton flavor.This remaining parts of the paper are organized as follows.Section II gives the general equation of quarkonium leptonic decay width including QCD correction.In section III, we give a brief review of the Bethe-Salpeter equation, and its instantaneous version, Salpeter equation.We then show the details of how to solve the full Salpeter equation and obtain the relativistic wave functions for vector meson in section IV.Section V shows our relativistic method how to calculate the decay constant with QCD correction.Finally, in section VI, we give the numerical results and some discussions.

II. DECAY WIDTH
The leptonic partial decay rate of a charmonium or a bottomonium nS vector state V is given by Γ where α em is the fine structure constant, e Q is the electric charge of the heavy quark Q in units of the electron charge, e Q = 2/3 for charm quark, e Q = −1/3 for bottom quark, M nS is the mass of quarkonium nS state, m ℓ is the mass of lepton, F V is the decay constant of a vector meson defined through the following matrix element of the electromagnetic current where P is quarkonium momentum, and ǫ the polarization vector.
In a non-relativistic limit and without QCD corrections, the well-known formula for the decay constant is where NR means the non-relativistic, Ψ V (0) is the non-relativistic meson wave function evaluated at the origin.In this non-relativistic limit, only one single radial wave function is existed, and a vector meson and its corresponding pseudoscalar have the same radial wave function and same decay constant.In a relativistic method, we can see that more than one radial wave functions have contribution to the vector decay constant.
In a relativistic method with QCD correction, the decay constant can be formulated as where F RE V is the vector decay constant calculated by a relativistic model but without QCD correction, while β V is the correction from QCD.In following part, we will focus on the calculating of the F RE V by a relativistic method, while for the perturbative QCD correction β V , we will choose the existing results in literature.

III. BETHE-SALPETER EQUATION AND SALPETER EQUATION
In this section, we briefly outline the Bethe-Salpeter (BS) equation [39], which is a relativistic dynamic equation describing two particle bound state, and its instantaneous version, Salpeter equation [40].The BS equation is provided by the Quantum Field Theory, for a meson which is a bound state of a quark labeled as 1 and anti-quark 2, the BS equation is read as [39] ( / p 1 − m 1 )χ P (q where χ P (q) is the relativistic wave function of the meson, V (P, k, q) is the interaction kernel between the quark and anti-quark, and p 1 , p 2 , m 1 , m 2 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, and we have the relations p 1 = α 1 P + q, p 2 = α 2 P − q, where α 1 = m 1 m 1 +m 2 and α 2 = m 2 m 1 +m 2 , so for quarkonium, α 1 = α 2 = 0.5.Since the BS equation is hard to be solved, and we only focus on double-heavy quarkonium, in this case the kernel between two heavy constituent quarks can be seen as an instantaneous one, so we will solve the instantaneous version of BS equation.Then it is convenience to divide the relative momentum q into two parts q µ = q µ + q µ ⊥ , where q µ ≡ (P • q/M 2 )P µ and q µ ⊥ ≡ q µ − q µ , M is the mass of the bound state, and we have P 2 = M 2 .Then we have two Lorentz invariant variables, q P = (P •q) M and q T = q 2 P − q 2 = −q 2 ⊥ , when P = 0, that is in the center-of-mass system of the meson, they turn to the usual component q 0 and | q|, and q ⊥ = (0, q).With these notations, the volume element of the relativistic momentum k can be written in an invariant form d 4 k = dk P k 2 T dsdφ, where ds = (k P q P − k • q)/(k T q T ) and φ is the azimuthal angle.When the kernel is instantaneous, in the center of mass system of the bound state, V (P, k, q) become to V (k ⊥ , q ⊥ , s), then we introduce the three dimensional wave function and the notation Then the BS Eq. ( 5) is rewritten as where S 1 (p 1 ) and S 2 (p 2 ) are the propagators of quark 1 and anti-quark 2, which can be decomposed as where we have defined the constituent quark energy ω i = m 2 i + q 2 T and the project operator M ω i ± (−1) i+1 (m i + / q ⊥ ) , i = 1, 2 for quark and anti-quark, respectively.
Using the project operators, we can divide the wave function into four parts with the definitions , where Ψ ++ P (q ⊥ ) and Ψ −− P (q ⊥ ) are called positive and negative energy wave functions, respectively.
After integrated over q P on both sides of Eq. ( 8) with contour integration, this means the instantaneous approach has been chosen, the BS Eq. ( 5) become to the famous Salpeter equation [40]: Or equally, the Salpeter equation can be written as four independent equations by using the project operators: Note that in literature, usually not the full Salpeter Eq. ( 11) is solved, or equally the four Eqs.(12)(13)(14)(15), but only the first positive wave function Eq. ( 12) is solved.There is a logical reason to make such approximation, because in its effective range, the numerical value of 12) is much smaller than that of the corresponding M + ω 1 + ω 2 in Eq. ( 13), thus result in dominant value of positive wave function Ψ ++ P (q ⊥ ), and the contribution of negative wave function Ψ −− P (q ⊥ ) is much smaller, which can be deleted safely.But we argue that, the contribution of Ψ −− P (q ⊥ ) can be deleted safely only after we already solved the four Eqs.(12)(13)(14)(15).Otherwise if we only consider and solve the Eq. ( 12) of positive wave function, we ignore not merely the contribution of negative wave function, but also the relativistic effects at the same time.The reason is simple, the number of eigenvalue equations will limit the number of the radial wave functions we obtained, and the coupled four Eqs.(12)(13)(14)(15) will provide us enough equations to give a relativistic wave function, we will show the details later.
The normalization condition for BS wave function is read as

IV. RELATIVISTIC WAVE FUNCTION AND THE KERNEL
Though the BS or the Salpeter equation is the relativistic dynamic equation describing two-particle bound state, but the equation itself cannot provide us the information of the wave function automatically, this means we have to provide a relativistic kinematic wave function form as input.
In literature, we are familiar with the non-relativistic wave function form for a 1 − vector where M, P and ǫ are the mass, momentum and polarization of the vector meson, respectively, q is the relative momentum between the quark and anti-quark.We note that there is only one unknown radial wave function ψ( q) which can be obtained numerically by solving the positive BS Eq. ( 12) or the non-relativistic Schrodinger equation.
The relative momentum q is related to the relative velocity v between quark and antiquark in the meson, the relation is q = m 1 m 2 m 1 +m 2 v, where m 1 and m 2 are the mass of quark 1 and anti-quark 2. So a relativistic wave function is must depend on the relative velocity v or momentum q, not only through the radial part ψ( q), because the radial part actually is ψ(| q|) or equally ψ( q 2 ).
To give the relativistic wave function form, we start from the J pc of a meson, because the J p or J pc is a good quantum number in any case, where J is total angular momentum, p and c are parity and charge conjugate parity of the meson, respectively.Parity transform change the momentum q = (q 0 , q) to q ′ = (q 0 , − q), so for a meson, after parity transform, its four-dimensional wave function χ P (q) become to p • γ 0 χ P ′ (q ′ )γ 0 .Charge conjugate transform change χ P (q) to c • Cχ T P (−q)C −1 , where C = γ 2 γ 0 is the charge conjugate transform operator, T is the transpose transform.Since the Salpeter equation is instantaneous, the input wave function Ψ P (q ⊥ ) is also instantaneous, in this case, the general wave function for a 1 − vector meson can be written as [41] Ψ Where we have 8 radial wave functions ψ i (q ⊥ ) = ψ i (| q|), i = 1 ∼ 8, their numerical values are unknown now, which obviously can not be obtained by solving only one equation, the positive Eq. ( 12), but will be obtained by solving the full Salpeter Eq. ( 11), actually the coupled Salpeter Eqs.(12)(13)(14)(15).One note that the general wave function do not include the terms with P • q, the reason is that, in the case of instantaneous, P • q = P • q ⊥ = 0.There are also no high order q ⊥ terms like q 2 ⊥ , q 3 ⊥ , q 4 ⊥ , etc, because the even power of q ⊥ can be absorbed into the radial part of ψ i (q ⊥ ), while the odd power q ⊥ terms can be changed to a lower power one, for example, q 3 ⊥ ψ ′ i (q ⊥ ) = q ⊥ ψ i (q ⊥ ).Compared with the non-relativistic wave function in Eq. ( 17), if we delete the terms with q except in the radial parts, then the wave function Eq. ( 18) become to M ǫ ⊥ ψ 5 (q ⊥ ) + ǫ ⊥ P ψ 6 (q ⊥ ), if we further set ψ 5 (q ⊥ ) = −ψ 6 (q ⊥ ) = ψ(q ⊥ ), then this wave function become to the non-relativistic one in Eq. (17).So the terms of ψ 1 , ψ 2 , ψ 3 , ψ 4 , ψ 7 and ψ 8 are all relativistic corrections in Eq. (18).Numerical values of radial wave functions will show the details and consist with this comment latter.
When we consider the charge conjugate parity, then the terms with ψ 2 (q ⊥ ) and ψ 7 (q ⊥ ) vanish because of positive charge conjugate parity c = +, so the general instantaneous wave function for a 1 −− quarkonium become to Before going on, we like to first give the interaction kernel between quark and anti-quark.
We know that from Quantum Chromodynamics, the strong interaction between quarks is given by exchange of gluon, so the basic kernel are a short-range γ µ ⊗ γ µ vector interaction − 4αs 3r and a long-range 1 ⊗ 1 linear confining scalar interaction λr suggested by Lattice QCD calculations [42].In Coulomb gauge and leading order, the kernel is the famous Cornell potential where λ is the string tension, V 0 is a free constant appearing in potential model to fit data, α s is the running coupling constant.In order to avoid infrared divergence and incorporate the screening effects in the potential [43], we add a factor e −αr , It is easy to show that when αr ≪ 1, Eq. ( 21) become to the Eq. ( 20).In the momentum space and the rest frame of the bound state, the potential is read as: where
Here α s ( q) is the running strong coupling constant with one loop QCD correction, and e = 2.71828.The constants λ, α, V 0 and Λ QCD are the parameters which characterize the potential, and N f = 3 for cc system, N f = 4 for b b system.
The readers may have a question why we choose a simple basic potential, not a relativistic one [42,44] including details of spin-independent potential and 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, not a relativistic one, otherwise we will meet the double counting problem.To explain this, let us show how to obtain the relativistic potential.The potential between a quark and an anti-quark is constructed by on-shell q q scattering amplitudes in the center-of-mass frame motivated by single gluon exchange, where the gluon propagator is in the Coulomb gauge, in this leading order, the basic non-relativistic vector potential − 4αs 3r can be obtained (usually in momentum space).To obtain the relativistic corrections of the potential, the on-shell Dirac spinors of the quark and anti-quark are expanded according to quantities like mass, momentum, etc, then the relativistic potential is obtained, at the same time, the relativistic effect from the free spinors (wave function for a bound state) is moved to the potential, then the corresponding wave function become non-relativistic one.
In our case, we have a relativistic wave function, then the potential is non-relativistic, otherwise both of them are relativistic, there is the double counting.So generally, a relativistic method, can be a relativistic wave function with a non-relativistic potential, or a non-relativistic wave function with a relativistic potential.In principle, a half-relativistic wave function with a half-relativistic potential is also permitted, but one has to be careful to avoid double counting.Usually, the method with a non-relativistic wave function but a relativistic potential, is good at calculating the mass spectrum of bound states, while the method with a relativistic wave function and a non-relativistic potential, not only good at calculating the mass spectrum as a eigenvalue problem, but also good at calculating transition amplitude.Now with the kernel Eq. ( 22) and relativistic wave function Eq. ( 18) or Eq. ( 19), we can solve the coupled Salpeter equations, Eqs.(12)(13)(14)(15).First, put Eq. ( 19) into Eq.( 14) and Eq. ( 15), take trace on both sides, then multiply a quantity with polarization vector on both sides, for example, q ⊥ • ǫ * or / ǫ * ⊥ • / P , etc, after using the completeness of the polarization vector, then we obtain the relations where we have use m 1 = m 2 for quarkonium state.Now we have only four independent unknown radial wave functions, ψ 3 (q ⊥ ), ψ 4 (q ⊥ ), ψ 5 (q ⊥ ), ψ 6 (q ⊥ ), whose numerical values can be obtained by solving the Eq. ( 12) and Eq. ( 13).Put wave function Eq. ( 19) into Eq.( 12) where we have ω 1 = ω 2 for quarkonium, and Since we have four coupled equations, then the numerical values of the four independent radial wave functions can be obtained, and as an eigenvalue problem, the mass spectra are obtained when solving the equations.Now the normalization condition Eq. ( 16) for the 1 −− wave functions is

V. THE DECAY CONSTANT IN SALPETER METHOD WITH QCD CORREC-TION
The relativistic decay constant F RE V of the Eq. ( 2) for a vector quarkonium in BS method can be calculated as where N c = 3 is the color number, T r is the trace operator.We note that, in calculating the decay constant, the Salpeter wave function Ψ P ( q), not only the positive wave function Ψ ++ P ( q) has contribution.For a vector quarkonium, with the wave function Ψ 1 −− P ( q), Eq. ( 19), we obtain the relativistic vector meson decay constant For the perturbative QCD correction, we choose the results in Ref. [22], where two-loop corrections have given both for charmonium and bottomonium.Then the decay constant for charmonium with two-loop QCD correction can be written as where n f = 3 for charmonium, and the value of α s is scale µ dependent.For the ground state J/Ψ, we will use µ = m c which is close to M J/Ψ /2, for excited state Ψ(nS), µ = M nS /2 is used.When we calculate the decay constant of bottomonium, the factor β V is read in Ref. [22] where n f = 4 for bottomonium is already used, and similarly, we choose for the ground bottomonium Υ, and for excited state Υ(nS), µ = M nS /2.

VI. RESULTS AND DISCUSSIONS
The parameters used in this paper: m b = 4.96 GeV, m c = 1.6 GeV, α = 0.06 GeV, Λ QCD = 0.21 GeV.In previous Letter [38], we chose different Λ QCD for cc and b b systems, since this parameter appeared only in the calculation of α s , which is the energy scale q (µ) dependent, so we choosing same Λ QCD for the two systems is better.At the same time, we choose different λ for the two systems.Then for cc system, λ = 0.23 GeV 2 and c = −0.249 GeV, for b b system, λ = 0.2 GeV 2 and c = −0.124GeV.
The mass spectra of vector charmonium and bottomonium systems are shown in Table I, the theoretical predictions are consist with experimental data.As an example, we show the figures of four different radial wave functions of J/ψ in Figure 1, where we plotted the dominant radial wave functions ψ 5 and ψ 6 , while for the minor ones, we plotted the q 2 ψ 3 /M 2 and q 2 ψ 4 /M 2 instead of the ψ 3 and ψ 4 , because they always appear in the former form in calculation.
As we described in Sec IV, wave function for the vector charmonium in Eq. ( 19), except the terms with radial wave functions ψ 5 and ψ 6 , others are all relativistic corrections, and their numerical values are much smaller than ψ 5 and ψ 6 , see Fig. 1.If we ignore all of them, and set ψ 6 = −ψ 5 , then the relativistic wave function Eq. ( 19) return to the non-relativistic wave function Eq. ( 17).From Fig. 1, it seems that the relativistic wave functions can be delete safely, but it is not true, because in Fig. 1, only in small | q| region, the wave function can be seen clearly.To see the details of wave function when | q| is large, in Figure 2, we plot the ratio of ψ 5 /( q 2 ψ 3 /M 2 ) for the J/ψ.We can see that, in | q| region, the value of ψ 5 is only several times larger than that of the ( q 2 ψ 3 /M 2 ), this means in a physical process, ψ 3 terms have sizable contributions in large q region, and in whole, we will obtain sizable relativistic contributions.
The results of ψ(nS) → ℓ + ℓ − are shown in Table II.Here in the first column 'NR',  the non-relativistic decay rates without QCD corrections are shown, that means when we calculate the decay constant in Eq. ( 29), only ψ 5 has contribution, and the contribution of ψ 3 is deleted.Second column 'Re' show the relativistic results including the contribution both of ψ 5 and ψ 3 in Eq. ( 29), but without QCD corrections.One can see that the relativistic results are deferent from the non-relativistic ones obviously in charmonium cases.To see this clearly, we add the third column '(NR-Re)/Re' in Table II to of the difference the non-relativistic result and the relativistic one over the later.
We can see that from the third column in Table II, for the J/ψ decay, the relativistic correction is about 22%, which is consist with our normal knowledge.While for excited states, the relativistic effects are much larger than the ground state case, for 2S, 3S, 4S and 5S states, the relativistic effects are about 34%, 41%, 52% and 62%, respectively, which are consist with our previous study on the semi-leptonic decays B + c → cc + ℓ + + ν ℓ , higher excited state has larger relativistic effect.This conclusion can also be obtained qualitatively from the diagrams of radial wave functions.First, we mention that relative momentum q relates to the relative velocity v between quark and antiquark in the meson, for charmonium and bottomonium, q = 0.5m Q v. Second, from Fig. 1, the J/ψ radial wave function, we can conclude that the dominant contribution come from the small q region, because in this region the numerical values of wave function are dominant large, and small contribution from the large q ( v) region means small relativistic correction.While for excited state, see figure 3 as an example, where we plot the radial wave functions for ψ 2S .We note that there is a node structure, before the node, the wave function is positive, after the node, the wave function become negative.The negative radial wave function will has negative contribution, and so large q ( v) region has sizable contribution, means larger relativistic contribution.
In Table II, the fourth and fifth columns show the relativistic result with QCD correction.
Where 'Re+α s ' means the relativistic result with the first order QCD correction, that is , and 'Re+α 2 s ' means the relativistic result with the two loop QCD correction, which include first order α s and second order α 2 s contributions, see Eq. ( 30).One can see that, with the first order QCD correction, our theoretical results become better comparing with experimental data, especially for J/ψ and ψ(2S).With the two loop QCD correction, higher excited states results become good, but even worse for the lower states.The reason is that we have a bad QCD convergence series in Eq. ( 30) because there is a large parameter (44.5 −0.14n f ) in α 2 s terms.For J/ψ, we choose µ = m c , then we obtain α s = 0.34 and the first order QCD correction factor β V (α s ) = 0.711, but the full two loop QCD factor β V (α 2 s ) = 0.195 result in a very small decay rate of J/ψ.When we consider highly excited state ψ(5S), we choose µ = M 5S /2, and obtained α s = 0.291, first order QCD factor β V (α s ) = 0.753, and two loop QCD factor β V (α 2 s ) = 0.374, which is much larger than J/ψ case.
In Table III, we show the results of bottomonium decays.One note that, the relativistic correction is sizable, for ground state Υ(1S), the relativistic correction is about 14%, for excited states Υ(2S − 5S), the relativistic corrections are from 20% to 28%, not very small, and can not be ignored.The unexpected large relativistic effect is special for the decay of Υ(nS) → ℓ + ℓ − , not a universal conclusion for a process when a bottomonium involved.
Because in this di-lepton decay, the amplitude is in proportion to the wave function, here 3M 2 ψ 3 ( q) , the power of the wave function is one, other process, for example, meson A semileptonic decay to meson B, the corresponding amplitude is in proportion to overlapping integral over two wave functions of initial and final states, d 3 q ϕ A • ϕ B , the power of wave function is two, so because of the behavior of wave functions (values are large in small q region), comparing with one wave function case, the product of two wave functions will suppress the large q region contribution, that is, suppress the contribution from relativistic correction.with the experimental data very well for all the bottomonium states.We also note that, for the Υ(1S) → e + e − , the PDG provide us two different data, it directly list Γ ee = 1.34 keV, but at the same time the PDG also list the branching ratio of Br = 2.38% [45], and the full width of Γ Υ(1S) = 54.02keV, then we get Γ ee = 1.29 keV for Υ(1S) → e + e − , which is exactly the same as our relativistic results without QCD correction.Same thing happened to the case of Υ(4S) → e + e − , the PDG list directly the Γ ee = 0.272 keV [45], but from the branching ratio given in PDG, we get Γ ee = 0.322 keV, we hope PDG update these data in near future.

One can see that in
The results of bottomonium decays including QCD corrections are also shown in the fourth and fifth column in Table III.In these cases, only the leading order QCD correction has sizable contribution, and the α 2 s term can be deleted.For Υ(1S), we choose µ = m b , then α s = 0.238, β V (α s ) = 0.798, while β V (α 2 s ) = 0.788.For Υ(5S), µ = M 5S /2 is chosen, then we has α s = 0.232, β V (α s ) = 0.803 and β V (α 2 s ) = 0.794, with these values, we got For the experimental data of 'PDG' [45], we have added the statistical and systematic uncertainties together.
τ τ ours 0.391 0.668 0.767 0.819 similar results in the fourth and fifth column, and the results are much different from the charmonium cases.
To check the lepton flavor university, we add the Table IV, where we show the values of R ψ nS τ τ and R Υ nS τ τ , their definitions are similar, for example, The deviation of the ratio R ψ nS τ τ from the lepton flavor universality will indicate the presence of new physics beyond the Standard Model.In Table IV, we only give one set of results, because no matter which value we choose, 'NR', 'Re', 'Re+α s ' or 'Re+α 2 s ', we got the same results, so the uncertainties from model are all cancelled.For charmonium, ratios of R ψ nS τ τ are much different from each other, while for R Υ nS τ τ , because the τ mass is much smaller than that of bottomonium, so we got almost the same value for all the R Υ nS τ τ .Our results of R Υ nS τ τ are consist with the existed experimental data within their error bars.But we note that, our center values are smaller than one, while experimental central data are larger than one.
To cancel the model dependence of our theoretical predictions, we add another three tables, Table V, Table VI and Table VII, where we show the ratios of Γ(ψ(nS) → e + e − )/Γ(J/ψ → e + e − ) and Γ(Υ(nS) → ℓ + ℓ − )/Γ(Υ(1S) → ℓ + ℓ − ), since J/ψ can not decay to τ + τ − , we also give the ratio Γ(ψ(nS) → τ + τ − )/Γ(ψ(2S) → τ + τ − ).In cases of Υ(nS), we obtained the same central values for e, µ and τ final states, so we only show the ratio Γ(Υ(nS) → ℓ + ℓ − )/Γ(Υ(1S) → ℓ + ℓ − ), and the calculation is based on the e + e − final states in Table VII.predictions, whether the pure relativistic results, or with QCD corrections, consist very well with experimental data.And for Γ(Υ(4S)) Γ(Υ(1S)) , our results consist with the data inside brackets.In conclusion, we studied the leptonic decays of heavy vector quarkonium.For the charmonium decays, bad consistency happened between theoretical predictions and experimental data, there is still long way to go for the theory as well as experiment before we get good matching.For the bottomonium di-lepton decays, our theoretical predictions consist very well with the data.
We found the relativistic corrections are large and important both for the decays ψ(nS) → ℓ + ℓ − and Υ(nS) → ℓ + ℓ − .The relativistic effects become larger for a higher excited states.

TABLE I :
[45] spectra of the S wave cc and b b vectors in unit of MeV.'Ex' is the experimental data from PDG[45].

TABLE II :
Decay rates of the ψ(nS) → ℓ + ℓ − in unit of keV.'NR' is the non-relativistic result where only ψ 5 is kept and without QCD correction, 'Re' is the relativistic result without QCD correction, 'Re+α s ' means the relativistic result with the first order (α s ) QCD correction, and 'Re+α 2 Table III, our relativistic results without QCD correction consist

TABLE III :
Decay rates of the Υ(nS) → ℓ + ℓ − in unit of keV.'NR' is the non-relativistic result where only ψ 5 is kept and without QCD correction, 'Re' is the relativistic result without QCD correction, 'Re+α s ' means the relativistic result with the first order (α s ) QCD correction, and 'Re+α 2