The charm/bottom quark mass from heavy quarkonium at N3LO

We determine the charm and bottom quark masses using the N3LO perturbative expression of the ground state (pseudoscalar) energy of the bottomonium, charmonium and Bc systems: the ηb, ηc and Bc 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 MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} masses reads m¯cm¯c=122333\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\overline{m}}_c\left({\overline{m}}_c\right)=1223(33) $$\end{document} MeV and m¯bm¯b=418637\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\overline{m}}_b\left({\overline{m}}_b\right)=4186(37) $$\end{document} MeV. We also extract a value of αs from a renormalon-free combination of the ηb, ηc and Bc masses: αs(Mz) = 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 [14], and later in [15,16], led to the realization [16] that using threshold masses [16][17][18][19][20] (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 [18,[21][22][23][24][25][26][27][28].
More recently, the B c spectrum has also reached N 3 LO precision [29]. 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 -1 -

JHEP09(2018)167
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) [30][31][32][33]. 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 [18].
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 mainly consider renormalon-free combinations, which are theoretically cleaner.

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. [34]) 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. [25,27]. 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.

Determination of m b
We extract the bottom quark mass from the condition: The determination of the bottom quark mass and the error analysis closely follows refs. [25,27] (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) [34] (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 [35]. 2 We determine the central value (and the renormalization scale associated error) using ν = 5 +3 −3 GeV. This scale is bigger than the scale 1 Here and throughout this paper we refer to η b,c (1S) as η b,c , Bc(1S) as Bc and η b,c (2S) as η b,c . 2 Note that there have been some small changes in αs and r3 with respect to refs. [25,27]. Indeed with the new value of r3 (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. used in refs. [25,27], 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.
Using eq. (2.1) in the RS and RS' approaches we extract, in MeV, respectively In eqs. 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 [25] 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. (2.7). In figure 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 1, we use the central value in eq. (2.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    in table 1 we explicitly incorporate the error from eq. (2.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 figure 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 m r = m b m c /(m b + m c ) for the B c . L E n are the non-Abelian Bethe logarithms, numerically determined in ref. [36] for l = 0 and in ref. [9] for l = 0 for low values of n.

JHEP09(2018)167
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 figure 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. (2.6).
We now compare our result with previous determinations of m b (m b ) from spectroscopy or from sum rules. We show such comparison in figure 2. We find that all determinations are perfectly consistent with each other. Determinations from sum rules with small number of momentums (low n sum rules) seem to give slightly smaller values of the mass than large n or spectroscopy determinations. We notice that the latter include the Coulomb resummation. Therefore, the physics they are sensitive to is different. Nevertheless, at this stage this difference is not statistically significant, and we are not in the position to push this discussion further. It is also worth noting that, at present, partial NNLL sum rules determinations produce slightly larger errors than other determinations. It would be interesting to see if the complete result, or its combination with the large n N 3 LO sum result, may lead to a more precise determination.
We now turn our discussion to the comparison between different spectroscopy determinations. Our value is slightlty smaller than the determinations in [26][27][28] (but well within each other's uncertanties). Compared with the determinations of Ayala et al. [27] and Mateu et al. [28], a part of the difference can be traced back to the fact that in the two latter references the Υ plays a predominant role in the fit (which leads to a larger value for the bottom mass than fits to the η b mass). In the first reference because the fit was only made to the Υ, in the second reference because, even if both the vector and the pseudoscalar were considered in the fit, their method weights the most the particle that produces the smaller errror, which in their paper is the vector. It is worth mentioning that if one looks at the individual fits in [28] to the η b mass one finds a smaller value of m b (m b ), quite close to ours. As at present all these determinations are consistent with each other, there is not a clear reason to prefer one versus the other. As we have already discussed before, the resummation of large logarithms may help to solve this problem. We hope to come back to this issue in the future. The other fact that produces a different central value of m b (m b ) is the choice of the factorization scale µ to fix the central value. In this paper we have chosen a larger value of the factorization scale than in [27,28]. We somewhat feel that a smaller value of the factorization scale leads to a region where the perturbative series does not behave well enough. In this respect the factorization scales we choose are closer to the factorization scales used in [26]. Actually, the difference with respect to the value obtained in [26] can be traced back to the fact that in that reference an average of the determination of the bottom mass from the Υ and the η b is made. Again pure determinations from the η b Figure 2. m b (m b ) determined from spectroscopy or from sum rules in the strict weak-coupling approximation. Chetyrkin et al. [37], HPQCD14 [38], and Dehnadi et al. [39] 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. Finite energy [40] or Laplace [41] sum rules have also been applied for low n. Penin et al. [42] and Beneke et al. [43,44] uses nonrelativistic (with Coulomb resummation) N 3 LO large n moments (the first reference uses a partial result). Pineda et al. [45] and Hoang et al. [46] use NNLO and (a partial) NNLL expression for nonrelativistic (with Coulomb resummation) large n moments. Kiyo et al. [26], Ayala et al. [27], Mateu et al. [28] and this paper use N 3 LO expressions for the heavy quarkonium masses. mass are closer to our numbers. It is also worth mentioning that the error of the m b (m b ) in [26] from the Υ or η b mass is quite similar (this is the reason why the value obtained in [26] is more or less the aritmetic average between the values of m b (m b ) obtained from the Υ and the η 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: (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. (2.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. 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: 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. We see that the perturbative expansions are reasonably convergent. If we compare eq. (2.18) with eq. (2.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 2.
Let us now study in more detail the convergence of the perturbative series and the reliability of the error estimate. In figure 3 we plot the theoretical prediction for M ηc and  M Bc − M η b /2 both in the RS and RS' as a function of ν truncating at different orders in the perturbative expansion. We emphasize that we use the same mass in the four cases: eq. (2.28) (and also eq. (2.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 figure 3. 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 figure 3 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. (2.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 .

JHEP09(2018)167
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. (2.33) is 1/2 of eq. (2.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 [47,48], the dimension six in [49] and the dimension eight (for l = 0) in [50]. 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 figure 3 covers the difference between the NNLO and N 3 LO results in the RS and in the RS' scheme for the whole scale variation in most cases. This difference is often used as an alternative determination of the error. Actually in the RS' scheme the difference between the different orders is quite small.
As a final check, we have performed the fit to the Ψ(1S) mass using the same setup. We find a shift of +41 MeV. This is a little bit less than one standard deviation of eq. (2.18), which we obtained from the η c . This is consistent with the error estimate. This difference is related with the hyperfine splitting. It has been argued that such energy splitting is rather sensitive to the resummation of large logarithms [51], which we do not incorporate in this paper (we hope to come back to this issue in the future). On top of that, the hyperfine splitting of the B c minus Υ/2 energy combination is expected to be much smaller than for charmonium. This reinforces in our choice of M Bc − M η b /2 as the optimal observable to determine the charm mass.
We now compare our result with previous determinations of m c (m c ) from spectroscopy or from sum rules. We show such comparison in figure 4. We find that all determinations are perfectly consistent with each other, but there is a small tension with determinations using low n sum rules, where our number is in the low range. In any case, the difference is hardly statistically significant. Low n sum rules determinations seems to give slightly larger values of the mass than large n or spectroscopy determinations. We notice that the latter include the Coulomb resummation. Therefore, the physics they are sensitive to is different. Nevertheless, at this stage this difference is no statistically significant, and we are not in the position to push this discussion further. It is also worth mentioning that, at present, there are no large n N 3 LO sum rules analyses (including the Coulomb resummation).
We now turn our discussion to the comparison between different spectroscopy determinations. The discussion follows parallel to the discussion we had for the bottom mass. Our value is slightlty smaller than the determinations in [26,28] (but well within each

JHEP09(2018)167
other's uncertanties). Compared with the determinations of Mateu et al. [28], a part of the difference can be traced back to the fact that in [28] the J/Ψ plays a dominant role in the fit (which leads to a larger value for the charm mass than fits to the η c mass). The reason is that, even if both the vector and the pseudoscalar were considered in the fit, their method weights the most the particle that produces the smaller errror, which in their paper is the vector. It is worth mentioning that if one looks to the individual fits in [28] to the η c mass one finds a smaller value of m c (m c ), quite close to ours. As at present all these determinations are consistent with each other, there is not a clear reason to prefer one versus the other. As we have already discussed before, the resummation of large logarithms may help to solve this problem. We hope to come back to this issue in the future. The other fact that produces a different central value is the choice of the factorization scale µ used to fix the central value of m c (m c ). In this paper we have chosen a larger value of the factorization scale than in [28]. We somewhat feel that a smaller value of the factorization scale leads to a region where the perturbative series does not behave well enough. In this respect the factorization scales we choose are closer to the factorization scales used in [26]. Actually, the difference with respect to the value obtained in [26] can be traced back to the fact that in that reference an average of the determination of the charm mass from the J/Ψ and the η c is made. Again pure determinations from the η c mass are closer to our numbers. It is also worth mentioning that the error of the m c (m c ) in [26] from the J/Ψ or η c mass is quite similar (this is the reason why the value obtained in [26] is more or less the aritmetic average between the values of m c (m c ) obtained from the J/Ψ and the η c ).

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). Figure 4. m c (m c ) determined from spectroscopy or from sum rules in the strict weak-coupling approximation. Chetyrkin et al. [52], HPQCD14 [38], and Dehnadi et al. [39] 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 (in [53] this was indeed assumed since threshold using quark-hadron duality). Finite energy [54] or Laplace [41] sum rules have also been applied for low n. Signer [55] uses NNLO and (a partial) NNLL expression for nonrelativistic (with Coulomb resummation) large n moments. JLQCD [56] uses N 3 LO (without Coulomb resummation) for large n. Kiyo et al. [26], Mateu et al. [28] and this paper use N 3 LO expressions for the heavy quarkonium masses.

JHEP09(2018)167
We can compare eq. (3.1) with the theoretical prediction using the masses obtained in section 2. We show this in figure 5 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 figure 5, 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: last equality of eqs. (3.3) and (3.5) we have symmetrized the error. The lower value in the variation in ν is motivated by the fact that, at smaller scales, it is not possible to find a value of α s that agrees with experiment. This hints that at smaller values of the scale, the observable does not behave as a (logarithmically modulated) polynomial in α s , with a small α s . The difference between the RS and RS' definitions is of the order of the error we obtain (unlike what happened for the mass determinations, where the difference was very small). 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. This is perfectly consistent with the world average albeit with larger errors. We show the plot with this value of α s in figure 6, where the yellow band is the error quoted in eq. (3.7). We see that it is quite conservative. We leave possible improvements for the future. We now move from the ground state to the more shaky domain of excitations. Our aim is to see how reliable possible predictions of excited states are. We consider only n = 2 bottomonium states, and energy differences among them and with the ground state. We show the plots in figure 7, where we work in the RS' scheme with ν f = 1 GeV. We see that even though they go in the right direction, the errors are large. The convergence is at most marginal. For the energy difference between the 1P and 1S state: M h b − M η b , the agreement with experiment is good for the N 3 LO prediction, albeit with large uncertainties, only looking at the scale variation one finds ∼ ±150 MeV. The M η b −M h b can be reproduced by the N 3 LO result at ν ∼ 1.5 GeV but the convergence is not good and the renormalization scale dependence is rather large. In both cases, the uncertainties of the perturbative expansion are large in comparison with the ultrasoft contribution.

Higher orders
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. (2.12), the ultrasoft correction): where the static potential will be approximated by a polynomial of order N + 1, n a n (ν; r) .

JHEP09(2018)167
In principle we would like to take N as large as possible (though we also want to explore the dependence on N ). In practice we take the static potential at most up to N=3, i.e. up to O(α 4 s ) including also the leading ultrasoft corrections. This is the order to which the coefficients a n are completely known: The O(α s ) term was computed in ref. [57], the O(α 2 s ) in refs. [58,59], the O(α 3 s ) logarithmic term in ref. [60], the light-flavour finite piece in ref. [61], and the pure gluonic finite piece in refs. [62,63].
We work in the following in the RS' scheme. Expressions for δm RS can be found in ref. [32]. The static potential we will use then reads Eq. (4.4) 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. (4.2). 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 Schrödinger 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. [64]).
(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. [64]. 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 α s .
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 (0) 10 for bottomonium in figure 8. 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 4 Strictly speaking the LO would be 2m b .
ν (GeV) Figure 9. Plot of the radius and v 2 of the B c in the RS' scheme with ν f = 1 GeV and ν r = ∞ (left panels), and with ν f = ν r = 0.7 GeV (right panels). They are obtained from solving eq. (4.1) for N = 0, 1, 2, 3. Continuous lines are evaluated with ν r = ∞ and dashed lines with ν r = 1 GeV.
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. [32]. Nevertheless, the B c system has never been considered before. Therefore, following ref. [32], we profit to get an estimate of the electromagnetic radius and the mean velocity of the B c . The analysis can be found in figure 9. 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. Once we have our new LO (the solution from the static potential), we can consider the incorporation of the relativistic and ultrasoft corrections. With the accuracy of this work, we only have to take the expectation value of δV where 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. (2.12). Overall the mass of the bound 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. (4.7) 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. [32,33], 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. [29], which uses results from refs. [65,66]. 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 with 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 well 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 figure 10, 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 fixedorder 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 figure 10 show a worse behavior. If we compute the expectation values of the relativistic potentials using the wave function obtained solving -19 -

JHEP09(2018)167
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 the matrix elements obtained solving eq. (4.1) 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 fixed-order. 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 section 3, we consider first the renormalon-free combination M h b − M η b . We plot the results of our analysis in figure 11. 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 P-wave 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 figure 12. 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. [32]). 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.  In figures 2 and 4 we have compared our numbers with alternative determinations from heavy quarkonium, either from spectroscopy or from sum rules, in the strict weakcoupling approximation. For the case of the bottom quark mass all determinations are consistent with each other. For the case of the charm quark mass our number is perfectly compatible with other determinations from spectroscopy, but there is a small tension with determinations using low n sum rules, where our number is in the low range. In any case, the difference is hardly statistically significant.
The consideration of the bottomonium, charmonium and B c together opens the possibility 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 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 M h b − M η b is compatible with experiment albeit the convergence is marginal. For M η b − M h b the situation 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.
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 . 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 renormalon-sensitive 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.

A The renormalon subtraction
We perform the renormalon subtraction in the RS and RS' schemes for all the energy levels and for different masses following the procedure in refs. [18,25]. We find that there is a non trivial mass dependence for different masses. The RS(RS') energy levels are: For different masses, following the notation in ref. [25] (see also this reference for extra details), we have where L ν = ln(nν/(2C F m r α s )) + H n+l , being H n the nth-harmonic number, andχ n are the coefficientsc n in eqs. (3.2b) and (3.2c) in ref. [25]. Note that in that reference the value of β 4 was not known and an estimate was used. Here we use the value obtained in [67,68]. This indeed changes the value of the RS/RS' masses.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.