Calculation of Quarkonium Spectrum and $m_b, m_c$ to Order $\alpha_s^4$

We include two loop, relativistic one loop and second order relativistic tree level corrections, plus leading nonperturbative contributions, to obtain a calculation of the lower states in the heavy quarkonium spectrum correct up to, and including, $O(\alpha_s^4)$ and leading $\Lambda^4/m^4$ terms. This allows us, in particular, to obtain a model independent determination of the pole masses of the $b, c$ quarks, $$m_b=5 015^{+110}_{-70} mev; m_c=1 884^{+222}_{-133} mev$$ to which correspond the $\bar{\hbox{MS}}$ masses, $$\bar{m}_b(\bar{m}_b^2)=4 453^{+50}_{-32} mev; \bar{m}_c(\bar{m}_c^2)=1 547^{+169}_{-102} mev.$$ The decay $\Gamma(\Upsilon\to e^+e^-)$ is found in agreement with experiment, $$\Gamma(\Upsilon\to e^+e^-)=1.135^{+0.27}_{-0.29} kev (\hbox{exp.}=1.320\pm0.04 kev),$$ and the hyperfine splitting is predicted to be $$M(\Upsilon)-M(\eta)=48.5^{+15.7}_{-12.2} mev.$$


Introduction
In recent years it has become possible to perform rigorous QCD analyses of heavy quarkonium systems, and this due to two reasons. First, radiative corrections have been calculated to increasing order of accuracy. The one loop corrections to the nonrelativistic (NR) spin-independent (SI) potential were calculated already in 1980 [1] . This was extended in refs. 2, 3 to the spin-dependent corrections, and in ref. 4 by including the velocity corrections to the SI part. Finally, the two loop nonrelativistic, spin independent corrections to the potential have been evaluated recently [5] .
Secondly, Leutwyler and Voloshin (ref. 6; see also ref. 7) have shown how to take into account, to leading order, nonperturbative (NP) effects, associated with the nonzero value of various condensates, of which the leading contribution is that of the gluon condensate α s G 2 . This has been implemented, together with the potential obtained with radiative corrections to one loop, in refs. 4, 8 where a study of bound statesbb with nl = 10, 20, 21 andcc sates with nl = 10 was given 1 . The analysis was extended in refs. 9, 7 with the inclusion of size effects and higher condensates.
The overall conclusion of these analyses is that pure QCD, without recourse to introducing phenomenological interactions, produces a good description with manageable errors of thebb ground state and, to a lesser extent, of the splitting M (Υ ) − M (η b ) and the decay Υ → e + e − . The description of the ground state of cc and of the excited states n = 2, l = 0, 1 ofbb was shown to be even less reliable: the corrections are large, in some cases larger than the nominally leading terms. Still, it was possible, by using the renormalization point µ as a free parameter, to get a fairly accurate description of all n = 2 states including tensor and LS splittings [8] .
In the present paper we extend this analysis by including the two loop corrections to the SI, NR potential recently calculated [5] adding also velocity corrections to certain one loop pieces to get a calculation accurate up to, and including, O(α 4 s ) corrections. By taking into account the leading nonperturbative terms, we also include in the analysis corrections of order Λ 4 /m 4 . The main improvement so obtained is that we get an extremely stable and precise determination of the ground state ofbb and, to a lesser extent,cc. If we invert the calculations we can deduce quark masses from the masses of the Υ, J/ψ particles. We then find, for the pole masses, (1.1b) The error includes the estimated theoretical error of the calculation, see the text. Note that (1.1) are very precise: as stated they are correct to order α 4 s , and leading, O(Λ 4 /m 4 ) nonperturbative effects. This is to be compared with estimates based on sum rules [10] which are only accurate to order α 2 s , or previous bound state calculations [4] , accurate only to third order in α s .

The Effective Potential
We follow the method of effective potentials of Gupta et al [3] , and the renormalization scheme of ref. 4. The Hamiltonian for quarkonium may then be split in the following form: and can (and will) be solved exactly. H 1 may be written as and log rµ r , (2.1d) Here the running coupling constant has to be taken to three loops. For the values of the constants entering above formulas, cf. the Appendix. A few words are due on Eqs. (2.1). First of all, they only take into account the perturbative part of the calculation; NP effects will be incorporated later. Secondly, it should be noted that H 1 contains a velocitydependent one loop piece, V s.rel . This is because the average velocity in a Coulombic potential is |v| ∼ α s , hence a calculation correct to order α 4 s requires tree level O(v 2 ) and one loop O(|v|) contributions. All these terms in H 1 may be treated as perturbations to first order, except V (L) . For this, the second order perturbative contribution is required as this also produces a correction of order α 4 s . A last comment concerns the renormalization scheme. We have followed ref. 4 in renormalizing α s in the MS scheme; but the mass m that appears in Eqs. (2.1) is the two loop pole mass. That is to say, it is defined by the equation, where S 2 (p / , m) is the quark propagator to two loops. One can relate m to the MS parameter, also to two loop accuracy, using the results of refs. 11: Here the δ V E nl may be easily evaluated with the formulas in the Appendix to ref. 4. We define generally the analogue of the Bohr radius, and then, (3.2e) We recall that constants are collected in the Appendix. For the vector states (Υ, Υ ′ , Υ ′′ ; J/ψ, ψ ′ , . . .) one has to add the hyperfine shift, at tree level, The calculation of the second order contribution of V E nl , is nontrivial, and may be found in the Appendix. We define and then one has, for the lowest states, In addition to this one has to consider the nonperturbative (NP) energy splittings. The dominant ones are associated with the gluon condensate and are [6] δ NP E nl = mǫ nl n 2 π α s G 2 na 2 4 = m ǫ nl n 6 π α s G 2 (mC Fαs ) 4 ; ǫ 10 = 1 872 1 275 , ǫ 20 = 2 102 1 326 , ǫ 21 = 9 929 9 945 .

(3.4)
Because α s G 2 ∼ Λ 4 , this is of order (Λ/m) 4 albeit with large coefficients: for all terms we have a fourth power of α s in the denominator, and for n > 1 the n 6 in the numerator of Eq. (3.4) grows very quickly out of hand. In fact it is the size of this term that limits the range of validity of our type of ab initio calculation.

Higher corrections
Besides the corrections reported in the previous subsection there are a few pieces of the higher order corrections that are known. First of all we have the relativistic, O(v 2 ) corrections to the one loop potential [4] . These produce corrections of higher order, α 5 s , but they are logically independent of the three loop ones that would produce terms of the same order but, presumably, smaller because of the extra 1/π 2 characteristic of radiative corrections. These corrections may be incorporated and then can be considered to give an indication of the error committed in neglecting higher order perturbative corrections. They produce, for the ground state, the energy shift (ref. 4; typos corrected in ref. 8), (3.5) Next we have higher order NP corrections. These include finite size corrections, estimated in ref. 7, and contributions of higher dimensional operators, some of which were evaluated in ref. 9. The last produce the shifts, and may be used to estimate the size of the higher order NP contributions. For the quark condensate the vacuum saturation approximation is assumed, and the value of κ is taken from ref. 10. For G 3 one takes the value 0.065 GeV 6 . Anyway, these quantities are poorly known.
It is important to realize that both (3.5), (3.6) should be taken as indications. With respect to the first, there is no guarantee that the coefficient of the three loop correction is not so large that it offsets the factors of 1/π; indeed, this already happens to two loops where the coefficient is large, b 1 ≃ 24. With respect to (3.6), and apart from the fact that it does not include all the higher dimensional operators (those associated with size corrections are neglected 2 ), it is clear that one cannot consider rigorously a contribution O(Λ 6 /m 6 ) so long as the radiative corrections to the O(Λ 4 /m 4 ) terms are not known. Nevertheless, we consider (3.5) and (3.6) as very useful for estimating the theoretical uncertainties of our calculation.

Numerical Results
Using the formulas deduced above one can evaluate the spectrum of heavy quarkonium systems. In principle one should take m, Λ, α s G 2 from other sources and predict the masses of the quarkonium states. In practice it is better to use the known masses of the lowest states (Υ and J/ψ) to evaluate the quark masses. The reason is that this produces by far the more precise evaluation available at present of these parameters, especially in the case of the b quark. The other parameters we take from independent sources. For the QCD parameter Λ we take, throughout this section, and for the gluon condensate, very poorly known, Another matter to be discussed is the choice of the renormalization point, µ. As our equations (3.2, 3) show, a natural value for this parameter is for states with the principal quantum number n, and this will be our choice. For states with n = 1 the results of the calculation will turn out to depend very little on the value of µ, provided it is reasonably close to (4.1c). Higher states are another matter; we will discuss our choices when we consider them.  Table I   In the estimate of the errors, the condition µ = 2/a is maintained satisfied when varying Λ while for the error due to the variation of µ the other parameters are kept fixed (i.e., one no more has then µ = 2/a). The dependence of m b on µ should be taken as an indication of the theoretical uncertainty of our calculation. To sestimate other theoretical errors in our evaluation we proceed as follows. We either calculate m b including the full O(v 2 ) corrections to one loop, Eqs. We consider that the best result is that of O(α 4 s ) reported in Table I, and take the difference with the quantities given in Eqs.  The values of α s (µ 2 ),α s (µ 2 ) corresponding to µ 2 = 7.019 GeV 2 are α s (µ 2 ) = 0.24 ,α s (µ 2 ) = 0.40.
The "theoretical" error coming from higher dimensional operators and higher order perturbative terms, Eqs. (4.2), is comfortably smaller that the errors due to the uncertainty on Λ, α s G 2 . We will henceforth omit these errors, so as not to double count them, and consider that the theoretical error is only that due to varying µ 2 by 25%. If we now compose all the errors quadratically, then we obtain the estimate reported in the Introduction, Eqs. (1.1).  The corrections here are fairly large, particularly the radiative correction [12]  and Note that, when varying Λ, α s G 2 , we have varied m b according to Eq. (4.3), but we have not varied m b when varying µ.
Higher order NP corrections due to the higher dimensional operators introduced in Eq. (3.6) are also known for the decay rate (see ref. 9). They read Size corrections, however, are not known now. δ N N P would produce a shift in the decay rate of ∼ 0.11 keV, smaller than the contribution of α s G 2 or the uncertainty caused by e.g. varying µ as in Eq. (4.3). We do not include δ N N P either in the evaluation or the error estimate.
The result for the decay is in reasonable agreement with experiment, Composing the errors we obtain the figures quoted in the Introduction, Eqs. (1.2, 3).

Higher states (n = 2) ofbb
The masses of the states with n = 2 will be next determined. As is clear from the expressions (3.  We only present the errors that follow from variation of the scale µ 2 by 25%; slightly smaller ones are produced by the errors of Λ, α s G 2 . We do not explicit this: because of the size of the errors in Eq. (4.7) there is no point in going for a more detailed error analysis.
Although they have decreased from th one loop evaluations (e. g. ref. 4),the errors are still fairly large here; within them, there is compatibility between theory and experiment. Agreement to a few MeV for both states is obtained if choosing µ 2 = 0.75/a ≃ 2.3 GeV 2 or keeping µ = 1/a and taking α s G 2 = 0.036 GeV 4 : unlike for the states with n = 1 we have now strong dependence of the results on the parameters of the calculation. This is due to the large size of the corrections, perturbative and (especially) nonperturbative. This last is made more apparent when considering contributions of higher dimensional operators [9] , which get completely out of hand for nl = 20 and are very large for nl = 21. In this context, it is satisfactory to realize that it is for this last state (21) for which agreement with experiment is best, and theoretical errors smaller.
We will not discuss here the spin and tensor splittings among the states with nl = 21. The inclusion of three loop corrections to the potential only affects their calculation in that the preferred value for m b will be different now, which is a minute effect compared with the uncertainties of the calculation: one should realize [8] that, while the NP corrections to the energy levels with principal quantum number n contain a coefficient n 6 /α 4 s , wave functions at the origin get a factor ∼ n 8 /α 6 s . This of course is what makes the calculation of M (Υ )−M (η b ) and the decay Υ → e + e − much less reliable than that of M (Υ ) (or, equivalently, m b ) and what makes the evaluation of tensor and spin splittings with n = 2 somewhat marginal. All one can do here is fit µ to the data; this is the procedure followed in ref. 8, and we have nothing new to report on this.  (4.9) As is obvious from these equations, the errors are now much larger than for the b quark case, but our determination of m c still competes in accuracy with those based on QCD sum rules.

Discussion
The calculations of this paper are rather straightforward; but there are a few points that merit further discussion. First of all, our values for the quark masses are somewhat larger than existing estimates based on sum rules; for the MS masses, of 100 ∼ 300 MeV, cf. ref. 10. In our opinion this is due to the influence of the terms of order α 3 s , α 4 s which we take into account, but which the sum rule evaluations, that stop at O(α 2 s ), do not. Thus we consider our estimates to be the more precise and reliable ones. Secondly, and from Table I, it may appear that the series is diverging: from the O(α 2 s ) to the O(α 3 s ) evaluation, m b increases by 106 MeV but from the last to the O(α 4 s ) the increase is of 157. Actually, convergence is reasonably good. The increase between O(α 3 s ) and O(α 4 s ) is due to two independent factors: inclusion of the two loop corrections to the potential, responsible for 46 MeV, and the relativistic corrections. Of these, 64 MeV for tree level corrections and 40 MeV for the mixed one loop-velocity correction. Each of the effects is small. Thus, if we included velocity corrections at every loop, the variation from zero to one to two loops would be diminishing. This is apparent if we compare the value obtained in ref. 4 (corrected for the increased values of Λ, α s G 2 we are using now) with our results, with a variation of only 60 MeV. One can see this more clearly if we arrange the calculation in increasing number of loops including, at every step, the pertinent relativistic corrections 4 or including all loop corrections, but in increasing order of the velocity corrections: Finally, we remark the satisfactory stability we now have against variation of the renormalization scale, µ. This stability of the results against changes of µ is made apparent by the fact that even multiplying or dividing the central value µ 2 = 7.019 by a factor of two only alters the central value of m b = 5.015 GeV by 98 MeV. The stability is due mostly to the inclusion of two loop effects; but attention should also be paid to the stabilizing influence of the NP corrections. These corrections are larger for larger µ, exactly the opposite to what happens to the perturbative corrections. One could even fix optimal values of µ as those where the combined perturbative-NP effects would show a minimal dependence on µ. This is essentially the procedure adopted in ref. 4. Here, and because of our much smaller dependence on µ, we need not have recourse to such methods.

Second order contribution
We now calculate δ This simple method gives correctly the coefficients of log aµ, log 2 aµ and misses the constant term by ∼ 10%.
and then δ