The bottom quark mass from the $\Upsilon(1S)$ system at NNNLO

We obtain an improved determination of the normalization constant of the first infrared renormalon of the pole mass (and the singlet static potential). For $N_f=3$ it reads $N_m=0.563(26)$. Charm quark effects in the bottom quark mass determination are carefully investigated. Finally, we determine the bottom quark mass using the NNNLO perturbative expression for the $\Upsilon(1S)$ mass. We work in the renormalon subtracted scheme, which allows us to control the divergence of the perturbation series due to pole mass renormalon. Our result for the ${\overline {\rm MS}}$ mass reads ${\overline m}_{b}({\overline m}_{b})=4201(43)$ MeV.

O(mα 5 s ). The use of effective field theory methods [11][12][13][14] was important for reaching this accuracy. One of the main motivations for this ongoing effort is the possibility to obtain accurate determinations of the bottom and (may be) charm quark masses by equating these theoretical expressions to the experimental values [4,7,[15][16][17][18][19][20][21]. This possibility is reinforced by the fact that the renormalon of the pole mass cancels with the renormalon of the static potential [22][23][24] making these energies mainly perturbative objects, and, therefore, ideal candidates for good determinations of the heavy quark masses. Taking advantage of this fact requires the use of the so-called threshold masses (see, for instance, [17,[24][25][26]), which explicitly take into account the cancellation between the pole mass and the static potential renormalon.
One of these analysis was made in Ref. [17]. In this reference the NNLO expression of the heavy quarkonium ground state mass was used, as well as some partial NNNLO effects (those obtained from the large-β 0 approximation [27,28], as well as the leading O(mα 5 s ) logarithms [5]) to obtain an accurate determination of the bottom mass. The charm quark was considered to be active and its mass approximated to zero. The threshold mass used was the so-called renormalon subtracted (RS) mass. This mass is defined such that the leading renormalon of the pole mass is explicitly subtracted. Therefore, it requires the knowledge of the normalization of the renormalon, which was also approximately computed in that reference.
There is a series of developments that motivate updating such analysis. One obvious improvement would be the incorporation of the charm quark mass effects. For the case of the heavy quarkonium mass (versus the pole mass) it was concluded that the charm quark decoupled and it was a good approximation to consider the theory with only three active massless flavours [18]. In this paper we argue that one should also use this approximation for the relation between the pole and the MS mass, producing much smaller shifts than if working with four active flavours. Therefore, the calculation should be redone accordingly.
It is also possible to improve the determination of the normalization of the pole mass (and the static potential) leading renormalon, N m . On the one hand the existence of the three-loop expression of the static potential allows the determination of N m to one order higher in the corresponding expansion. On the other hand, recent analysis in lattice simulations [29] suggest that the direct determination of N m from the last known coefficients of the perturbative series may actually produce more accurate results than previous estimates, which were obtained using the Borel transform of the perturbative series as their key quantity. Such improved value would have an immediate impact in heavy quark physics in general, and in the determination of the heavy quark mass from the heavy quarkonium spectrum in particular.
With respect to the latter, the complete NNNLO correction to the perturbative expression of the Υ(1S) mass is now known. By including the complete NNNLO expression we can study this term without scheme ambiguity and assess its impact. Even more important, at this order ultrasoft effects appear for the first time. There is the worry that physics at the ultrasoft scale can not be computed in perturbation theory. The reason is that the natural scale associated to those degrees of freedom is of order mα 2 s , which, up to numerical factors, is a low scale. Yet, there have been some analysis where the ultrasoft scale has been treated in perturbation theory. For instance, in Ref [30] the perturbative expression of the static potential (which includes ultrasoft effects) was compared with lattice simulations. Also in Refs. [16,18] it has been argued that the nonperturbative effects are small for the low states of heavy quarkonium. On the other hand two existing analyses that incorporate the NNNLO expression yield bigger values of m b [7,21]. 1 In our analysis we would like to quantify the real impact of these corrections, as it is important to know what pure perturbation theory has to say before asking for non-perturbative corrections. Finally, the existence of the complete NNNLO result allows to compare with the large-β 0 estimate of the NNNLO correction, and see how reliable such approximation is.
In Sec. II we study the corrections to the pole mass and the static potential due to the charm quark. In Sec. III we present the calculation of the normalization constant (residue parameter) N m of the leading infrared renormalon: in Subsec. III A from the static quark potential V (r) and in Subsec. III B from the ratio m q /m q . In addition, in Subsec. III C we elaborate over these determinations and obtain an improved estimate of N m and the coefficients r 3 and r 4 of the perturbative expansion in α s of the mass ratio m q /m q . In 1 Ref. [7] uses an estimate for the three-loop static potential.
Sec. IV we extract the bottom quark mass from the energy of the quarkonium ground state in the RS and RS' scheme defined in Ref. [17] and later in the text. In the conclusions we summarize the results obtained.

II. CHARM EFFECTS IN THE POLE MASS AND STATIC POTENTIAL
In this section we present the perturbation expansions of the pole mass and static potential with special emphasis as to how to incorporate charm quark effects.

A. Charm quark effects in the pole mass
In this subsection we assume that we have N l massless quarks, one active massive quark with mass m c , and a (non-active) heavy quark with mass m b (such that m b > m c ). Therefore, we have a total of N f = N l + 1 active quarks. This is the situation relevant for the bottom quark (where N f = 4 and N l = 3).
The pole mass m b and the MS mass m b ≡ m b (µ = m b ) of the quark b are related by the following equality where the coefficients r n have been evaluated with N f active massless quarks. The "+" stands for the fact that a is also evaluated with N f (massless) active quarks: a + (µ) = a(µ; N f ) ≡ α s (µ; N f )/π. The coefficients R 0 , r 1 , and r 2 were obtained in Refs. [31], [32], [33,34], respectively: The value of r 3 is unknown (except for the N 3 l [35] and N 2 l [36] dependence). Therefore, the value of r 3 will be estimated, see Table II in Sec. III C where we use the notation c j ≡ β j /β 0 for j ≥ 1. β 3 was computed in Refs. [37,38].
The sum in Eq. (1) can be reexpressed in terms of a + (µ) at an arbitrary renormalization scale µ: where (r L m (µ) = ln(µ 2 /m 2 b ), and we maintain, for simplicity, the notation r j ≡ r j (m b ). Finite-mass charm effects are incorporated in which vanishes in the m c → 0 limit. The first term of the series is known [32] and reads where (see also [18]) 5 The exact expression of δm (2) (c,+) was obtained in Ref. [39] and it will be considered later. On the other hand O(a 4 + ) terms or higher are unknown. It has been noticed in Ref. [40] that δm (1) (c,+) is mainly determined by the infrared behaviour of the loop integral, which is saturated to a large extent by virtualities of order ∼ m c . In this approximation we have for m c = 1.27 GeV (to be compared with the exact result δm (1) (c,+) = 1.8058 MeV with m b = 4.2 GeV). This shows that already at this order δm (+) c is dominated by the infrared behaviour of the loop integral. We expect that this will be even more so at higher loops. On the other hand the infrared behaviour of δm (+) c can be related with the infrared behaviour of the static potential. For the static potential the charm mass dependence is known with two loop accuracy [41]. In Ref. [28] this observation was used to obtain the infrared behaviour of δm (2) (c,+) , ie. the linear behavior of δm (2) (c,+) : where b 2 = 1.12; f 2 = 0.47; f 1 = ln(A/b 2 )/ ln(f 2 /b 2 ); b 1 = ln(A/f 2 )/ ln(b 2 /f 2 ) and ln(A) = 13ζ (3) 19 + 161 228 − ln(2). The coefficients b 2 , f 2 were obtained from an approximate numerical fit to [41]. The last term in Eq. (9) comes from the fact that we are using the MS charm mass (otherwise it could be absorbed in Eq. (8)). In any case it is small compared with the rest of the coefficient. Eq. (9) is then approximated by The exact analytic expression of δm (2) (c,+) is extremely lengthy. An accurate approximated numerical form can be found in Ref. [39], which is enough for our purposes.
The difference between δm and δm (2,lin) (c,+) is due to the approximations involved in obtaining Eq. (10) (see the discussion in Ref. [39]). Therefore, we take Eq. (12) as the exact expression for the linear approximation. We observe that the linear approximation represents a quite good approximation of the exact result.
We can now compare the size of δm  Ref. [42], at high orders in perturbation theory the charm quark decouples. The reason is that at order n, the natural scale of the loop integral is me −n , which for n large enough becomes smaller than m c . Therefore, the charm mass acts as an infrared cutoff killing the low energy contributions to the integral of the fourth flavour that would otherwise produce the factorial behaviour. Thus, working with N f active flavours produces spurious contributions that deteriorate the convergence of δm c . This problem can be solved by decoupling the charm quark by expanding a + in powers of a − ≡ a(µ; N l ) (see App. A for details). The relation between the pole and the MS mass now reads where (r and we have absorbed the effects of the decoupling of S in δm c , which now reads where δm (i) (c,dec.) are generated by this decoupling and read If we put numbers we obtain δm We observe that the series is now convergent, and the strong cancellation between δm . After the cancellation, the O(a 3 ) charm-mass effect is clearly negligible compared with other uncertainties. Note also that, even though δm (2,lin) (c,+) reproduces quite well the magnitude of δm (2) (c,+) , it does not well enough to get an accurate value of the NLO correction after the cancellation. Therefore, the linear approximation could only be used to get the order of magnitude of the NLO effect (once the cancellation has been incorporated in the computation). We show this effect in Fig. 1 by comparing the exact NLO (solid line) correction with the linear approximation of the NLO correction (dotted line). Note also that the precision required is such that δm

B. Charm quark effects in the static potential
In this subsection we directly work with N l massless active quarks (as motivated by the analysis of Ref. [18]). The effects of the charm quark are included as an explicit correction to the potential.
The perturbation expansion of the QCD q-q static singlet potential is known with high accuracy. Its O(a 2 ) contribution was obtained in Ref. [1], the O(a 3 ) in Refs. [3,43], the O(a 4 ) logarithmic term in Ref. [44], the O(a 4 ) light-flavour finite piece in Ref. [8], and the O(a 4 ) pure gluonic finite piece in Refs. [9,10]. In momentum space the potential reads where L = ln(µ 2 /|k| 2 ) and µ is the renormalization scale. The coefficients a 1 , a 2 and a 3 read where the coefficients in a 3 are The terms involving powers of L in Eq.
Besides, at O(a 4 ) there is a factorization scale dependence that can not be absorbed in α s .
We have singled out this contribution, it is proportional to and depends on the infrared cutoff µ f , which cuts out ultrasoft (us) degrees of freedom The existence of the infrared divergent terms at ∼ a 4 in the static Wilson loop was first pointed out in Ref. [45].
The three-dimensional Fourier transformation of Eq. (18) gives the perturbation expansion of the static potential in position space where the notation l = ln(µr exp(γ E )) is used, with γ E being the Euler constant (γ E = 0.5772...).
The leading charm quark correction to the potential is the following Its effect will be quite tiny. Therefore, we have only incorporated Eq. (25) in our final evaluations and have not considered any other subdominant effects in the charm mass.

III. LEADING RENORMALON OF THE POLE MASS AND THE SINGLET STATIC POTENTIAL
The determination of the normalization constant of the leading infrared renormalon of the pole mass (and the singlet static potential) is an essential ingredient for the RS scheme defined in Ref. [17]. Therefore, in this section we want to improve over previous determinations of this quantity.
It is clear from the discussion of Sec. II that the charm quark decouples at large orders in perturbation theory (in practice this happens at rather low orders). Therefore, we work in the theory with N l (N l = 3 for the bottom case) active (masless flavours), and all the coefficients in this section should be understood with N f = N l .
The leading asymptotic behaviour of the perturbation series of m b (see Eq. (14)) is determined by the leading infrared renormalon ambiguity of m b . This ambiguity δm b is renormalization scale and scheme independent and is a QCD scale with the dimension of energy; therefore, it must be proportional to the QCD scale Λ QCD : δm b = const × Λ QCD [46]. This scale, written in terms of a(µ) and of the renormalization scale µ, has the form where The renormalon ambiguities give us information on the asymptotic behaviour of the perturbation expansion. This information is easily encoded in the Borel transform of S, This function has renormalon singularities at u = 1/2, 3/2, 2, . . . , −1, −2, . . . [25,47,48], and likely also at u = +1 [49]. Except for the normalization, the leading infrared renormalon ambiguity of m b , δm b = const × Λ QCD , completely determines [46] the behaviour of the nearest singularity to the origin where κ is a µ-independent constant and 'BI' denotes the Borel-integrated expression for S. 2 We then have B S (u; µ) where N m is the residue (normalization constant) parameter of the renormalon.c 1 was first computed in Ref. [46], andc 2 in Refs. [17,48]. We also give a value for c 3 using the estimate for the MS scheme coefficient c 4 = β 4 /β 0 obtained in Ref. [51] by Padé-related methods The coefficients h N (µ) of the analytic part (31) where we recall that r 0 = c 0 = 1. The sum in Eq. (32) introduces O(1/N ) corrections to the leading asymptotic behaviour. The numbers c s entering the sum in Eq. (32), are given by Eqs. (27) and are known for s ≤ 3. 3 Therefore, by default, we truncate the sum in (32) at s = 3. This introduces an error of order O(1/N 4 ) for the asymptotic behaviour (we will typically take the difference between truncating the sum at s = 2 or s = 3 to check the quality of the approximation). We also set h N = 0 since they yield (in comparison) exponentially suppressed terms. Overall, we approximate the asymptotic behaviour of r N by the following equality: We can have a similar discussion for the static potential. The Borel transformation of the dimensionless potential (− 3 4π )rV (r) in (24) is given by: where v j is the coefficient at the power a j (µ) in the expansion (24) (v 1 = a 1 /4 + 2β 0 l, etc.).
The asymptotic behaviour of v N is equal to the behaviour of r N except for the normalization: The coefficients d N are analogous to the coefficients h N for the pole mass. We will set them equal to zero for the same reason, as they yield exponentially suppressed terms in comparison. We also truncate the sum to the first known terms (s ≤ 3). Therefore, we approximate the asymptotic behaviour of v N by the following equality:

A. Determination of N V
In this subsection we compute N V using two methods that we name A) and B).
The method A) uses the idea of Refs. [53,54] of, instead of working with Eq. (35), using an associated function that kills the leading singularity in the Borel plane. This idea was first applied to the static potential (and the pole mass) in Ref. [17] and also used in [30,55,56].
One uses which is defined such that the leading singularity at u = 1/2 of B V is eliminated. Its evaluation at u = 1/2 gives In this paper we carefully study the case when the sum is truncated at N = 3. Truncating the infinity sum to its first orders produces some remaining scale dependence. In particular, from N = 3 on a dependence on the ultrasoft factorization scale µ f appears, which we set equal to µ. We plot the scale dependence of −N V /2 (the relevant quantity to be compared with N m ) in Fig. 3 for N = 0, 1, 2, 3 for N l = 3. We observe a nicely convergent pattern, specially for the difference between the N = 2 and N = 3 computation. Note also that we expect the results to be better for µr ∼ 1, since ln(µr) terms are not large. If we compare method A) and B), we observe that method B) yields a more convergent series and a milder scale dependence for N = 3 than method A). Therefore, we fix the central As expected, determinations with x ≡ µr ∼ 1 yield the best results, since this minimizes possible large ln(µr) terms.
The factorization scale µ f is not related with the leading infrared renormalon. Eliminating this contribution altogether allows us to measure the quality of our error estimate. We plot the determination of N V if we completely eliminate the ultrasoft term in Eq. (24) in Figs. 4 and 5. We observe that the associated shift is much smaller than the scale variation of the NNNLO curve, specially for the N l = 0 case.
The fact that method B) yields a more convergent (and stable) series was also clearly observed in Ref. [29], where N m was determined from the perturbative computation of the self-energy of a static source to O(α 20 s ) (compare Fig. 12 with Fig. 14 in this reference).  (40). We also plot the NNNLO curve without the US term (thin solid line).

B. Determination of N m from the pole mass
In this subsection we obtain N m from the perturbation expansion of the pole mass using the methods A) and B) described in the previous section.

From method A) one obtains N m from
where Truncating the sum to N = 2, we recover the results of Refs. [17,19,56]. Unlike in the previous section, we can not go to one order higher since the O(a 4 ) term of the pole-MS mass relation is not known. Therefore, we do not dwell further on this method.
With method B) we determine N m by dividing the exactly known coefficients r N (cf. Eq.  Therefore, the predictions in this section are generically less precise.

C. Final determination of N m
In the situation where 1/r Λ QCD , one can do the matching between NRQCD and pNRQCD in perturbation theory, and 2m q + V (r) can be understood as an observable up to O(r 2 Λ 3 QCD , Λ 2 QCD /m) renormalon (and/or nonperturbative) contributions (see the discussion in Ref. [17]). This implies that the leading infrared renormalon of the singlet static potencial must cancel with the leading renormalon of twice the pole mass, so that the following relation between N m and N V holds:   In Sec. III B we have determined N m using what we have named method A) and B). The results from method A) are nothing but those obtained in Ref. [17], where the estimated uncertainty was of around 10%. The results from method B) are new, and summarized in Fig. 6. Both methods use the perturbation expansion of the pole mass to O(a 3 ), which is 4 One may think of other observables, the perturbation expansions of which are dominated by the pole mass renormalon. One of those is the anomalous magnetic moment of the heavy quark. In Ref. [57] it was shown (see also Ref. [58]) that, (1 + κ)/m is renormalon free. More precisely, the leading renormalon of 1 + κ cancels with the leading renormalon of m and one can write 1+κ s is free of the u = 1/2 renormalon. Therefore, we are in the situation wherem(1 + κ) = m(1 + C κ ) has the pole mass (renormalon) but modulated by a nontrivial Wilson coefficient. This changes the 1/N corrections of the asymptotic expression. We have studied this quantity and obtained a number for the normalization of the renormalon compatible with ours though less precise. one order less than for the static potential. Actually, these fits typically yield a stronger scale dependence than for the case of the static potential. Therefore, we will not dwell further with these determinations.
In Sec. III A, we have obtained two new determinations of N V using the perturbation expansion of the static potential to O(a 4 ), one order more than for the case of the pole mass. Therefore, we expect our new determinations to yield more accurate results. We have found this is specially so for method B of Sec. III A. Since N V and N m are related by Eq. (43), this gives a determination of N m , which we take as the most precise and display in Table I as our final numbers. To illustrate this, in Figs. 7 and 8, we compare the NNNLO evaluations using method B) (of the static potential) with method A) (of the static potential). We also include the NNLO evaluations using methods A) and B) (of the pole mass). Around x ∼ 1 all of them agree within one standard deviation. This signals that alternative errors estimates would give similar numbers. We see how the NNNLO evaluation using method B yields the more stable result under scale variations. For the N l = 0 case we can compare with the value obtained in Ref. [59]. We agree within one standard deviation. This is quite rewarding as these numbers have been obtained with completely different methods.
It is also interesting to study the N l dependence of N m . The u = 1/2 infrared renormalon should disappear when N l → ∞, as the theory is not asymptotically free anymore. We plot N m as a function of N l (we fix x = 1) using our preferred method (method B) from the static potential in Fig. 9. We observe how N m tends to zero in the range of N l ∈ (12, 23), a range of values for which one could expect a conformal window. This shows the disappearance of u = 1/2 infrared renormalon for N l > 12. In the range N l ∈ (25, 40) the evaluation of N m with method B) is unstable because the asymptotic expression of r 3 : r asym 3 (see Eq. 33) have a couple of zeros in this range (due to the beta coefficients, therefore it is very sensitive to subleading corrections and the truncation), which produces divergences for the theoretical expression that we use to determine N m . This is nothing but the reflection of the fact that we are in a transition region before we reach the behavior expected for N l → ∞. In this limit we expect, not only the disappearance of the u = 1/2 infrared renormalon, but its transformation into a ultraviolet u = −1/2 renormalon (so that the perturbative series is sign alternating), for which the normalization can be computed in the large N l limit [47]: N We have also done some fit-play of the N l dependence of N m for small N l using a polynomial function: N m (N l ) = N m (0)+d 1 N l +d 2 N 2 l +· · · . We observe that subsequent coefficients d n get smaller as we increase n for small N l , with the leading coefficient, d 1 , of order ∼ −10 −2 (see also Table I).
Since we have a reliable determination of N m , we can obtain a prediction for the high order coefficients of the perturbation expansions of the pole mass and static potential. We explicitly show them in Table II for  We now compare our predictions with earlier estimates. The quality of large-β 0 predictions is worse [35]. This is to be expected as they do not incorporate the right large-N    (29) asymptotic behaviour. Dispersion-like analyses [60] also seem to have problems to capture the right asymptotic as they yield significantly smaller numbers than the ones obtained here.
In Ref. [17] (see also Ref. [55] for N l = 0) the NNLO prediction from method B) was used, which is around one sigma away from our new number. In App. C of Ref. [21], a variant of this method using Padé approximant was worked out and the number obtained was quite similar. There has also been a recent prediction of r 3 , made in Ref. [61] by demanding stability of the perturbation expansion of the heavy quarkonium energy (in the static limit) to the next order. Our numbers are bigger than his for small N l . Note though that for large N l , we get similar numbers. This points to a different N l dependence, which in our case is more pronounced. Finally, for N l = 0 we can also compare with Eq. (13) of Ref. [59]. Their value for r 3 is in agreement with ours within one standard deviation. This is quite remarkable as that method is completely different, based on lattice simulations, and, therefore, with where m is the pole mass (m q ), and we use the notations of Eqs. (27)- (29). Therefore, we have the following explicit expression for m RS : where δm RS is the residual mass (we recall that c 0 = 1): Note that we work in the theory with three active flavours only, as the charm decouples at large orders in perturbation theory (which is the regime δm RS deals with).
Equation (45) is still formal. In practice, one rewrites m in terms of m using Eq. (1) and reexpands the perturbation series in Eq. (46) around the same coupling a − (µ), at fixed but otherwise arbitrary scale µ: and this leads to a relation analogous to Eqs. (47) m where h N (ν f ; µ) in Eq. (49b) are obtained by expanding a − (ν f ) in Eq. (49a) in powers of a − (µ). The explicit relation between the two Borel transforms is The perturbation expansion of the Υ(1S) mass is presently known up to O(m b a 5 ) where µ is the renormalization scale, m b is the pole mass of the bottom quark and The numerical expressions of the coefficients K i,j (N f ) and K 3,0,j are given for reference in Appendix B.
As we have already discussed throughout the paper, it is compulsory to implement the cancellation of the leading infrared renormalon (u = 1/2) in the above perturbation series to get a convergent series. We do so by working in the RS scheme. In practice this means to rewrite m b in terms of m b,RS in Eq. (51). The resulting expression reads M Υ(1S) : where we denoted In the expression (53) for M Υ(1S) , the terms of the same order (ν f /m b )a n and a n+1 were combined in common brackets [. . .], in order to account for the renormalon cancellation.
If using the RS' mass in our approach instead, the above expressions are valid without changes, except that m RS → m RS and K 0 → 0 (and: h 0 (µ) → 4/3).
We note that we take N l = 3 active flavours and that the expression above does not incorporate yet the charm quark effects. The leading one is due to the potential Eq. (25) and reads (see, for instance, [62]) whereρ = 3m c /(2m b πa − ). It produces a shift of order 1 MeV, completely negligible in comparison with other uncertainties. This has also been stressed in Ref. [18]. Nevertheless, it got obscured because specific numbers were given with N f = 4.
We now investigate the dependence of our results on the theoretical and experimental parameters. To estimate the errors, we vary µ, ν f , α s and N m as follows: µ = 2.5 +1.5 −1 GeV, ν f = 2 ± 1 GeV, α s (M z ) = 0.1184(7) [63] (with decoupling at 4.2 GeV and 1.27 GeV for the bottom and charm MS masses, respectively) and N m = 0.563 (26). For the RS scheme, we obtain the following result 56 For the RS' scheme, we obtain the result (with the same variation of the parameters) m RS (2 GeV) = (4201 + 189 + 36 where the -0 of the last equality is accidental for the specific scale chosen (see Fig. 10). 5 Here and in the following, in the determination of m MS (m MS ) ≡ m(m) ≡ m, we have used our estimate of the four-loop relation. 6 Note that the scale dependence of m is the one associated to the fit to m RS . In the RS' we obtain m RS (2 GeV) = (4206 + 476 + 60 + 18 + 1) MeV .
We observe a nicely convergent perturbative series in the relation between the RS(RS') masses and the MS mass. In the perturbative relation between the Υ(1S) mass and the RS(RS') masses we also observe a convergent series except when we consider the difference between the NNLO and NNNLO results. In Figs. 10 [5]). Specially problematic is the appearance of the ultrasoft scale, since it is potentially a rather small scale. Whether this scale can be treated within perturbation theory can only be elucidated by higher order computations, as well as by analyses using renormalization group techniques. Without such analyses it is not possible to unambiguously set the scale of the NNNLO contribution, and any estimate of higher order effects due to the ultrasoft contributions will suffer from some scheme dependence. Nevertheless, we can observe some general trends. The ultrasoft logarithmic dominated terms [5,64] yield a positive contribution to the NNNLO term. This contribution would be even bigger if we set the scale of one of the powers of α s at the ultrasoft scale (as the effective field theory suggests). This would go in the direction of making the NNNLO result (and also the value of the bottom mass) smaller and improve convergence. Nevertheless, without a renormalization group analysis we can not make this discussion more quantitative. Finally, without a better control of the ultrasoft effects it would be premature to consider nonperturbative corrections, which we neglect in this analysis. As for the large-β 0 approximation we observe that they give numbers in the right ballpark, yet one should keep in mind that this comparison will depend on the constant term in the ultrasoft logarithm. 7 To this analysis one has to add charm effects. The leading correction to the heavy quarkonium spectrum can be found in Eq. (55). It produces a negligible correction ∼ −0.6 MeV. Therefore, it barely changes our determination above. The corrections to the relation between m RS and m, Eqs. shows a better convergence, as expected.

V. CONCLUSIONS
In this paper we have considered different improvements over the analysis made in Ref. [17].
First, we have studied whether (or when) the charm quark decouples in the perturbative relation between the pole and the MS mass. For the case of the heavy quarkonium mass (versus the pole mass) it was seen that the charm quark decoupled and it was a good approximation to consider the theory with only three active massless flavours [18]. In this eventual agreement will deteriorate at higher orders. = 0.563 (26) .
This has an immediate impact in the determination of the heavy quark mass from the heavy quarkonium ground state mass but it is also applicable to other observables in heavy quark physics. This improvement is twofold. On the one hand the existence of the three-loop expression of the static potential allows us to determine the normalization to one order higher in the corresponding expansion. On the other hand, we obtain the normalization directly from the last known coefficient of the perturbation expansion. This leads to a more stable result compared with previous approaches, as it has already been observed in lattice simulations [29] for the case of the self-energy of a static quark.
using the renormalon subraction scheme. In this analysis we have worked with three active flavours, as motivated by the previous discussion. We have also used the updated value of N m obtained in Eq. (66). By including the complete NNNLO expression we can study this term without scheme ambiguity. At this order ultrasoft effects appear for the first time.
Consistent with this fact we observe that the NNNLO result is more scale dependent than the NNLO one, and the convergence of the perturbative series deteriorates somehow. It remains to be analyzed whether a renormalization group analysis could reduce the scale dependence of the result and improve the convergence. Yet, the magnitude of the NNNLO 8 As they suffer from different systematics, we can consider combining this result with N m N l =0 = 0.620 (35) [59] and obtain an even more accurate value: N m N l =0 = 0.608 (22).
correction is small, so that the final number is quite close to the value obtained by the NNLO analysis made in Ref. [17]. On the other hand, our determination is considerably smaller than the NNNLO determinations in Refs. [7,21]. [7] worked with N f = 4 and an estimate for a 3 but did not implement explicitly the cancellation of the renormalon.
In Ref. [21] the MS scheme and another related scheme were used, the calculation was performed at NNNLO for the Υ(1S) mass with N f = 3 but N f = 4 was used in the relation for m/m, and a lower soft renormalization scale µ ≈ 2 GeV was used; the (strong) m c effects were not under control at NNNLO (due to the use of N f = 4), besides producing a small mismatch in the renormalon cancellation. Our number is somewhat in the middle between two recent determinations from bottomonium NR sum rules [65,66]. Both of them worked with N f = 4 active massless quarks, though the latter reference estimated the shift produced by the finite mass charm quark effects. [65] performed a partial NNLL computation. The difference with the (also partial) NNLL determination in Ref. [67] stems from extra/different terms incorporated in the analysis. [66] performs a partial NNNLO computation. Note also that our number is similar to the number obtained from a lattice HQET determination [68], and not very far away from the numbers obtained using low-n or finite-energy bottomonium sum rules [69,70], or a lattice determination using NRQCD [71]. Those have complete different systematics. This agreement may indicate that nonperturbative corrections are indeed small, as advocated in Refs. [16,18].
where the coefficients z 1 account for the N f = 4 → 3 quark threshold effects and the