Bottom Quark Mass with Calibrated Uncertainty

We determine the bottom quark mass $\hat{m}_b$ from QCD sum rules of moments of the vector current correlator calculated in perturbative QCD to ${\cal O} (\hat\alpha_s^3)$. Our approach is based on the mutual consistency across a set of moments where experimental data are required for the resonance contributions only. Additional experimental information from the continuum region can then be used for stability tests and to assess the theoretical uncertainty. We find $\hat{m}_b(\hat{m}_b) = (4180.2 \pm 7.9)$ MeV for $\hat\alpha_s(M_Z) = 0.1182$.


Introduction
Highly precise values of the charm and bottom quark masses can be obtained in QCD perturbation theory, because they are sufficiently large to suppress non-perturbative effects. The object of interest is the vector current correlation function, which can be studied experimentally in a clean environment in electron-positron annihilation. Furthermore, by considering moments of the correlator one arrives at theoretically most accessible inclusive observables, which -at least in the case of the vector current -offers excellent perturbative convergence even in the context of the charm quark mass, m c , where the strong coupling, α s (m c ) ∼ 0.4, is not all that small. By the specific method [1] reviewed in Sect. 2, which is a concrete implementation of the general QCD sum rule idea [2,3], we were able to determine m c with a controlled theory uncertainty [4], competitive even with the results from lattice gauge theory simulations [5], and in excellent agreement with them [6].
In the case of the bottom quark mass, m b , lattice gauge theory faces an impediment, as the strong interaction scale of O(m ρ ) differs significantly from m b itself. By contrast, this separation of scales is a virtue in any approach effectively utilizing the operator product expansion (OPE). Together with the smaller value of the strong coupling, α s (m b ) ∼ 0.23, this turns the QCD sum rule approach into the method of choice to determine m b .
On the other hand, much less experimental information is available on the bottom quark current correlator compared to that of the charm quark. This is because the bb electroproduction cross section is by more than an order of magnitude smaller than non-bb quark production, so that B tagging is needed in order to determine the exclusive cross section. Furthermore, while formally the domain of the dispersion integration extends to infinite energy, the experimentally scanned kinematical region for bottom meson pair production does not exceed √ s ≈ 11.2 GeV, leaving a roughly four times smaller window in relative comparison to open charm production. Fortunately, this problem can be solved by considering higher moments, which in contrast to the charm case [4,[7][8][9] is a viable option for m b .
The essential feature of our approach (Sect. 2) is that the masses and electronic decay widths of the low-lying Υ resonances provide sufficient experimental knowledge to determine m b , as long as the 0 th moment is considered alongside the more standard positive-n moments. We may then use the limited experimental information from the continuum region that is available to test the stability of our results in Sect. 3 as a function of the moment number, and to control (in fact over-constrain) the theoretical uncertainty (see Sect. 4). We present our conclusions and a comparison with other approaches in Sect. 5.

Moment sum rules
The transverse part of the correlation functionΠ q (t) (quantities marked with a caret are defined in the MS renormalization scheme) of two heavy quark vector currents obeys the subtracted dispersion relation [10], where R q (s) = 12πImΠ q (s), and wherem q =m q (m q ) is the heavy quark mass. Taking derivatives in the limit t → 0, one obtains the moments [2,3,11], A 0 th moment [1] can also be defined, provided the limit t → ∞ and the integration over ImΠ(s) at s → ∞ is regularized by properly chosen subtractionsΠ ∞ (−t) and R ∞ q (s) [1,4]. M 0 is then obtained from the dispersion relation for the differenceΠ q (−t) −Π q (0) where the (unphysical) constantΠ q (0) is subtracted. The subtraction and the explicit sum rule for M 0 will be given below.  (7). The values quoted with an uncertainty are taken from Ref. [22].
The left-hand sides of Eqs. (2) and (3) can be calculated in perturbative QCD (pQCD) order by order in the strong couplingα s (m q ) as a function ofm q . On the right-hand side one can use the optical theorem to relate R q (s) to the cross section for heavy quark production in e + e − annihilation. It can be split into a contribution from a small number of narrow resonances below the heavy quark production threshold and a continuum contribution above, One possible method to determinem q is thus to combine data, where available, for the evaluation of the integrals on the right-hand side of Eqs. (2) and (3) with predictions from pQCD at large s where there are no data. This approach has been followed, for example, in Refs. [7,9,[12][13][14]. A certain amount of modeling is necessary since experimental information about R q (s) is restricted to relatively small energies.
Here, we will choose a different strategy. The idea is to describe the continuum region above the heavy quark production threshold on average only, not having to rely on local quark-hadron duality. We follow Refs. [1,4] and use the simple ansatz , where 3Q 2 q λ q 1 (s) is the zero-mass limit of R q (s) and s : . M is taken as the mass of the lightest pseudoscalar meson, i.e. in the case of the bottom quark, M = M B ± = 5.27934 GeV [6]. This ansatz guarantees a smooth transition between the onset of the heavy quark production threshold at 2M and pQCD at large s. Since we only need to consider moments, fine details of the ansatz are not very important. However, we will also investigate variations of our ansatz where the resonances above the threshold 4M 2 , Υ(4S), Υ(5S), and Υ(6S), are explicitly added to the expression (5).
The two unknowns, namely the heavy quark massm q (m q ), and the single free parameter in Eq. (5), λ q 3 , will be determined from Eq. (3) and one of the Eqs. (2). The other moments are then fixed and can be used to check the consistency of the approach [4]. Thus, besides the value of M , only the masses and electronic decay widths of the low-lying resonances are needed as the experimental input to extractm q (m q ). The quark mass and λ q 3 can, in principle, be determined from any combination of two moments. But only including the 0 th moment provides the leverage to sufficiently break the correlation between λ q 3 fromm q . We now give the explicit expressions needed for our numerical evaluation. From now on we particularize to the bottom quark case, in which we may neglect higher-dimensional operators in the OPE, such as from the gluon condensate 1 . Perturbative QCD predictions for the positive moments can be cast into the form, withĈ The coefficientsĈ n are known [15][16][17][18][19] up to O(α 3 s ) for n ≤ 3, and up to O(α 2 s ) for the rest [20,21]. The numerical values required for our analysis are collected in Table 1.
Since we evaluate the moments up to O(α 3 s ), we use the predictions for n > 3 provided in Ref. [22], inducing the uncertainties shown in the  Table 2: Resonance data [6] used in the analysis. The uncertainties from the resonance masses are negligible. The first three resonances are below the continuum threshold and define R res b (s), while the higher ones will be needed later when we evaluate the theoretical uncertainties.
Taylor expansion at q 2 = 0, the threshold expansion at q 2 = 4m 2 q , and the high-energy expansion at q 2 → ∞. The reconstruction is analytic and systematic, and is controlled by an error function which becomes smaller as more terms in the expansions are known. Once the correlator is reconstructed, one can calculate the moments in Eq. (2). An overview of our theory errors for the moments up to n = 7 is shown in Fig. 1. An alternative prediction of these coefficients in Ref. [23] has quoted errors that are smaller by one order of magnitude or more, leading to a total error in the extractedm b about an MeV smaller, but we believe that the conservative approach of Ref. [22] is a better reflection of the corresponding error.
We shall approximate the contributions of the narrow resonances by δ-functions, where the masses M R and electronic widths Γ e R [6] are listed in Table 2. The values of the running fine structure constant at the resonance are also given in the table 2 .
Finally, we need the regularized expression for M 0 , which requires to subtract the zero- While it is known to order O(α 4 s ), we need only the third-order expression [25], , and n b = 5 is the total number of active flavors. Using the results of Refs. [26,27], the sum rule for M 0 defined in Eq. (3) reads explicitly, whereα s =α s (m b ). The third-order coefficient A 3 is available in numerical form [4,23,28], In the last line of Eq. (10) we show the numerical values for n b = 5. The onset of the continuum is at 2M , the pseudoscalar threshold. The lower integration limit in the subtraction term involving λ b 1 (s) is, in principle, arbitrary, but is set tom 2 b in concordance with the choice to evaluateα s on the right-hand side of Eq. (10) at scalem b .
For both the theoretical predictions of the moments and the contributions from resonances and continuum to the sum rules one has to assess the uncertainties. To assign a truncation error to the pQCD prediction of the moments we follow the method proposed in Refs. [1,4] and consider the largest group theoretical factor in the next un-calculated perturbative order, . Alternatively, the dependence on the renormalization scale is often used to estimate theory errors, where, for example, in Refs. [9,13] the scale was varied between 5 and 15 GeV. Our prescription is more conservative, as has already been observed in our previous analysis [4] of the charm quark mass. In order to determine the error from the continuum contribution we proceed as follows. First, we choose a pair of moments (M 0 , M n ) from whichm b (m b ) and λ b 3 are determined. Then we input this value ofm b (m b ) into Eq. (5) and integrate with the weight corresponding to the 0 th moment as in Eq. (3), but with the energy integration range restricted to the threshold region, 2M B ≤ √ s ≤ 11.20 GeV. As this is a function of λ b 3 , we can adjust its value to coincide with the corresponding integral over the experimentally determined threshold region (see Sect. 4) yielding an experimental value, denoted λ b,exp 3 . In the final step, we use λ b,exp 3 in the n th moment sum rule to re-calculatem b (m b ), and treat the difference between these twom b (m b ) values as an additional uncertainty. It serves as a control of the error component associated with the entire methodology which we will denote by λ b . For example, neglected non-perturbative contributions to the moments such as from condensates or from residual duality violations would become visible in the comparisons of the values λ b 3 from the theoretical moments with λ b,exp

Numerical results and determination ofm b
We have analyzed the determination ofm b (m b ) from different pairs of moments and using different prescriptions to include resonances on top of the continuum. The results are shown in figures and tables in this section. We find that the largest source of uncertainty is from the continuum contribution. Indeed, the values of λ b 3 derived from the mutual consistency of the moments deviate from λ exp 3 determined from data if none of the resonances above threshold are taken into account explicitly. The lower moments are more sensitive to the continuum region, and this deviation indicates that the simple ansatz using only R cont q (s) from Eq. (5) does not capture the strong onset of the cross section for energies just above the threshold for open bottom production. As a consequence, stable results are not reached for lower moments. However, the stability improves greatly with the inclusion of the Υ(4S) and Υ(5S) states. We parametrize them as Gamma distributions, where and α and β are chosen such that the peak location M R and the second derivative coincide with those of a relativistic Breit-Wigner distribution with widthΓ R , We use the peak positions M R and total width Γ R of the resonances as given in Ref. [6] and collected above in Table 2.
In order to understand why the inclusion of resonances above threshold lead to an improved determination of the b quark mass, we provide in Fig. 2 a graphical account of the landscape of R b (s) above threshold. The upper plot shows that the continuum alone does not describe the data in the energy range of the Υ(4S), Υ(5S), and Υ(6S) resonances and the pQCD limit is reached only when √ s is above threshold by an amount of the order of the b quark mass, i.e. far above the energy range where data are available. It is therefore not a suprise that with the continuum ansatz alone one cannot obtain stable solutions from the set of sum rules. The second row of plots in Fig. 2 shows how the global description of data for R b (s) can be improved by the inclusion of Gamma distributions, Eq. (13), for the Υ(4S) and Υ(5S) resonances. If we use the total decay widths Γ R in Eq. (15) as given by the PDG [6] the local description of the data is still not good; however, the moments, i.e. integrals over R b (s) can be matched. To see this more clearly one can exploit the fact that moments do not change even if the total widths are significantly increased (which we denote byΓ R ) if one aims at a better visual representation of the local behavior of the data, as done for the right plot of the middle row of Fig. 2. Here a good description of the data on average is clearly visible. The lower row of plots in Fig. 2 shows other possible choices, namely to add only one resonance, the Υ(4S), or three resonances, Υ(4S), Υ(5S) and Υ(6S), on top of the continuum. The first (latter) choice would lead to an underestimate (overestimate) of moments in the region above threshold. As a consequence, these choices would lead to solutions for λ b 3 from the set of sum rules in disagreement with λ b,exp  We therefore determine the two free parameters, λ b 3 andm b (m b ) from pairs of sum rules using the continuum ansatz where we include the Υ(4S) and Υ(5S) as described above. The   , determined from different pairs of moments as described in the text, where the Υ(4S) and Υ(5S) resonances have been added explicitly to the ansatz in Eq (5). Below the double line: breakdown of the uncertainties inm b (m b ) followed by the total errors. The dependence onα s is shown in the next-to-last line, where the minus sign indicates thatm b decreases whenα s is increased relative to the reference valueα s (M Z ) = 0.1182. The last line contains the uncertainty induced from ∆α s (M Z ) = ±0.0016, i.e. the error obtained from the global fit to electroweak precision data [6]. See also Fig. 3 for a graphical representation of these results. results are summarized in Table 3 and Fig. 3, including the breakdown of the uncertainties from the different sources as discussed before. Results for other options, (i) where we do not include resonances on top of the continuum, (ii) where we include only the Υ(4S), (iii) or where we include additionally the Υ(6S) parametrized as a Gamma distribution as well, are presented in Fig. 4. The shift of the value for the b quark mass induced by these different options is small; for example including three resonances above threshold,m b (m b ) would be reduced by 1.3 MeV, i.e. by much less than our error estimate. The most stable result and smallest overall uncertainty is obtained with our default option in Fig. 3.
A summary of the best determination in each scenario is shown in Table 4.   Table 3), and the uncertainty induced by ∆α s (M Z ) = ±0.0016 is shown in purple.

Experimental moments
Our determination ofm b (m b ) described above does not rely on the details of experimental data for R b except resonance parameters. However, a comparison with data for R b allows us to calibrate the uncertainty of them b (m b ) determination. As described above, this is done by calculating moments from data and extracting an experimental value for λ b,exp 3 which can be compared with the value of λ b 3 obtained from the consistency relations for moments. In this section we present the details of our determination of λ b,exp 3 . We take data from the BaBar Collaboration [34]. These data cover the range of energies between √ s = 10.54 and 11.20 GeV (cf. Fig. 5). Data from the Belle Collaboration [35] will be used to obtain a cross-check, but they cover too short a range in energies to be useful for a calculation of moments for our purpose.

Data and corrections
The published experimental data for continuum heavy quark production must be corrected for vacuum polarisation and QED radiative effects before they can be used in our analysis.   Corrections due to vacuum polarisation can be taken into account by substituting the value for α em used in the experimental work by the running fine structure constant, α em ( √ s). Since the variation of α em ( √ s) in the considered energy range is very small, we take it to be constant and use (α em (0)/α em (M R )) 2 = 0.93, (see Table 2). This factor should be multiplied with the measured R b ratio.
BaBar experimental data are available for energies above the open bottom threshold. In this energy range, the radiative tails from the Υ(1S), Υ(2S), Υ(3S) resonances contribute. The required corrections are provided by BaBar in supplementary material to Ref. [34] and are easily subtracted from the data.
To remove initial-state radiative (ISR) effects from the continuum data after subtracting radiative tails from the resonances, we use the prescription following Refs. [29,30] (see also Refs. [9,14]). The measured R ratio,R, is given by a convolution, of the true R ratio with the radiator function G(z, s) describing QED corrections. G(z, s) is taken from Ref. [14] and includes next-to-next-to-leading order contributions. The lower integration limit of the integral in Eq. (17) should start at the onset of the continuum region, which we fix at z 0 = s 0 /s with s 0 = (10.54 GeV) 2 . The true R ratio must be determined by inverting (i.e. unfolding) Eq. (17). This can be done iteratively imposing the boundary condition R(s 0 ) = 0. This condition is automatically satisfied by the BaBar data after subtraction of the radiative tails. The BaBar data corrected for vacuum polarization, radiative tails and ISR is shown in Fig. 6 (red points). We also show the uncorrected data (blue points), which are the same as shown in Fig. 5. BaBar data contain an outlier at √ s = 10.86 GeV (not shown in Fig. 6). At this energy, there are two different experimental measurements, separated by only ∆ √ s = 0.0005 GeV which disagree among themselves. Instead of removing this point, as has been suggested in Ref. [9], we take the average of the two points and ascribe, as an error, the difference of the two measured R values. We have checked that either option, removing the outlier or averaging with the close-by point, translates into a tiny difference for the experimental moments.

Numerical results for moments
Experimental moments are calculated as numerical integrals over the ISR corrected R values, using the trapezoidal rule. We collect our results in Table 5    is obtained by direct integration over corrected data. The first error is due to the uncorrelated statistical and systematic uncertainties of the data, while the second is the correlated one. The third and fourth columns show the moments calculated from our ansatz for R b (s) with λ b,exp 3 = 0.82(20) (column 3) and λ b 3 = 1.53 (column 4) as input. In both cases,m b = 4.1802 GeV was used. The last column collects the experimental moments when BaBar data is used without any kind of correction or subtraction.
following the prescription given by the BaBar Collaboration [34]. For comparison, we show in Table 5 also the moments calculated from R b (s), but using the value λ b 3 = 1.53. This value was obtained in our preferred scenario where the Υ(4S) and Υ(5S) resonances are included on top of the continuum and using the pair of moments (M 0 , M 7 ). In the last column of Table 5 we also show moments calculated from uncorrected data. One can see that ISR corrections are indeed very small and do not introduce an additional source of uncertainties.
In Table 6 we compare our determination of moments with those from Ref. [9] and Ref. [14]. To do so, we have to adjust the energy range correspondingly. For both references the lower limit of the energy range was chosen at √ s = 10.62 GeV. The upper integration limit was √ s = 11.20 GeV in Ref. [9] and √ s = 11.24 GeV in Ref. [14]. We also follow Refs. [9,14] and subtract the Υ(4S) resonance, which is parameterized by a Breit-Wigner distribution, as well as its radiative tail. Above √ s = 11.20 GeV, we use our ansatz to extrapolate up to 11.24 GeV. As can be seen from Table 6, we find good agreement with both references.
We have used the BaBar data since it covers an energy range large enough to extract a reliable description of the continuum region. The Belle Collaboration [35] also provides a measurement of R b (s), but only the narrow energy range between √ s = 10.620 and 11.047 GeV is covered with the first three experimental points quite disconnected from the fine-scan around the Υ(5S) and Υ(6S) resonances, i.e. 10.754 GeV ≤ √ s ≤ 11.047 GeV. If we use Belle data we find that this short energy range contributes 0.198 (7) to the 0 th experimental moment, to be compared with 0.172(5) from BaBar data for the same energy region. These results are compatible at the 3σ level only. Such a difference could by attributed to the different treatment of QED radiative effects of the narrow resonances in the case of Belle data. For our calculation of the 0 th moment we have used Belle data corrected for vacuum polarisation effects, but without subtracting radiative tails. The Belle Collaboration does not provide the corresponding information. If we had used the radiative tail provided by BaBar, we would find 0.167(7) for the 0 th moment. This would bring the values of the 0 th moment calculated from Belle or from BaBar data in very good agreement.
In our previous analysis of data for charm quark production [4] we found that using data in an energy range between 3.7 and 5 GeV, extended by using pQCD above, can lead to a consistent picture and a reliable determination of the charm quark mass. A correspondingly large energy window for the bottom quark would cover energies up to 15    one unit of the heavy quark mass above threshold. Unfortunately, data are available only for √ s ≤ 11.2 GeV. The treatment of the energy range 11.2 GeV ≤ √ s 15 GeV requires special care and may lead to additional uncertainties. In Ref. [14] it was argued that the gap above √ s = 11.2 GeV should be described by pQCD. However, this introduces a discontinuity with the experimental data. In Ref. [9] a smooth polynomial fit was used instead. We opt for using our own ansatz , which approaches pQCD only for s → ∞. In Table 7 we compare two possible options: calculating moments from data in the window 2M B ≤ √ s ≤ 11.2 GeV, combined either with pQCD or our ansatz for R b (s), both up to 15 GeV. If one uses data and pQCD above √ s = 11.2 GeV, i.e. a prescription with a discontinuity, one obtains moments which are larger than for the case where we use our ansatz with a smooth √ s-dependence. Correspondingly, this will result in smaller values for the bottom mass, as found in Ref. [14].
Differences in the b quark mass determination between Refs. [9,14] and our work can thus be traced to a different prescription for including contributions to the moments from an energy range where no data are available. One could argue that this should be considered as an additional systematic error for the b quark mass. Experimental data covering the energy range between 11 and 15 GeV are definitely needed to ultimately solve this issue. Until such data will become available, we believe that a description of the unknown part of R b (s) with a smooth function is preferable over one with a discontinuity.  Figure 7 where we group the results in two sets and in chronological order within each set. The first set is based on phenomenological approaches, extractingm b (m b ) by comparing theory predictions with data. This includes other results based on relativistic sum rules, Chetyrkin (2009) [14], Bodenstein (2012) [36], and Dehnadi (2015) [9], shown as red diamonds. The methodology of these publications is closest to our own approach. Our higher value form b (m b ) can be traced to the treatment of the intermediate energy behavior where our method approaches the perturbative regime of QCD at higher energies, as discussed in detail above. We also display results based on non-relativstic sum rules (orange squares), Laschka (2011) [37], Penin (2014) [38], Beneke (2015) [39], Kiyo (2016) [40], and Peset (2018) [41], as well as on other sum rule methods (green stars), Lucha (2013) [42] and Narison (2020) [43]. The bottom quark mass (brown stars) was also determined as a byproduct in a global fit to inclusive semileptonic B-meson decays to obtain the CKM matrix element V cb , Alberti (2015) [44], and from an analysis of deep inelastic scattering data at HERA compared with perturbative QCD calculations, Abramowicz (2018) [45].
The results in the lower part of Fig. 7 (blue points) are lattice QCD calculations. They are based on an improved non-relativistic QCD action, Lee (2014) [46], on Heavy Quark Effective Theory non-perturbatively matched to QCD, Bernardoni (2014) [47], on using time-moments of the vector current-current correlator, Colquhoun (2015) [48], as well as the MILC highly improved staggered quark ensembles with four flavors of dynamical quarks, Bazavov (2018) [49].
We also show the average of the 2021 FLAG Review [50] for N f = 2 + 1 + 1.