Bottom Quark Mass from $\Upsilon$ Sum Rules to ${\cal O}(\alpha_s^3)$

We use the ${\cal O}(\alpha_s^3)$ approximation of the heavy-quark vacuum polarization function in the threshold region to determine the bottom quark mass from nonrelativistic $\Upsilon$ sum rules. We find very good stability and convergence of the perturbative series for the bottom quark mass in $\rm\overline{MS}$ renormalization scheme. Our final result is $\overline{m}_b(\overline{m}_b)=4.169\pm 0.008_{th}\pm 0.002_{\alpha_s}\pm 0.002_{exp}$.


Introduction
The bottom quark mass m b is a fundamental parameter of the Standard Model of particle interactions. It is an essential input parameter for the analysis of the B-meson decays and CKM quark mixing matrix as well as the Higgs boson decay rates and branching ratios. The precise value of the bottom quark mass is also crucial for testing possible extensions of the Standard Model such as grand unified theories. Thus determination of the bottom quark mass with the best possible precision is an important problem of particle phenomenology. A unique tool for such a determination is given by the analysis of the family of Υ resonances within quantum chromodynamics. Direct application of perturbative QCD to the description of the heavy-quarkonium properties such as the resonance mass and width suffers from sizable long-distance nonperturbative effects [1,2] resulting in large uncertainties even if high-order approximation is available [3][4][5][6][7]. The sum rules approach [8,9] suggests an elegant solution of the problem. It relates the moments of the spectral density saturated with the contribution of Υ resonances to the derivatives of the heavy quark vacuum polarization function in a deep Euclidean region, which can be reliably computed in perturbation theory. Thus the approach provides a model independent determination of the bottom quark mass entirely based on the first principles of QCD. The low-moment or "relativistic" and the high-moment or "nonrelativistic" sum rules require essentially different experimental and theoretical input and can be considered as complimentary methods. Both methods have been extensively applied to the bottom quark mass determination. The most recent analysis of the lowmoment sum rules includes the third-order corrections in the strong coupling constant α s to the leading-order result [10][11][12]. At the same time the high-moment sum rules have been evaluated only through the next-to-next-to-leading order (NNLO) [13][14][15][16][17][18][19][20][21] though the effect of higher order logarithmically enhanced terms have been considered [22,23]. In this paper we present the complete O(α 3 s ) corrections to the heavy quarkonium parameters required for the N 3 LO analysis of the nonrelativistic Υ sum rules and apply the result to the determination of the bottom quark mass. In the next section we outline the main concept of the nonrelativistic Υ sum rules and describe the perturbative approximation. The numerical analysis is given in Section 3. In Section 4 our new estimate of the bottom quark mass is compared to the existing high-order results.

Υ sum rules
We consider the vacuum polarization function Π(q 2 ) defined through the two-point vacuum correlator of the heavy-quark electromagnetic current j µ =bγ µ b q µ q ν − g µν q 2 Π(q 2 ) = i d d x e iqx 0|T j µ (x)j ν (0)|0 . (2.1) Its nth moment is given by the normalized derivative where s = q 2 and R(s) = 12πImΠ(s + i ) is the spectral density. The long-distance nonperturbative contribution to Eq. (2.2) is parametrically suppressed as Λ 4 QCD /m 4 b and the moments can be reliably computed in perturbative QCD [9]. On the other hand the optical theorem relates the spectral density to the experimentally measured cross section of bb hadron production in electron-positron annihilation where where α(2m b ) is the running QED coupling constant, M Υ(mS) (Γ Υ(mS)→l + l − ) is the resonance mass (leptonic width), and ellipsis stand for the nonresonant contribution. The sum rules then read M exp n = M th n , (2.5) where the theoretical moment M th n is evaluated with the perturbative QCD approximation for the spectral density. For small n the moments get sizable contribution from relativistic region above bb pair threshold where the uncertainty of the measured cross section is relatively large. For large n the experimental moments are saturated with the contribution of the lowest Υ resonances measured with very high accuracy. Moreover for large n the moments (2.4) are very sensitive to the value of m b , which significantly reduces the uncertainty of the extracted bottom quark mass. The number of a phenomenologically relevant moment is limited only by the magnitude of nonperturbative contribution which grows as n 3 and becomes sizable for n > ∼ 20 [14]. The perturbative description of the high moments, however, is nontrivial. For large n the theoretical moments are saturated with the nonrelativistic threshold region where the heavy quark velocity v is of order 1/ √ n [13]. This results in enhancement of the Coulomb effects which are characterized by the expansion parameter α s /v ∼ √ nα s rather then α s . For n > ∼ 1/α 2 s ∼ 10 the Coulomb terms are not suppressed and have to be resummed to all orders. Thus for high moments the QCD perturbation theory should be build up about the nonrelativistic Coulomb solution instead of the free heavy-quark approximation used for analysis of the low moments. In the next section we outline how this can be done systematically within the nonrelativistic effective field theory framework.

Effective theory approach to nonrelativistic sum rules
In the threshold region the heavy-quark velocity v is a small parameter. An expansion in v may be performed directly in the QCD Lagrangian by using the concept of effective field theory [24]. The relevant modes are characterized by the hard (k 0 , k ∼ m q ), the soft (k 0 , k ∼ m q v), the potential (k 0 ∼ m q v 2 , k ∼ m q v), and the ultrasoft (k 0 , k ∼ m q v 2 ) scaling of energy k 0 and three-momentum k in respect to the heavy-quark mass m q . Integrating out the hard modes matches QCD onto non-relativistic QCD (NRQCD) [25]. By subsequent integrating out the soft modes and potential gluons one obtains the effective theory of potential NRQCD (pNRQCD) [26][27][28][29], which contains potential heavy quarks and ultrasoft gluons as dynamical fields relevant for the description of the heavy-quark threshold dynamics. In this theory the propagation of a color-singlet quark-antiquark pair is described by the Green function G s (r, r ; E). In the leading-order Coulomb approximation the Green function satisfies the Schrödinger equation with the Hamiltonian where r = |r|, m q the heavy-quark pole mass, and C F = (N 2 c −1)/(2N c ), N c = 3. The leadingorder Green function gets corrections due to the high-order terms in the nonrelativistic and perturbative expansion of the effective Hamiltonian as well as due to the multipole interaction of the potential quarks to the ultrasoft gluons. In the effective theory the electromagnetic current is represented by a series of operators composed of the nonrelativistic quark and antiquark two-component Pauli spinor fields ψ and χ where the matching coefficients c v = 1 + O(α s ) and d v = 1 + O(α s ) are the series in α s . The threshold behaviour of the heavy-quark polarization function is given by the following pNRQCD expression where E = q 2 − 2m q ∼ v 2 m q and only the component of total spin one is kept in the Green function as dictated by the form of the production current. The corrections to the leading order Coulomb approximation for the polarization function (2.9) can be computed in pNRQCD as a series in α s and v ∼ α s according to the effective theory power counting, which gives a systematic perturbative expansion for the high moments. The explicit result for the polarization function up to the NNLO can be found in Refs. [18,19]. The spectral representation of the color-singlet Green function includes an infinite number of Coulomb-like bound state poles where the ellipsis stand for the continuum contribution, and E n (ψ n (r)) is the energy (wave function) of the bound state with spin S = 1 and orbital angular momentum l = 0. In the Coulomb approximation they read Thus, the effective theory expression for the moments can be written in the following form where and R(s) for s > 4m 2 b is determined by the imaginary part of Eq. (2.9). The quantities E n and Z n determine the perturbative QCD predictions for the Υ(nS) resonance mass and leptonic width (2.14) The nonperturbative corrections to the binding energy and the width in Eq. (2.14) are parametrically as large as Λ 4 QCD /(α 6 s m 4 b ), being significantly enhanced in comparison to the moments [1,2]. They grow rapidly with the principal quantum number and result in a large theoretical uncertainty even in the case of the Υ(1S) state. On the other hand for high moments the contribution of the sum to Eq. (2.12) significantly exceeds the one of the integral. The evaluation of the perturbative corrections to the parameters E n and Z n is therefore crucial for high-order analysis of the nonrelativistic sum rules. In the next section we present the complete O(α 3 s ) result for these quantities.
2.2 Heavy quarkonium mass and leptonic width to O(α 3 s ) We parameterize the perturbative series for the resonance mass and leptonic width as follows is the MS renormalized coupling with n l light-quark flavors and e (0) n = z (0) n = 1. The first two coefficients of the series for the binding energy are well known [30] and listed in the Appendix A. The third-order term can be decomposed according to the powers of the logarithm where L n = ln (nµ/α s C F m q ) and L = ln (µ/m q ). The full analytical expression of the coefficients in Eq. (2.16) for arbitrary n [31,32] is too cumbersome and we only present their values for the six lowest states. For the ground state they read [4] where ζ(z) is the Riemann zeta-function, ζ(3) = 1.20206 . . .. The last coefficient incorporates the three-loop contribution to the static potential [33,34]. For n = 2, . . . , 6 the coefficients are listed in the Appendix B. The corrections to the leptonic width up to O(α 2 s ) can be read off the results [18,19,35,36] and are given in the Appendix A. The third-order term for general n has only been evaluated in the logarithmic approximation [5,6,37]. Recently the total third-order correction for n = 1 has been published in a numerical form [7]. Below we present the result for the excited states. As it follows from Eq. (2.15), the O(α 3 s ) contribution consists of: (i) Interference of the leading and next-to-leading order corrections from each of three factors in Eq. (2.13).
(ii) The one-loop correction to the matching coefficient d v [38] and recently completed three-loop corrections to matching coefficient c v [39][40][41].
(iii) The corrections to the wave function [31,32,42,43] due to the N 3 LO operators in the effective Hamiltonian [33,34,44,45] and due to the multiple iterations of the next-toleading and NNLO operators.
(iv) The ultrasoft correction to the wave function [46] Combining all the above contribution we get the complete result for the third-order term in the following form For the ground state we obtain the following values of the n-dependent coefficients where the error in δ (0) z is due to the numerical evaluation of c v [41]. The corresponding expressions for n = 2, . . . , 6 are listed in the Appendix B.

Numerical analysis
Now we are in a position to apply the method described in the previous section for the determination of the bottom quark mass. We use the parameters of the six Υ resonances listed  (12) 10.876 (11) 11.019(8) Γ Υ(nS)→e + e − (keV) 0.272 (29) 0.31 (7) 0.130(30) Table 1. Experimental values of the Υ-resonance masses and leptonic widths [47].
in Table 1 as the experimental input. The 5th and 6th resonances actually lie above the Bmeson production threshold and their contribution does not represent the total experimental spectral density in this region. We use this contribution only to estimate the experimental uncertainty of our result. For the QED running coupling we adopt the value α(2m b ) = 1.036 α [48], where α = 1/137.036 is the fine structure constant.
On the theory side we use the complete O(α 3 s ) expression for the contribution of the perturbative heavy-quarkonium bound states below the bb threshold given in Section 2.2 for n l = 4. The continuum contribution from the above-threshold region is strongly suppressed for high moments and high accuracy of the theoretical approximation there is not mandatory. Thus, without introducing significant error we emulate the N 3 LO spectral density for s > 4m 2 b by rescaling the NNLO result [16] To estimate the error introduced by this approximation we multiply the total continuum contribution to the theoretical moments by 1/2 < ρ < 2, which corresponds to the variation of its absolute value by factor four. The sequence of the values of the bottom quark pole mass m b extracted from the sum rules order by order in perturbation theory does not converge well. This is expected since the pole mass is widely believed not to be a good parameter of perturbative QCD due to infrared renormalon contribution (see [49] for a review and [50] for a recent high-order analysis). The perturbative behavior of the "short-distance" mass parameter m b (µ) defined in MS renormalization scheme is supposed to be much better. Therefore we convert the extracted pole mass value into m b (m b ) according to the relation where the coefficients r (n) have been evaluated up to n = 3 [51,52]. To achieve the cancellation of factorially growing terms associated with the infrared renormalon one has to correlate perturbative approximations for the mass relation and the sum rules in such a way that the series (3.2) is truncated at one order higher than the series (2.15), i.e. the one-loop mass relation is used with the Coulomb approximation for the moments, and so on [53,54]. Our analysis therefore requires the four-loop coefficient r (4) . Since the exact value of this coefficient is not yet available we use the renormalon-based estimate [55] r (4) ren ≈ 1346 . This method reproduces the value of the three-loop coefficient r (3) with the precision better than one part in a thousand. We take the difference between Eq. (3.3) and the large-β 0 prediction r An important issue of the numerical analysis of high moments is whether the factor 1/(2m b + E m ) 2n+1 in Eq. (2.12) should be expanded about 1/(2m b + E C m ) 2n+1 . Formally all the terms of this expansion beyond the N 3 LO are suppressed according to the √ n ∼ 1/α s power counting and can be neglected in our approximation. Such an expansion, however, merely violates the cancellation of the large perturbative corrections related to the infrared renormalon in the series for the binding energy (2.15) and for the pole mass (3.2) [18,20,21,[53][54][55]. This spoils the convergence of the resulting series for m b since the factorial growth of the coefficients e For the numerical estimates we use the moments in the interval 10 > ∼ n > ∼ 20, where the nonrelativistic perturbative approximation for the spectral density is valid, the nonperturbative effects are under control, and the result is almost insensitive to the continuum contribution to the theoretical moments. The renormalization scale is varied in a physically motivated interval between the soft scale µ ∼ α s m b and the hard scales µ ∼ m b . We take α s (M z ) = 0.1184 ± 0.0007 [47] as an input and run it down to µ = m b with the four-loop beta-function [56]. To convert α s (m b ) into α s (µ) used in our numerical analysis we correlate the order of the renormalization group evolution of the strong coupling constant with the perturbative expansion for the sum rules so that the one-loop running is used in the leading approximation, an so on. To ensure the renormalon cancellation we reexpress Eq. (3.2) through α s (µ) and use the same renormalization scale both for the sum rules and the mass relation.
A stable perturbative result is achieved for n > ∼ 10 and µ > ∼ 3.5 GeV, see Fig. 2. We take m b (m b ) = 4.194 as a central value of our estimate. It corresponds to n = 15 in the center of the allowed interval and to the renormalization scale µ = m b , which belongs to the stability plateau and provides nα 2 s ≈ 0.8 in agreement with the power counting rules. The uncertainty budget of our estimate is summarized in Table 2. The experimental part of the uncertainty ∆ exp accounts for both the error bars in the measured values of the resonance mass and width, and the contribution from the region above the B-meson production threshold. We estimate the latter by the size of the 5th and 6th Υ-resonance contribution to the experimental moments. The uncertainty ∆ αs corresponds to the error in the input value of α s (M Z ). The quantities ∆ ρ and ∆ r (4) account for the approximate character of the Eqs. (3.1, 3.3). The value of ∆ ρ is given by a half of the variation of the result with the parameter ρ changing from    1/2 to 2. In the same way ∆ n is given by a half of the variation of the result with the moment number spanning the interval 10 < n < 20. We take one half of the third-order correction as an estimate of the uncertainty ∆ p.t. introduced by truncation of the perturbative expansion. This procedure can be verified through the N 3 LO. Indeed, the numerical series for the bottom quark mass reads where the third-order correction amounts only about a quarter of the second-order term. Each individual contribution to the total uncertainty is computed with all other parameters fixed at their central values and at the normalization scale µ = m b . We refrain from using the scale dependence for the uncertainty estimate since this procedure strongly depends on the low boundary of the allowed scale variation which cannot be unambiguously defined. We therefore restrict the analysis only to the large stability region where the variation of the result due to the change of the scale is much smaller than our estimate of perturbative uncertainty given above. The choice of the hard renormalization scale may look ambiguous since for large n the soft scale m b / √ n and the ultrasoft scale m b /n are involved. As a consequence the coefficients of the series get contributions enhanced by logarithm of a scale ratio proportional to ln n, which may affect the convergence of the perturbative expansion. The logarithmic terms cannot be completely eliminated by adjusting the renormalization scale of α s but can be resummed to all orders through the effective theory renormalization group [22,23,57]. However, for our choice of the moments ln n < 3, i.e. the asymptotic regime is not yet reached. In fact for such n the logarithmic terms do not saturate the coefficients of the series numerically and do not have to be distinguished from the nonlogarithmic contributions. This justifies our choice of the renormalization scale dictated solely by the convergence of the perturbation theory.
The moments get also a contribution from the nonperturbative scale Λ QCD . Within the operator product expansion it is given by a series in nΛ 2 QCD /m 2 b . The leading nonperturbative contribution to the high moments due to the gluon condensate turns out to be numerically suppressed [14]. For n = 15 and αs π G 2 ≈ 0.012 GeV 4 the corresponding correction to the bottom quark mass is about −0.8 MeV. We take this value as the nonperturbative uncertainty ∆ n.p. of our result. Though the nonperturbative uncertainty rapidly increases with the moment number, the perturbative one is suppressed for higher moments and their sum does not significantly change over the whole range of n considered in the paper. One may be concerned that for large n the ultrasoft scale m b /n approaches Λ QCD , which questions the perturbative treatment of the ultrasoft contribution. However, the result of Ref. [14] ensures that even for n = 20 the moments do not get a sizable contribution from the gluonic field fluctuations at the scale Λ QCD and the perturbative description of the moments is justified.

Charm mass effect
So far we considered the charm quark to be massless. Since its mass m c is not much smaller than the soft scale, the charm quark mass effect on the effective potential and the bound state dynamics may not be negligible [58,59]. This effect has been analyzed in the context of Υ sum rules through the NNLO [60]. By using the results of Ref. [60] for the charm quark mass correction to the moments and to the mass relation, for m c (m c ) ≈ 1.3 GeV [11] and a set of input parameters similar to the one adopted in this section we get a negative shift of approximately 25 MeV in the extracted value of m b (m b ) with an estimated error ∆ mc = 5 MeV. We incorporate this correction in the analysis. Our final prediction for the bottom quark mass reads m b (m b ) = 4.169 ± 0.008 th ± 0.002 αs ± 0.002 exp , (3.5) where the theoretical error corresponds to ∆ ρ , ∆ r (4) , ∆ n , ∆ p.t. , ∆ n.p. , and ∆ mc added up in quadrature.

Υ(1S) mass and leptonic width
Though significant nonperturbative effects are expected in the QCD analysis of the Υresonance mass and width, it is instructive to figure out how perturbative QCD results [4,7] reproduce the experimental data for the lowest Υ(1S) state, where the nonperturbative contribution is minimal. From the result of the previous section we get the following numerical series for the ground state mass and width for m c = 0 in agreement with [4,7]. In Fig. 3 we plot M Υ(1S) and Γ Υ(1S)→e + e − evaluated according to Eqs. The O(α 3 s ) approximation is in a rather good agreement with the experimental data for the resonance mass (width). The difference of about 60 MeV (0.3 keV) quantitatively agrees with an estimate of the nonperturbative contribution of the gluon condensate [61]. The inclusion of the charm mass effects increases the difference between the perturbative result and the measured Υ(1S) mass to approximately 90 MeV since the reduction of the binding energy [60] does not fully compensate the negative corrections to the input value of m b (3.5). The interpretation of the nonperturbative contribution to Eqs. (3.6, 3.7) in this case is more ambiguous [7] but still consistent with the gluon condensate. This validates our estimate of the nonperturbative correction to the sum rules, which is also based on the gluon condensate contribution. It is interesting that for large n the nonperturbative correction to Eq. (3.5) is of the order of n 2 Λ 4 QCD /m 4 b and for n ∼ 1/α 2 s is not parametrically suppressed in comparison to the first equation of (2.14). Nevertheless, numerically the gluon condensate contribution to the sum rules result for m b is almost two orders of magnitude smaller than the nonperturbative correction to the resonance mass, i.e. even for the large n the moments are significantly less sensitive to the long-distance phenomena.

Method
Approximation m b (m b ) (GeV) Ref. [ Table 3. The results of the bottom quark mass determination from Υ family properties beyond the NNLO. In the last line all the errors given in Table. 2 are added up in quadrature.

Summary and discussion
The main result of this work is the new value of the bottom quark mass (3.5) from the O(α 3 s ) analysis of the nonrelativistic Υ-sum rules. In Table 3 we confront Eq. (3.5) with the existing results of the bottom quark mass determination from Υ phenomenology beyond the NNLO of perturbation theory. In Refs. [4,32] the bottom quark mass has been obtained from the O(α 3 s ) approximation for M Υ(1S) . A relatively large value of the bottom quark mass reported in Ref. [4] is due to the choice of the "physical" soft renormalization scale µ = 2.7 GeV natural for the bound state dynamics. However, the perturbative expansion becomes unstable at such a low scale (cf. Fig. 3(a)). For µ = 4.20 GeV the analysis [4] gives the value m b (m b ) = 4.22 ± 0.07 consistent with Refs. [21,32] and Eq. (3.5).
The high-moment sum rules have been considered in a context of the effective theory renormalization group in Refs. [22,23]. Both analyses involve partial resummation of the next-to-next-to-leading logarithms (NNLL) of the form α m+2 s ln m α s for all m. 1 The result of Ref. [22] agrees with Eq. (3.5) within the error bars. Though the renormalization group resummation improves the behavior of the perturbative expansion especially at a low renormalization scale, the logarithmic terms do not dominate the perturbative series (2.15) and cannot be used for a precise quantitative estimate of the third-order corrections. As a consequence, the theoretical error of the NNLL approximation is significantly larger than the one of the N 3 LO result.
The most accurate value of the bottom quark mass up to date has been reported in Ref. [11] and is obtained from the relativistic sum rules at O(α 3 s ). The central value given in Table 3 corresponds to n = 2. The error estimate includes ±10 MeV due to uncertainty of the experimentally measured cross section, ±12 MeV due to the input value α s (M Z ) = 0.1189 ± 0.002, and the theoretical uncertainty ±3 estimated by the variation of the renormalization scale. The overall accuracy of the relativistic sum rules is comparable to Eq. (3.5) for a given interval of α s (M Z ). However for the low moments the experimental error clearly dominates the theoretical one. Our result is in a perfect agreement with the analysis [11]. The amazing agreement of two approaches based on significantly different theoretical and experimental input boosts our confidence in the result and, in particular, in the uncertainty assessment.