1 Introduction

As a fundamental parameter of the standard model (SM) of particle physics that enters the calculation of many physical observables, it is important to determine the charm quark mass as accurately as possible and by using several different methods. Moreover, it is one of the central goals of future high-energy lepton colliders to measure the charm quark Yukawa coupling to sub-percent precision in order to test the mass-coupling relation, as predicted in the minimal SM with a single Higgs doublet. Thus, the precision in the charm quark mass should at least match and ideally surpass that in the charm Yukawa coupling.

Most recent determinations [1] of the \(\overline{\mathrm{MS}}\) charm quark mass, \(\hat{m}_c(\mu = \hat{m}_c)\), including all of the most precise ones, rely on either lattice simulations or variants of QCD sum rules [2,3,4]. Both approaches are based on rigorous field theoretical principles and can be systematically improved. Nevertheless, the resulting uncertainties are dominated by theory errors which are notoriously difficult to estimate and often subject of vigorous debate.

The goal of this article is then to revisit the relativistic sum rule method with emphasis on the evaluation of the uncertainty and specifically to show that the overall error may also be constrained within our approach. To this end we will adopt a strategy first considered in Ref. [5] where the only exploited experimental information are the masses and electronic decay widths of the narrow resonances in the sub-continuum charm region, \(J/\Psi \) and \(\Psi (2S)\). Consistency between different QCD sum rules will be seen to suffice to constrain the continuum of charm pair production with good precision. For this procedure to work it is crucial to include alongside the first and higher moment sum rules also the zeroth moment, as the latter exhibits enhanced sensitivity to the continuum. This is in contrast to other analyses which do not use the zeroth moment. Comparison with existing data on the R-ratio for hadronic relative to leptonic final states in \(e^+ e^-\) annihilation will then serve as a control, providing an independent error estimate which we interpret as the error on the method and add it as an additional error contribution. In this way, we can show that the overall precision in \(\hat{m}_c\) from relativistic sum rules is at the sub-percent level. Compared with the previous work [5], we push the calculation one order higher and use more recent and more precise experimental data including updated values for some input parameters, as the strong coupling \(\alpha _s\) and the gluon condensate.

Specifically, there are two sources of uncertainty which are particularly difficult to quantify, and therefore occasionally not even included in the error budget. One are condensate contributions within the operator product expansion (OPE) beyond the leading gluon condensate term. The other issue is quark–hadron duality: Formally, the sum rule is based on the optical theorem, where the real part of the correlator is evaluated theoretically with quark and gluon degrees of freedom, while the imaginary part involves as dual degrees of freedom the hadrons observed in experiments. While such a global change of variables is highly non-trivial, it is innocuous as long as one does not mix these descriptions on either side of the sum rule. As far as the evaluation of the imaginary part is concerned, one is forced, however, to switch at some specific point in the squared energy s from experimental data to perturbative QCD since in practice data are necessarily constrained to a finite region while the upper integration limit is unbounded. In this way one has to rely on local quark–hadron duality which—at least as a matter of principle—is unjustified. Still, as long as s is large enough, one introduces little additional uncertainty. Now, the largest value of \(\sqrt{s}\) for which data of the total hadronic cross section in \(e^+ e^-\) annihilation are available (see the data points in Fig. 5 in Sect. 3) is 5 GeV, and beyond \(\sqrt{s} = 4.6\) GeV data points are scarce and have very large errors. Around this energy there is still considerable fluctuations in the measured cross section, shedding some doubt on the applicability of local quark–hadron duality even in practice. One of the features of our work is that it merely relies on quark–hadron duality in a finite region, namely between the \(\psi (2S)\) resonance and the continuum threshold. While this is still not rigorously justified, it should largely mitigate the aforementioned problem. Furthermore, by using continuum data only as a calibration of the uncertainty, we control the error associated with the necessary deviation from strict global quark–hadron duality. Our approach for the uncertainty calibration of the charm mass is conceptually new since it allows us to estimate the effect from correlated errors. This should lead to a more reliable uncertainty estimate than is possible in other approaches.

The paper is organized as follows. In the next section we review the formalism following Ref. [5] and describe how self-consistency between moments of the current correlator allows us to constrain the continuum region and determine a precise value of the charm quark mass. This section will already contain our main result, \(\hat{m}_c(\hat{m}_c) = 1272 \pm 8\) MeV for \(\hat{\alpha }_s(M_Z) = 0.1182\). In Sect. 3 we add a detailed discussion of the influence of the continuum region using experimental data and confirm the validity of the result of Sect. 2. Section 4 offers details of a more general fit procedure where parameter uncertainties can be taken into account in a more systematic way, and we also include a comparison with earlier results. We summarize in Sect. 6.

2 Formalism and charm-mass determination

We consider the transverse part of the correlator \(\hat{\Pi }_q(t)\) of two heavy-quark vector currents where the caret indicates \(\overline{\mathrm{MS}}\) subtraction. \(\hat{\Pi }_q(t)\) can be calculated in perturbative QCD (pQCD) order by order and obeys the subtracted dispersion relation [6]

$$\begin{aligned} 12 \pi ^2 \frac{\hat{\Pi }_q (0) - \hat{\Pi }_q (-t)}{t} = \int _{4 \hat{m}_q^2}^\infty \frac{\mathrm{d} s}{s} \frac{R_q(s)}{s + t}, \end{aligned}$$
(1)

where \(R_q(s) = 12 \pi \text{ Im } \hat{\Pi }_q(s)\), and \(\hat{m}_q = \hat{m}_q (\hat{m}_q)\) is the mass of the heavy quark. In the limit \(t\rightarrow 0\), Eq. (1) coincides with the first moment, \(\mathcal{M}_1\), of \(\hat{\Pi }_q(t)\). In general, there is a sum rule for each higher moment as well [2,3,4],

$$\begin{aligned} \mathcal{M}_n := \left. \frac{12\pi ^2}{n !} \frac{\mathrm{d}^n}{\mathrm{d} t^n} \hat{\Pi }_q(t) \right| _{t=0} = \int _{4 \hat{m}_q^2}^\infty \frac{\mathrm{d} s}{s^{n+1}} R_q(s). \end{aligned}$$
(2)

Taking the opposite limit of Eq. (1), i.e. \(t \rightarrow \infty \), multiplying with t and properly regularizing, one can also define a sum rule for the \(0^{th}\) moment [5]. \(\mathcal{M}_0\) is then obtained from the dispersion relation for the difference \(\hat{\Pi }_q (0) - \hat{\Pi }_q (-t)\) where the (unphysical) constant \(\hat{\Pi }_q (0)\) is subtracted, tantamount to the definition of higher moments as derivatives of \(\hat{\Pi }_q (-t)\). At a given order in pQCD, the required regularization can be obtained by subtracting the zero-mass limit of \(R_{q}(s)\), which we write as \(3 Q_q^2 \lambda _1^q(s)\) with \(Q_q\) the quark charge. \(\lambda _1^q(s)\) is known up to \(\mathcal{O}(\hat{\alpha }_s^4)\), but we will only need the third-order expression [7],

$$\begin{aligned} \lambda ^q_1(s)= & {} 1 + \frac{\hat{\alpha }_s}{\pi } + \frac{\hat{\alpha }_s^2}{\pi ^2} \left[ \frac{365}{24} - 11 \zeta (3) + n_q \left( \frac{2}{3} \zeta (3) - \frac{11}{12} \right) \right] \nonumber \\&+\; \frac{\hat{\alpha }_s^3}{\pi ^3} \left[ \frac{87029}{288} - \frac{121}{8} \zeta (2) - \frac{1103}{4} \zeta (3) + \frac{275}{6}\zeta (5) \right. \nonumber \\&+\; \left. n_q \left( - \frac{7847}{216} + \frac{11}{6} \zeta (2) + \frac{262}{9} \zeta (3) - \frac{25}{9}\zeta (5) \right) \right. \nonumber \\&+\;\left. n_q^2 \left( \frac{151}{162} - \frac{\zeta (2)}{18} - \frac{19}{27} \zeta (3) \right) \right] , \end{aligned}$$
(3)

where \(\hat{\alpha }_s = \hat{\alpha }_s (\sqrt{s})\), \(n_q = n_l + 1\) and \(n_l\) is the number of light flavors (taken as massless), i.e., quarks with masses below the heavy quark under consideration.

By the optical theorem, \(R_{q}(s)\) can be related to the measurable cross section for heavy-quark production in \(e^+e^-\) annihilation. Below the threshold for continuum heavy-flavor production, the cross section is determined by a small number of narrow Breit–Wigner resonances [2], approximated by \(\delta \)-functions,

$$\begin{aligned} R_q^\mathrm{Res}(s) = \frac{9\pi }{\alpha _\mathrm{em}^2(M_R)} M_R \Gamma _R^e\delta (s-M_R^2). \end{aligned}$$
(4)

The masses \(M_R\) and electronic widths \(\Gamma ^e_R\) of the resonances [8] needed for the determination of the charm mass are listed in Table 1 and \(\alpha _\mathrm{em}(M_R)\) is the running fine structure constant at the resonance.Footnote 1

Table 1 Resonance data [8] used in the analysis. The uncertainties from the resonance masses are negligible for our purpose

We invoke global quark–hadron duality, and we assume that continuum production can be described on average by the simple ansatz [5],

$$\begin{aligned} R_q^\mathrm{cont}(s) = 3 Q^2_q \lambda ^q_1 (s) \sqrt{1 - \frac{4\, \hat{m}_q^2 (2 M)}{s^\prime }} \left[ 1 + \lambda ^q_3 \frac{2\, \hat{m}_q^2(2 M)}{s^\prime } \right] ,\nonumber \\ \end{aligned}$$
(5)

where \(s' := s + 4(\hat{m}_q^2(2M) - M^2)\), and M is taken as the mass of the lightest pseudoscalar charmed meson, i.e., \(M = M_{D^0} = 1864.84\) MeV [8]. \(\lambda ^q_3\) is a free parameter to be determined. \(R_q^\mathrm{cont}(s)\) interpolates smoothly between the threshold and the onset of open heavy-quark pair production and coincides asymptotically with the prediction of pQCD for massless quarks.

Using the results of Refs. [10, 11] and \(R_q(s) = \sum \nolimits _\mathrm{resonances} R_q^\mathrm{Res}(s) + R_q^\mathrm{cont}(s)\), the sum rule for \(\mathcal{M}_0\) obtained from Eq. (1) in the limit \(t \rightarrow \infty \) and dividing by \(3Q_q^2\) reads explicitly

$$\begin{aligned}&\sum \limits _\mathrm{resonances} \frac{9\pi \Gamma ^e_R}{3 Q_q^2 M_R \alpha _{em}^2 (M_R)} + \int \limits _{4 M^2}^\infty \frac{\mathrm{d} s}{s} \left( \frac{R_q^\mathrm{cont}(s)}{3 Q_q^2} -\lambda _1^q(s)\right) \nonumber \\&\qquad -\;\int \limits _{\hat{m}_q^2}^{4M^2} \frac{\mathrm{d} s}{s} \lambda ^q_1 (s)= - \frac{5}{3} + \frac{\hat{\alpha }_s}{\pi }\left[ 4 \zeta (3) - \frac{7}{2} \right] \nonumber \\&\qquad +\;\left( \frac{\hat{\alpha }_s}{\pi }\right) ^2 \left[ \frac{2429}{48} \zeta (3) - \frac{25}{3} \zeta (5) - \frac{2543}{48}\right. \nonumber \\&\qquad \left. +\;n_q \left( \frac{677}{216} - \frac{19}{9} \zeta (3) \right) \right] + \left( \frac{\hat{\alpha }_s}{\pi }\right) ^3 A_3, \nonumber \\&\quad = -1.667 + 1.308\, \frac{\hat{\alpha }_s}{\pi }+ 1.595 \left( \frac{\hat{\alpha }_s}{\pi }\right) ^2 - 8.427 \left( \frac{\hat{\alpha }_s}{\pi }\right) ^3,\nonumber \\ \end{aligned}$$
(6)

where \(\hat{\alpha }_s = \hat{\alpha }_s(\hat{m}_q)\) and the third-order coefficient \(A_3\) is available in numerical form [12, 13],

$$\begin{aligned} A_3 = -9.863 + 0.399 \, n_q - 0.010 \, n_q^2\ . \end{aligned}$$
(7)

In the last line of Eq. (6) we show numerical values for the case of charm quarks, i.e. \(n_q=4\), exhibiting a marked breakdown of convergence of the perturbative expansion. Notice that the continuum \(R_q^\mathrm{cont}(s)\) contributes with the lower integration limit \(4M^2\), while the subtraction term \(\lambda _1^q(s)\) is integrated starting from \(\hat{m}_q^2\). We reiterate that this definition of \(\mathcal{M}_0\) does not involve the unphysical constant \(\hat{\Pi }_q (0)\) separately, but only the difference \(\hat{\Pi }_q (0) - \hat{\Pi }_q (-t)\).

Table 2 The coefficients \(C_n^{(i)}\) for the perturbative expansion of the QCD moments, Eqs. (8, 9), as used in this work. The first line for \(n=\)1 is taken from Ref. [14], the coefficients for \(n=2\), 3 from [18] and the last ones for \(n=4\), 5 from Ref. [21]

Equation (6) contains two unknowns, the quark mass \(\hat{m}_q(\hat{m}_q)\) and the parameter \(\lambda _3^q\) entering in our prescription for \(R_q^\mathrm{cont}(s)\). The zeroth sum rule is the most sensitive to the continuum region and therefore to \(\lambda _3^q\), while the other moments mostly determine the quark mass.

Theory predictions for the higher moments in perturbative QCD can be cast into the form

$$\begin{aligned} \mathcal{M}_n^\mathrm{pQCD} = \frac{9}{4} Q_q^2 \left( \frac{1}{2\hat{m}_q(\hat{m}_q)} \right) ^{2n} \hat{C}_n \end{aligned}$$
(8)

with

$$\begin{aligned} \hat{C}_n= & {} C_n^{(0)} + \left( \frac{\hat{\alpha }_s}{\pi }+ \frac{3 Q_q^2 \alpha _\mathrm{em}}{4 \pi }\right) C_n^{(1)} + \left( \frac{\hat{\alpha }_s}{\pi }\right) ^2C_n^{(2)} \nonumber \\&+\;\left( \frac{\hat{\alpha }_s}{\pi }\right) ^3C_n^{(3)} + \mathcal{O}(\hat{\alpha }_s^4). \end{aligned}$$
(9)

The \(\hat{C}_n\) are known up to \(\mathcal{O}(\hat{\alpha }_s^3)\) for \(n \le 3\) [14,15,16,17,18], and up to \(\mathcal{O}(\hat{\alpha }_s^2)\) for the rest [19, 20]. For the convenience of the reader we collect the numerical values of the coefficients \(C_n^{(i)}\) required for our analysis in Table 2. Since we evaluate the moments up to \(\mathcal{O} (\hat{\alpha }_s^3)\) we use the predictions for \(n>3\) provided in Ref. [21]. The QED contribution to the QCD moments proportional to \(\alpha _\mathrm{em} / \pi \) is very smallFootnote 2 and it is sufficient to include the first-order term.

In general, vacuum expectation values of higher-dimensional operators in the OPE contribute to the moments of the current correlator as well. These condensates may be important for a high-precision determination of heavy-quark masses, in particular in the case of the charm quark. The leading term involves the dimension-4 gluon condensate [2],

$$\begin{aligned} \mathcal{M}^\mathrm{cond}_n = \frac{12 \pi ^2 Q_q^2}{(4 \hat{m}_q^2)^{n+2}} \left\langle \frac{\hat{\alpha }_s}{\pi } G^2 \right\rangle \, a_n \left( 1 + \frac{\hat{\alpha }_s(\hat{m}_q)}{\pi } b_n\right) . \end{aligned}$$
(10)

The coefficients \(a_n\) and \(b_n\) can be found in Refs. [22,23,24] and are collected in Table 3. We will use \(C_G^\mathrm{exp} := \langle \frac{\hat{\alpha }_s}{\pi } G^2 \rangle ^\mathrm{exp} = 0.005\) GeV\(^4\) as our central value and assume a 100% uncertainty, \(\Delta _G = 0.005\) GeV\(^4\), taken from the recent analysis [25].

Table 3 The coefficients \(a_n\) and \(b_n\) needed to calculate the contribution from the dimension-4 gluon condensate to the QCD moments, Eq. (10). The numbers are taken from Refs. [22,23,24]
Table 4 Results for the lowest moments, \(\mathcal{M}_n\), defined in Eq. (6) for \(n = 0\) (multiplied by \(3Q_q^2\)) and Eq. (2) for \(n \ge 1\) and including the contribution from the gluon condensate, Eq. (10) (multiplied by \(10^n\, \text{ GeV }^{2n}\)). The first error for the continuum part is obtained from shifting the central value \(\lambda _3^c = 1.23\) (determined from the zeroth and the second moment) to \(\lambda _3^{c,\mathrm{exp}} = 1.34(17)\) (determined from fixing the zeroth moment by experimental data; see the next section). The second error is propagated from the experimental uncertainties of data in the threshold region, \(\Delta \lambda _3^{c,\mathrm{exp}} = \pm 0.17\). The third error is due to the uncertainty of the gluon condensate. The errors for the sum of the resonance and continuum parts (denoted ‘Total’ in column 4) combine all errors in quadrature. The last column shows the theoretical predictions using Eq. (8) for \(\hat{m}_c (\hat{m}_c) = 1.272\) GeV and \(\hat{\alpha }_s(M_Z) = 0.1182\) [26] with the truncation error determined from Eq. (12)

In the following we will use this formalism for a determination of the charm quark mass and specify the notation accordingly, using the index c instead of q where appropriate. One can use

$$\begin{aligned} \mathcal{M}_n^\mathrm{pQCD} + \mathcal{M}^\mathrm{cond}_n= & {} \int _{4 \hat{m}_c^2}^\infty \frac{\mathrm{d} s}{s^{n+1}} R_c(s) \nonumber \\= & {} \int _{4 \hat{m}_c^2}^\infty \frac{\mathrm{d} s}{s^{n+1}} [R_c^\mathrm{Res}(s) + R_c^\mathrm{cont}(s)] \end{aligned}$$
(11)

for two different n to determine values for the heavy-quark mass \(\hat{m}_c(\hat{m}_c)\) and the constant \(\lambda ^c_3\). The other moments are then fixed and can be used to check the consistency of our approach. No experimental data other than the resonance parameters in Table 1 are necessary. Numerical results choosing \(n=0\) and \(n=2\) in Eq. (11) are shown in Table 4. We show separately the contributions from the narrow resonances (column 2) and from the continuum part (column 3), evaluated using \(R_c^\mathrm{cont}(s)\) in Eq. (5), including the condensate contribution from Eq. (10). The sum of these two contributions (column 4) can be compared with the theory prediction (column 5) from Eq. (8). The second column in Table 4 accounts for the narrow vector resonances below the heavy-quark threshold, in the charm sector the first two charmonium resonances, \(J/\Psi \) and \(\Psi (2S)\). The higher resonances are included in the continuum. The errors for the resonance contributions shown in Table 4 are exclusively determined by the uncertainty of the electronic widths from Table 1, taken to be uncorrelated.Footnote 3

In order to determine an error for the continuum contributions (column 3) we proceed in the following way: we compare experimental data with the zeroth moment in the restricted energy range of the threshold region, \(2 M_{D^0} \le \sqrt{s} \le 4.8\) GeV, to obtain an experimental value for \(\lambda _3^c\), denoted \(\lambda _3^{c,\mathrm{exp}}\) (details will be described below in Sect. 3). Here we fix \(\hat{m}_c(\hat{m}_c) = 1.272\) GeV. Then we can also determine an error, \(\Delta \lambda _3^{c,\mathrm{exp}}\), from the experimental uncertainty of the data in this threshold region. The shift in the moments resulting from the different values for \(\lambda _3^{c}\) (either from two moments combined with resonance data only, or from the comparison of the zeroth moment with continuum data in the threshold region) turns out to be small. Strictly speaking this shift is a one-sided error, but to be conservative we include it as an additional error in the results of Table 4. Note that non-perturbative effects not included in our expressions for the moments, as for example omitted higher-dimensional condensates, would become visible in the comparison of the values of \(\lambda _3^{c,\mathrm{exp}}\) determined from data with \(\lambda _3^{c}\) obtained from moments. Since we include this difference as a contribution to the uncertainty budget of the charm mass, we do not include an additional uncertainty due to the truncation of the OPE. More details of our determination of \(\lambda _3^{c,\mathrm{exp}}\) are given in the next section.

Finally, we assign a truncation error to the theory prediction of the moments. Following the method proposed in Ref. [5] we consider the largest group theoretical factor in the next uncalculated perturbative order as a way to estimate these errors,

$$\begin{aligned} \Delta \mathcal{M}_n^{(i)} = \pm Q_q^2 N_C C_F C_A^{i-1} \left( \frac{\hat{\alpha }_s (\hat{m}_q)}{\pi }\right) ^i \left( \frac{1}{2 \hat{m}_q(\hat{m}_q)} \right) ^{2n}, \nonumber \\ \end{aligned}$$
(12)

(\(N_C = C_A = 3\), \(C_F = 4/3\)). At order \(\mathcal{O}(\hat{\alpha }_s^4)\), this corresponds to an uncertainty of \(\pm 48 (\hat{\alpha }_s/\pi )^4\) for \(\hat{C}_n^{(4)}\) in Eq. (8). Applying this prescription to the known \(\mathcal{O}(\hat{\alpha }_s^3)\) terms, we observe that the resulting errors are quite conservative: the error estimate from this approach for \(\hat{m}_c(\hat{m}_c)\) is more conservative by a factor given in the last column of Table 5 compared to the alternative to base the error estimate on the last known term in the perturbative expansion. We have convinced ourselves that this carries over to the errors for the charm mass: the above prescription is more conservative than the often used alternative approach to vary the renormalization scale within a conventional factor of four.

For the moments with \(n>3\) (which, however, will not enter into our final result) taken from Ref. [21] we have to include additional uncertainties specific to the method used to obtain predictions for \(\mathcal{M}_n\). These errors are small, but they are included for completeness. An overview of our theory errors for the moments up to \(n=5\) is shown in Fig. 1.

Table 5 Ratios of the truncation errors \(\Delta \mathcal{M}_n^{(k)}\) from Eq. (12) and the known contribution to the moments \(\mathcal{M}_n^{(k)}\) at order \(O(\hat{\alpha }_s^k)\) (see Eq. (8)). In the last column, \(\hat{m}_c (\hat{m}_c) = 1.272\) GeV has been used. The ratios are often much greater than unity, showing that our estimated truncation errors are conservative
Fig. 1
figure 1

Theoretical moments, Eq. (8), multiplied by \(10^n\, \text{ GeV }^{2n}\), for the charm quark using \(\hat{m}_c(\hat{m}_c) = 1.272\) GeV at different orders of \(\hat{\alpha }_s\). Blue squares show results of Eq. (8) up to \(\mathcal{O}(\hat{\alpha }_s)\), red circles up to \(\mathcal{O}(\hat{\alpha }_s^2)\), and green triangles up to \(\mathcal{O}(\hat{\alpha }_s^3)\). The error bars show the truncation errors given in Eq. (12) at the given order

Table 6 Values and breakdown of the uncertainties of \(\hat{m}_c(\hat{m}_c)\) (in MeV) and \(\lambda _3^c\) determined from different pairs of the moments. The line denoted ‘Total’ gives the quadratic sum of the errors from \(\lambda _3^c\), the resonances, the gluon condensate and the pQCD truncation. Numerical values for the uncertainties from \(\hat{\alpha }_s\) and \(C_G = \langle \frac{\hat{\alpha }_s}{\pi } G^2 \rangle \) (in units of GeV\(^4\)) are shown in separate lines. In the line labeled ‘Electroweak fit’ we use \(\Delta \hat{\alpha }_s(M_z) = 0.0016\) [26], in the last line denoted ‘Lattice’, we use \(\Delta \hat{\alpha }_s(M_z) = 0.0012\) [8]. See Fig. 2 for a graphical representation

The charm mass and the continuum parameter \(\lambda _3^c\) can, in principle, be determined from any combination of two moments. The zeroth moment, however, provides the highest sensitivity. The results for combinations of the zeroth with one higher moment are summarized in Table 6, including the breakdown of the uncertainties due to their different sources. We include the difference between the two possibilities to determine \(\lambda _3^c\) as described above as an estimate of the error due to the method, in particular the truncation of the OPE. In Fig. 3 we also show the error budget at different orders of \(\hat{\alpha }_s\) as obtained from the combination of the moments \(\mathcal{M}_0\) and \(\mathcal{M}_2\). It is remarkable that the parametric uncertainty due to the input value of \(\hat{\alpha }_s(M_Z)\) increases at the last step from \(O(\hat{\alpha }_s^2)\) to \(O(\hat{\alpha }_s^3)\) (see the last, purple error bars), while the central value for \(\hat{m}_c (\hat{m}_c)\) remains stable. A closer look at the underlying formulas reveals that this is due to the bad convergence of the right-hand side of Eq. (6). It is therefore unlikely that knowledge of the next terms in the perturbative series will lead to a reduction of the error on the charm mass. Our final result for the charm quark mass,

$$\begin{aligned} \hat{m}_c(\hat{m}_c) = (1272 + 2616 \Delta \hat{\alpha }_s \pm 8)~\mathrm{MeV}, \end{aligned}$$
(13)

is based on \((\mathcal{M}_0, \mathcal{M}_2)\) where we explicitly exhibit its dependence on \(\hat{\alpha }_s\) relative to our adopted central value,Footnote 4 i.e., \(\Delta \hat{\alpha }_s = \hat{\alpha }_s(M_z) - 0.1182\).

3 Continuum data

Our final result for the charm quark mass given above in Eq. (13) is obtained from an analysis of sum rule moments without using experimental data from the continuum region. Instead, the continuum is described by the parameter \(\lambda _3^c\), and determined simultaneously with \(\hat{m}_c\) within our formalism. In this section we show that Eq. (5) indeed reproduces the experimental moments well (even though our ansatz must obviously fail to describe the underlying cross-section data locally).

Fig. 2
figure 2

\(\hat{m}_c(\hat{m}_c)\) using different combinations of the moments and the corresponding uncertainties. Blue represents the full error, red is the one from the resonance region, green from the truncation errors of the theoretical moments, cyan from the error in \(\lambda _3^c\) (which is the symmetrized combination of the shift and the experimental error on \(\lambda _3^c\)), orange from the gluon condensate and purple from the uncertainty induced by \(\Delta \hat{\alpha }_s(M_Z) = 0.0016\). Notice that all determinations are mutually consistent within our error estimates

Fig. 3
figure 3

Error budget for \(\hat{m}_c(\hat{m}_c)\) determined from the combination of \(\mathcal{M}_0\) and \(\mathcal{M}_2\) at different orders in \(\hat{\alpha }_s\). The color coding is as in Fig. 2

In order to do so, we have to add the background from light-quark contributions [6, 23, 34, 35], \(R_\mathrm{background}(s) = R_\mathrm{uds}(s) + R_\mathrm{uds(cb)}(s) + R_\mathrm{singlet}(s) + R_\mathrm{QED}(s)\) to our model for the charm continuum. The first part, \(R_\mathrm{uds}(s)\), is given by Eq. (3) with \(n_q=3\). Contributions from virtual or secondary heavy quarks, \(R_\mathrm{uds(cb)}(s)\), appear first at \(\mathcal{O}(\hat{\alpha }_s^2)\) and are known up to \(\mathcal{O}(\hat{\alpha }_s^3)\) in approximate form [34]. In addition, at order \(\mathcal{O}(\hat{\alpha }_s^3)\) one has to take into account the singlet contribution \(R_\mathrm{singlet}(s)\) [34] which is found to be numerically small, \(-0.55 (\hat{\alpha }_s/\pi )^3\) [23], but included for completeness. Finally, also corrections due to QED have been taken into account through \(R_\mathrm{QED}(s)\) [6].

In Fig. 4 we show the resulting R ratio, \(R(s) = R_c^\mathrm{cont}(s) + R_\mathrm{background}(s)\), compared with experimental \(e^+ e^-\) annihilation data for final state hadrons. The data are compiled from Crystal Ball [28], BES [29,30,31,32], and CLEO [33]. The inset zooms into the region above threshold in the energy range of the \(\Psi (3770)\) resonance, \(2 M_{D^0} \le \sqrt{s} \le 3.83\) GeV, with data points from BES [30, 31]. The full red curve in this figure is obtained from our parametrization of \(R_c^\mathrm{cont}(s)\) in Eq. (5), with \(\lambda _3^c = 1.23\) and \(\hat{m}_c(\hat{m}_c)\) in Eq. (13).

Fig. 4
figure 4

Data for the ratio R(s) for \(e^+e^- \rightarrow \) hadrons in the charm threshold region: Crystal Ball CB86 (green) [28]; BES00, 02, 06, 09 (black, blue, cyan, and red) [29,30,31,32], and CLEO09 (orange) [33]. The full (red) curve is our parametrization of \(R_c^\mathrm{cont}(s)\) with \(\lambda _3^c = 1.23\) and \(\hat{m}_c(\hat{m}_c) = 1.272\) GeV. The inner plot is a zoom into the energy range \(2 M_{D^0} \le \sqrt{s} \le 3.83\) GeV

Fig. 5
figure 5

Same as Fig. 4 at different orders of \(\hat{\alpha }_s\). For the full (red) curve we have used \(\lambda _3^c\) determined from the pair of the moments \(\mathcal{M}_0\) and \(\mathcal{M}_2\) at the respective order indicated in the plots with errors as described in the text, while for the blue band we have used \(\lambda _3^{c,\mathrm{exp}}\) determined from data and their uncertainties at each order in \(\hat{\alpha }_s\)

As described in the previous section, we can use data directly to obtain an experimental value \(\lambda _3^{c, \mathrm{exp}}\). A comparison of our ansatz Eq. (5) using either this experimental parameter or the one determined via pairs of the moments as explained above is shown in Fig. 5. Here we compare the resulting parametrizations of the charm continuum with data based on a determination of either \(\lambda _3^c\) or \(\lambda _3^{c, \mathrm{exp}}\) performed at different orders of pQCD. The red bands in Fig. 5 use \(\lambda _3^c\) and \(\hat{m}_c(\hat{m}_c)\) obtained at each order in \(\hat{\alpha }_s\) while the blue bands use the same value for \(\hat{m}_c(\hat{m}_c)\), but \(\lambda _3^{c,\mathrm{exp}}\) instead, obtained using Eq. (16) at each order in \(\hat{\alpha }_s\). We observe that the differences when going from \(O(\hat{\alpha }_s)\) to \(O(\hat{\alpha }_s^2)\) and from \(O(\hat{\alpha }_s^2)\) to \(O(\hat{\alpha }_s^3)\) are of very similar size. This can be traced back to the slow convergence of the coefficients in \(\lambda _1^c(s)\) of Eq. (3). It is interesting to note that this does not lead to large changes in the values of \(\hat{m}_c(\hat{m}_c)\) as is evident from Fig. 3. Note that we do not assign an extra error to \(R_\mathrm{background}(s)\) since any uncertainty in this quantity is relegated to \(\Delta \lambda _3^{c,\mathrm{exp}}\).

We now have a closer look at the determination of the moments from data in the energy range \(2 M_{D^0} \le \sqrt{s} \le 4.8\) GeV. After light-quark background subtraction, we calculate the contributions to the moments in the considered energy range by numerical integration over the data. In Table 7 we show the results for the zeroth, first, and second moments with data taken from Refs. [28, 30, 31, 33]. We display the separate contributions from each individual collaboration split into five different energy intervals. These intervals have been chosen in such a way that none of the published data had to be split up into smaller segments than originally reported by the collaborations. The first errors are statistical and the second systematic. The statistical errors are uncorrelated. The systematic errors are taken uncorrelated among different collaborations, but completely correlated for data from the same collaboration. The last part in Table 7 (labeled ‘Total’) shows the average of the contributions from the different collaborations in each energy interval. The averaging procedure follows Ref. [36]. This prescription takes the relative weights of statistical and systematic errors from each collaboration into account. The breakdown of the total error \(\Delta _\mathrm{tot}\) into statistical (\(\Delta _\mathrm{stat}\)) and systematic (\(\Delta _\mathrm{sys}\)) contributions can then be calculated in terms of the statistical (\(\Delta \mathrm{C}^i_\mathrm{stat}\)) and systematic (\(\Delta \mathrm{C}^i_\mathrm{sys}\)) errors of collaboration i contributing to the given interval as [36]

$$\begin{aligned} (\Delta _\mathrm{stat})^2= & {} \sum _i (\Delta \mathrm{C}^i_\mathrm{stat})^2 \left( \frac{\Delta _\mathrm{tot}}{\sqrt{(\Delta \mathrm{C}^i_\mathrm{stat})^2 + (\Delta \mathrm{C}^i_\mathrm{sys})^2}} \right) ^4,\nonumber \\ \end{aligned}$$
(14)
$$\begin{aligned} (\Delta _\mathrm{sys})^2= & {} \sum _i (\Delta \mathrm{C}^i_\mathrm{sys})^2 \left( \frac{\Delta _\mathrm{tot}}{\sqrt{(\Delta \mathrm{C}^i_\mathrm{stat})^2 + (\Delta \mathrm{C}^i_\mathrm{sys})^2}} \right) ^4.\nonumber \\ \end{aligned}$$
(15)
Table 7 Contributions to the charm moments (\(\times 10^{n}\, \text{ GeV }^{2n}\)) from different energy intervals (given in GeV) and experimental collaborations (CB86 [28], BES02 [30], BES06 [31], CLEO09 [33]). The first errors are statistical and the second systematic. The last lines show the averaged results
Table 8 Contributions to the charm moments (\(\times 10^{n}\, \text{ GeV }^{2n}\)) from the energy range \(2 M_{D^0} \le \sqrt{s} \le 4.8\) GeV. For the results in the column labeled ‘Data’, light-quark contributions have been subtracted using the pQCD prediction at order \(\mathcal{O} (\hat{\alpha }_s^3)\). These entries are obtained from the data displayed in Table 7, taking into account the correlation of systematic errors within each experiment. The third column uses \(\hat{m}_c = 1.272\, \mathrm{GeV}\) and \(\lambda _3^{c,\mathrm{exp}}\) determined by the zeroth experimental moment (see text for details). The last column shows the theoretical prediction for the moments using \(\hat{m}_c = 1.272\, \mathrm{GeV}\) and \(\lambda ^c_3 = 1.23\), not including condensates

The final averaged result for the energy range between \(2 M_{D^0}\) and \(\sqrt{s} = 4.8\) GeV is given in Table 8 up to the fifth moment. The errors shown there are the combined statistical and systematic ones. The third column shows the theoretical moments, again restricted to the region \(2 M_{D^0} \le \sqrt{s} \le 4.8\) GeV, using \(\hat{m}_c = 1.272\) GeV as input, and choosing the parameter \(\lambda _3^c\) to match the zeroth experimental moment,

$$\begin{aligned} \int _{(2M_{D^0})^2}^{(4.8\mathrm{\, GeV})^2} \frac{\mathrm{d} s}{s} R_c^\mathrm{cont}(s) \Big |_{\hat{m}_c = 1.272 \, \mathrm{GeV}} = \mathcal{M}_0^\mathrm{Data} = 0.6367(195),\nonumber \\ \end{aligned}$$
(16)

which gives the value \(\lambda _3^{c,\mathrm{exp}} = 1.34(17)\). Finally, the fourth column shows the theoretical calculation of the moments in the threshold region for \(\hat{m}_c = 1.272\) GeV and using the previous value \(\lambda _3^c = 1.23\) found above. The agreement between these three columns is remarkable and shows that the data in this restricted energy range is well described by \(R_c^\mathrm{cont}(s)\) in Eq. (5).

In the analysis of the previous section only the two narrow resonances below the heavy-quark threshold, i.e., the \(J/\Psi \) and \(\Psi (2S)\), have been treated separately, while the other charmonium states were included in the continuum. We now study the effect of treating also the \(\Psi (3770)\) resonance separately and include it based on the narrow-width approximation.

Table 9 Contributions to the charm moments (\(\times 10^{n}\, \text{ GeV }^{2n}\)) from the energy range \(2 M_{D^0} \le \sqrt{s} \le 3.83\) GeV. For the results in the column labeled ‘Data’, light-quark contributions have been subtracted using the pQCD prediction at order \(\mathcal{O} (\hat{\alpha }_s^3)\). ‘BW’ refers to the Breit–Wigner expression in Eq. (17). The last column shows the theoretical prediction for the moment using \(\hat{m}_c = 1.272\, \mathrm{GeV}\) and \(\lambda ^c_3 = 1.23\), not including condensates

Table 9 shows a comparison of charm moments in the closer vicinity of the \(\Psi (3770)\) resonance. In the column labeled ‘Data’, the moments have been calculated as above, i.e. directly from data, but now restricted to the energy range \(2 M_{D^0} \le \sqrt{s} \le 3.83\) GeV (see the inner plot of Fig. 4) again after subtraction of the light-quark contribution. The upper limit of this energy interval is chosen to cover the \(\Psi (3770)\) resonance completely. In the third column (labeled \(\Psi (3770)\)) we display results obtained assuming that the resonance is not narrow enough to justify using Eq. (4). Instead, we use the full Breit–Wigner form,

$$\begin{aligned} R^\mathrm{\Psi (3770)}_\mathrm{BW} = \frac{9 \Gamma ^e_R \Gamma _R}{\alpha _\mathrm{em}^2(M_R^2)} \, \frac{M_R^2}{(s-M_R^2)^2+ \Gamma _R^2 M_R^2} \end{aligned}$$
(17)

with \(M_{\Psi (3770)} = 3.773\) GeV, \(\Gamma _{\Psi (3770)}^e = 0.262(18)\) keV, and \(\Gamma _R = 27.2 \pm 1.0\) MeV for the total width of the \(\Psi (3770)\) [8]. Even though one cannot assume that the \(\Psi (3770)\) resonance saturates the charm cross section in this narrow energy range, a Breit–Wigner description agrees with the prediction based on Eq. (5) when \(R_c^\mathrm{cont}(s)\) is used with \(\lambda _3^c = 1.23\).

Another way to estimate the impact of the \(\Psi (3770)\) region on the charm quark mass is to include the \(\Psi (3770)\) into our parameterization for \(R_c^\mathrm{Res}(s)\) as a \(\delta \)-function (4) and re-evaluate \(\hat{m}_c(\hat{m}_c)\) and \(\lambda ^c_3\) using the pair of the moments \((\mathcal{M}_0,\mathcal{M}_2)\). We find shifts of \(\Delta \hat{m}_c(\hat{m}_c) = -2.4\) MeV and \(\Delta \lambda ^c_3 = -0.12\). This procedure adds more weight at lower \(\sqrt{s}\); therefore the second moment will increase relative to the zeroth moment. As a compensation, the charm mass will decrease. The uncertainties associated with the \(\Psi (3770)\) resonance parameters are negligible.

Table 10 Contributions to the charm moments (\(\times 10^{n}\, \text{ GeV }^{2n}\)) from the energy range \(2 M_{D^0} \le \sqrt{s} \le 4.8\) GeV evaluated after taking into account a \(2\,\%\) shift of the normalization of the pQCD prediction to fit the data, to be compared with Table 8. In this case we obtained \(\hat{m}_c(\hat{m}_c) = 1.272\) GeV and \(\lambda ^c_3 = 1.15(16)\) after matching with the zeroth experimental moment

One may also raise the question whether the specific prescription for the subtraction of the light-quark background in the charm sub-threshold region should be modified, or whether an additional error should be assigned to this prescription. To check the impact of this effect, we have fitted the pQCD prediction including a free normalization factor to the data between \(\sqrt{s} = 2\) GeV and \(2M_D\). We found a normalization factor of \(1.02 \pm 0.01\) and the moments changed to the values given in Table 10. The charm mass is shifted by less than 1 MeV if this modified light-quark subtraction is used instead of the prescription described above.

4 Alternative fits

The approach described in Sect. 2 requires one to make a choice of a pair of the moments from which to determine a value for the charm mass. The results displayed in Table 6 (see the line labelled ‘Total’) show that the pair \((\mathcal{M}_0,\mathcal{M}_2)\) leads to the smallest error (not taking into account the uncertainty from \(\hat{\alpha }_s\)). A more careful examination reveals that different sources of the errors contribute with different weights. While the pair \((\mathcal{M}_0,\mathcal{M}_1)\) is the least sensitive to the gluon condensate and most strongly affected by the uncertainty of the experimental data from the continuum region, the pair \((\mathcal{M}_0, \mathcal{M}_3)\) has little sensitivity to the continuum data and is more strongly affected by the truncation of the OPE. In this section we investigate the possibility to perform a fit using also more than two moments. In addition, the approach that we discuss now will allow us to take into account uncertainties on the parameters entering the moments in a more systematic way. This will confirm the robustness of our main result in Eq. (13).

Table 11 Fit results using the first three (four) moments with different scenarios for the correlation between the moments (see text)

We define a \(\chi ^2\)-function by

$$\begin{aligned} \chi ^2 = \frac{1}{2} \sum _{n,m} \left( \mathcal{M}_n - \mathcal{M}_n^\mathrm{pQCD} \right) \left( \mathcal{C}^{-1}\right) ^{nm} \left( \mathcal{M}_m - \mathcal{M}_m^\mathrm{pQCD} \right) + \chi ^2_c \nonumber \\ \end{aligned}$$
(18)

where the \(\mathcal{M}_n\) are the sum rule expressions for the moments defined in Eq. (2) for \(n>0\), and Eq. (6) multiplyed by \(3Q_c^2\) for \(n=0\). The corresponding predictions from perturbative QCD, \(\mathcal{M}_n^\mathrm{pQCD}\), were defined in Eq. (8). We will consider different fit scenarios where we include different sets of the moments in the sum of Eq. (18). The coefficients \(\left( \mathcal{C}^{-1}\right) ^{nm}\) are the elements of a correlation matrix which we choose as

$$\begin{aligned} \mathcal{C} = \frac{1}{2} \sum _{n,m} \rho ^{|n-m|} \Delta \mathcal{M}_n^{(3)} \Delta \mathcal{M}_m^{(3)} \end{aligned}$$
(19)

with the truncation error \(\Delta \mathcal{M}_m^{(3)}\) defined in Eq. (12). The factor \(\rho ^{|n-m|}\) serves as an ansatz to include correlations between different moments, where \(\rho \) itself is a free parameter. In addition, we will also consider cases with vanishing correlations. In particular, the correlation between the zeroth and the other moments may be much smaller than the correlations between the higher moments (especially adjacent ones). The term \(\chi ^2_c\) in Eq. (18) imposes constraints on the variation of parameters like the electronic widths in Table 1, the strong coupling constant, \(\hat{\alpha }_s(M_z)^\mathrm{exp} = 0.1182 \pm 0.0016\) [8], and the gluon condensate, \(C_G^\mathrm{exp} = 0.005 \pm 0.005\) GeV\(^4\) [25],

$$\begin{aligned} \chi ^2_c= & {} \left[ \frac{\Gamma ^{e}_{J/\Psi } - \Gamma ^{e,\mathrm{exp}}_{J/\Psi }}{\Delta \Gamma ^{e}_{J/\Psi }}\right] ^2 + \left[ \frac{\Gamma ^{e}_{\Psi (2S)} - \Gamma ^{e,\mathrm{exp}}_{\Psi (2S)}}{\Delta \Gamma ^{e}_{\Psi (2S)}}\right] ^2\\&+\;\left[ \frac{\hat{\alpha }_s(M_z) - \hat{\alpha }_s(M_z)^\mathrm{exp}}{\Delta \hat{\alpha }_s(M_z)}\right] ^2 + \left[ \frac{C_G - C_G^\mathrm{exp}}{\Delta _G}\right] ^2. \end{aligned}$$

Thus, \(\chi ^2\) is a function of six parameters, \(\hat{m}_c(\hat{m}_c)\), \(\lambda _3^c\) and \(\hat{\alpha }_s\) for the continuum part of the moments, the electronic widths \(\Gamma ^e_{J/\Psi }\) and \(\Gamma ^e_{\Psi (2S)}\) for the resonance contributions, and \(C_G\). In addition there is the correlation parameter \(\rho \).

We now determine all or a subset of the parameters by minimizing \(\chi ^2\). First, we observe that without correlations between the moments, the minimal \(\chi ^2\) is small. This is likely an indication that the theoretical moments are indeed correlated. Taking into account correlations as described above, we now determine the correlation parameter \(\rho \) from the condition that \(\chi ^2_\mathrm{min}\) be equal to the number of degrees of freedom for each considered fit scenario. Then we determine the allowed parameter ranges by solving the equation \(\chi ^2 = \chi ^2_\mathrm{min} + 1\).

In Table 11 we collect the results of three typical fit scenarios based on the first three (four) moments. Columns 3,4, and 5 are results of six-parameter fits. In the first case (column 3) we assume that all three moments are correlated as described above. In the second case (column 4) we assume that only the first and the second moments are correlated. For the results in column 5 we have also included the moment \(\mathcal{M}_3\). In this case the uncertainty is slightly smaller, but we have more confidence in lower moments which are less sensitive to the condensates. Moreover, in this scenario the correlation is larger. The errors are obtained by projecting the contour \(\chi ^2 = \chi ^2_\mathrm{min} + 1\) onto each of the six parameters considered. We observe that the results for the charm mass as well as for the continuum parameter \(\lambda _3^c\) are stable. The other four parameters are shifted away from their experimental values only slightly and stay well within their uncertainties. We have performed additional fits in modified scenarios, e.g., using the first three and four moments without correlations. We observe only small shifts in the fitted charm mass of not more than 4 MeV, which is well within our quoted uncertainty. Also fits with only two moments are stable and the results for \(\hat{m}_c\) and \(\lambda _3^c\) agree within errors. A few typical examples are shown in Table 12 where no correlation is assumed.

Table 12 Fit results using pairs of the moments without correlation
Table 13 Error budget for the fit scenario \(\mathcal{M}_0,~(\mathcal{M}_1, \mathcal{M}_2)_\rho \) from one-parameter fits for the parameter displayed in each line and using the other parameters fixed at their values given in column 4 of Table 11. The last column compares this with our main result determined from the pair \((\mathcal{M}_0, \mathcal{M}_2)\) in Table 6

Table 13 contains the error budget for the uncertainty of the charm mass in the fit scenario using the moments for \(n=0\), 1, and 2 with a correlation between \(\mathcal{M}_1\) and \(\mathcal{M}_2\). The errors are from one-parameter fits with all other parameters fixed at their best fit values, given in the fourth column of Table 11.

5 Comparison with previous work

Our result agrees within errors with the previous analysis [5] where \(\hat{m}_c(\hat{m}_c) = 1.289^{+0.040}_{-0.045}\) GeV was found. We traced the moderate numerical difference to the following individual shifts:

  • The resonance parameters have been updated in recent years. In particular, the value of \(\Gamma ^e_{J/\Psi }\) changed from 5.26(37) keV to 5.55(14) keV. The larger electronic widths increase the resonance contribution to the higher moments more than to the zeroth moment. As a consequence of the new values of \(\Gamma ^e_R\) the charm quark mass is shifted by \(-12\) MeV.

  • Both the central value and the uncertainty of \(\hat{\alpha }_s(M_z)\) have changed. In Ref. [5], \(\hat{\alpha }_s(M_z) = 0.1221\) was used, while in the current work we use \(\hat{\alpha }_s(M_z) = 0.1182\). For \(\hat{m}_c(\hat{m}_c)\) this leads to a shift of \(-10\) MeV.

  • In the current work, we use a non-zero central value for the contribution from the gluon condensate. \(C_G = 0.005\) GeV\(^4\) with \(\hat{\alpha }_s(M_z) = 0.1221\) induces a shift of the charm quark mass of \(-2\) MeV.

  • Terms proportional to \(\pi ^2\) (i.e., \(\zeta (2)\)) in the definition of \(\lambda _1^c(s)\) in [5] are in fact not present. Correcting this leads to a shift of \(+3\) MeV for \(\hat{m}_c(\hat{m}_c)\).

  • Progress in theoretical calculations of higher-order terms allows us to include effects of \(\mathcal{O}(\hat{\alpha }_s^3)\) in the perturbative coefficients of the moments, in the definition of \(\lambda _1^c(s)\), and in the running of \(\hat{\alpha }_s\). The combined effect of all \(\mathcal{O}(\hat{\alpha }_s^3)\) terms induces a shift of \(+4\) MeV.

  • There are far more data for the cross section of \(e^+ e^-\) annihilation in the charm threshold region, as well as in the region below the charm threshold. These new data affect only the calibration of the uncertainty, however.

  • The inclusion of the singlet piece in the light-quark background subtraction affects the charm mass only at the sub-MeV level.

  • Shifts induced by possible variations in \(\hat{m}_b(\hat{m}_b)\) amount to at most a few hundred keV.

Table 14 Comparison between Eq. (13) and the result of Ref. [5] where \(\hat{\alpha }_s(M_z) = 0.1221\) was used and the breakdown of the different sources of errors. See the text for details

We summarize the various shifts in the charm quark mass in Table 14 where we also include the corresponding shifts in the continuum parameter \(\lambda _3^c\).

Fig. 6
figure 6

Recent charm quark mass determinations (see text for details and references)

We also find agreement within errors with other charm quark mass determinations in the literature. The result of Ref. [24], \(\hat{m}_c(\hat{m}_c)=1.279(13)\) GeV, obtained from the first moment and using \(\hat{\alpha }_s(M_z) = 0.1189\) and \(C_G = 0.006(12)\) GeV\(^4\), agrees very well with the value we find from the pair \((\mathcal{M}_0, \mathcal{M}_1)\) in Table 6. Reference [38] used the second moment and a modification of the kernel of the sum rule in Eq. (1) to obtain \(\hat{m}_c(\hat{m}_c) = 1.278(9)\) GeV using as input \(\hat{\alpha }_s(M_z) = 0.1184(7)\) and \(C_G = 0.01(1)\) GeV\(^4\). Reference [39] is based on higher moments extracting \(C_G = (0.022 \pm 0.004)\) GeV\(^4\), which then yields \(\hat{m}_c(\hat{m}_c) = 1.261(16)\) GeV. While the charm production cross section is described by a set of six narrow-width resonances, combined with a continuum from perturbative QCD, the difference between Ref. [39] and our result can be explained to a large extent by the different values of the gluon condensate. A more recent analysis [40] extracted the charm quark mass from the first moment sum rule using experimental data from a previous analysis [41]. We can reproduce their final result, \(\hat{m}_c(\hat{m}_c) = 1.288(20)\) GeV, if we use their input values, \(\hat{\alpha }_s(M_z) = 0.1184(21)\) and \(C_G = 0.006(12)\) GeV\(^4\), and the electronic widths of the \(\Psi (2S)\) resonance before its recent update. The authors of Ref. [41] have also discussed the Z-boson contribution to the QCD moments and found it to be very small. This contribution is neglected in our work.

Beyond the phenomenological charm quark mass determinations discussed so far, there are also precise lattice QCD calculations. Reference [42] reports \(\hat{m}_c(\hat{m}_c) = 1.2715(95)\) GeV with \(\hat{\alpha }_s(M_z) = 0.11822(74)\) from a study of pseudoscalar–pseudoscalar correlators. Reference [43], using Möbius domain-wall fermions to describe both light and charm quarks, obtained the value \(\hat{m}_c(\hat{m}_c) = 1.2871(106)\)GeV for \(\hat{\alpha }_s(M_z) = 0.1177(25)\). A value for the gluon condensate was also fitted and found to be compatible with zero. Finally, Ref. [44] used the Highly Improved Staggered Quark action and moments of the pseudoscalar–charmonium correlator, and found \(\hat{m}_c(\hat{m}_c) = 1.267(11)\) GeV and \(\hat{\alpha }_s(M_z) = 0.1162(75)\). These results are displayed in Fig. 6 together with our final result for \(\hat{\alpha }_s(M_Z) = 0.1182(16)\) [8].

6 Conclusions

We presented a determination of the charm quark mass from QCD sum rules with a careful determination of the uncertainty. Our analysis is based on [5], but goes beyond that previous work in several respects: first, on the theoretical side, we were now in a position to include terms of one order higher than previously. Secondly, more precise experimental data are available now, including values for the strong coupling \(\alpha _s\) and for the gluon condensate. In addition our new approach for the estimate of uncertainties allows us to take into account the effect from error correlations. This should lead to a more reliable uncertainty estimate than is possible in other approaches.

In contrast to other analyses we included the zeroth moment in our strategy and exploited the consistency between different moments. The usefulness of the zeroth moment is not its sensitivity to \(\hat{m}_c\) itself (which is mild), but to the continuum contribution to which it is very sensitive. While the coefficient \(A_3\) entering the zeroth moment through Eq. (6) is uncomfortably large, we fully account for this by our conservative error estimate in Eq. (12). It is interesting to note, however, that the dispersive calculation of the value of the electromagnetic coupling \(\alpha _\mathrm{em}\) at the Z-pole (an essential input into the calculation of the W boson mass, as well as for the electroweak program at LEP and any future Z factory), effectively relies on precisely the zeroth moment. One usually neglects making any reference to what value of \(\hat{m}_c\) the charm quark contribution to \(\alpha _\mathrm{em}(M_Z)\) would correspond to and whether that value is compatible with other determinations. It is this relation that shows poor convergence. Our approach amounts to determining the error on \(\hat{m}_c\) by demanding compatibility with some higher moment. Our value of \(\hat{m}_c\) in Eq. (13) can then be used along the lines of Ref. [37], where perturbative QCD was used to directly compute \(\alpha _\mathrm{em}(M_Z)\) in the \(\overline{\mathrm{MS}}\) scheme (with \(\hat{m}_c\) as an external input).

Note that at this level of accuracy a subtle issue arises: we have implicitly defined \(\hat{m}_c(\hat{m}_c)\) as the \(\overline{\mathrm{MS}}\) mass w.r.t. both QCD and QED, as can be seen e.g. from Eq. (9). Alternatively adopting a hybrid definition, \(\hat{m}_c(\hat{m}_c)^\mathrm{(pole-QED)}\), which assumes the on-shell mass w.r.t. QED so that the QED treatment of quarks mirrors more closely that of leptons, one has

$$\begin{aligned} \hat{m}_c(\hat{m}_c)^\mathrm{(pole-QED)}= & {} \hat{m}_c(\hat{m}_c) \left[ 1+ Q_c^2 \frac{\alpha _{em}(m_c)}{\pi } + \mathcal{O}(\alpha _{em}^2) \right] \nonumber \\= & {} (1274 + 2616 \Delta \hat{\alpha }_s \pm 8)~\text{ MeV }. \end{aligned}$$
(20)

The largest contribution to the total uncertainty in the charm quark mass is the truncation error. Since the coefficients of the perturbative expansion of the moments exhibit no compelling evidence for a convergent behavior, we consider it unlikely that knowledge of the next terms in the perturbative series will lead to a reduction of the error on the charm mass.

In future work we will apply our approach also to the case of the bottom quark. We believe that we should be able to justify an error estimate for \(\hat{m}_b(\hat{m}_b)\) of about 15 MeV. This follows from the fact that our determination of the charm mass scales with the difference between half of the lowest-mass resonance, \(M_{J/\Psi }/2\), and the quark mass \(\hat{m}_c(\hat{m}_c)\). Our uncertainty of \(\Delta \hat{m}_c(\hat{m}_c) = 8\) MeV corresponds to \(2.9\,\%\) of this difference. Translated to the case of the b quark, using \(M_{\Upsilon } = 9.460\) GeV and \(\hat{m}_b(\hat{m}_b) \sim 4.200\) GeV, this would correspond to 15 MeV.