Bottom and Charm Mass Determinations with a Convergence Test

We present new determinations of the MS-bar charm quark mass using relativistic QCD sum rules at O(alpha_s^3) from the moments of the vector and the pseudoscalar current correlators. We use available experimental measurements from e+e- collisions and lattice simulation results, respectively. Our analysis of the theoretical uncertainties is based on different implementations of the perturbative series and on independent variations of the renormalization scales for the mass and the strong coupling. Taking into account the resulting set of series to estimate perturbative uncertainties is crucial, since some ways to treat the perturbative expansion can exhibit extraordinarily small scale dependence when the two scales are set equal. As an additional refinement, we address the issue that double scale variation could overestimate the perturbative uncertainties. We supplement the analysis with a test that quantifies the convergence rate of each perturbative series by a single number. We find that this convergence test allows to determine an overall and average convergence rate that is characteristic for the series expansions of each moment, and to discard those series for which the convergence rate is significantly worse. We obtain mc(mc) = 1.288 +- 0.020 GeV from the vector correlator. The method is also applied to the extraction of the MS-bar bottom quark mass from the vector correlator. We compute the experimental moments including a modeling uncertainty associated to the continuum region where no data is available. We obtain mb(mb) = 4.176 +- 0.023 GeV.


Introduction
Precise and reliable determinations of the charm and bottom quark masses are an important input for a number of theoretical predictions, such as Higgs branching ratios to charm and bottom quarks or for the corresponding Yukawa couplings [1,2]. They also affect the theoretical predictions of radiative and inclusive B decays, as well as rare kaon decays. For example, the inclusive semileptonic decay rate of B mesons depends on the fifth power of the bottom quark mass. These weak decays provide crucial methods to determine elements of the CKM matrix, which in turn are important for testing the validity of the Standard Model, as well as for indirect searches of new physics. In this context, having a reliable estimate of uncertainties for the quark masses is as important as knowing their precise values [3]. Due to confinement quark masses are not physical observables. Rather, they are scheme-dependent parameters of the QCD Lagrangian which have to be determined from quantities that strongly depend on them.
One of the most precise tools to determine the charm and bottom quark masses is the QCD sum rule method, where weighted averages of the normalized cross section R e + e − → qq +X , with q = c, b, can be related to moments of the quark vector current correlator Π V [4,5]: , j µ (x) =q(x)γ µ q(x) , Here Q q is the quark electric charge and √ s = q 2 is the e + e − center-of-mass energy. Given that the integration over the experimental R-ratio extends from the quark pair threshold up to infinity but experimental measurements only exist for energies up to around 11 GeV, one relies on using theory input for energies above that scale (which we call the "continuum" region). For the charm moments, the combination of all available measurements is actually sufficient to render the experimental moments essentially independent of uncertainties one may assign to the theory input for the continuum region [6]. For the bottom moments, the dependence on the continuum theory input is very large, and the dependence of the low-n experimental moments on unavoidable assumptions about the continuum uncertainty can be the most important component of the error budget, see e.g. [7]. In fact, the use of the first moment M V 1 to determine the bottom mass appears to be excluded until more experimental data becomes available for higher energies.
Alternatively one can also consider moments of the pseudoscalar current correlator to extract the heavy quark masses. Experimental information on the pseudoscalar correlator Π P is not available in a form useful for quark mass determinations, but for the charm quark very precise lattice calculations have become available recently [8]. For the pseudoscalar correlator it turns out that the first two Taylor coefficients in the small-q 2 expansion need to be regularized and defined in a given scheme, and that the third term (which we will denote by M P 0 ) is hardly sensitive to m q . We adopt the definitions Π P (s) = i dx e iqx 0 |T j P (x)j P (0)| 0 , j P (x) = 2 m q iq(x)γ 5 q(x) , where the explicit mass factor in the definition of the pseudoscalar current ensures it is renormalization-scheme independent.
For small values of n such that m q /n Λ QCD , the theoretical moments for the vector and pseudoscalar correlators can be computed in the framework of the OPE, i.e. as an expansion in vacuum matrix elements involving operators of increasing dimension [4,5]. The leading matrix element corresponds to the perturbative QCD computations, which greatly dominates the series. Nonperturbative corrections are parametrized by vacuum condensates, and we find that even the leading correction, given by the gluon condensate, has a very small effect for low n, particularly so for the bottom correlator. For moments at low values of n, it is mandatory to employ a short-distance mass scheme such as MS [9], which renders the quark mass m q (µ m ) dependent on its renormalization scale µ m , similar to the strong coupling α s (µ α ), which depends on µ α . This method of determining the heavy quark masses with high precision is frequently called relativistic charmonium/bottomonium sum rules.
The most recent determinations of the MS charm mass from charmonium sum rules for the vector correlator [6,26,27] obtain very accurate results, but differ in the way they estimate theoretical uncertainties, and also in the computation of the moments from experimental data. Concerning the estimate of the perturbative uncertainties, Ref. [6] obtained 19 MeV compared to 1 and 2 MeV obtained in Refs. [26] and [27], respectively. The discrepancy arises from two differences. First, in Refs. [26,27] the renormalization scales µ m , and µ α were set equal, while in Ref. [6] is was argued that they should be varied independently. Second, in Refs. [26,27] the lowest renormalization scale was chosen to be 2 GeV, while in Ref. [6] variations down to the charm mass value were considered. For the case of the pseudoscalar moments, the HPQCD collaboration has made a number of very accurate predictions for charm and bottom masses [8,28,29], the last of which has the smallest uncertainty claimed so far for the charm mass. In all these analyses the renormalization scales µ m , and µ α are set equal when estimating the truncation uncertainty. A detailed discussion on the estimates of theoretical and experimental uncertainties can be found in Secs. 3 and 10 of this article (see also Ref. [6]).
Similarly, bottomonium sum rules have been used to determine the bottom mass from low-n moments. To the best of our knowledge, the most recent and precise determinations are from Refs. [27,30]. These two analyses estimate their perturbative uncertainties in the same way as the corresponding charm mass extractions from the same collaborations [26,27]. Furthermore, when it comes to compute the experimental moments, they use theoretical input at O(α 3 s ) with perturbative uncertainties to model the high-energy region (continuum region) of the spectrum. As we discuss in this work, similar caveats as for their charm analyses can be argued to also affect their bottom quark results.
In this work we revisit the charmonium sum rules for the vector correlator, refining our perturbative error estimate from Ref. [6] by incorporating a convergence test. The convergence test addresses the issue that the independent variation of µ m and µ α together with the relatively large value of the α s close to the charm mass scale, might lead to an overestimate of the perturbative uncertainty. The convergence test allows to quantify the convergence property of each perturbative series with a single parameter and to discard series for which the convergence is substantially worse then for the rest of the series. We show that this procedure is meaningful, since the complete set of series for the moments shows a strongly peaked distribution in these convergence values, which allows to define an overall convergence for the set of perturbative series. This leads to a reduction of the perturbative uncertainty from 19 to 14 MeV, and the corresponding result for the MS charm mass supersedes the main result given in Ref. [6]. We also apply this improved method of estimating theory uncertainties to obtain a new charm mass determination from the pseudoscalar correlator, and to extract the bottom mass from the vector correlator. For the latter, we compute the bottom experimental moments by combining contributions from narrow resonances, experimental data taken in the continuum, and a theoretical model for the continuum region. We carefully study the assignment of adequate uncertainties to this last contribution, to make sure that the the model dependence is reduced to an acceptable level.
This paper is organized as follows: In Sec. 2 we summarize the theoretical framework introduced in [6], and adapt it to cover the case of the pseudoscalar moments. We also introduce the ratios of moments and the perturbative expansions associated to them. Sec. 3 contains a brief summary of the results obtained in [6], and the discussion is extended to the case of the pseudoscalar correlator and the bottom mass. In Sec. 4 we introduce the convergence test, and discuss how it allows to identify and discard series with a bad convergence. Sec. 5 contains a discussion on the lattice simulation results we use for our analysis. In Sec. 6 we present our computation of the experimental moments for the bottom correlator. The results are compared to previous determinations in Sec. 7. The computation of the ratio of experimental moments is presented in Sec. 8. Our final results for the quark masses are given in Sec. 9. The results are compared to previous charm and bottom mass analyses in Sec. 10. We present our conclusions in Sec. 11. In Appendix A the numerical values of the perturbative coefficients that enter into our analysis and are not yet provided by Ref. [6] are collected for the convenience of the reader.

Perturbative Contribution
The moments of the vector and pseudoscalar current correlators are defined in Eqs. (1.1) and (1.3), respectively. In the OPE framework they are dominated by the perturbative contribution (that is, a partonic computation), which exhibits a nonlinear dependence on the mass. Within perturbation theory one can decide to manipulate the series expansion to get a more linear dependence on the mass. Conceptually there is no preference. As advocated in our previous analyses [6], one might consider various versions of the expansion to reliably estimate the perturbative uncertainties. Four types of expansion were suggested in Ref. [6], which we briefly review below.

(a) Standard fixed-order expansion
We write the perturbative vacuum polarization function as where X = V, P for vector and pseudoscalar currents, respectively. Note that for notation reasons, in Eqs. (2.1, 2.5, 2.7, 2.8) we use Π P (q 2 ) = P (q 2 ), where P is the twicesubtracted pseudoscalar correlator defined in Eq. (1.3). Here Q q is the quark electric charge with q = c, b, and n f = 4, 5 for charm and bottom, respectively. In full generality, the perturbative momentsM n can be expressed as the following sum: This is the standard fixed-order expression for the perturbative moments. The numerical values for the [C V (n f = 5)] a,b n,i coefficients are given in Table 8 (for the vector current with n f = 4, these coefficients can be found in Table 1 of Ref. [6]). Likewise, [C P (n f = 4)] a,b n,i are collected in a numerical form in Table 11.
The expression in Eq. (2.2) is the common way to write the perturbative series of the moments. However, as noted in Ref. [6], the nonlinear dependence on m q of the standard fixed-order expansion has the disadvantage that for charm (bottom) quarks there are frequently no solutions for the mass in the sum rule mass determination, for moments higher than the first (second), for some set of values of the renormalization scales.

(b) Linearized expansion
One can linearize the the fixed-order form expansion of Eq. (2.2) with respect to the exponent of the quark mass pre-factor by taking the 2n-th root. This choice is e.g. made in Ref. [28], and in general one can write: 3) The coefficients [C V (n f = 5)] a,b n,i and [C P (n f = 4)] a,b n,i are given in Tables 9 and 12, respectively (for n f = 4 the coefficients for the vector current can be found in Table 2 of Ref. [6]). Even though relation (2.3) still exhibits some nonlinear dependence on m q through perturbative logarithms, we find that it always has a solution for the quark mass.

(c) Iterative linearized expansion
For the expansion methods (a) and (b) shown in Eqs. (2.2) and (2.3), one solves for the quark masses m c,b (µ m ) numerically keeping the exact mass dependence on the theory side of the equation. Alternatively, one can solve for m c,b (µ m ) iteratively order by order, which is perturbatively equivalent to the exact numerical solution, but gives different numerical results. The method consists of inserting the lower order values for m c,b (µ m ) in the higher order perturbative coefficients, and re-expanding consistently. This method has been explained in detail in Sec. 2.1(c) of Ref. [6], and we only quote the final results here: where the numerical value of the coefficients [C V (n f = 5)] a,b n,i and [C P (n f = 4)] a,b n,i are collected in Tables 10 and 13, and the values for the vector current with n f = 4 can be found in Table 3 of Ref. [6].
By construction, the iterative expansion always has a solution for the quark mass. Accordingly, potential biases on the numerical analysis related to any possible nonlinear dependence are eliminated.

(d) Contour-improved expansion
For the expansions (a), (b) and (c) the moments and the quark masses are computed for a fixed value of the renormalization scale µ α . Using the analytic properties of the vacuum polarization function, one can rewrite the fixed-order moments as integrals in the complex plane. This opens the possibility of making µ α dependent on the integration variable, s (1,0) One possible integration path in the complexs-plane for the computation of the contour-improved moments.
in analogy to the contour-improved methods used for τ -decays (see e.g. Refs. [31][32][33][34][35][36]). Therefore we define the contour-improved moments [37] as (see Fig. 1), and we employ the following path-dependent µ c α , first used in Ref. [37] (µ q α ) 2 (s, . (2.6) It weights in a different way the threshold versus the high energy parts of the spectrum. It was shown in Ref. [6] that the resulting momentsM X,C n can be obtained analytically from the Taylor expansion around s = 0 of the vacuum polarization function using an s-dependent µ c α (s, m 2 q ): This trick works because α s (µ c α (s, m 2 q )) has the same cut as the fixed-order expression for Π X . Other choices could spoil this property. Expanding the analytic expression for M X,C n on α s at a given finite order, one recovers the fixed-order momentsM X n . This shows that the dependence on the contour is only residual and represents an effect from higher order terms from beyond the order one employs for the calculation.
The contour-improved moments have a residual sensitivity to the value of the vacuum polarization function at zero momentum transfer. 1 For the case of the vector correlator this value depends on the UV-subtraction scheme and corresponds to Π(0) =M V 0 . For the case of the pseudoscalar correlator,M P 0 is scheme-independent, since P (q 2 ) already includes two UV subtractions. However, one could as well define a three-times-subtracted pseudoscalar correlator of the form P (q 2 ) = P (q 2 ) − P (0). Slightly abusing notation, we denote P as the "on-shell" scheme for P (q 2 ), and the twice subtracted (original) definition as the MS scheme for P (q 2 ). Using the OS scheme with Π X (0) = 0 for either vector or pseudoscalar correlator, we find that the first moment for the contour-improved expansion gives exactly the first fixed-order moment,M X,C 1 =M X 1 . Thus, in order to implement a non-trivial modification, and following Ref. [6], we employ the MS scheme for Π V (0) defined for µ = m q (m q ), and the twice-subtracted expression for P (q 2 ). Generically it can be written as The numerical values for the coefficients [C X ] a,b 0,i are collected in Table 7 for the vector correlator with 5 flavors and the pseudoscalar correlator with 4 flavors. In Table 4 or Ref. [6] one finds the the numerical values of [C V (n f = 4)] a,b 0,i .

Gluon Condensate Contribution
We estimate nonperturbative power corrections by including the gluon condensate contribution. The gluon condensate is a dimension-4 matrix element and gives the leading power correction in the OPE for the moments [38,39] Here the ellipses represent higher-order power corrections of the OPE involving condensates with dimensions bigger than 4. The Wilson coefficients of the gluon condensate corrections are known to O(α s ) accuracy [25]. Following Ref. [40], we express the Wilson coefficient of the gluon condensate in terms of the pole mass, since in this way the correction is numerically more stable for higher moments. However, as we did in Ref. [6], we still write the pole mass in terms of the MS quark mass at one loop. The resulting expression reads We use the renormalization group invariant (RGI) scheme for the gluon condensate [41]. The numerical value of the [a V (n f = 5)] a n and [a P (n f = 4)] a n coefficients are collected in Table 6. The values for [a V (n f = 4)] a n can be found in Table 5 of Ref. [6]. For methods (b) and (c) one can obtain the gluon condensate contribution by performing simple algebra operations and re-expansions in α

Ratios of Moments
An alternative set of observables which are also highly sensitive to the quark masses are the ratios of consecutive moments. To that end we define R X n (n f ) ≡ M X n+1 (n f )/M X n (n f ). Such ratios are proportional to the inverse square of the quark mass for any value of n. Their perturbative series can be expressed as an expansion in powers of α . Their computation is trivial, as one only needs to take the ratio of the two consecutive theoretical moments and re-expand as a series in α (n f ) s . We call this the standard fixedorder expansion, analogous to the expansion (a) of Sec. 2.1. The numerical expressions for the [R X (n f )] i,j a,b coefficients for the vector correlator with n f = 4, 5 are given in Table 14, and for the pseudoscalar correlator with n f = 4 in Table 15. As for the regular moments, we find that the nonlinear dependence of R X n on the quark mass sometimes causes that there is no numerical solution for m q .
By taking the square root of the ratio of two consecutive moments one gets a linear dependence on the inverse of the quark mass. The corresponding theoretical expression is obtained by re-expanding the perturbative expansion of R X n (n f ) as a series in powers of α (n f ) s . Thus we obtain an expression of the form of Eq. (2.3) with the replacement . This is referred to as the linearized expansion, in analogy to the expansion (b) of Sec. 2.1. The numerical values for the [R X (n f )] i,j a,b coefficients are collected for the vector correlator with n f = 4, 5 in Table 16, and for the pseudoscalar correlator with n f = 4 in Table 17.
Finally, one can use R X n (n f ) to solve for m q (µ m ) in an iterative way, exactly as explained in Sec. 2.1 for expansion (c). The theoretical expression is analogous to Eq.
We collect the numerical values for the [R X (n f )] i,j a,b coefficients, in Tabs. 18 and 19 for the vector (n f = 4, 5) and pseudoscalar (n f = 4) correlators, respectively. We call this the iterative linearized expansion.
One cannot implement a contour-improved expression for the ratios of moments, as they cannot be computed as the contour integral of a correlator. For the ratios of moments, in any of the three expansions, one can include non-perturbative corrections in the form of a gluon condensate OPE term, just using Eq. (2.10) and performing simple algebra operations and re-expansions in α (n f ) s and G 2 .

Previous Results and Scale Variations
In a number of recent low-n sum rule analyses [8,17,[26][27][28][29][30]44], which determined the charm and bottom quark masses with very small uncertainties using O(α 3 s ) theoretical computations for the moments [21][22][23][24], the theory uncertainties from the truncation of the perturbative series have been estimated with the scale setting µ m = µ α based on just one type of expansion, which was either the fixed-order [expansion (a)] for the vector correlator, or the linearized [expansion (b)] for the pseudoscalar correlator. In Ref. [6] we analyzed the perturbative series for the moments M V 1,2,3,4 of the charm vector correlator at O(α 3 s ) using an alternative way to estimate the perturbative uncertainties, based on the four different expansion methods (a)-(d), as explained in Sec. 2.1. We also focused on the question whether renormalization scale variation restricted to µ m = µ α leads to a compatible estimate of the perturbative uncertainties. From our analysis we found: • The extractions for the MS charm mass using the expansions (a)-(d) with correlated variations of µ m and µ α (e.g. µ m = µ α ) for the vector correlator can lead to very small scale variations, which are not compatible with each other. 2 Moreover, for some expansions also the results from the different orders can be incompatible to each other. It was therefore concluded that using correlated scale variation and one type of expansion can lead to an underestimate of the perturbative uncertainty.
• Uncorrelated (i.e. independent) variation of µ m and µ α leads to charm mass extractions with perturbative uncertainty estimates that are in general larger, but fully compatible among the expansions (a)-(d) and for the different orders. It was therefore concluded that µ m and µ α should be varied independently to obtain a reliable estimate of the perturbative uncertainty.
• The size of the charm mass perturbative uncertainty has a significant dependence on the value of the lower bound of the range of the scale variation. The choice of the upper bound has a marginal impact.
• The pattern of size of the correlated scale variations for the different expansions can be traced back to the form of the contours of constant charm mass in the µ mµ α plane, which happens to be located along the diagonal µ m ∼ µ α for expansions (a) and (b), but roughly orthogonal for expansions (c) and (d), see Fig. 6 of Ref. [6].
For example, 3 in Ref. [27] method (a) has been used for M V 1 with µ m = µ α varied between 2 and 4 GeV, quoting a perturbative error estimate of 2 MeV. For the expansion methods 4 (a)-(d) we obtain for the same scale variation 1.2893 ± 0.0007, 1.2904 ± 0.0004, 1.2963 ± 0.0045 and 1.3009 ± 0.0020 GeV, respectively, for m c (m c ), which are inconsistent. This can be compared to the corresponding results using independent variations as suggested in Ref. [6]. Using 2 GeV ≤ µ α , µ m ≤ 4 GeV we obtain 1.291 ± 0.003, 1.291 ± 0.003, 1.296 ± 0.005 and 1.302 ± 0.003 GeV, respectively, for expansions (a)-(d). These results are not consistent either. It was furthermore argued in Ref. [6] that an adequate variation range should include the charm mass itself (after all, that is the scale that governs the series), motivated by the range 2 m c ± m c around the pair production threshold. Thus, adopting independent scale variation in the range m c (m c ) ≤ µ m , µ α ≤ 4 GeV one obtains 1.287 ± 0.018, 1.287 ± 0.015, 1.282 ± 0.019 and 1.291 ± 0.014 GeV respectively. The results show consistency and demonstrate the strong dependence on the lower bound of the renormalization scale variation. The outcome is illustrated graphically in Fig. 4(a), and the order-by-order dependence in Fig. 1 of Ref. [45]. In Ref. [6] we also explored scale setting in which µ m was fixed to m c (m c ) and only µ α was varied. The outcome is shown in Figs. 4 and 5 of that reference. The contour lines in the µ mµ α plane for the mass extraction from the first moment of the vector correlator for all methods are shown in Fig. 6 of Ref. [6]. The final result quoted in Ref. [6], using α s (m Z ) = 0.1184 ± 0.0021, was 3 The size of the scale variations quoted in this paragraph applies to mc(mc) as well as to mc(3 GeV), and all numerical results are obtained at O(α 3 s ). We also stress that there are no perturbative instabilities concerning the use of the RGE down to the scale mc(mc). 4 To compute the charm and bottom masses in this section we use M V, exp 1 = 0.2121 GeV −2 for the charm vector correlator, result obtained in Ref. [6], M P, latt 1 = 0.1402 GeV −2 for the charm pseudoscalar correlator, from Ref. [8], and our own computation M V, exp 2 = 2.858 × 10 −5 GeV −4 for the bottom vector correlator, see Sec. 6. We also use αs(mZ ) = 0.1184.    We have repeated this analysis for the first moment of the pseudoscalar correlator M P 1 . Ref. [8] uses method (b) with the same scale variation as Ref. [27], quoting 4 MeV for the truncation error. For methods (a)-(d) and using 2 GeV ≤ µ m = µ α ≤ 4 GeV we obtain 1.276 ± 0.003, 1.277 ± 0.004, 1.275 ± 0.005 and 1.297 ± 0.004 GeV, respectively. For independent double scale variation between 2 and 4 GeV we obtain, 1.276 ± 0.013, 1.277±0.012, 1.271±0.012 and 1.294±0.012 GeV, and if we use m c (m c ) as the lower bound to we obtain 1.260±0.039, 1.267±0.037 GeV, 1.259±0.041 and 1.272±0.034. These results are displayed graphically in Fig. 4(b). The contour lines for the mass extraction from the first moment of the pseudoscalar correlator for all methods are shown in Fig. 2. We see that the results show a qualitative agreement with the situation for the vector current, but at a level of perturbative scale variations that are in general roughly larger by a factor of two.
A similar study can be performed for the extraction of the bottom mass m b (m b ) from the second moment of the vector correlator M V 2 . Ref. [27] uses the fixed-order expansion [method (a)] and correlated scale variation between 5 GeV ≤ µ m = µ α ≤ 15 GeV, quoting a perturbative error of 3 MeV. For the same variation we obtain 4.1691 ± 0.0005, 4.1682 ± 0.0015, 4.1729 ± 0.0034 and 4.1702 ± 0.0044 GeV for methods (a)-(d), respectively. As in the charm case the results are not consistent, but the variations of the results have a much smaller size, as is expected from the fact that for the bottom the renormalization scales are much larger. For independent variation between the same values we get 4.174 ± 0.008, 4.173 ± 0.006, 4.172 ± 0.006 and 4.177 ± 0.013 GeV. Finally, if the lower limit of the double variation starts at m b (m b ) we find 4.170 ± 0.011, 4.173 ± 0.011, 4.167 ± 0.011 and 4.176 ± 0.0015 GeV. These results are collected in Fig. 4(c). The corresponding m b (m b ) contours in the µ mµ α plane are shown in Fig. 3. As for the charm case, we find fully consistent results for the independent scale variation and using m b (m b ) as the lower bound.
We have also studied the ratio of the first over the second moments for the three cases, and observe a very similar pattern. We do not provide a detailed discussion in the text, but display the outcome graphically in Figs. 4(d) to 4(f).

Convergence Test
At this point it is useful to consider the perturbative series for all choices of µ α and µ m as different perturbative expansions, which can have different convergence properties. To estimate the perturbative uncertainties one analyzes the outcome of this set of (truncated) series. While the uncorrelated scale variation certainly is a conservative method, one possible concern is that it might lead to an overestimate of the size of the perturbative error. For instance, this might arise for a non-vanishing value of ln(µ m /µ α ) in connection with sizeable values of α s (µ α ) for µ α close to the charm mass scale, which might artificially spoil the convergence of the expansion. One possible resolution might be to simply reduce the range of scale variation (such as increasing the lower bound). However, this does not resolve the issue, since the resulting smaller variation merely represents a matter of choice. Furthermore, there is in general no guarantee that the series which are left have a better convergence despite the fact that the overall scale variation might become reduced. Preferably, the issue should be fixed from inherent properties of the perturbative series themselves. It is possible to address this issue by supplementing the uncorrelated scale variation method with a convergence test constraint, which we explain in the following.
We implement a finite-order version of the root convergence test. Let us recall that in mathematics, the root test (also known as Cauchy's radical test) states that for a series of terms a n , S[a] = n a n , if the quantity V ∞ defined as is smaller (bigger) than 1, the sum is absolutely convergent (divergent). If V ∞ approaches 1 from above then the series is still divergent, otherwise the test is not conclusive. In Eq. (4.1) lim sup stands for the superior limit, which essentially means that in case of oscillating series, one takes the maximum value of the oscillation. In the context of our analysis with truncated series the relevant property is that a smaller V ∞ implies a better convergent series. For the different expansion methods we use, it is simplest to apply the method directly to the sequence of quark masses that are extracted order by order, rewriting the results as a series expansion. Since we only know a finite number of coefficients of the perturbative series, we need to adapt the test. We now introduce V c and proceed as follows: 5 (a) For each pair of renormalization scales (µ m , µ α ) we define the convergence parameter V c from the charm mass series m c (m c ) = m (0) + δm (1) + δm (2) + δm (3) resulting from the extractions at O(α 0,1,2,3 s ): The resulting distribution for V c values can be conveniently cast as a histogram, and the resulting distribution is a measure for the overall convergence of the perturbative expansion being employed. We apply the convergence analysis to the region If the distribution is peaked around the average V c it has a well-defined convergence.
Hence discarding series with V c V c (particularly if they significantly enlarge the estimate of the perturbative error) is justified. , the average V c values are significantly larger than for the vector correlator, indicating that the overall perturbative convergence for the pseudoscalar moment is still excellent but worse than for the vector moment. This means that the vector correlator method is superior, and we expect that the perturbative uncertainty in the charm mass from the pseudoscalar is larger. This expectation is indeed confirmed as we discussed in Sec. 3, see also Sec. 9. Fig. 6(b) shows that the effect of discarding the series with the worst convergence is very similar to that of the vector correlator. For our determination of the bottom mass we use the second moment M V 2 (see Sec. 6 for a discussion on why we discard the first moment), and employ uncorrelated scale variations in the range m b (m b ) ≤ µ m , µ α ≤ 15 GeV. Fig. 5(c) shows the corresponding histograms, and we find that the convergence test yields V c double = (0.13, 0.11, 0.12, 0.15) for expansions (a)-(d) [for the correlated variation with scales set equal and 5 GeV ≤ µ α = µ m ≤ 15 GeV we find V c corr = (0.13, 0.09, 0.13, 0.15)]. As expected, the averages for the bottom are much smaller than for the charm. We further find that discarding series with the highest V c values only has minor effects on the perturbative error estimate for fractions up to 5%, see Fig 6(c). This is a confirmation that the series for bottom moments overall are more stable, which is again expected from the fact that perturbation theory should work better for the bottom than for the lighter charm.
The behavior of the ratios of moments is very similar as that for regular moments, as can be seen in Figs the ratios of moments is in general terms a bit worse than that of regular moments, except for the linearized iterative method of the pseudoscalar ratios.
In our final numerical analyses we discard 3% of the series with the worst V c values. As can be seen from Fig. 5, this only affects series with V c values much larger than the average values for the whole set of series. It is our intention to keep the fraction of discarded series as small as possible, since it is our aim to remove only series with convergence properties that are obviously much worse than those of the bulk of the series. We call this procedure trimming in the following. As we see in Figs

Lattice Simulation Data
The pseudoscalar current is not realized in nature in a way which is useful to compute the moments of the corresponding correlator from experimental data. Results for the moments of the pseudoscalar current correlator can, however, be obtained from simulations on the lattice. The strategy of these numerical simulations is to tune the lattice parameters (such as bare coupling constant and masses) to a selected number of observables (e.g. the energy splitting of Υ resonances). Once this tuning is performed, the lattice action is fully specified and no further changes are implemented. The tuned lattice action can then be used to perform all sorts of predictions, moments of correlators among them. Lattice simulations have to face a number of challenges, which usually translate into sizeable uncertainties. Among those are the extrapolations to the infinite volume and the zero lattice spacing (the latter being much harder), as well as the extrapolation to physical quark masses. On top of these systematic uncertainties, there are also statistical ones, which are related to the finite sampling used to perform the path integrations. On the other hand there are also concerns on the type of lattice action which is being used for the fermions. According to Ref. [8], the moments of the pseudoscalar correlator are least affected by systematic uncertainties, and so HPQCD has focused on those for their subsequent high-precision analyses.
To the best of our knowledge, HPQCD is the only lattice collaboration which has published results on QCD correlators. They have used the staggered-quarks lattice action, and MILC configurations for gluons. These results have been used to determine the charm mass and the strong coupling constant [8,28,29] with high accuracy, as well as the bottom mass [28,29], with smaller precision. We will use the simulation results as given in Ref. [8], even if the results quoted in of [28,29] are a bit more precise. The reason for this choice is that while [8] makes a straightforward extrapolation to the continuum, which is independent of the charm mass and α s extractions, in [28,29] the fit for the quark masses, the strong coupling constant and the extrapolation to the continuum is performed all at once, in a single fit. Furthermore that fits contains a lot of priors for the parameters one is interested in fitting. In any case, as we have seen, the charm mass extraction from the pseudoscalar correlator is dominated by perturbative uncertainties, as a result of the bad convergence of the series expansion for its moments.
Ref. [8] provides simulation results for the so-called reduced moments R k , which are collected in their Table II. The index k takes only even values, and starts with the value k = 4, which is fairly insensitive to the charm mass. Hence the lowest moment we consider is R 6 . Reduced moments are defined as (up to a global power) the full moment divided by the tree-level result. By taking this ratio, the authors of Ref. [8] claim that large cancellations between systematic errors take place. The reduced moments are scaleless, and the massdimension that one obviously needs to determinate the charm quark mass is regained by dividing with the mass of the η c pseudoscalar particle. Thus one can easily translate the reduced moments into the more familiar correlator moments M P n with the following relation: where the C P coefficients correspond to the tree-level terms of the standard fixed-order expansion of Eq. (2.2). Although the experimental value for η c is 2.9836(7) GeV, we use the value m ηc = 2.980 GeV given in Ref. [8], in order to ease comparison with that analysis. In Ref. [29] the value m ηc = 2.9863 (27) GeV is used. It is claimed that (as for the lattice action) it has no QED effects, and the error accounts for cc annihilation. Using the quote in Ref. [29] changes M P 1 by 0.4% and the effect on the charm mass is of the order of 2 MeV. The uncertainty in the η c mass has no effect on the M P n errors. In Table 1 we quote the lattice simulation results written as regular moments M P n .

Computation of the Experimental Moments for the Bottom Correlator
In this section we present our computation of the moments for the bottom vector current correlator from experimental data. These are made of three distinct contributions: the narrow resonances below threshold, the region of broader resonances, explored by BABAR [46], and the continuum region, where no data has been taken and some modeling is required.
The BABAR data has to be corrected for initial-state radiation. In the continuum region we use perturbation theory as a model for missing experimental data, and assign a conservative uncertainty by comparing the pQCD prediction with BABAR data in the region with measurements with the highest energy. Our determination of the experimental bottom moments differs from Ref. [27] in the way we model the uncertainties for the hadronic cross section in the continuum region, plus other minor differences in the contributions from the narrow widths and the threshold region, see discussion in Sec. 7. We also provide the correlation matrix among different moments, which cannot be found in the literature. Therefore, even though there are some similarities with the computations outlined in [27], we find it justified to discuss our computation of the experimental moments in some detail. We note that our results for the moments of the bottom vector current correlator have already been used in the analyses of high-n moments of Ref. [47], in the context of nonrelativistic large-n sum rules.

Narrow Resonances
The contribution of resonances below the open bottom threshold √ s = 10.62 GeV includes Υ(1S) up to Υ(4S). We use the narrow width approximation to compute their contribution to the experimental moments, finding (6.1) The masses and electronic widths of these four resonances are taken from the PDG [48], and the values of the effective electromagnetic coupling constant evaluated at the Υ masses are taken from Ref. [44]. This information is collected in Table 2. We have also checked that if one uses a Breit-Wigner instead of the narrow width approximation the results change by an amount well within the error due to the uncertainty in the electronic width.
In analogy to what we found in our study of the charm moments [6], the effect of the mass uncertainty in the moments is negligible. Therefore one only needs to consider the experimental uncertainty in the electronic widths. There is no information on the correlation between the measurements of these widths. The PDG averages of the electronic partial widths for the first three resonances is dominated by the CLEO measurement [49]. 7 Therefore we take the approach that half of the width's uncertainty (in quadrature) is uncorrelated (therefore mainly of statistical origin), whereas the other half is correlated among the various resonances (therefore coming from common systematics in the measurements).

Threshold Region
The region between the open bottom threshold and the experimental measurement of the R b -ratio at the highest energy, 10.62 GeV ≤ √ s ≤ 11.2062 GeV, is referred to as the threshold region. The region above the last experimental measurement will be collectively (26)  denoted as the continuum region. The first experimental data close to the B meson threshold were taken by the CLEO [51][52][53] and CUSB [54] collaborations. The measurements at each c.m. energy have a 6% systematic uncertainty. More recently the BABAR collaboration [46] has measured the R b -ratio in the energy region between 10.54 GeV and 11.20 GeV, with significantly higher statistics and better control of systematic uncertainties (of the order of 3%). These measurements are taken in small energy bins, densely populating the threshold region. The BABAR data supersedes the older data of CLEO and CUSB, and it has already been used in Refs. [27,30,40], in which the bottom mass was also determined. This BABAR data for the R b -ratio has not been corrected for initial-state radiation effects. Moreover, the effect of the Υ(4S) resonance has not been subtracted, 8 so we have performed the subtraction ourselves, using the Breit-Wigner approximation and using for the total width the PDG value Γ 4S = 20.5 MeV: For the subtraction of the Υ(4S) resonance and the correction for the initial state radiation we take an approach similar to Ref. [27].

Subtraction of the Υ(4S) Radiative Tail
The contribution to be subtracted from the BABAR data is the ISR-distorted tail of the Υ(4S), which reaches to energies above its mass. The cross section R and the ISR-distorted cross-sectionR are related by a convolution relation which can be used to determine the ISR effects on the Υ(4S) resonance given in Eq. (6.2). Here the lower integration bound is z 0 = (10 GeV) 2 /s. This value is not fully fixed by theoretical arguments, and it is chosen such that it excludes the narrow resonances, but keeps the major part of the Υ(4S) line shape. The radiator function G is given as [55,56] where the specific form of β, F and the two δ's can be found in Eq. (7) of [27]. Note that the function G(z, s) is divergent as z → 1, but since 0 < β − 1 < −1, it is integrable. The divergent behavior is absent in G, which in the limit z → 1 reduces to After subtracting the radiative tail of the Υ(4S) we find that to a good approximation the cross section vanishes for energies below 10.62 GeV. Therefore we add an additional point to our BABAR dataset: R b (10.62 GeV) = 0 and take R b = 0 for energies below 10.62 GeV. Since the subtracted cross section does not exactly vanish between 10.5408 GeV and 10.62 GeV, we take the (small) contribution of the subtracted cross section in that region to the moments as an additional source of systematic correlated uncertainty.

Deconvolution of Initial-State Radiation
After subtraction of the radiative tails, the BABAR threshold data are corrected for ISR. The inversion of the convolution in Eq. (6.3), can be carried out in an iterative way [27]. Defining δG(z, s) = G(z, s) − δ(1 − z) one can use a successive series of approximations where we denoted the j-th approximation of R(s) as R j (s) and use as starting point R 0 (s) =R(s), the BABAR data. In Eq. (6.6) we take z 0 = (10.62 GeV) 2 /s, using as a starting point the energy value for which the cross section vanishes after the subtraction of the radiative tails. To isolate the singularity at the higher endpoint one can perform a subtraction at z = 1, resulting in: We use the trapezoidal rule to evaluate the integration on the discrete set of experimental data measurements labeled by the index i. Changing the integration variable from z to energy we find where we have used R j 1 = R(10.62 GeV) ≡ 0 for all iterations. After applying the procedure as many times as necessary to obtain a stable solution, one obtains the ISRcorrected cross section. Among the experimental measurements one finds two data points GeV. It turns out the fact that they lie very close makes the iterative procedure unstable. Therefore we drop the latter point from our analysis. In Fig. 8 we show the BABAR after the subtraction of all radiative tails, before (red) and after (blue) the ISR deconvolution.

Determination of the Unfolding Error Matrix
The BABAR collaboration splits the experimental uncertainties into statistical, systematic uncorrelated, and systematic correlated. We add the two former in quadrature to obtain the total uncorrelated uncertainty uncor and rename the latter as the total correlated uncertainty cor . The removal of the radiative tails of the Υ mesons has no effect on these uncertainties. Therefore, the correlation matrix for the BABAR data after the subtraction of the radiative tails, before it is corrected for ISR effects, can be written as One needs to compute a new correlation matrix after each iteration. In this way we determine the unfolding error matrix.
The master formula in Eq. (6.8) can be cast in a matrix form as follows: where R j i is to be thought as the i-th component of the column vector R j , and G ik represents the (i, k)-component of a matrix G. Here R j i depends only on the initial value R 0 i and the result of the previous iteration R j−1 i . The G ik do not depend on R j k or the iteration step j. Therefore, for the error propagation one uses both of them j-independent. We will denote with M j j the correlation matrix among the entries of the vector R j for a given iteration j. We also find it convenient to introduce the correlation matrix among R j and R 0 , referred to as M j 0 . Finally we use the notation M 0 j ≡ (M j 0 ) T . We find for the correlation matrix after j iterations: were the elements of the matrix M 0 0 are given in Eq. (6.9). We find that after five iterations the result has converged already to a level well below the experimental uncertainties.

Contribution of the Threshold Region
After having corrected BABAR data for ISR effects, we use the trapezoidal rule for integrating the threshold region between 10.62 GeV and 11.20 GeV: where R i has been already ISR corrected and N is the number of data points. We have added the boundary condition point R 1 = R(10.62) = 0. From Eq. (6.13) one can compute the correlation matrix among M thr n for various n values, using the unfolding matrix among the R i computed in Sec. 6.2.3.

Continuum Region
For the determination of the experimental moments from the region above 11.2 GeV we use pQCD (which has essentially negligible errors) supplemented by a modeling uncertainty. Comparing pQCD and re-binned data (with 25 MeV bins) in the region between 11.06 GeV and 11.2 GeV we find a 4% discrepancy concerning the central values. This can be seen in Fig. 9. Given that the relative discrepancy between experiment and pQCD for R b at the Z-pole is about 0.3% [57], we adopt a relative modeling error that decreases linearly from 4% at 11.2 GeV to 0.3% at m Z , and stays constant for energies larger than m Z . This uncertainty makes up for 96.9% of the total error for the first moment M V 1 (which has an total 2.45% relative error), and 86.27% of the second moment M V 2 (which has a total 1.83% relative error). Note that if we would adopt a constant 4% error for all energies above R b Figure 9. Comparison of ISR-corrected BABAR data in the continuum region (black dots with error bars) with pQCD (purple band). Since the experimental measurements are very dense in this region, we perform an average rebinning, shown as the red band.
11.2 GeV, this continuum uncertainty would make up for 97.25% of the total error for the first moment M V 1 (from a total 2.59% relative error), and 86.57% of the second moment M V 2 (from a total 1.86% relative error). The difference is small because contributions from higher energies are suppressed. Following our procedure in Ref. [6] we consider this uncertainty as fully correlated for the various moments, but without any correlation to the narrow resonances or the threshold region.
The perturbative QCD theoretical expression which we use to determine this contribution includes the non-singlet massless quark cross section supplemented with bottom mass corrections up to O(m 4 b )/s 2 ). It reads: [58][59][60][61][62] where  Table 3. Results for our computations of the experimental moments. The second column collects the contribution from the first four Υ resonances (using the narrow width approximation). The third to fifth columns show the contributions from the threshold (using ISR-corrected BABAR data) and continuum (using pQCD as a model for the lack of data) regions, and the total moment determinations, respectively. The two numbers quoted in parentheses correspond to the uncorrelated and correlated experimental uncertainties, respectively. All numbers are given in units of We use the initial conditions m b (m b ) = 4.2 GeV and α s (m Z ) = 0.118. Therefore, for the continuum region we use the following expression, with s 0 = (11.2062 GeV) 2 . Here γ is the auxiliary variable used to parametrize our uncertainty, which we consider as 100% correlated among the various moments. The related entries of the correlation matrix are trivially computed as

Final Results for the Experimental Moments
The full result for the experimental moments is obtained by summing up all the portions described before, We determine two correlation matrices among the first four moments. One of them comes from the various uncorrelated uncertainties, whereas the other encodes the systematic uncertainties. We denote them as the correlated and uncorrelated correlation matrices, respectively. These are computed by summing up the respective individual matrices from each region, and in the same way as we did for our charm analysis [6], we assume there is no region-to-region correlation. We find: where the (n, m) entry of each matrix is given in units of 10 −2(n+m+1) GeV −2(n+m) , and the total correlation matrix is the sum of C exp uc and C exp cor . The contribution of each region to the final experimental moments and the corresponding uncertainties are presented in Table 3.

Comparison to other Determinations of the Experimental Moments
In this section we compare our result for the experimental moments for the bottom vector current correlator with previous determinations. These are collected in Table 4. 9 The most relevant comparison is between the second and third columns, where the most recent data on the narrow resonances and the BABAR continuum data are used. For the contributions from the narrow resonances we have a perfect agreement with [27], although slightly larger errors. For the threshold region our results are slightly bigger, and our uncertainties are almost identical. The main difference between these two determinations is the estimate of the uncertainties coming from the continuum region, where the pQCD prediction for the R b -ratio is employed. Whereas we adopt the more conservative approach described in Sec. 6.3, Ref. [27] employs only the perturbative uncertainties related to the purple band in Fig. 9. In Ref. [40] the same collaboration presents a more critical analysis of their errors. In particular they observe that the last experimental measurement of BABAR, after being corrected for ISR, disagrees with the pQCD prediction at the 20% level (way outside the corresponding uncertainties). 10 To resolve this discrepancy they assume two possible scenarios: a) pQCD starts being reliable at energies above 13 GeV (therefore the authors interpolate between the last experimental point and pQCD at 13 GeV); b) BABAR systematic errors have been underestimated (therefore the central values of the experimental measurements are rescaled by a factor of 1.21). Ref. [40] quotes the values of the experimental moments and the resulting values for the bottom mass for these two scenarios. Since the effect of this difference on the bottom mass determination from M V 2 is only slightly higher than the size of the other uncertainties (that is, the uncertainties of the theoretical moments plus the other experimental errors) added together, it is argued that this issue can be ignored. We disagree with this argument, since the rest of the uncertainties are 100% correlated in each of the scenarios and drop out when comparing to one another. Therefore this shift 9 In the case of Ref. [7], we reconstruct the experimental moments from their Table 3, where the moments are split in several different contributions. For the reconstructed uncertainty, we take one half of the error of the narrow resonances correlated to each other, and the other half as uncorrelated. The errors from patches where theory input is used are taken as fully correlated to one another. The total narrow-resonance error, and the total "theory-patch" error are added in quadrature to get the final uncertainty. 10 From our own computation of the ISR-corrected R b -ratio, we do not observe such a large deviation between the last data point and the pQCD prediction, see Fig. 9.
must be taken as an additional source of error on the experimental moments (and indeed would then dominate the corresponding total error). The additional error (to be added in quadrature to the one in round brackets) is quoted in square brackets in the second column of Table 4. It amounts to an additional error of 30, 18, 11 and 7 MeV for m b (m b ) extracted from moments M V 1 to M V 4 , respectively. Refs. [7,17,44] have used the older CLEO and CUSB experimental measurements, resulting in relatively large uncertainties. In Ref. [44] the CLEO measurements are divided by a factor of 1.28, and an error of 10% is assigned. It is argued that this procedure is necessary to reconcile old and new CLEO measurements, as well as to improve the agreement with pQCD predictions. Ref. [7] uses values for the Υ-states electronic partial widths given by the PDG 2002, which have larger uncertainties. This makes their determination of the experimental moments rather imprecise.
Concerning the continuum region where no measurements exist, while some previous analyses have taken a less conservative approach than ours, in Ref. [63] a much more conservative approach is adopted. In this region they consider the R b -ratio as constant with a 66% uncertainty. In Ref. [7] also a more conservative approach is adopted. Between 11.1 and 12 GeV O(α 2 s ) pQCD errors are used, which are larger than 10%; for energies above 12 GeV a global 10% correlated error is assigned.  Table 4. Comparison of our results for the experimental moments of the bottom vector correlator (2nd column) to previous determinations (3rd to 5th columns). The 2nd and 3rd columns use BABAR data from Ref. [46], while 4th and 5th use older data from Refs. [51,52]. The 3rd and 4th columns use perturbative uncertainties in the continuum region, while 2nd and 5th use a more conservative estimate based on the agreement of data and pQCD. In the 2nd column, we quote in square brackets our own estimate of an additional systematic error from the considerations made in Ref. [40]. All numbers are given in units of 10 −(2n+1) GeV −2n .  correlations is not provided in Ref. [8]. Therefore we make the simplest possible assumption, which is that the moments are fully uncorrelated. This will most certainly overestimate the uncertainties for the ratios of moments, but given that we are anyway dominated by perturbative uncertainties, our approach appears justified. We collect our results for the computation of the ratios of experimental moments in Table. 5. Readers interested in the full correlation matrix among them can send a request to the authors.

Results
In this section we present the final results for our analyses at O(α 3 s ). We take method (c) (linearized iterative expansion) as our default expansion. For the estimate of the perturbative uncertainty, we perform double scale variation in the ranges m c (m c ) ≤ µ α , µ m ≤ 4 GeV for charm (either correlator), and m b (m b ) ≤ µ α , µ m ≤ 15 GeV for bottom, and we discard 3% of the series with the worst convergence (that is, with highest values of the V c convergence parameter). For the charm mass determinations (either vector or pseudoscalar correlator) we use the first moment as our default, given that it is theoretically more reliable than the higher moments. For the analysis of the bottom mass from regular moments, we use M V 2 as our default, since it is less afflicted by systematic experimental errors than the first moment, and is nevertheless theoretically sound. For the charm and the bottom mass analyses we also examine the ratio of the second over the first moment as a cross check and validation of the results from regular moments. The results for the experimental moments are collected in: the last column of Table 9 in Ref. [6] (charm vector correlator regular moments); the last column of Table 3 (bottom vector correlator regular moments); Table 1 (lattice regular moments); and Table 5 (all ratios of moments).
We also analyzed higher (and also lower for the case of bottom) moments and ratios of moments for all correlator and quark species. Since, as already discussed, fixed-order and contour-improved higher moments are afflicted by a nonlinear dependence on the quark mass, we only consider the linearized and iterative methods for this analysis. In any case, since higher moments have a larger sensitivity to infrared effects and are therefore theoretically less sound, the analysis involving higher moments mainly aims at providing cross checks. The results are collected in a graphical form in Fig. 11, and the numerical results can be obtained from the authors upon request.
Our final determinations include nonperturbative effects through the gluon condensate including its Wilson coefficients at order O(α s ). Furthermore, we assign as a conservative estimate of the nonperturbative uncertainty twice the shift caused by including the gluon condensate. In any case, this error is very small, particularly for the bottom mass determination. One source of uncertainty which we have not discussed so far is that coming from the strong coupling constant. Although the world average α s (m Z ) = 0.1185 ± 0.006 has a very small error, see Ref. [48], one cannot ignore the fact that it is fairly dominated by lattice determinations, e.g. [29]. Furthermore, there are other precise determinations with lower central values and in disagreement with the world average from event-shapes [64][65][66][67] and DIS [68]. A review on recent α s determinations can be found in Refs. [69][70][71]. Therefore, in analogy with Ref. [6], we perform our analyses for several values of α s (m Z ) between 0.113 and 0.119, and provide the central values and perturbative errors as (approximate) linear functions of α s (m Z ). The other uncertainties are essentially α s -independent, so we just provide the average. We also quote quark mass results for α s taken from the world average: where we adopt an uncertainty 3.5 times larger than the current world average [48]. We note that in Ref. [6] we have taken α s (m Z ) = 0.1184 ± 0.0021 as an input which causes only tiny sub-MeV differences in the quark masses. We refrain ourselves from presenting the α s dependence of the higher-moment result, which the reader can get from the authors upon request.
For the numerical analyses that we carry out in this article we have created two completely independent codes: one using Mathematica [72] and another using Fortran [73], which is much faster and suitable for parallelized runs on computer clusters. The two codes agree for the extracted quark masses at the level of 10 −6 .

Results for the Charm Mass from the Vector Correlator
For the analysis using the first moment of the charm vector correlator we use the experimental value quoted in Eq. For the ratio of the second over the first moment of the vector correlator we use as the experimental input R V, exp 1 = (6.969 ± 0.032 stat ± 0.059 syst ) × 10 −2 GeV −2 , which yields the following result for the charm mass: m c (m c ) = 1.271 ± (0.003) stat ± (0.005) syst ± (0.016) pert (9.4) ± (0.004) αs ± (0.004) GG GeV .
With correlated variation 2 GeV ≤ µ α = µ m ≤ 4 GeV we get 1.258 ± (0.005) pert and 1.279 ± (0.007) pert for methods (a) and (c), respectively. In this case the scale variations are a factor of 2 to 3 smaller than our perturbative error estimate. 11 The α s dependence, which can be seen in Figs. 12(c) and 12(d), has the form: We observe that the central value for the ratios of moments is 17 MeV smaller than for the first moment analysis, but fully compatible within theoretical uncertainties. Furthermore, the dependence on α s of the central value obtained from the regular moment analysis is larger, which translates into a corresponding larger error due to the uncertainty in α s . Both determinations have very similar perturbative uncertainties for any value of α s . We also see that the charm mass from the ratio of moments has smaller experimental uncertainties, as a result of cancellations between correlated errors. Moreover, the charm mass result from R V 1 has a nonperturbative error twice as large as that from M V 1 . The two charm mass results from the first moment and the moment ratio are compared graphically in Fig. 15(a).

Results for the Charm Mass from the Pseudoscalar Correlator
For the analysis of the first moment of the charm pseudoscalar correlator we employ M P, latt 1 = (0.1402 ± 0.0020 latt ) GeV −2 [8], which yields the following charm mass determination: m c (m c ) = 1.267 ± (0.008) lat ± (0.035) pert ± (0.019) αs ± (0.002) GG GeV . which is also displayed in Figs. 13(a) and 13(b). As expected, the perturbative error is much larger than for the vector correlator, and has a stronger dependence on α s . We see that the central value has a much stronger dependence on α s as well, which again translates into a large error due to the uncertainty in the strong coupling. The central value is 21 MeV lower than Eq. (9.2), but fully compatible within errors [see Fig 15(a)]. The nonperturbative uncertainties are identical to the vector current case.
For the ratio of second over the first moment of the pseudoscalar correlator we use R P, latt 1 = (0.0971 ± 0.0032 latt ) GeV −2 . We find for the charm mass m c (m c ) = 1.266 ± (0.020) latt ± (0.018) pert ± (0.006) αs ± (0.002) GG GeV . The central values for both M P 1 and R P 1 are almost identical, but their α s dependence is not: the latter is much smaller (even smaller than for M V 1 , but larger than for R V 1 ). Note that the lattice error is larger for the ratio since we made the very conservative assumption that they are fully uncorrelated. This is because the correlation matrix for various lattice moments is unknown. The perturbative error reduces by a factor of two for any value of α s when using the ratio, but we have checked that this only happens for the iterative expansion. On the other hand, the α s dependence of the perturbative uncertainty is smaller for the regular moment determination. The nonperturbative errors are identical.
All charm determinations are illustrated graphically in Fig. 15(a), where in red we show our preferred determination from the first moment of the vector correlator. The perturbative error is 30% smaller than for the charm vector correlator analysis, as a result of the smaller value of α s at the scales close to the bottom mass. This is consistent with our discussion on the convergence properties of perturbation series for the bottom quark carried out in Sec. 4. The total error is dominated by the experimental systematic uncertainty, which in turn mainly comes from the continuum region where one relies on modeling in the absence of any experimental measurements. The nonperturbative error is completely negligible. This is expected since it is suppressed by two powers of the bottom mass. Using the correlated scale variation 5 GeV ≤ µ α = µ m ≤ 15 GeV for methods Although the central value for the ratio analysis is 7 MeV higher, this has no significance given the size of the uncertainties. The dependence of the central value on α s is three times smaller for the ratio analysis. The perturbative error and its α s dependence are roughly the same for the ratio and the single moment analysis. Moreover, the two experimental errors are very similar. This is because, even though there is some cancellation of correlated errors in the ratio, a significant part of the huge systematic error of the first moment remains uncanceled. A graphical illustration of the two bottom mass determinations is shown in Fig. 15(b). Both combined uncertainties and central values are rather similar, and we adopt the result from the second moment (in red) as our default result.

Comparison to other Determinations
In this section we make a comparison to previous analyses of our updated charm mass determination from the vector correlator, our new results for the charm mass from the pseudoscalar current correlator and of our bottom mass determination. We restrict our discussion to determinations which use QCD sum rules with infinite as well as finite energy range for the vector or pseudoscalar current correlators, and including relativistic and nonrelativistic versions of the sum rules. We do not cover charm mass determinations from DIS or bottom mass determinations from jets (which are in any case rather imprecise), as well as determinations which are based on the mass of bound states (B mesons or quarkonia) or B decays.
In Figs. 16(a) and 16(b) we present in a graphical form a compilation of recent sum rule determinations of the charm and bottom masses, respectively. We have labeled them from top to bottom with numbers from 1 to 14. We note that comparing these results, one has to keep in mind that different analyses in general employed different values and uncertainties for the strong coupling. Only the analyses in Refs. [6,27,37,47,74] and ours have provided the dependence of their results on the value of α s (m Z ).

Charm Mass
Let us first focus our attention on the charm mass, Fig. 16(a). Within each color, determinations are ordered according to publication date. In red (determinations 12 to 14) we show the results of our collaboration: 12 and 13 for the vector correlator, the former (dashed) corresponding to Ref. [6] without trimming procedure, and the latter (solid) corresponding to this work, which includes the trimming procedure. Determinations 12 to 14 are the only analyses using uncorrelated scale variation. Determination 6 (gray) [37] sets µ m = m c (m c ), and all the other analyses have set µ m = µ α . Determinations in blue (1 [29], 2 [28] and 3 [8]) were performed by the HPQCD collaboration, which employ method (b) for the pseudoscalar correlator moments used for the mass determination in their lattice analyses. Only 1 to 3 and 14 use pseudoscalar moments, while all the other analyses use the vector correlator. Among those 7 [26], 8 [27], 9 [44] and 10 [17] use data in the threshold region only up to 4.8 GeV; analysis 6 uses two patches of data in the threshold region, one from threshold to 4.7 GeV, and another between 7.2 and 11 GeV; analyses 12 and 13 use all available data (see Ref. [6] for the complete bibliographic information on charm data). Analysis 7 (orange) corresponds to weighted finite-energy sum rules. They employ a kernel which enhances the sensitivity to the charm mass, and at the same time reduces the sensitivity to the continuum region. Green color analyses, collected in 4 [43,75] and 5 [74], apply other kinds of sum rules. Analysis 5 uses a finite energy sum rule similar to 7, but the kernel makes the sensitivity to the charm mass quite small. On the other hand, the two determinations of 4 use shifted moments, ratios of shifted moments, and exponential sum rules, and consider only the contributions from the first 6 vector resonances in the narrow width approximation, and the pQCD prediction for the continuum. The lower analysis of 4 [75] includes contributions from condensates up to dimension 6; the higher [43] includes condensates up to dimension 8. In purple (11 [76]) we show the only analysis which uses large-n moments for the charm mass fits, employing NRQCD methods to sum up large logs and threshold enhanced perturbative corrections, supplemented with fixed order predictions for the formally power suppressed terms. This analysis uses contributions from narrow resonances, plus a crude model for the threshold and continuum patches, for which a conservative uncertainty is assigned. We note, however, that this analysis might be questioned since perturbative NRQCD is in general not applicable for the charmonium states.
Our new vector correlator result agrees well with the world average, having a similar uncertainty. Our result is fully compatible with the other determinations shown in Fig. 16(a). As mentioned already before, we disagree with the small perturbative uncertainties related to the scale variations of the vector and/or pseudoscalar moments adopted in analyses 1 to 3 and 7 to 10.

Bottom Mass
Let us now turn our attention to the bottom mass results, see Fig. 16(b). The coloring and chronological conventions are analogous to Fig. 16(a), and we try to keep a similar ordering. We show three nonrelativistic determinations (11 [63], 12 [47] and 13 [77]) in purple; O(α 2 s ) fixed-order analyses are shown in gray (5 [7]), black (10 [78]) and green (4 [79]); finite energy sum rules also based on fixed-order appear in orange (6 [30]) and green (4); there are two lattice analyses in blue, collected in 1 [28,29]. Analyses 3, 6, 7 and 11 to 13 use the new BABAR data, whereas the others use the older CLEO and CUSB data. Analyses 4 and 11 include only the contributions of the first six vector resonances. Analyses 4 and 5 use older measurements of the electronic width for the narrow resonances. Analyses 3, 4, 6 to 9 use pQCD in the high-energy spectrum for the experimental moments. The theoretical treatment of the bottom mass analyses in red, gray, black, blue and green are in complete analogy to their charm mass analyses: 3 and 4 for bottom correspond to 4 and 5 for charm, respectively; 1 in bottom corresponds to 1 and 2 for charm.
The upper analysis of 1 [29] uses a computation of m b /m c and their result for m c to make a prediction for the bottom mass. The analysis 10 uses a combination of M V 6 with the infinite momentum transfer moment, both in fixed-order, in order to constrain the continuum region. They only use experimental information on narrow resonances, and model the rest of the spectrum with theory predictions. Finally, they make the following scale choice: µ α = µ m = m b (m b ) and estimate the truncation error from an ansatz for the O(α 3 s ) term. 12 Analyses 12 and 13 use large-n moments and NRQCD methods for their theoretical moments. Analysis 12 uses NRQCD fixed-order perturbation theory at N 3 LO (which accounts for the summation of the Coulomb singularities) and 13 uses renormalization group improved perturbation theory in the framework of vNRQCD at N 2 LL order. Both analyses employ low-scale short-distance masses to avoid ambiguities related to the pole mass renormalon. Analysis 11 also uses N 3 LO NRQCD fixed-order input, but is incomplete concerning the contributions from the continuum region in the theoretical moments. Moreover, they extract the pole mass which is then converted to the MS scheme.
Our result 14 is in full agreement with the world average, having a slightly smaller uncertainty. It also agrees with the other analyses shown, with slightly smaller or comparable uncertainties. We disagree with the small perturbative uncertainties related to scale variations quoted in 6 to 10.

Conclusions
In this work we have determined the MS charm and bottom quark masses from quarkonium sum rules in the framework of the OPE, using O(α 3 s ) perturbative computations, plus nonperturbative effects from the gluon condensate including its Wilson coefficient at O(α s ). For the determination of the perturbative uncertainties we independently varied the renormalization scales of the strong coupling and the quark masses, in order to account for the variations due to different possible types of α s -expansions, as suggested earlier in Ref. [6].
In order to avoid a possible overestimate of the perturbative uncertainties, coming from the double scale variation in connection with a low scale of α s and resulting in badly convergent series, we have re-examined the charm mass determination from charmonium sum rules (vector correlator) supplementing the analysis with a convergence test. The convergence test is based on Cauchy's radical test, which is adapted to the situation in which only a few terms of the series are known, and quantifies the convergence rate of each series by the parameter V c . We find that the distribution of the convergence parameter V c coming from the complete set of series peaks around its mean value, and allows to quantify the overall convergence rate of the set of series for each moment in a meaningful way. This justifies discarding (or "trimming") series with values of V c much larger than the average. For our analysis we discard 3% of the series having the highest V c values, which results in a reduction of the perturbative uncertainty for the MS charm mass m c (m c ) from 19 MeV in [6]  where all sources of uncertainty have been added in quadrature. This result supersedes our corresponding earlier result from Ref. [6], which was 1.282 ± 0.024 GeV. This makes it clear that the trimming procedure discards series which produce small values of the charm mass.
We have applied the same method of theory uncertainty estimate to analyze the HPQCD lattice simulation results for the pseudoscalar correlator. Our convergence test signals that the pseudoscalar moments have far worse convergence than the corresponding vector ones. This translates into an uncertainties of 35 MeV due to the truncation of the perturbative series and the error in α s (roughly twice as big as for the vector determination). In contrast, using correlated scale variation (e.g. setting the scales in the mass and the strong coupling equal) the scale variation can be smaller by a factor of 8. Our new determination from the first moment (again being the most reliable theoretical prediction) of the pseudoscalar correlator reads m c (m c ) = 1.267±0.041 GeV, where again all individual errors have been added in quadrature. The combined total error is twice as big as for the vector correlator, and therefore we consider it as a validation of Eq. (11.1) in connection with the convergence test. The result is in sharp contrast with the analyses carried out by the HPQCD collaboration [8,28,29], where perturbative uncertainties of 4 MeV are claimed. We have checked that, as for the vector correlator, for the different possible types of α s -expansion the correlated variation in general leads to a bad order-by-order convergence of the charm mass determination.
The second important result of this work is the determination of the bottom quark mass from the vector correlator. We have reanalyzed the experimental moments by combining experimental measurements of the first four narrow resonances, the threshold region covered by BABAR, and a theoretical model for the continuum. This theoretical model is pQCD, to which we assign a 4% systematic uncertainty which decreases linearly to reach 0.3% at the Z-pole, and stays constant at 0.3% for higher energies. Our treatment is motivated by the discrepancy between experimental measurements and theoretical predictions in the energy range between 11.0 and 11.2 GeV and at the Z-pole. This results into a large error for the first moment, and therefore we choose the second moment (which is theoretically as clean as the first one for the case of the bottom quark) for our final analysis, giving a total experimental uncertainty of 18 MeV. Our treatment of the experimental continuum uncertainty is in contrast to Ref. [27], where instead the very small perturbative QCD uncertainties (less than 1%) are used, claiming an experimental uncertainties of 6 MeV. In the light of the analysis carried out here, supported by the observations made in Ref. [40], we believe this is not justified. Our convergence test reveals that, as expected for the heavier bottom quark, the perturbative series converge faster than for the charm quark. Correspondingly, the perturbative and α s uncertainties are ∼ 30% smaller than those for charm. Taking correlated scale variation as used in Refs. [27,30] the perturbative error estimate can shrink up to a factor of 20. We also find that correlated variation leads to incompatible results for the different types of α s -expansions. Our final result for the bottom mass from the second moment, with all errors added in quadrature, reads: were the total error is fairly dominated by the systematic error, which comes from the continuum region of the spectrum. Our uncertainty is very similar to the one obtained by the HPQCD analysis, but 30% larger than the 16 MeV claimed by [27]. Our central value is 4 MeV larger than the latter. This almost head-on agreement is a result of two effects that push in opposite directions: smaller value of the second experimental moment, and different perturbative analysis. Curiously enough, a similar accidental cancellation was observed for the charm mass in [6].
In order to further validate the results discussed above, we have also analyzed the ratios of consecutive moments of each one of the three correlators as alternative observables. In all cases the results from the moment ratios agree very well the regular moment analyses.

A Numerical Values for the Perturbative Coefficients
In this appendix we succinctly collect the numerical values for all of the coefficients appearing in the perturbative series, such that our analysis can be reproduced. We organize these values in tables, each of them corresponding to a different equation.  Table 6. Numerical values for the coefficients of Eq. (2.10) for the vector correlator with n f = 5 (first two columns), and for the pseudoscalar correlator with n f = 4 (last two columns). (Gluon condensate contribution).  Table 7. Numerical values of the coefficients for Eq. (2.8) for Π V (0) in the MS scheme and n f = 5, and P (q 2 = 0) for n f = 4.  Table 8. Numerical values of the coefficients for Eq. (2.2) for the vector current with n f = 5.
(Standard fixed-order expansion).  Table 9. Numerical values of the coefficients for Eq. (2.3) for the vector current with n f = 5.
(Linearized expansion).   Table 11. Numerical values of the coefficients for Eq. (2.2) for the pseudoscalar correlator with n f = 4. (Standard fixed-order expansion).        Table 15. Numerical values for the coefficients of the linearized expansion of the ratios of pseudoscalar moments with n f = 4.     Table 17. Numerical values for the coefficients of the linearized expansion of the ratios of pseudoscalar moments with n f = 4.     Table 19. Numerical values for the coefficients of the iterative linearized expansion of the ratios of pseudoscalar moments with n f = 4.