The charm/bottom quark mass from heavy quarkonium at N$^3$LO

We determine the charm and bottom quark masses using the N$^3$LO perturbative expression of the ground state (pseudoscalar) energy of the bottomonium, charmonium and $B_c$ systems: the $\eta_b$, $\eta_c$ and $B_c$ masses. We work in the renormalon subtracted scheme, which allows us to control the divergence of the perturbation series due to the pole mass renormalon. Our result for the $\overline{\rm MS}$ masses reads ${\overline m}_{c}({\overline m}_{c})=1223(33)$ MeV and ${\overline m}_{b}({\overline m}_{b})=4186(37)$ MeV. We also extract a value of $\alpha_s$ from a renormalon-free combination of the $\eta_b$, $\eta_c$ and $B_c$ masses: $\alpha_s(M_z)=0.1195(53)$. We explore the applicability of the weak coupling approximation to bottomonium $n=2$ states. Finally, we consider an alternative computational scheme that treats the static potential exactly and study its convergence properties.

The cancellation of the leading renormalon of the pole mass and the static potential, first found in [12], led to the realization [13] that using threshold masses [13][14][15][16][17] (which explicitly implement the cancellation of the renormalon in heavy quarkonium observables) improves the convergence of the perturbative series. This has made these very precise computations useful not only for academical purposes but also for phenomenological applications. As a consequence, determinations of the bottom, or bottom and charm quark masses have been obtained using these new results [15,[18][19][20][21][22][23][24].
More recently, the B c spectrum has also reached N 3 LO precision [25]. This theoretical expression has not yet been used for phenomenology and confronted with experiment. One can use the B c spectrum to determine the heavy quark masses. In principle this system is more perturbative than charmonium. Hence, it should lead to a more accurate determination of the charm quark mass. We aim to do so in this paper. In particular, we consider specific energy combinations that are more suitable for clean theoretical analyses. This results in more accurate determinations of the charm quark mass and also allows an independent determination of α s .
The other main motivation of this paper is to study the feasibility of an alternative computational scheme that reorganizes the perturbative expansion of the above analyses. This scheme is characterized by solving the Schroedinger equation including the static potential exactly (to the order it is known). This incorporates formally subleading terms in the leading order (LO) solution. On the other hand the relativistic corrections to the spectrum are included perturbatively. This working scheme performs a partial resummation of higher order effects. This may accelerate the convergence of the perturbative series. This is indeed the effect seen in the cases where it has been applied (spectrum and decays) [26][27][28][29]. This scheme naturally leads to the organization of the computation in powers of v, the relative velocity of the heavy quark in the bound state.
Irrespectively if working with strict fixed-order perturbation theory or with an improved perturbation scheme, one has to implement the renormalon cancellation in the computation.
In this paper we use the RS threshold mass as the expansion parameter [15].
Finally, we explore the applicability of the weak-coupling approximation to the first excited state (n = 2) of heavy quarkonium, both in the strict weak-coupling approximation and in the alternative perturbation scheme. For this analysis we also consider renormalonfree combinations, which are theoretically cleaner.

II. GROUND STATE HEAVY QUARKONIUM ENERGY AT N 3 LO
In this section we determine the bottom and charm quark masses from the experimental values (as quoted from Ref. [30]) of the ground state masses of the bottomonium, charmonium and B c , and the corresponding theoretical expressions to N 3 LO. Whereas the bottomonium and charmonium spectra have been used before, the use of the B c spectrum is novel.
An analysis of the Υ(1S) mass to N 3 LO in the strict weak-coupling expansion, and using the RS threshold mass, has been done in Refs. [22,23]. Here we will follow a similar methodology but considering instead the η b . 1 The underlying reason for using pseudoscalar masses in this paper is that the experimental information of the B c spectrum is incomplete.
For the ground state, only the pseudoscalar mass is known. Therefore, the experimental set of masses we will use are the η b , η c and B c masses.

A. Determination of m b
We extract the bottom quark mass from the condition: (1) The determination of the bottom quark mass and the error analysis closely follows Refs. [22,23] (see those references for extra details and notation). We vary the parameters in the following way: ν f = 2 +1 −1 GeV, α s (M z ) = 0.1184 (12) [30] (with decoupling at m b = 4.2 GeV and at m c = 1.27 GeV; the specific location of the decoupling plays a marginal role in the determination of the bottom quark mass), N m = 0.5626(260), and (4/3)r 3 (m b ; N l ) = 1698.59 ± 1.74 [31]. 2 We determine the central value (and the renormalization scale associated error) using ν = 5 +3 −3 GeV. This scale is bigger than the scale used in Refs. [22,23], which was ν = 2.5 +1.5 −1 GeV. This is motivated by inspecting the scale dependence of the observable. In the RS scheme, ν ∼ 5 GeV is roughly the scale where the function is extremal. The RS' scheme suggests even higher scales. For scales smaller than 2 GeV the result shows a strong scale dependence. This motivates the minimum that we take for the ν, and we take the maximum symmetrically, thus fixing the scale variation.
1 Here and throughout this paper we refer to η b,c (1S) as η b,c , B c (1S) as B c and η b,c (2S) as η b,c . 2 Note that there have been some small changes in α s and r 3 with respect to Refs. [22,23]. Indeed with the new value of r 3 (the four loop coefficient relating the pole and the MS mass) the associated error is negligible: ∼ 0 MeV in all cases. Therefore, we will not explicitly display it in the following.
Using Eq. (1) in the RS and RS' approaches we extract, in MeV, respectively In Eqs. (3) and (5), m b has been determined from the RS masses at the scale ν = 2.5 GeV.
Setting ν = 5 GeV would not change much the value: bigger scales produce small changes in the value of the mass but smaller scales produce larger variations. In both cases reasonable renormalization scale variations in the relation between the RS masses and m b are well inside the uncertainties of our determination of m b . Overall, the uncertainties in m b are dominated by the variation of the renormalization scale ν in the fit.
Taking the average of Eqs.
(3) and (5) we obtain where we have rounded the ± variation of each parameter to the maximum (for the combined error we take the largest of both determinations) and added them in quadrature. To this number we add the leading charm corrections following the analysis of [22]. It is summarized by a shift in the final value of ∼ +2 MeV. Our final number reads For illustration, we show the value of the different orders in the perturbation series for the central value parameters: m b,RS (2 GeV) = (4185 + 145 + 58 Overall, we see a convergent pattern. Let us now study in more detail the convergence of the perturbative series, and the reliability of the error estimate in Eq. (7). In Fig. 1 we study the convergence of the perturbative expansion of M η b and its dependence on the factorization scale ν, which is the major source of uncertainty. We do so for both the RS and RS' schemes (even though the differences between both schemes are large at low orders they become quite small at N 3 LO).   In Table I, we use the central value in Eq. (6) to determine the RS masses (for latter use we also give the values for ν f = 0.7, 1 GeV). The difference between the values in the last column of Table I and the values in Eqs. (2) and (4) is completely marginal. In the values displayed in Table I we explicitly incorporate the error from Eq. (6) to show the sensitivity to the error of m b (m b ). As expected, and needed for a precise determination of the heavy quark mass, we find such sensitivity to be very large (see the yellow bands in Fig. 1).
In general, we observe that the convergence is quite reasonable irrespectively of the scheme (with ν f = 1 GeV the RS' scheme does not converge that well, though the final value of the bottom quark mass is very similar); except for the last (N 3 LO) term, where the factorization scale dependence becomes stronger. This could signal the importance of ultrasoft effects. In order to quantify their importance we have also plotted the N 3 LO prediction of M η b minus δE U S 10 , the pure ultrasoft contribution to the energy, which for a general quantum number reads where m r = m b /2 for bottomonium, m r = m c /2 for charmonium, and for the B c . L E n are the non-Abelian Bethe logarithms, numerically determined in Ref. [32] for l = 0 and in Ref. [7] for l = 0 for low values of n.
We observe that for scales larger than 2 GeV, δE U S 10 is small and well inside uncertainties. It is also comforting that the yellow error band in Fig. 1 (associated to the uncertainty of m b ) covers the difference between the NNLO and N 3 LO results in the RS and in the RS' scheme for the whole scale variation. This difference is often used as an alternative determination of the error. Actually, in the RS' scheme the distance between the different orders is quite small, which may ask for future more aggressive determinations of the bottom mass. We refrain from doing so until a more clear understanding of the renormalization scale dependence of the observable at low scales and of the role played by the ultrasoft scale is achieved.
As a final check, we have performed the fit to the Υ(1S) mass using the same setup. We find a shift of almost +20 MeV, well within the uncertainties of Eq. (6).

B. Determination of m c
We first extract the charm quark mass from the charmonium system using and follow the analysis made in the previous section adapted to this case. The central value determination of m c,RS/RS and the error estimates are computed using the following values for the parameters: ν = 2.5 +1.5 −1.0 GeV, ν f = 1 +0.5 −0.3 GeV, α s (M z ) = 0.1184 (12), N m = 0.5626(260), and (4/3)r 3 (m c ; N l ) = 1698.59 ± 1.74. The scale ν ∼ 2.5 GeV is chosen after inspecting the renormalization scale dependence of the observable, since around this value M (th) ηc is extremal in the RS scheme. We also observe that for scales smaller than 1.5 GeV there is a strong scale dependence. This fixes the scale variation we take. We also take a smaller value of ν f , more natural for charmonium.
Using Eq. (13) in RS and RS' approaches we extract, in MeV, respectively where m c is fixed from the RS masses at the scale ν = 1.5 GeV (considering ν = 2.5 or ν = 1 GeV produces small changes in m c ). Nicely enough, the uncertainty associated to the scheme dependence is quite small, i.e. the values obtained in Eqs. (15) and (17) where we have rounded the ± variation of each parameter to the maximum (for the combined error we take the largest of both determinations) and added them in quadrature.
We see that the perturbative expansions are reasonably convergent.
Before we study in more detail the numbers obtained above and the error analysis, we turn to our alternative determination of the charm mass. We aim to use the experimental value of M (exp) Bc . Nevertheless, such observable is strongly dependent on the bottom quark mass. To eliminate such dependence it is better to consider the following observable instead:  Table I that follow from Eq. (7). Overall, using Eq. (23) in RS and RS' approaches we extract, in MeV, respectively where m c is fixed from the RS masses at the scale ν = 1.5 GeV (again considering ν = 2.5 or ν = 1 GeV produces small changes in m c ). Nicely enough the uncertainty associated to the bottom quark mass is quite small and so is the scheme dependence, i.e. the values obtained in Eqs. (25) and (27)  Taking the average of the RS and RS' determinations of m c we obtain For illustration, we show the value of the different orders in the perturbation series for the central value parameters: parameters: We see that the perturbative expansions are reasonably convergent.
If we compare Eq. (18) with Eq. (28) the central values are almost identical. The determination from the B c has a slightly smaller error, even though this is strongly dependent on the scale variation we choose to take for the observables. Moreover, on theoretical grounds, the transfer energy in the B c (and η b ) is bigger than for charmonium. Therefore, we take this determination as more solid and choose it as our final determination. The RS and RS' masses determined from this value (and the m c associated error) are collected in Table II.  Let us now study in more detail the convergence of the perturbative series and the reliability of the error estimate. In Fig. 2   the perturbative expansion. We emphasize that we use the same mass in the four cases: Eq. (28) (and also Eq. (7) for M Bc − M η b /2). In the four cases the N 3 LO predictions are in perfect accordance with experiment (for not too small scales). We also show the sensitivity of the results to variations of the charm mass in the yellow band of Fig. 2. As expected, and needed for a precise determination of the heavy quark mass, we find such sensitivity to be very large. In the plots in Fig. 2 we observe that the convergence is quite reasonable irrespectively of the scheme.
Ultrasoft effects enter at O(α 5 s ). As we did for the bottomonium, in order to quantify their importance (or whether one can even compute their effects at weak coupling), we have also plotted the N 3 LO prediction minus the pure ultrasoft contribution in Eq. (12).
We observe that for the scales where our prediction is more robust (above 1.5-2 GeV) such contribution is small and well inside uncertainties. Actually, the fact that we obtain basically the same value for m c from M Bc − M η b /2 or M ηc , hints that possible nonperturbative effects are strongly constrained. If ultrasoft effects cannot be computed in perturbation theory one can roughly estimate their size (for scaling purposes) to be of order Λ 3 QCD r 2 . We would then expect to obtain different values for the charm mass from the two fits, since for the former the energy shift is and for the latter These two quantities are, in principle, different. Still, we basically find the same value for m c . This is compatible with assuming that either: 1) the ultrasoft scale can be dealt with in perturbation theory, 2) the nonperturbative effects are small, or 3) numerically, Eq. (33) is 1/2 of Eq. (34). In the scenario 1), the leading nonperturbative effects can be written as an expansion in terms proportional to local condensates of higher and higher dimensionality.
The dimension four corrections (proportional to the gluon condensate) were computed in [33,34], the dimension six in [35] and the dimension eight (for l = 0) in [36]. We refrain from using those corrections in this paper until a clear signal of the associated asymptotic behavior of the perturbative expansion related with these condensates is seen.
It is also comforting that the yellow error band associated to the uncertainty of m c in As a final check, we have performed the fit to the Ψ(1S) mass using the same setup. We find a shift of +14 MeV, well within the uncertainties of Eq. (28).

III. RENORMALON-FREE COMBINATIONS AND DETERMINATION OF α s
It is interesting to consider renormalon-free combinations. Particularly compelling is the combination On the one hand we know its experimental value. On the other hand, for this observable, the leading renormalon ambiguity has cancelled 3 . A possible u = 1 renormalon ambiguity associated to the kinetic term also exactly cancels in this combination. Therefore, this observable is (potentially) particularly good for the analysis of ultrasoft effects, as any noise in the signal associated to these renormalons disappears. We emphasize that the determinations of the bottom and charm mass in the previous section assume that perturbative computations are reliable at the ultrasoft scale, or at least that possible nonperturbative effects are small. Roughly, one can assume the nonperturbative effects of ∆ RF to scale as (note though that this estimate looses the ultrasoft logarithms) This combination of expectation values is not zero (though some cancellation is expected).
We can compare Eq. (35) with the theoretical prediction using the masses obtained in Sec. II. We show this in Fig. 3 for the RS and RS' schemes. We observe that the theory prediction is consistent with experiment within the expected uncertainties. We take this plot as an indication that we can use perturbation theory for the ground state of the three systems. Looking at Fig. 3, the incorporation of the perturbative ultrasoft contribution indeed seems to improve the agreement with data, even though the effect is small and can be perfectly masked by the uncertainties.
Indeed, assuming the weak-coupling approximation, ∆ RF eliminates most of the mass dependence (both on the bottom and the charm quark masses). This makes this observable specially sensitive to α s . If we use it to fit α s we obtain: where we take ν = 3 +2 −1 GeV, ν f = 1 GeV, and the masses from Tables I and II. In the last equality of Eqs. (37) and (39)  The difference between the NNLO and N 3 LO evaluation also gives an estimate for the error associated to higher order corrections. Again in this case it is of the order of the error. Note also that for the mass determinations our final error was bigger (and much bigger in some cases) than the difference between the N 3 LO and the NNLO evaluation. Therefore, for the final error we take the minimum-maximum error of each determination in Eqs. (37) and (39): This is perfectly consistent with the world average albeit with larger errors. We show the plot with this value of α s in Fig. 4, where the yellow band is the error quoted in Eq. (41).
We see that it is quite conservative. We leave possible improvements for the future.  In the previous sections we have seen that a strict (fixed-order) weak-coupling computation works well for the ground state of bottomonium, charmonium and the B c . We have also explored the n = 2 excitation for bottomonium. The convergence in this case was, at most, marginal. We now study a computational scheme that reorganizes the perturbative expansion such that it performs a selective sum of higher order corrections. We want to test if such scheme could improve/accelerate the convergence. In this method we incorporate the static potential exactly (to a given order) in the leading order Hamiltonian (the explicit ν dependence of the static potential appears at N 3 LO order and partially cancels with the explicit ν dependence of Eq. (12), the ultrasoft correction): where the static potential will be approximated by a polynomial of order N + 1, n a n (ν; r) .
We work in the following in the RS' scheme. Expressions for δm RS can be found in Ref. [28]. The static potential we will use then reads Eq. (45) encodes all the possible relevant limits: (a). The case ν r = ∞, ν f = 0 is nothing but the on-shell static potential at fixed order, i.e. Eq. (43). Note that the N = 0 case reduces to a standard computation with a Coulomb potential, for which we can compare with analytic results for the matrix elements. We use this fact to check our numerical solutions of the Schroedinger equation.
(b). The case ν r = ∞ (with finite non-zero ν f ) is nothing but adding an r-independent constant to the static potential (see the discussion in Ref. [44]).
(c). The case ν r =finite (and, for consistency, ν r ≥ ν f ). We expect this case to improve over the previous results, as it incorporates the correct (logarithmically modulated) short distance behavior of the potential. This has to be done with care in order not to spoil the renormalon cancellation. For this purpose it is compulsory to keep a finite, non-vanishing, ν f , otherwise the renormalon cancellation is not achieved order by order in N , as it was discussed in detail in Ref. [44]. We have explored the effect of different values of ν f in our analysis. Large values of ν f imply a large infrared cutoff. In this way our scheme becomes closer to a MS-like scheme. Such schemes still achieve renormalon cancellation, yet they jeopardize the power counting, as the residual mass δm RS does not count as mv 2 . As a consequence, consecutive terms of the perturbative series become bigger. Therefore, we prefer values of ν f as low as possible, with the constraints that one should still obtain the renormalon cancellation, and that it is still possible to perform the expansion in powers of nl correctly incorporates the N N LO corrections to the spectrum associated to the static potentials. It also includes higher order corrections (those generated by the iteration of the static potential). In order for this computational scheme to make sense, it first requires that the N → ∞ converges, or at least that the error is small compared with the relativistic correction. This is indeed so for the bottomonium ground state. We show the result of the computation of E 10 for bottomonium in Fig. 6. We see that the convergence is good. There is an error associated to the computation, which we estimate as the difference between the last two terms in the expansion. Setting ν r = 1 GeV makes the result more stable under scale variations (it can even be made more stable if we take ν f = ν r = 0.7 GeV). Typically, the ν r = ∞ and the ν r = 1 GeV evaluation agree in the 1-2 GeV ν range. Moreover, the differences with the strict fixed-order computations are not large. This allows us to take 4 E (0) nl as the (more or less well defined) LO on top of which one can incorporate relativistic and/or ultrasoft corrections. A similar analysis for B c , η c and bottomonium n = 2 with ν f = 1 GeV does not show such convergent behavior, though it does with ν f = 0.7 GeV.
If we instead consider renormalon-free observables, they typically show a better behavior. Therefore, these are the only combinations we consider in this paper and postpone a more detailed study of renormalon-dependent observables to future work. In this paper we usually work with ν r = 1 or ∞ GeV but explore in some cases ν r = 0.7 GeV. An earlier discussion, including also the charmonium ground state, can already be found in Ref. [28]. Nevertheless, the B c system has never been considered before. Therefore, following Ref. [28], we profit to get an estimate of the electromagnetic radius and the mean velocity of the B c . The analysis can be found in Fig. 7. We find The radius of the ground state B c is bigger than the radius of the ground state of bottomonium, and so is the mean velocity of the constituents. Overall, the convergence is good.
stands for the relativistic potential (V s is the total singlet potential) that contributes up to N 3 LO and also add the ultrasoft correction from Eq. (12). Overall the mass of the bound states reads where E (0) nl counts as v 2 , (0) n, l|δV |n, l (0) counts as v 4 (including also v 4 α s corrections) and δE U S nl as v 5 . Eq. (48) is numerically correct with N 3 LO precision and incorporates extra subleading terms (albeit in an incomplete way).
This computational scheme resums a subset of subleading corrections in the hope that they would account for the bulk of such subleading terms. This could be achievable if the higher order corrections that we infer from our knowledge of the static potential are responsible of the leading corrections. This computational scheme has already been applied in Refs. [28,29], where indeed an improved description of experiment was observed (see these references for extra details).
The expressions we use for the relativistic potential (valid also in the unequal mass case) are taken from Ref. [25], which uses results from Refs. [45,46]. We can use any of the bases for the potentials presented in that paper, which were referred as: Wilson, onshell, Coulomb or Feynman. At strict N 3 LO they all yield the same result. Since the computational scheme we implement in this section partially resums higher orders some dependence on the basis of potentials shows up. We have checked that, for the set of bases we consider, the dependence is quite small.
The computation of the relativistic corrections opens new issues compared with the static potential. In the case of the static potential the natural scale is ν ∼ 1/r, except in the O(α 4 s ) term where also the ultrasoft scale ν us appears. The case of the relativistic potentials is quite different. They are much more dependent on the hard, and above all, the ultrasoft scale.
Moreover, in order for the computation with the static potential to be a more or less reasonable approximation we need to have at least three or more terms (also important is the resummation of soft logarithms). For the case of the relativistic potentials, we have at most two terms. This, together with a much stronger scale dependence can make that inefficiencies of the description of the relativistic potentials get amplified when computing the expectation values. If, for instance, we consider M η b , we find that the relativistic corrections obtained with this computational scheme and from the strict fixed-order computation are very different. For the former the relativistic corrections are much larger than for the latter. This is different to the case of the static potential, where both approaches show comparatively small differences. This could be due to the fact that, for the relativistic corrections, we do not have enough terms (two at most) or that the unknown scale dependence of higher orders is very important and we are not describing it properly enough. Actually, preliminary computations suggest that the resummation of logarithms is indeed important. We postpone this discussion for future research.
To gain further insight we study in detail M Bc − M η b /2 − M ηc /2, which is a renormalonfree observable. We show our results in Fig. 8, where we plot our predictions for ν r = 1 GeV or ∞. We observe that the LO solution is nicely convergent in N in both cases, and also quite similar to the strict fixed-order result. This is quite remarkable as the individual contributions to the static potential do not converge that well for the case of the η b , and certainly less for the case of the η c and the B c (nor they individually agree with the strict fixed-order computation, in particular in the last two cases). Overall, the static potential gives an energy shift ∼ 60 MeV. For both values of ν r the difference with the strict fixed-order computation and among themselves is quite small for scales bigger than 1.5 GeV, for scales below the result is quite sensitive to the partial incorporation of higher order corrections.
Setting ν r = 1 GeV makes the result quite scale independent.
The relativistic corrections in Fig. 8 show a worse behavior. If we compute the expectation values of the relativistic potentials using the wave function obtained solving the Schrodinger equation with the N=3 static potential, the result is different from the result obtained from the strict fixed-order computation by around ∼ 20 MeV, as we can see in the plot. One can basically recover the strict fixed-order computation by evaluating the relativistic correction with N = 1 for the leading relativistic correction and with N = 0 for the subleading relativistic corrections. Setting a smaller value of ν r ∼ 0.7 GeV further deteriorates the agreement with experiment. We do not have a clear explanation for this fact. As we already mentioned before, we conjecture that it is fundamental to properly account for the different scales, as they may significantly affect the strength of the relativistic correction.
Nevertheless, it is also worth to mention that we are talking about a difference with the strict fixed-order result of order ∼ 20 MeV. This is much smaller than the differences we found for the relativistic corrections to the energies of each individual state. This strong cancellation may indicate that effects such as accounting for the different scales or truncating the series of the relativistic potentials produce smaller errors for the observable ∆ RF . Actually, ∼ 20 MeV is of the order of the difference between the RS and RS' computation at strict N 3 LO fixedorder. Finally, we explore the effect of the ultrasoft term δE U S 10 . We find it is comparatively small for scales bigger than 1.5 GeV.
We now turn to n = 2 bottomonium states. Following the discussion in Sec. III, we consider first the renormalon-free combination M h b −M η b . We plot the results of our analysis in Fig. 9. The behavior of the static potential is convergent for this energy combination, though less than in the previous case (higher order corrections could still give non-negligible contributions). This is to be expected since, for the first time, we consider the n = 2 states where the weak coupling approximation should be worse (still, n = 2 Pwave states behave typically better than 2S states). The static potential gives an energy shift ∼ 400 MeV. Setting ν r = 1 GeV does not make the result significantly more scale independent for the static potential (in this respect ν r = 0.7 GeV produces flatter curves). Indeed, the difference between both computations is small for the static potential and also not very different from the strict fixed-order result. The incorporation of the relativistic corrections (and the ultrasoft correction) improves the agreement with experiment. It also increases the dependence on the scale, more for the case with ν r = ∞ than for the case ν r = 1 GeV (a flatter curve is indeed obtained for a smaller value of ν r = 0.7 GeV, as expected). The case ν r = ∞ is also quite similar to the strict fixed-order computation. We can estimate the size of the shift of the O(v 4 ) to be of order 200 MeV. The difference with experiment is of order 100 MeV, which is consistent with expected uncertainties.
Finally, we consider the renormalon-free combination M η b − M h b . We show the analysis in Fig. 10. The behavior of the static potential is not really convergent for this energy combination. This is to be expected since it is unclear whether the weak-coupling expansion is valid for 2S states (somewhat it depends on the observable, as we saw in Ref. [28]). The scale dependence on ν of this observable is large, and changing ν r does not improve the scale dependence. The difference with the strict fixed-order computation is not large. We cannot really give quantitative estimates. On the other hand, it is fair to say that, had we considered M η b − M η b instead, the general convergence pattern would have been better. The fact that we have considered a more fine-tuned observable has amplified the relative errors.  [47], HPQCD14 [48], and Dehnadi et al. [49] use N 3 LO low n moments equated to the corresponding experimental moments evaluated using the available experimental information (masses and decays) supplemented with lattice information and/or assuming perturbation theory at high energies. Penin et al. [50] and Beneke [51] uses nonrelativistic (with Coulomb resummation) N 3 LO large n moments (the first reference uses a partial result). Pineda et al. [52] and Hoang et al. [53] use NNLO and (a partial) NNLL expression for nonrelativistic (with Coulomb resummation) large n moments. Kiyo [54], Ayala [23], Mateu [24] and this paper use N 3 LO expressions for the heavy quarkonium masses. of using the renormalon-free energy combination: M Bc − M η b /2 − M ηc /2, which is weakly dependent on the heavy quark masses but shows a strong dependence on α s . We use this observable to obtain a determination of α s : This number is perfectly compatible with the world average, albeit with larger errors.
We have also explored the applicability of the weak coupling approximation to n = 2 states. We have limited the analysis to bottomonium. The N 3 LO prediction for is somewhat worse: the experimental value is reproduced by the N 3 LO prediction at around ν ∼ 1.5 GeV but the convergence is not good and the renormalization scale dependence is rather large. Signer '08, NR moments, NNLO+NNLL * Large n sum rules: FIG. 12: m c (m c ) determined from spectroscopy or from sum rules in the strict weak-coupling approximation. Chetyrkin et al. [55], HPQCD14 [48], and Dehnadi et al. [49] use N 3 LO low n moments equated to the corresponding experimental moments evaluated using the available experimental information (masses and decays) supplemented with lattice information and/or assuming perturbation theory at high energies. Signer [56] uses NNLO and (a partial) NNLL expression for nonrelativistic (with Coulomb resummation) large n moments. JLQCD [57] uses N 3 LO (without Coulomb resummation) for large n. Kiyo [54], Mateu [24] and this paper use N 3 LO expressions for the heavy quarkonium masses.
Finally, we have also studied an alternative computational scheme that reorganizes the perturbative expansion by treating the static potential (truncated to order α N +1 s ) exactly.
We have observed that the solution coming from the static potential converges (as we increase N ) for the bottomonium ground state mass for ν f = 1 GeV. Incorporating the resummation of soft logarithms (ν = 1/r for 1/r < ν r ) further diminishes the scale dependence. A similar picture holds for the B c and charmonium if we set ν f = 0.7 GeV (also, but with less convergence, for the P-wave n = 2 bottomonium state). Nevertheless, the relativistic corrections do not generally show a good behavior. The difference with the corresponding relativistic strict fixed-order computation is large. However, if we turn to renormalon-free energy combinations the situation improves. First, the convergence of the static potential is better, specially for the combinations M Bc − M η b /2 − M ηc /2 and M h b − M η b (for M η b − M h b the improvement is marginal). The relativistic corrections also show a better behavior. The theoretical prediction for M h b − M η b shows a better agreement with experiment after the inclusion of the relativistic corrections. Nevertheless, this is not really so for M Bc − M η b /2 − M ηc /2, even though the degree of cancellation between very large individual relativistic corrections is still quite remarkable. At this stage we cannot make firm statements about the behavior of the relativistic corrections. We can speculate though that including more terms in the perturbative expansion and/or performing a renormalization group analysis could improve their behavior. On top of that, the different behavior between renormalonsensitive and renormalon-free observables (for which the renormalon cancellation is achieved from the start) may indicate that small inefficiencies in the renormalon cancellation using threshold masses get amplified in this computational scheme: affecting direct determinations of the heavy quarkonium masses, but largely canceling for renormalon-free observables. We hope to come back to all these issues in the near future.