The bottom-quark mass from non-relativistic sum rules at NNNLO

We determine the mass of the bottom quark from high moments of the bottom production cross section in e+ e- annihilation, which are dominated by the threshold region. On the theory side next-to-next-to-next-to-leading order (NNNLO) calculations both for the resonances and the continuum cross section are used for the first time. We find mb,PS(2 GeV) = 4.532 +0.013 -0.039 GeV for the potential-subtracted mass and mb,MSbar(mb,MSbar) = 4.193 +0.022 -0.035 GeV for the MSbar bottom-quark mass.


Introduction
The bottom-quark mass is one of only a handful of fundamental QCD parameters, and thus its precise knowledge is of considerable interest by itself. Furthermore there are also phenomenological applications which benefit from a small uncertainty in the value of the bottom-quark mass. Examples include Higgs decays to bottom quarks and decays of B mesons.
Sum rules provide a well-established method for the determination of heavy-quark masses [1,2]. Considering the normalised bottom production cross section where s is the square of the e + e − center-of-mass energy, we can define its moments M n as .
(1.2) Π b denotes the contribution from bottom quarks to the vacuum polarisation function. C is an arbitrary closed contour that encloses the origin and does not cross the branch cut, i.e. the part of the positive real half-axis where R b (s) = 12π Im Π b (s+iǫ) > 0. Quark-hadron duality now permits to determine the bottom-quark mass by equating moments obtained from the measured hadronic cross section with moments calculated from derivatives of the theoretically predicted polarisation function, which can themselves be obtained from the theoretically predicted partonic cross section by the above dispersion relation (1.2).
For small values of n it is most convenient to calculate moments directly from the derivatives of the polarisation function in conventional perturbation theory. As n increases the theory uncertainty grows and for n 10 the perturbation series shows no signs of convergence anymore (cf. [3,4]).
For larger values of n the integral over R b in equation (1.2) is increasingly dominated by small kinetic energies E = √ s − 2m b ∼ m b /n [5,6], where m b is the bottom-quark mass. Thus, the quarks are non-relativistic with a small velocity v ∼ 1/ √ n ≪ 1. At the same time the strong coupling constant α s is of the same order as the quark velocity. For this reason terms scaling with powers of α s /v, which originate from Coulomb interaction, have to be summed to all orders. This is achieved in the effective theory of potential nonrelativistic QCD (PNRQCD) [7]. In summary, for large-n, "non-relativistic" moments the power counting for the cross section up to NNNLO is given by (1. 3) The admissible values of n are limited by the requirement that the ultrasoft scale m b /n remains larger than the intrinsic strong interaction scale of a few hundred MeV [8].
In this work we present the first determination of the bottom-quark mass from large-n sum rules at next-to-next-to-next-to-leading order (NNNLO). The corresponding NNLO calculations have been done about fifteen years ago [8][9][10]. A partial NNNLO result that uses NNNLO accuracy for the bound-state contribution and NNLO accuracy for the continuum contribution to the moments appeared recently [11]. The outline of this paper is as follows. First, in section 2, we derive the values of the relevant moments from experimental data. We proceed in section 3 with an outline of the corresponding PNRQCD calculation and a discussion of suitable mass schemes. In section 4 we summarise our results for R b and for the bottom-quark mass and compare them to other recent high-precision sum rule analyses. For the comparison we consider the works by Chetyrkin et al. [12], Hoang et al. [13], and Penin and Zerf [11]. We conclude in section 5. Since we are aiming at high precision, we include the effects of the non-zero charm-quark mass. The details of this computation are given in four appendices.

Experimental moments
For sufficiently large values of n the experimental moments M exp n are dominated by Υ bound states. In the narrow-width approximation for the resonances Υ(1S) to Υ(4S) we obtain where we take the current PDG values [14] for the masses M Υ(N S) and the leptonic widths Γ Υ(N S)→l + l − of the Υ resonances. In the energy region of interest we approximate the running QED coupling by a constant value, which is given in terms of the fine structure constant by α(2m b ) ≈ 1.036 α [15]. For the remaining continuum integral we use experimental data [16] corrected for initial state radiation between √ s cont = 10.6178 GeV and 11.2062 GeV (see [12] for details). In the absence of data for higher center-of-mass energies we assume R b to stay roughly constant with R b = 0.3 ± 0.2. For n = 6, 10, 15 the unknown high-energy part constitutes approximately 53%, 35%, 21% of the total continuum contribution, respectively. The large uncertainty of this part is therefore not expected to be important. The resulting values and uncertainties for the experimental moments are shown in table 1.

Theory moments 3.1 The cross section in PNRQCD
The normalised cross section for bottom production at NNNLO in PNRQCD has the form where j i are the spatial components of the relativistic currentbγ µ b, and ψ (χ) the non-relativistic quark (antiquark) spinors. Finally, G(E) is the non-relativistic current correlator Below threshold, i.e. for E < 0, G(E) develops infinitely many poles which can be interpreted as S-wave, spin-1 bb bound states. In the vicinity of such a bound state G(E) behaves as where ψ N (0) is the bound state wave function at the origin while E N is the binding energy. Splitting off the bound-state contribution we can thus write the moments as with the residues Note that Z N is a physical quantity, while the wave function ψ N (0) and the matching coefficients separately depend on the factorisation scale and scheme that separates shortand long-distances. R b (3.1) and Z N (3.6) are understood to be strictly expanded up to NNNLO according to the PNRQCD power counting (1.3). For example, terms of order α 4 s and higher that are generated through the product of lower-order terms are to be dropped. Contrarily, we will leave the factor (2m b + E N ) −(2n+1) in equation (3.5) unexpanded for reasons discussed in section 3.3. Not expanding this factor is equivalent to resumming the poles back into the polarisation function (cf. [18]), i.e. replacing in the contour integral that defines the moments (1.2).

Technical implementation
We require the following ingredients for the NNNLO analysis of non-relativistic moments: • The hard Wilson coefficient c v at order α 3 s [19] and d v at order α s [20].
For an extensive overview of all third-order corrections see [17,25]. Note that the results in [25,33] are calculated for the case of top quarks, that is in the complex energy plane for finite imaginary part Γ corresponding to the top-quark width. The application of these results to bottom quarks requires to take the limit Γ → 0, which is non-trivial, since the analytic expressions cannot be straightforwardly evaluated numerically for vanishing width. In [33] the Γ = 0 result has been constructed for the ultrasoft contribution by extrapolation. Evaluating the third-order potential corrections given in [25] for real energy E requires a substantial amount of extra work, which we briefly discuss in the following.
The higher-order potential corrections to the Green function are expressed in terms of nested harmonic sums, sums over gamma and polygamma functions, and generalised hypergeometric functions [23,25]. The complex variable λ = (α s C F /2) −m/(E + iΓ) appears for example in the argument of (poly-)gamma functions or as one of the parameters of the hypergeometric functions. In our application, we have to evaluate the Green function for positive values of the energy E, starting at E = 0. For vanishing width Γ, λ tends to +i∞ as E tends to zero. Thus, we have to ensure that all expressions are well-defined in this limit and that their numerical evaluation is possible. In most cases it is possible to express the correction to the Green function in terms of harmonic sums. This is in particular the case for most of the generalised hypergeometric functions, which can be treated as described in appendix A.1 of [34]. The harmonic sums can then be analytically continued with the methods of [35,36]. However, in some cases the correction to the Green function is expressed in terms of single or even double sums which could not be expressed as nested harmonic sums. For such sums it was often necessary to truncate the summation at some (λ-dependent) value and construct suitable asymptotic expansions to approximate the remainder. In all cases we have checked that the numerical precision is sufficient for our extraction of the bottom-quark mass, such that the numerical uncertainty can be neglected.
A relatively simple example of this procedure is given by the sum which appears in the insertion of the Darwin term. Here ψ is the logarithmic derivative of the gamma function and ψ (1) the first derivative of ψ. The sum converges only slowly when λ is large, which makes the numerical evaluation difficult. Therefore, we introduce a cut-off Λ for the summation and explicitly sum all terms up to this cut-off. Choosing Λ to be much larger than |λ|, we can approximate the remainder by expanding the summand in the limit k → ∞. Note that for ψ(k − λ) this is not simply an expansion in |λ|/k, but rather a double expansion for k ≫ 1 and k ≫ |λ|. In the first step the entire argument of the ψ function is considered large and the terms of the resulting asymptotic series are further expanded for |λ|/k → 0 in the second step, yielding (3.9) After expanding the summand in (3.8), the sum over k from Λ + 1 to infinity can be evaluated in terms of Hurwitz zeta functions (in more complicated cases we also encounter derivatives of this function). The first three terms are where ζ(s, a) = ∞ k=0 (k + a) −s . Due to the nature of the double expansion of ψ(k − λ), higher order terms in this series can have the same powers of λ as lower order ones. However, these terms are still suppressed by additional powers of Λ, since the first argument of ζ increases with each term and ζ(s, Λ) with s > 1 behaves like 1/Λ s−1 as Λ tends to infinity. In practice, we fix for each sum an appropriate expansion depth N max for the remainder. This makes it possible to speed-up the calculation by pre-computing the expansion for arbitrary values of the cut-off. The latter is then determined at runtime as Λ = max(Λ 0 , f |λ|), where Λ 0 is a lower limit for the cut-off and f is a positive integer. Λ 0 and f are again chosen separately for each sum. In the above example we use N max = 15, Λ 0 = 100, and f = 2. By varying N max and f , we can check the numerical stability of the procedure.

Mass schemes
It is well-known that for quark masses the pole scheme is ambiguous due to its sensitivity to infrared renormalons [37][38][39]. Numerous short-distance mass schemes have been developed to cure this shortcoming [40][41][42][43]. In this work we consider the potential-subtracted (PS) mass introduced in [41] and the mass in the MS scheme. They are related to the pole mass via perturbative series of the form For the cancellation of leading infrared renormalons we have to take into account the first correction term δm M 0 already at leading order and one additional correction term for each further order in the PNRQCD expansion.
In the PS scheme the leading subtraction term is given by where we choose the subtraction scale as µ f = 2 GeV. The choice of the renormalisation scale µ is discussed in section 4.1. The higher-order terms up to δm PS 3 required for our NNNLO analysis can be found in [23].
In the MS scheme the higher-order corrections in relation (3.11) are only known up to i = 2 [44,45]. 1 It is, however, expected that the correction terms δm MS i are dominated by the leading infrared renormalon already at relatively low orders [46,47]. On this basis an approximation was constructed in [43]. The deviation from the known result is about 10% for i = 1 and less than 1% for i = 2. We employ this approximation for i = 3, which corresponds to setting the correction term δm MS s in the notation of [43]. To estimate the resulting uncertainty we vary the value ofr 3 by 10%. This range encompasses the valuer 3 = 13.5972 obtained from the large-β 0 approximation [46,47].
When eliminating the pole mass in the cross section formulae in favour of a shortdistance mass the question arises whether the resulting expression should be expanded in the correction terms δm M i , i ≥ 1. For the MS scheme such an expansion is not sensible in the threshold region [41]. This can for instance be seen by considering the factor 1 Note that in our notation δm MS i denotes the (i + 1)-loop correction.
(2m b + E N ) −(2n+1) in the resonance contribution to the theory moments (3.5). On the one hand, the correction of order i to the energy levels E N is parametrically of the same order as δm MS i+1 according to the PNRQCD power counting α s ∼ v ≪ 1. On the other hand, renormalon cancellation only occurs between the correction of order i to the binding energies and δm MS i , so an expansion consistent with our power counting would spoil this cancellation.
In the PS scheme both expanding and not expanding the PS-pole mass relation (3.11) appear to be viable options. However, as discussed in [25], an expansion induces unphysical behaviour near threshold in the continuum cross section. Therefore we will not expand the cross section in δm PS i , which corresponds to the PS-shift (PSS) prescription of [25]. At the same time this again implies that the entire factor (2m b +E N ) −(2n+1) must remain unexpanded to ensure the cancellation of leading infrared renormalons. In practice, this means that in both schemes, for given m M b and order N k LO, we first compute m b from (3.11) truncating the sum at i = k and add the bound state energy evaluated to the same order.

Expansion in the kinetic energy
Up to now, in (3.1), (3.5) and (3.6) an overall factor of 1/s stemming from the relativistic polarisation function has been expanded for E ≪ m b . Since we chose not to expand the factor 1/s n+1 in the definition (1.2) of the moments, we may also contemplate to keep this factor unexpanded. The corresponding expressions, denoted by a tilde, take the following forms: Formally, the difference to (3.1), (3.5) and (3.6) is of higher order (NNNNLO). These higher-order corrections are numerically non-negligible, possibly due to subleading renormalon contributions. For the Nth resonance 1/s assumes the form (2m b + E N ) −2 , where renormalon contributions cancel between the binding energy and the mass. If we now expand this factor for E N ≪ m b and, according to the PNRQCD power counting, discard contributions to E N that are beyond NLO the renormalon cancellation will again be spoilt. Since the exponent now is not of order n as was the case for the 1/s n+1 factor, the generated renormalon ambiguity is only of order Λ QCD /m b and therefore subleading. Since this is not the only source of sub-leading renormalon contributions, in order to decide which prescription for the expansion of 1/s is the preferred one it would be necessary to analyse carefully how all of these contributions can be cancelled in the determination of the moments. To our knowledge such an analysis has not been performed yet. We will find in section 4.3 that not expanding the factor 1/s in the polarisation function improves the consistency of mass values extracted from different moments. In the following, we will therefore mainly concentrate on the "unexpanded" moments defined by (3.15). We then check that the difference to the "expanded" approach is compatible with our estimate of the perturbative uncertainty.

Charm-quark mass effects
The ingredients for the NNNLO cross section described in section 3.1 all assume the light quarks to be massless, which is well justified as long as the light-quark masses are smaller than the physical scales that appear in the problem. However the actual charmquark mass is of the order of the soft scale (inverse Bohr radius) m b α s ∼ m b v of the bb system near threshold and therefore has to be included in a consistent treatment. At NNLO the effects of a non-zero charm-quark mass on the moments have been discussed thoroughly in [48], where they were found to lead to a sizeable shift of around −30 MeV in the extracted MS bottom-quark mass.
Since the factor 1/s n+1 in the moments is expanded non-relativistically in [48] the results cannot be included in our expressions (3.5) and (3.15). In the following we discuss in which steps in the computation of the cross section the effects of a non-zero charm-quark mass are relevant and determine the missing contributions up to NNLO. The NNNLO charm-mass effects are unknown at the time of this writing and their determination is clearly outside the scope of this work.
The computation of the cross section proceeds in three separate steps, which are the hard matching, the soft matching and the computation of the spectral function in the resulting effective theory PNRQCD, as was discussed in detail in [17]. Since m c is considered to be soft the results of the hard matching must be analytic in m 2 c /m 2 b ∼ α 2 s and charm-mass effects due to the insertion of a charm loop into a gluon line scale as α 2 s m 2 c /m 2 b ∼ α 4 s compared to α s for the gluon line without charm-loop insertion. To NNNLO the hard matching procedure is therefore unaffected.
In the soft matching procedure the charm quark is integrated out. The only part of the resulting PNRQCD Lagrangian that is affected at NNLO is the Coulomb potential, which can be split in two parts where V mc is defined such that it vanishes for m c → 0. 2 It has been computed to two loops in [49,50]; convenient representations are given in [48]. In order to offer a selfcontained discussion we quote these results in appendix A. As discussed in [17] the soft matching of the external vector current with massless light quarks is trivial to all orders, because the respective integrals are scaleless. The introduction of the charm-quark mass as a soft scale means that starting at two loops this is no longer true. However the corresponding matching coefficient must be trivial in the limit m c → 0 and since the only other scale relevant for the external current is m b , the correction must contain at least one power of m c /m b , which again is beyond NNLO. 3 The spectral function receives contributions from the charm-quark mass through insertions of the Coulomb potential. At NLO only the single insertion of V (1) mc contribute. Our results for the charm corrections to the energy levels E N and the wave functions |ψ N (0)| 2 can be found in appendix C. We find numerical agreement in those parts that are available in the literature [48,51].
Since, for the reasons discussed in section 3.3, we use the PS or MS mass instead of the pole mass in the cross section as well as the moments, the charm-mass effects also have to be considered in the relation between different mass schemes. In the conversion between the pole and the MS scheme, these effects are known at order α 3 s [52]. For our analysis we have also computed the corresponding corrections to the relation between the PS and the pole mass. These results are summarised in appendix B.
Anticipating our numerical analysis in section 4 let us now discuss the impact of a non-zero charm-quark mass. We parametrise the corrections in terms of the MS charmquark mass at our overall renormalisation scale µ, which we obtain via 4-loop running from the initial value m MS c (3 GeV) = 0.986 GeV [12]. 4 It should be noted that effects from a non-zero charm-quark mass are expected to be large for quantities that have large infrared sensitivity. The reason for this is that the charm-loop correction effectively acts as an infrared cut-off on the virtuality of the gluon line into which it is inserted. Table 2 illustrates the numerical impact of a non-zero charm-quark mass on the first energy level and the corresponding wave function at the origin and demonstrates the cancellations in the transition to a short-distance mass scheme. Let us now discuss the effect on the final value of the MS bottom-quark mass extracted from the tenth moment. If, at first, we only consider the corrections from single insertions to the binding energies and the wave functions at NNLO (NLO) the resulting mass value receives a shift of +16 MeV (+4 MeV). In line with our previous discussion we expect a considerable compensation from the charm effects in the relation between the infrared-sensitive pole mass and the PS mass. In fact, including also these corrections the mass shift reduces to around +0.5 MeV (−0.5 MeV). The cancellation of infrared contributions can also be observed analytically. Expanding all corrections in the limit of a small charm-quark mass, the linear terms in the expansion cancel exactly between the corrections to the energy levels and wave functions on the one side and the relation between PS and pole mass on the other side (cf. [48]). The charm mass corrections to the relation between the PS and the MS scheme are relatively small; adding them to the analysis leads to a total shift of −3.5 MeV (−2 MeV) in the MS bottom-quark mass.
9.50 23.21 -6.82 -19.50 0.0335 0.0600 Table 2: Contributions of a non-zero charm-quark mass to the binding energy E 1,mc and wave function of the Υ(1S) resonance. The corrections to the wave function are given by |ψ , and f once a short-distance mass scheme is used, we also show results for the charm-mass effects in the relation between the PS and pole mass with Up to now, we have neglected charm corrections to double insertions and to the continuum cross section. In agreement with [48] we find that the charm contributions from double insertions are suppressed compared to the single insertions, causing an additional mass shift of only +0.5 MeV. Since the overall continuum contribution to the tenth moment is already small we expect charm effects in the cross section above threshold to be negligible. Indeed we find an extra shift of less than 0.1 MeV from the NLO corrections listed in appendix D.
In the light of these findings we will use the computationally expensive charm corrections to the continuum cross section and the double insertions only to extract our central value and neglect them in the error analysis. To account for the unknown NNNLO charm mass corrections we will assign an uncertainty of 100% to the known NNLO effects.
Concerning the size of the charm corrections there is a large discrepancy between the final mass shift of −3 MeV in our analysis and about −30 MeV in [48]. We find that differences between the two approaches, such as the choice of renormalisation scales, expansion of factors 1/s, and different values for the charm-quark mass, cannot account for this disparity. It seems suspicious that [48] claims large effects of around 50% for n = 20. Such high moments are dominated by the first resonance. By the definition of the 1S mass, however, charm effects cancel completely in the combination 2m Thus, only the charm corrections to the residue Z 1 contribute. These effects are independent of n, and, according to our findings, rather moderate in size (cf. table 2). 4 Numerical analysis 4

.1 Choice of the renormalisation scale
Since the moments receive contributions from several distinct physical regions it is important to choose an adequate value for the overall renormalisation scale µ. A priori "natural" options include the hard scale µ ∼ m b , the soft scale µ ∼ 2m b / √ n and the ultrasoft scale µ ∼ m b /n. In this work we do not consider the possibility of summing logarithmic effects ln √ n related to the presence of different scales by renormalisation group methods (see [13] in the present context), since ln √ 10 is not a particularly large number. Figure 1 shows the scale dependence of the tenth moment for µ between 2 GeV and 10 GeV at different orders. There is clearly no convergence of the perturbation series for µ 3 GeV. 5 We therefore adopt µ = m PS b as our central scale and vary µ between 3 GeV and 10 GeV to estimate the uncertainty. Note that there is no overlap in the moment values at NLO and NNLO over the entire scale range, which is one of the motivations to perform the third-order calculation. The figure shows that the NNNLO curve lies within the interval determined by the NNLO scale variation.

Comparison to the fixed-order continuum
The present work is the first to include the continuum cross section with NNNLO accuracy (heavily relying on the input from [25,33] as described in section 3.2), which is also the most complicated part of the NNNLO calculation. Since sum rule determinations of the bottom quark mass rely on quark-hadron duality, the continuum is conceptually a crucial ingredient in the analysis.
The summation of factors α s /v implicit in the PNRQCD calculation of the correlator G(E) is only necessary close to threshold. It is therefore reasonable to expect that there is a region with α s ≪ v ≪ 1 where the continuum cross section can be calculated reliably both in PNRQCD (requiring v ≪ 1) and conventional perturbation theory (requiring v ≫ α s ). In figure 2 we show the respective predictions adopting the pole-mass scheme with m b = 5 GeV, and without expanding the factor 1/s in the polarisation function (see section 3.4). For the fixed-order curves we have used the analytically known result up to order α s [54] and Padé approximation [55][56][57] at orders α 2 s and α 3 s . At NLO and NNLO (upper panel) there is apparently a good agreement between the two theories for v ∼ 0.4. The NNNLO curve in PNRQCD (lower panel), however, lies significantly below the NNLO and the fixed-order curves and shows a very strong dependence on the renormalisation scale. 6 For smaller scales µ we even observe unphysical negative values for R b .
One reason for the considerable difference between fixed-order QCD and PNRQCD at NNNLO may be the limited information from the threshold region used for the Padé approximation. In particular, the behaviour of the NNNLO fixed-order curve for v → 0 is by construction determined by the naïve expansion of the NNLO PNRQCD polarisation function to order α 3 s [56,57]. Although they are suppressed by an additional factor of v, corrections from NNNLO PNRQCD could alter this behaviour considerably as the NNNLO PNRQCD result differs significantly from the NNLO one.
One example for such a big missing correction is given by the product of the thirdorder correction to the hard Wilson coefficient and the leading-order G(E). Parametrically this correction is of order v α 3 s . Due to the numerically large factors involved it still causes a shift of ∆R b ∼ −0.2 at v = 0.4. At higher orders in α s there will be further sizeable contributions from this product, e.g. an additional negative shift of −0.2 at order α 4 s v 0 . While it may be possible to reconcile the discrepancy between the relativistic fixedorder and PNRQCD predictions there still remains the problem that the convergence of the continuum cross section in resummed non-relativistic perturbation theory appears to be quite poor. Neither the difference between consecutive orders nor the scale uncertainty shrink when considering higher-order corrections to the continuum cross section. In fact, at NNNLO the scale uncertainty is much larger than at any lower order.
In our analysis, however, we are of course not interested in the continuum cross section itself but in non-relativistic moments. These are expected to receive their dominant contribution from the resonances. Figure 3 shows the continuum contribution to the tenth moment at different orders in perturbation theory as a function of the renormalisation scale 3 GeV < µ < 10 GeV. As expected we observe that the continuum contribution at NNNLO is rather small, amounting to less than 5% at our central scale µ = m b and about 15% at µ = 10 GeV. Furthermore, the continuum contribution reduces both the distance between consecutive orders and the scale dependence at each order. We conclude that the seemingly problematic behaviour of the continuum cross section does not impede the extraction of the bottom-quark mass from large-n moments at NNNLO.

Determination of the bottom-quark mass
We are now in a position to determine the bottom-quark mass by requiring M th n = M exp n for moments with n ≈ 10. We first eliminate the pole mass in favour of the PS mass m PS b (µ f ), and then convert the resulting mass to the MS scheme. Irrespective of the order of the moments we always perform the conversion at order α 4 s . As a part of our error analysis we will also first convert to the MS mass m MS b (µ) at an intermediate scale µ, which we keep separate from the overall renormalisation scale µ, and then use renormalisation group evolution to find the scale-invariant mass m MS b (m MS b ). To estimate the error of the resulting mass value for a given moment M n we consider the following sources of uncertainties.
• Experimental uncertainty. We add in quadrature the errors induced by uncertainties in the Υ masses, the leptonic widths, the BaBar data directly above the resonances [16], and our estimate 0.1 ≤ R b ≤ 0.5 in the high-energy region.
• Unknown higher-order corrections. As detailed in section 4.1 we choose the extracted PS mass as our central renormalisation scale and vary µ between 3 GeV and 10 GeV.
• QED effects. In our analysis, we include the leading correction in the conversion to the MS mass scheme and the QED Coulomb potential. Assuming α ∼ α 2 s , these effects are formally of NLO. In practice, we find a shift of less than 1 MeV in the extracted quark mass.
• Number of theoretical resonances. In practice we only take into account the contribution from the first six resonances in the formulae (3.5), (3.15) for the theoretical moments and estimate the resulting uncertainty from the difference to considering only four resonances.
• Scheme conversions. When extracting the PS mass, we vary µ f between 1 GeV and 3 GeV. In the conversion to the MS scheme, we vary the intermediate scale µ between 3 GeV and 10 GeV and estimate the error from the unknown value of the conversion coefficient δ MS 3 as described in section 3.3. All scheme conversion errors are added in quadrature.
• Charm mass effects. As detailed in section 3.5 we assign a 100% uncertainty to these corrections.
• Non-perturbative effects. To estimate the order of magnitude, we follow [6] (see also [58]) and consider the leading contribution to the operator product expansion (OPE) from the gluon condensate. It is assumed that higher-dimensional condensates and higher-order corrections to the Wilson coefficients can be neglected. Note that the OPE is only valid for m b v 2 ∼ m b /n ≫ Λ QCD . The need to avoid an uncontrolled, non-perturbative ultrasoft contribution limits the admissible values of n. Under these assumptions we find a negligible effect of less than 1 MeV on the value of the bottom mass.
The smallness of the gluon condensate contribution to moments of order n = 10 is surprising, since the corresponding contribution to the Υ(1S) is large (though uncertain, see the recent discussion in [28]) and the high moments are already completely determined by the Υ(1S). While a certain amount of cancellation of non-perturbative effects is expected in the moments due to quark-hadron duality, the degree of cancellation (about one part in 500 for n = 10) is somewhat puzzling. We therefore advocate that the estimate of non-perturbative effects from the gluon condensate correction is considered with some caution and refrain from using it to limit the allowed values of n. More details on the gluon condensate contribution are given in appendix E.
Our results for the MS masses obtained from moments M n with 6 ≤ n ≤ 15 are shown in figure 4. There appears to be a good agreement of the mass values extracted from moments with n ≈ 10 and a reasonable convergence of the perturbation series. For smaller values of n the behaviour is considerably worse as the non-relativistic approximation becomes less reliable. As can be seen from figure 5, the mass values at NNNLO decrease due to the small continuum cross section near threshold (cf. section 4.2). The behaviour is worse if "expanded" moments M n are considered (lower panel).
As anticipated in section 3.4 there is a considerable difference between the mass values extracted from "unexpanded" moments M th n and "expanded" moments M th n . For the In addition to the uncertainties discussed above we added the differences to the central mass values obtained from M th 8 and M th 12 to our error estimate. The corresponding term is marked by the label (n). Our value for the MS mass is in good agreement with determinations from relativistic (small n) sum rules at order α 3 s [12] and approximate NNNLO non-relativistic sum-rules [11]. A more detailed comparison is given in section 4.4.
An alternative method is to forego the PS scheme and directly use relation (3.11) to eliminate the pole mass in favour of the MS mass. For this "direct" extraction of the In the upper panel we show the resulting mass values for "unexpanded" theoretical moments defined according to (3.15); in the lower panel we have used "expanded" moments (3.5). To facilitate the comparison to existing results we include the mass values obtained from a fixed-order determination [12] for n ≤ 4.  Table 3: Bottom-quark masses obtained from different sum rule analyses. In the last column experimental and theoretical errors were added in quadrature.

Comparison with previous works
As can be seen from table 3, recent sum rule determinations of the bottom-quark mass are in reasonably good agreement. In the following we summarise the differences between the listed analyses.

Chetyrkin et al. (CKMMMSS) [12]
CKMMMSS consider relativistic moments M n , n ≤ 4 in fixed-order perturbation theory to order α 3 s . Our experimental input largely corresponds to the one used by CKMMMSS. However, as low moments are much more sensitive to the experimental high-energy continuum, CKMMMSS rely on the prediction from perturbation theory above 11.24 GeV and use a linear interpolation between this region and the last data point at 11.2062 GeV. CKMMMSS choose the same scale for the strong coupling and the MS mass. Since we vary the two scales independently, the perturbative error estimates are not directly comparable. We expect however that an independent scale variation generally yields a more conservative error estimate. Conversely, if we set µ ≡ µ in our analysis we obtain m MS b (m MS b ) = 4.192 +0.015 −0.031 GeV in place of (4.2) with a significantly smaller estimated uncertainty. The more conservative error estimate using uncorrelated coupling and mass renormalisation scale variations that should be compared to our result (4.2) or (4.3) is unfortunately not available.

Hoang, Ruiz-Femenía and Stahlhofen (HRS) [13]
HRS perform a renormalisation group improved analysis of the moments M n with 6 ≤ n ≤ 14 and a default value of n = 10 in the framework of vNRQCD [59]. In addition to the usual terms scaling as α s √ n also logarithms α s ln √ n are summed, achieving NNLO + (partial) NNLL accuracy. Like CKMMMSS, HRS use the perturbative QCD prediction in place of experimental data for the high-energy continuum, but assign an uncertainty of 10% to it. We agree with the experimental moments used by HRS within the errors. HRS first determine the so-called 1S mass, which is then converted to the MS scheme. The analysis of HRS does not include charm-mass effects, which are currently unknown in their renormalisation group improved framework. HRS assume an uncertainty of 15 MeV in the conversion to the MS mass. This is comparable to our estimate, albeit slightly lower. Like HRS we find that the dependence on instead of µ = m b , but use the hard scale for the vector current matching coefficients. Furthermore HRS neglect charm mass effects and use the MS-pole mass relation at order α 3 s . We find that it is mainly the scale choice which is responsible for the difference between the NNLO results given by HRS and ours. The remaining differences in the analyses have only moderate numerical effects, though their precise size depends on the order in which they are implemented. As a net result, when we adapt our analysis to the one of HRS, we reproduce their NNLO value up to a negligible difference of 2 MeV. We note that the estimate of the perturbative uncertainty at NNLO based on the scale variation of HRS is about twice as large as the one used by us. Given that the convergence of the successive approximations from NNLO to NNNLO at the scale µ = m b is very good (see figure 1), while low scales generally seem to lead to a break-down of non-relativistic perturbation theory [23], we conclude that our error estimate is sufficiently conservative.
The final result of HRS refers to the (partial) NNLL calculation. The resummation of logarithms reduces the scale dependence compared to the NNLO result and therefore compensates part of the large positive shift of the bottom quark mass introduced by choosing a lower scale, leading to the bottom quark mass quoted in table 3.

Penin and Zerf (PZ) [11]
PZ also use non-relativistic sum rules, but choose even higher moments with a central value of n = 15. Such moments are rather insensitive to the experimental high-energy continuum, which in their work is neglected, but potentially introduce an unspecified systematic uncertainty from ultrasoft effects that may already be in the non-perturbative regime (further discussion in appendix E). Choosing n = 15 instead of n = 10 increases the resulting mass value by approximately 12 MeV.
On the theory side, PZ employ the complete NNNLO PNRQCD prediction for the resonances, and an estimate of R b = ρ Z NNNLO 1 /Z NNLO 1 R NNLO , 0.5 ≤ ρ ≤ 2, for the continuum to account for the (then) unknown NNNLO contribution. Our calculation shows that the true result lies somewhat below the band spanned by the variation of the auxiliary parameter ρ. Using rescaled NNLO instead of NNNLO for the continuum leads to an increase of about 5 MeV in the mass extracted from M 15 . Not expanding the factor 1/s in the polarisation function (cf.

Conclusions
We have presented the first complete NNNLO determination of the bottom-quark mass from non-relativistic sum rules. We find a mass of m PS b (2 GeV) = 4.532 +0.013 −0.039 in the PS scheme, which corresponds to an MS mass of m MS b (m MS b ) = 4.190 +0.026 −0.037 . Compared to previous NNLO analyses of non-relativistic moments we observe a significantly reduced uncertainty. This reduction is mostly due to the choice of a higher central scale µ = m PS b , which follows from better insight into the convergence of successive approximations which are now known up to NNNLO. In spite of poor behaviour of the continuum cross section we observe that the NNNLO moments are stable under scale variation, and moments with different n ≈ 10 are in good agreement with each other. Our results agree reasonably well with other recent determinations of the bottom-quark mass from the inclusive e + e − → bb cross section, including NNNLO fixed-order analyses based on relativistic sum-rules. Conservative uncertainty estimates of the bottom-quark MS mass have now reached the ± (25-30) MeV range.

A Charm effects in the Coulomb potential
To determine the effects of a non-zero charm-quark mass we split the Coulomb potential into two partsṼ =Ṽ massless +Ṽ mc , (A.1) whereṼ massless refers to the Coulomb potential for a massless charm quark, so that the charm correctionṼ mc vanishes for m c = 0. The corrections to the Coulomb potential due to the charm-quark mass have been determined at NNLO in [49,50]. We use the following dispersion relation representations from [48], which are very convenient for computations. In momentum space the potential is given bỹ The parameters c 1 , c 2 , d 1 , d 2 are adopted from [48] as , (A.10) The potential in configuration space can be obtained by a Fourier transformation. We obtain In (A.13) we have corrected some typos in Eq. (30) of [48]. This also affects an integral representation in [48] that should read However, to our understanding this does not affect other parts of [48].

B Charm effects in the relation between PS and pole mass
The PS mass at a subtraction scale µ f is defined by [41] Defining the charm corrections to the PS-pole relation For our analysis we set µ c to our overall renormalisation scale µ.

C Charm corrections to bound-state energies and wave functions
Writing the bound-state energies and wave functions in the form we can again split the higher-order terms e (i) N,massless are the respective corrections for the case of massless charm quarks. A general strategy for computing bound state corrections is described in detail in [25]. We refrain from repeating this discussion here. For i = 1 (NLO), the mass corrections originate from a single potential insertion into the Coulomb Green function. where the subscripts {m c } and {m c , m c } denote single and double insertions of the potential V mc , respectively. {m c , massless} denotes the contribution from the mixed double insertion of one potential V mc and one potential V massless . In the following we list our results for these corrections.

C.1 Single insertions
The corrections to the energy levels and wave functions from single insertions of V mc read e (1) , (C.14) The integration contour is chosen such that 0 < Re(u) < 1 on the real axis. 7 Expressing the charm mass corrections in terms of the MS mass m MS c (µ c ) changes the NNLO contributions as follows.
We checked that our results for the binding energy of the Υ(1S) resonance up to NNLO agree numerically with [48]. Furthermore, the energy levels e

C.2 Massive double insertion
For double insertions of the potential V mc we obtain the following corrections to the binding energies and the wave functions: N,mc f (1) The coefficients h (j) mc for j = 0, 1, 2 are 23) and ξ N and I (1) (u) as defined in (C.11), (C.14). We again chose a contour with 0 < Re(u) < 1. It should be noted that during numerical evaluation large cancellations arise between the summands in (C.20) so that typically high-precision arithmetic is required.

C.3 Mixed double insertion
Our result for the mixed double insertion of V mc and V massless is where the coefficients h Here S i (n) = n k=1 k −i denote the generalised harmonic numbers of order i.

D NLO charm effects in the non-relativistic current correlator
The perturbative expansion of G(E) as defined in (3.3) can be written as Splitting off contributions from a non-zero charm-quark mass we define such that δ i G mc (E) vanishes for m c → 0. At NLO charm mass effects enter through single insertions of the potential (A.13). Following [25] we split the correction δ 1 G mc (E) into a part A with no additional Coulomb exchange between the quark and the anti-quark and a part B which resums all ladder diagrams with at least one Coulomb exchange. In contrast to the cases considered in [25] splitting the correction into two parts is not strictly necessary, but still helps to elucidate the structure of the result.
Part A corresponds to a two loop diagram, that will be computed in momentum space. This diagram has no ultraviolet divergence and can be calculated directly in d = 3 dimensions. We obtain and the contour must be chosen such that 0 < Re(u) < 1/2 on the real axis. For part B we perform the calculation in position space. The result can be written as a two-dimensional Mellin-Barnes integral with I (1) (u) as defined in (C.14) and In the integral over the variable u the contour should be fixed such that 0 < Re(u) < 1 on the real axis. Note that by choosing the contour to the right of the pole at u = 0 we have accounted for the subtraction term in the potential (A.13). The notation Γ * (−w−1) denotes that the contour should be chosen to the right of the first pole at w = −1. For positive E it can be fixed parallel to the imaginary axis with −1 < Re(w) < 0, since the real part of λ vanishes and left and right poles are separated. For positive integer values N of λ Eq. (D.5) contains poles in N − λ due to the pinching of the contour in the complex w plane by left and right poles. These poles determine the charm mass corrections to the resonances as obtained in appendix C.

E Gluon condensate correction
The gluon condensate correction to the PNRQCD Green function is given by [58] δ G 2 G(E) = − 0|πα s G 2 |0 18 d 3 r d 3 r ′ (r · r ′ ) G (0) (0, r; E)G (8) (r, r ′ ; E)G (0) (r ′ , 0; E), where the superscript (0), (8) refers to the colour-singlet and colour-octet Coulomb Green function, respectively (cf. [17]). Proceeding as in [25], we find the representation s! H G 2 (s) 2 (s + 3)!(s + 2 + λ/8) , where λ is defined in (D.6), and The sum in (E.2) can be evaluated in terms of polygamma functions. Expanding (E.2) around the bound-state poles at λ = n determines the condensate correction to the S-wave energy levels E N and wave functions at the origin, |ψ N (0)| 2 [61,62]. Eq. (E.2) can be integrated in the complex energy plane to yield the condensate contribution to the moments. It is worth noting that splitting this contribution into a resonance and continuum contribution is ill-defined, since the two are separately divergent. For the resonances the divergence arises in the sum over N, since the condensate correction rises too rapidly with principal quantum number N. For the continuum, the integral over energy is too singular at E = 0. This reminds us that the moment calculation really refers to the calculation of high derivatives of Π b q 2 at q 2 = 0, and assumes the validity of the operator product expansion (OPE) for this quantity, for which the split into resonances and continuum contributions is artificial.
It is well-defined to split the condensate contribution into the one from the Υ(1S) resonance and the rest. The result normalised to the NNNLO theoretical moments (without the condensate contribution) is given in table 4 together with the corresponding split-up of the perturbative contribution to the moment, and the experimental moment.
We first note that the theoretically computed, perturbative Υ(1S) contribution to the moments is very close to the experimental one. Both are large, increasing from 80% for the 10th moment to more than 95% for n = 24. Next, we observe that the gluon condensate correction to the Υ(1S) contribution to the moment is extremely large. Taken at face value, it exceeds the perturbative moment by a factor of about two. While there is a large ambiguity in the absolute size of the gluon condensate contribution (mainly related to the choice of scale in α s in the expression for K), as discussed in [28], the enormous cancellation between the contribution to the Υ(1S) resonance, δ G 2 M Υ(1S) n / M pert n , and the remainder, δ G 2 M rest n / M pert n , is independent of this size. For n = 8 it is effective to one part in 1000, and it remains at the 1% level even for very high moments n = 24. As a consequence, the total gluon condensate correction (last row in the table) remains very small, around 2%, for n = 24 when the ultrasoft scale m b /n is clearly in the non-perturbative regime.