---
author:
- |
    Bahman Dehnadi ${}^a$\
    E-mail:
- |
    Andre H. Hoang ${}^{b,c}$\
    E-mail:
- |
    Vicent Mateu ${}^{c,d,e}$\
    E-mail:
- |
    S. Mohammad Zebarjad ${}^a$\
    E-mail:\
    ${}^a$ Shiraz University, Physics Department, Shiraz 71454, Iran.\
    ${}^b$ University of Vienna, Faculty of Physics, Boltzmanngasse 5, A-1090 Vienna, Austria.\
    ${}^c$ Max-Planck-Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, D-80805 München, Germany.\
    ${}^d$ Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139.\
    ${}^e$ Instituto de Física Corpuscular, UVEG - Consejo Superior de Investigaciones Científicas, Apartado de Correos 22085, E-46071, Valencia, Spain.
bibliography:
- 'charm2.bib'
title: 'Charm Mass Determination from QCD Charmonium Sum Rules at Order $\alpha_s^3$'
---

Introduction {#sectionintroduction}
============

Accurate determinations of the charm quark mass are an important ingredient in the prediction of inclusive and radiative $B$ decays or exclusive kaon decays such as $K\to\pi\nu\bar{\nu}$. Since these decays are instruments to either measure CKM matrix elements or to search for new physics effects, appropriate and realistic estimates of the uncertainties are also an important element of these analyses [@Antonelli:2009ws].

One of the most powerful methods to determine the charm quark mass is based on sum rules for the charm-anticharm production rate in $e^{+}e^{-}$ annihilation [@Novikov:1977dq]. Here, moments of the correlation function of two charm vector currents at zero momentum transfer $$\begin{aligned}
\label{momentdef1}
M_{n}^{{\rm th}} & = & \dfrac{12\pi^2 Q_c^2}{n!}\,\dfrac{{\rm d}}{{\rm
    d}q^{2n}}\left.\Pi(q^{2})\right|_{q^{2}=0}\,,\\
\left(g_{\mu\nu}q^{2}-q_{\mu}q_{\nu}\right)\,\Pi(q^{2}) 
& = & 
-\, i\int\mathrm{d}x\, e^{iqx}\left\langle \,0\left|T\,
    j_{\mu}(x)j_{\nu}(0)\right|0\,\right\rangle 
\,,\nonumber \\[2mm]
j^{\mu}(x) 
& = & 
\bar{\psi}(x)\gamma^{\mu}\psi(x)
\,,\nonumber \end{aligned}$$ $Q_c$ being the charm quark electric charge, can be related to weighted integrals of the normalized charm cross section $$\begin{aligned}
\label{momentdef2}
M_{n} & = & 
\int\dfrac{{\rm d}s}{s^{n+1}}R_{e^{+}e^{-}\to\, c\bar{c}\,+X}(s)\,,\\
R_{e^{+}e^{-}\to\, c\bar{c}\,+X}(s) & = & \dfrac{\sigma_{e^{+}e^{-}\to\, c\bar{c}\,+X}(s)}{\sigma_{e^{+}e^{-}\to\,\mu^{+}\mu^{-}}(s)}\,,\nonumber \end{aligned}$$ which can be obtained from experiments. For small values of $n$ such that $m_{c}/n\gtrsim\Lambda_{{\rm QCD}}$ the theoretical moments $M_{n}^{{\rm th}}$ can be computed in an operator product expansion (OPE) where the dominant part is provided by perturbative QCD supplemented by small vacuum condensates that parametrize nonperturbative effects [@Shifman:1978bx; @Shifman:1978by]. The leading gluon condensate power correction term has a surprisingly small numerical effect and is essentially negligible for the numerical analysis as long as $n$ is small.

This allows to determine the charm mass in a short distance scheme such as $\overline{{\rm MS}}$ to high precision. This method to determine the $\overline{{\rm MS}}$ charm mass is frequently called charmonium sum rules. For the theoretical moments the perturbative part of the OPE is known at ${\mathcal O}(\alpha_{s}^{0})$ and ${\mathcal O}(\alpha_{s})$ for any value of $n$ [@Kallen:1955fb]. At ${\mathcal O}(\alpha_{s}^{2})$ the first 30 moments are known [@Boughezal:2006uu; @Maier:2007yn], and to ${\mathcal O}(\alpha_{s}^{3})$ for $n=1$ [@Chetyrkin:2006xg; @Boughezal:2006px], $n=2$ [@Maier:2008he], and $n=3$ [@Maier:2009fz]. Higher moments at ${\cal O}(\alpha_s^3)$ have been determined by a semianalytical procedure [@Hoang:2008qy; @Kiyo:2009gb] (see also [@Greynat:2010kx]). The Wilson coefficient of the gluon condensate contribution is known to ${\mathcal O}(\alpha_{s})$ [@Broadhurst:1994qj]. On the experimental side the total hadronic cross section in $e^{+}e^{-}$ annihilation is known from various experimental measurements for c.m. energies up to $10.538\,$GeV. None of the experimental analyses actually ranges over the entire energy region between the charmonium region and $10.538$ GeV, but different analyses overlapping in energy exist such that energies up to $10.538\,$GeV are completely covered [@Bai:1999pk; @Bai:2001ct; @Ablikim:2004ck; @Ablikim:2006aj; @Ablikim:2006mb; @:2009jsa; @Osterheld:1986hw; @Edwards:1990pc; @Ammar:1997sk; @Besson:1984bd; @:2007qwa; @CroninHennessy:2008yi; @Blinov:1993fw; @Criegee:1981qx; @Siegrist:1976br; @Rapidis:1977cv; @Abrams:1979cx; @Siegrist:1981zp].[^1] Interestingly, to the best of our knowledge, the complete set of all available experimental data on the hadronic cross section has never been used in previous charmonium sum rule analyses to determine the experimental moments. Rather, sum rule analyses have relied heavily on theoretical input using different approaches to determine the corresponding “experimental error” and intrinsically leading to a sizable modeling uncertainty for energy regions below $10.538\,$GeV for low values of $n$ [@Hoang:2004xm].

The most recent charmonium sum rule analysis based on Eqs. (\[momentdef1\]) and (\[momentdef2\]), carried out by Kühn et al. [@Chetyrkin:2009fv; @Kuhn:2007vp] using input from perturbative QCD (pQCD) at ${\mathcal O}(\alpha_{s}^{3})$ for the perturbative contribution, obtained $\overline{m}_{c}(\overline{m}_{c})=
1279\pm(2)_{{\rm pert}}\pm(9)_{{\rm exp}}
\pm(9)_{\alpha_{s}}\pm(1)_{{\rm \left\langle GG\right\rangle}}\,$ MeV where the first error is the perturbative uncertainty and the second is the experimental one. The third and the fourth uncertainties come from $\alpha_s$ and the gluon condensate correction, respectively. To our knowledge this result, the outcome of similar analyses in Ref. [@Chetyrkin:2006xg] and by Boughezal, Czakon and Schutzmeier [@Boughezal:2006px][^2], and a closely related analysis based on lattice results instead of data for pseudoscalar moments [@Allison:2008xk; @McNeile:2010ji] represent the analyses with the highest precision achieved so far in the literature. If confirmed, any further investigations and attempts concerning a more precise charm quark $\overline{{\rm MS}}$ mass would likely be irrelevant for any foreseeable future.

We therefore find it warranted to reexamine the charmonium sum rule analysis with special attention on the way how perturbative and experimental uncertainties have been treated in Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp]. A closer look into their analysis reveals that the quoted perturbative uncertainty results from a specific way to arrange the $\alpha_s$ expansion for the charm mass extractions and, in addition, by setting the $\overline{{\rm MS}}$ renormalization scales in $\alpha_{s}$ and in the charm mass (which we call $\mu_{\alpha}$ and $\mu_{m}$, respectively) equal to each other (i.e., they use $\mu_{\alpha}=\mu_{m}$). Moreover, concerning the experimental moments, only data up to $\sqrt{s}=4.8\,$GeV from the BES experiments [@Bai:2001ct; @Ablikim:2006mb] were used, while for $\sqrt{s}>4.8$ GeV perturbative QCD predictions were employed. Conceptually this approach is somewhat related to the method of finite energy sum rules (see e.g. Ref. [@Penarrocha:2001ig]), which we, however, do not discuss in this work. While this approach might be justified to estimate the overall nominal contribution for the experimental moments from $\sqrt{s}>4.8$ GeV, since perturbative QCD predictions describe quite well the measured total hadronic cross section outside the resonance regions, it is certainly not a suitable method to determine the experimental uncertainty. Since the region $\sqrt{s}>4.8$ GeV constitutes about $30\%$ of the first moment $M_{1}$, which is theoretically most reliable, this approach contains a significant intrinsic model dependence that cannot be quantified unambiguously.

In this work we reexamine the charmonium sum rules analysis for low values of $n$ using the latest ${\mathcal O}(\alpha_{s}^{3})$ perturbative results, and we implement improvements which concern the two issues just mentioned:

1.  We analyze several different types of perturbative expansions and examine in detail how the result for the $\overline{{\rm MS}}$ charm mass depends on independent choices of $\mu_{\alpha}$ and $\mu_{m}$. We show in particular that the interplay of the perturbative expansion and the scale setting $\mu_{\alpha}=\mu_{m}$ used in previous ${\mathcal O}(\alpha_s^3)$ analyses leads to sizable cancellations of the dependence on $\mu_{\alpha}$ and $\mu_{m}$ that in the light of our new analysis has to be considered as accidental. As the outcome of our analysis we estimate the current ${\mathcal O}(\alpha_{s}^{3})$ perturbative error as around $20$ MeV, which is an order of magnitude larger than that of Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp; @Chetyrkin:2006xg; @Boughezal:2006px].

2.  Using a clustering method [@D'Agostini:1993uj; @Takeuchi:1995xe; @Hagiwara:2003da] to combine correlated data from many different experimental measurements we show that the $e^{+}e^{-}$ total hadronic cross section relevant for the charmonium sum rules can be determined with a complete coverage of center of mass energies above the $J/\psi$ and $\psi'$ resonances up to $10.538\,$GeV. Conservatively estimated modeling uncertainties coming from the energy range above $10.538\,$GeV then only lead to an insignificant contribution to the total uncertainty of the experimental moments.

This paper is organized as follows: In Sec. \[sectiontheory\] we introduce the theoretical framework and review the current status of perturbative computations. We also show various equivalent ways of arranging the perturbative series in $\alpha_s$ for the charm mass. Finally we discuss how to properly estimate theoretical uncertainties due to the truncation of the perturbative series. In Sec. \[sectiondata\] we present all the experimental information that goes into our analysis. We discuss a clustering fit procedure that allows to combine data from different experiments accounting for their correlation and show the results. In Sec. \[sectionanalysis\] we carry out the numerical charm mass analysis concentrating on the first moment using arbitrary values of $\alpha_s$, and we present our final charm mass result. In Appendix A we present more details on the outcome of our clustering fit procedure for the charm , and in Appendix B we prove the equivalence of different versions of $\chi^2$ functions when auxiliary fit parameters are employed.

Theoretical Input {#sectiontheory}
=================

Perturbative Contribution {#subsectionperturbative}
-------------------------

The moments of the vector current correlator are defined in Eq. (\[momentdef1\]). Their perturbative contribution in the framework of the OPE has a non-linear dependence on the charm quark mass. Thus in principle no conceptual preference can be imposed on any of the possible perturbative series that arises when solving for the charm mass. As a consequence, different versions of the expansion should be considered to obtain reliable estimates of the perturbative uncertainty. As indicated in Sec. \[sectionintroduction\] we use in the following $\mu_\alpha$ as the renormalization scale in $\alpha_s$ and $\mu_m$ as the renormalization scale in the $\overline{\rm MS}$ charm quark mass $\overline m_c$.

[**(a) Standard fixed-order expansion**]{}\
Writing the perturbative vacuum polarization function as $$\begin{aligned}
\label{Mnpertfixedorder1}
\Pi^{\rm pert}(q^2,\alpha_s(\mu_\alpha), \overline m_c(\mu_m), 
\mu_\alpha, \mu_m) 
\, = \, \dfrac{1}{12\pi^2 Q_c^2}\sum_{n=0}^\infty q^{2n} M_n^{\rm pert}
\,,\end{aligned}$$ we have for the perturbative moments $M_n^{\rm pert}$ $$\begin{aligned}
\label{Mnpertfixedorder2}
M_n^{\rm pert} & = &
\frac{1}{(4\overline m^2_c(\mu_m))^{n}}
\sum_{i,a,b} \left(\frac{\alpha_s(\mu_\alpha)}{\pi}\right)^i
C^{a,b}_{n,i}\,\ln^a\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_m}\right)   
\ln^b\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_\alpha}\right).\end{aligned}$$ This is the standard fixed-order expression for the perturbative moments. At ${\cal O}(\alpha_s^3)$ the coefficients $C^{0,0}_{n,3}$ were recently determined for $n=1$ [@Chetyrkin:2006xg; @Boughezal:2006px], $n=2$ [@Maier:2008he], $n=3$ [@Maier:2009fz] and higher [@Hoang:2008qy; @Kiyo:2009gb; @Greynat:2010kx]. We refer to Ref. [@Boughezal:2006uu; @Maier:2007yn] for the coefficients at ${\cal O}(\alpha_s^2)$. We have summarized the numerical expressions for the $C^{a,b}_{n,i}$ coefficients of the first four moments in Tab. \[tabcfixedorder\].

The standard fixed-order expansion in Eq. (\[Mnpertfixedorder2\]) is the common way to represent the perturbative moments. However, written in this form the non-linear dependence on $\overline m_c$ does for some values of the experimental moments and the renormalization scales not yield numerical solutions for $\overline m_c$.[^3]

[**(b) Linearized expansion**]{}\
Concerning the charm mass dependence, a more linear way to organize the perturbative expansion is to take the root of Eq. (\[Mnpertfixedorder2\]): $$\begin{aligned}
\label{Mnpertlinearized1}
\Big(M_n^{\rm th, pert}\Big)^{1/2n} & = &
\frac{1}{2\overline{m}_c(\mu_m)} \,\sum_{i,a,b}\left(\frac{\alpha_s(\mu_\alpha)}{\pi}\right)^i \tilde C_{n,i}^{a,b}
\ln^a\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_m}\right)   \ln^b\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_\alpha}\right)
,\end{aligned}$$ or equivalently $$\begin{aligned}
\label{Mnpertlinearized2}
\overline m_c(\mu_m) & = & \frac{1}{2\Big(M_n^{\rm th, pert}\Big)^{1/2n}} \sum_{i,a,b}\,\left(\frac{\alpha_s(\mu_\alpha)}{\pi}\right)^i\tilde C_{n,i}^{a,b}
\ln^a\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_m}\right)   \ln^b\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_\alpha}\right)
.\end{aligned}$$ The coefficients $\tilde C_{n,i}^{a,b}$ using again $\mu_\alpha$ for the renormalization scale in $\alpha_s$ and $\mu_m$ for the renormalization scale in the $\overline{\rm MS}$ charm mass are given in Tab. \[tabctildefixedorder\]. Although relation (\[Mnpertlinearized2\]) involves a non-linear dependence on ${\overline m}_c$, we find that it always has a numerical solution.

[**(c) Iterative linearized expansion**]{}\
For the standard and the linearized expansions in Eqs. (\[Mnpertfixedorder2\]) and (\[Mnpertlinearized1\]) one searches for numerical solutions of the charm mass $\overline m_c(\mu_m)$ keeping the exact mass dependence on the respective RHS at each order in $\alpha_s$. An alternative way that is consistent within perturbation theory is to solve for $\overline
m_c(\mu_m)$ iteratively order-by-order supplementing appropriate lower order values for $\overline m_c(\mu_m)$ in higher order perturbative coefficients. To be more explicit, we describe the method in the following. As the basis for the iterative expansion carried out in our analysis we use the linearized expansion of Eq. (\[Mnpertlinearized2\]).

In the first step we determine $\overline m_c(\mu_m)$ employing the tree-level relation $$\begin{aligned}
\label{Mnpertiterative1}
\overline m_c^{(0)}(\mu_m) =
\frac{1}{2\Big(M_n^{\rm th, pert}\Big)^{1/2n}} \,\tilde C_{n,0}^{0,0}
\,,\end{aligned}$$ giving the tree-level charm mass $\overline m_c^{(0)}(\mu_m)$. In the next step one employs the relation $$\begin{aligned}
\label{Mnpertiterative2}
\overline m_c^{(1)}(\mu_m) &=&
\frac{1}{2\Big(M_n^{\rm th, pert}\Big)^{1/2n}}\,\Biggr\{\,
\tilde C_{n,0}^{0,0} \, + \,
\frac{\alpha_s(\mu_\alpha)}{\pi}\bigg[\tilde C_{n,1}^{0,0}+\tilde C_{n,1}^{1,0}
\ln\left(\frac{\overline{m}_c^{(0)\,2}(\mu_m)}{\mu^2_m}\right)\Bigg]\Biggr\}
\,,\quad\end{aligned}$$ to determine the ${\cal O}(\alpha_s)$ charm mass $\overline m_c^{(1)}(\mu_m)$. In the ${\cal O}(\alpha_s)$ terms on the RHS of Eq. (\[Mnpertiterative2\]) the tree-level charm mass $\overline m_c^{(0)}$ is used, which is consistent to ${\cal O}(\alpha_s)$. At ${\cal O}(\alpha_s^2)$ for the determination of $\overline m_c^{(2)}(\mu_m)$ one uses $\overline m_c^{(0)}$ for the ${\cal
  O}(\alpha_s^2)$ coefficient and $\overline m_c^{(1)}(\mu_m)$ for the ${\cal
  O}(\alpha_s)$ correction, which in the strict $\alpha_s$ expansion yields $$\begin{aligned}
\label{Mnpertiterative3}
\overline m_c^{(2)} & = &
\frac{1}{2\Big(M_n^{\rm th, pert}\Big)^{1/2n}}\,\Biggr\{\,
\tilde C_{n,0}^{0,0} \, +\,
\frac{\alpha_s(\mu_\alpha)}{\pi}\Bigg[\tilde C_{n,1}^{0,0}+\tilde C_{n,1}^{1,0}
\ln\left(\frac{\overline{m}_c^{(0)\,2}(\mu_m)}{\mu^2_m}\right)\Bigg]+\nonumber\\
&&\left(\frac{\alpha_s(\mu_\alpha)}{\pi}\right)^2\Bigg[2\,\dfrac{\tilde C_{n,1}^{1,0}\,\tilde C_{n,1}^{0,0}}{\tilde C_{n,0}^{0,0}}+
2\,\dfrac{(\tilde C_{n,1}^{1,0})^2}{\tilde C_{n,0}^{0,0}}\ln\left(\frac{\overline{m}_c^{(0)\,2}(\mu_m)}{\mu^2_m}\right)+\nonumber\\
&&\sum_{a,b}\,\tilde C_{n,2}^{a,b}
\ln^a\left(\frac{\overline{m}_c^{(0)\,2}(\mu_m)}{\mu^2_m}\right)   \ln^b\left(\frac{\overline{m}_c^{(0)\,2}(\mu_m)}{\mu^2_\alpha}\right)\Bigg]\Biggr\}
\,.\end{aligned}$$ Here the second line contains the derivative of the ${\mathcal O}(\alpha_s)$ terms with respect to the charm mass. The determination of the ${\cal O}(\alpha_s^3)$ charm mass $\overline 
m_c^{(3)}(\mu_m)$ is then carried out in an analogous way involving the second (first) derivative with respect to the mass in the ${\cal O}(\alpha_s)$ (${\cal
  O}(\alpha_s^2)$) correction and using again $\overline m_c^{(0)}(\mu_m)$ for the ${\cal O}(\alpha_s^3)$ coefficient.

In general we can write the iterative expansion as follows: $${\overline m}_c(\mu_m)\,=\,{\overline m}_c^{(0)}
\sum_{i,a,b}\left(\frac{\alpha_s(\mu_\alpha)}{\pi}\right)^i \hat C_{n,i}^{a,b}\,
\ln^a\left(\frac{{\overline m}_c^{(0)\,2}}{\mu^2_m}\right)   
\ln^b\left(\frac{{\overline m}_c^{(0)\,2}}{\mu^2_\alpha}\right),
\label{eq:iterative-general}$$ where the numerical value of the coefficients $\hat C_{n,i}^{a,b}$ are collected in Tab. \[tab:chat\].

The iterative way to treat the perturbative series for the charm mass has the advantage that solving for the charm mass involves equations that are strictly linear in the charm mass at any order of the $\alpha_s$ expansion and thus always have solutions. In this way any possible influence on the analysis arising from a non-linear dependence is eliminated.

[**(d) Contour improved expansion**]{}\
For the expansion methods (a)-(c) the moments and the charm quark mass are computed for a fixed choice of the renormalization scale $\mu_\alpha$ in the strong coupling $\alpha_s$. In analogy to the contour improved methods used for (see e.g. Refs. [@LeDiberder:1992te; @Pivovarov:1991rh; @Braaten:1991qm; @Narison:1988ni; @Braaten:1988ea; @Braaten:1988hc]) one can employ a path-dependent $\mu_\alpha$ in the contour integration that defines the perturbative moments [@Hoang:2004xm], see Fig. \[figcontour\], $$\begin{aligned}
\label{Mnpertcontour1}
M_n^{\rm c, pert} & = &
\frac{6\pi Q_c^2}{i}\,\int_c\,\frac{{\rm d}s}{s^{n+1}}
\Pi(q^2, \alpha_s(\mu_\alpha^c(s,\overline m_c^2)), \overline{m}_c(\mu_m),
\mu_\alpha^c(s,\overline m_c^2), \mu_m)
\,.\end{aligned}$$

![ Path of integration in the complex for the computation of the moments. \[figcontour\] ](contour.eps){width="30.00000%"}

Due to the independence of the moments on $\mu_\alpha$ and since no large logarithms are being generated anywhere for a path with distance of order $\overline m_c$ from the cut on the real axis, this method is a viable alternative to carry out the perturbative expansion. The different orders in the expansion of the contour improved moments $M_n^{\rm c, pert}$ are generated from the fixed-order $\alpha_s$ expansion of the vacuum polarization function $\Pi$ in Eq. (\[Mnpertcontour1\]).

A useful path-dependent choice for $\mu_\alpha^c$ is given by [@Hoang:2004xm] $$\begin{aligned}
\label{mualphacontour}
(\mu_\alpha^c)^2(s,\overline m_c^2) & = &
\mu_\alpha^2\,\bigg(\,1-\frac{s}{4\overline m_c^2(\mu_m)}\,\bigg) 
\,,\end{aligned}$$ which implements a modified weighting of threshold versus high energy contributions. It is straightforward to prove that the resulting moments $M_n^{\rm c,
  pert}$ can be obtained from the small-$q^2$ expansion of the perturbative vacuum polarization function using $\mu_\alpha^c$ as the renormalization scale of $\alpha_s$, $$\begin{aligned}
\label{Mnpertcontour2}
\Pi^{\overline {\rm MS}}\Big(q^2,
\alpha_s(\mu_\alpha^c(q^2,\overline m_c^2)),
\overline m_c(\mu_m), 
\mu_\alpha^c(q^2,\overline m_c^2),
\mu_m\Big) & = &
\sum\limits_{n=0}^\infty \, q^{2n}\,M_n^{\rm c, pert}
\,.\end{aligned}$$ From Eq. (\[Mnpertcontour2\]) we can see that the $M_n^{\rm c, pert}$ can be derived from the expressions for the fixed-order moments $M_m^{\rm pert}$ with $m\le n$ given in Eq. (\[Mnpertfixedorder2\]).[^4] They also depend on the QCD and its derivatives, which arise in the small-$q^2$ expansion of $\alpha_s(\mu^c_\alpha(q^2,\overline m_c^2))$. Note that the that has to be employed must be exactly the same that is used for the contour integration in Eq. (\[Mnpertcontour1\]). Expanding the dependence of the on $\alpha_s$ strictly in fixed-order one recovers the fixed-order moments $M_n^{\rm pert}$. So the dependence of the contour improved moments $M_n^{\rm c, pert}$ on the fixed-order moments $M_m^{\rm pert}$ with $m<n$ is only residual due to the truncation of the $\alpha_s$ series, representing yet another alternative parametrization of higher order perturbative corrections. The contour improved moments do in particular have a residual dependence on the UV-subtraction scheme for the vacuum polarization function, i.e. on $\Pi(0)=M_0^{\rm pert}$. Using the “on-shell” scheme with $\Pi(0)=0$ one finds that $M_1^{\rm c,pert} = M_1^{\rm pert}$. For our analysis we employ the $\overline{\rm MS}$ scheme for $\Pi(0)$ defined for $\mu={\overline m}_c({\overline m}_c)$. Expressed in terms of $\alpha_s(\mu_\alpha)$ and ${\overline m}_c(\mu_m)$ it has the form [@Chetyrkin:2006xg] $$\begin{aligned}
\label{Pi0msbar}
\Pi^{\overline {\rm MS}}(0) & = &
\sum_{i,a,b} \left(\frac{\alpha_s(\mu_\alpha)}{\pi}\right)^i
C^{a,b}_{0,i}\,\ln^a\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_m}\right)   \ln^b\left(\frac{\overline{m}^2_c(\mu_m)}{\mu^2_\alpha}\right)
\,.\end{aligned}$$ The numerical values for the coefficients $C^{a,b}_{0,i}$ can be found in Tab. \[tabcPi0\]. Using contour improved moments to determine $\overline m_c$ also involves non-linear relations, which implies that in some cases there is no solution. Again this can happen for $n=>1$ and for $\mu_m\sim 3$ GeV at any order.

Gluon Condensate Contribution {#subsectioncondensate}
-----------------------------

The dominant subleading contribution in the OPE for the theory moments $M_n^{\rm
  th}$ is arising from the gluon condensate. Writing $$\begin{aligned}
\label{MnOPE1} 
M_n^{\rm th} & = & M_n^{\rm pert} + \Delta M_n^{\langle G^2\rangle}\,+\,\ldots
\,,\end{aligned}$$ and neglecting any higher order power-suppressed condensate contributions of the OPE we have [@Novikov:1977dq; @Baikov:1993kc; @Broadhurst:1994qj] $$\label{DeltaMncondensate1}
\Delta M_n^{\langle G^2\rangle}
 = \dfrac{1}{(4\overline{m}_c^2(\mu_m))^{n+2}}\Big\langle\frac{\alpha_s}{\pi} G^2\Big\rangle_{\rm RGI}  \left\{ a^{0,0}_{n}+\dfrac{\alpha_{s}(\mu_\alpha)}{\pi}\left[a^{1,0}_{n}+a^{1,1}_n \log\dfrac{\overline{m}_c^2(\mu_m)}{\mu_m^{2}}\right]\right\}\,.$$ We have written Eq. (\[DeltaMncondensate1\]) using the renormalization group invariant (RGI) scheme for the gluon condensate [@Narison:1983kn] and the $\overline{\rm
  MS}$ charm quark mass. In Ref. [@Chetyrkin:2010ic] it is suggested that the Wilson coefficient of the gluon condensate should be expressed in terms of the pole rather than the $\overline{\mbox{MS}}$ mass based on the observation that the pole mass leads to a condensate correction that is numerically quite stable for higher moments. We disagree with this prescription due to the infrared renormalon issue in the pole mass and since it is expected that the condensate corrections for the moment increases with $n$. On the other hand, for a sum rule analysis of the first moment with $n=1$ this issue is absolutely irrelevant due to the small size of the gluon condensate contribution compared to the level of the remaining perturbative uncertainties. Numerical values for the coefficients $a_n^{i,j}$ are given in Tab. \[tabgluoncondensate\] for $n=1,2,3,4$. For the RGI gluon condensate we adopt [@Ioffe:2005ym] $$\begin{aligned}
\label{condensatevalue1}
\Big\langle\frac{\alpha_s}{\pi} G^2\Big\rangle_{\rm RGI} 
& = & 0.006\pm0.012\;\mathrm{GeV}^4\,.\end{aligned}$$ The overall contribution of the gluon condensate correction in Eq. (\[DeltaMncondensate1\]) in the charm quark mass analysis is quite small. Its contribution to the moments amounts to around $0.2\%$, $0.6\%$, $2\%$, and $3\%$ for the first four moments, respectively. For $n=1$ it leads to a correction in the $\overline{\rm MS}$ charm quark mass at the level of $1$ MeV and is an order of magnitude smaller than our perturbative uncertainty. We therefore ignore the condensates correction for the discussion of the perturbative uncertainties in Sec. \[subsectionmcerror\]. Its contribution is, however, included in the final charm mass results presented in Sec. \[sectionanalysis\].

Running Coupling and Mass {#subsectionrunning}
-------------------------

The analysis of the charmonium sum rules naturally involves renormalization scales around the charm mass, $\mu\sim {\overline m}_c\sim 1.3$ GeV, which are close to the limits of a perturbative treatment. In fact, parametrically, the typical scale relevant for the perturbative computation of the moment $M_n^{\rm th}$ is of order $\mu\sim {\overline m}_c/n$ (see e.g. Ref. [@Hoang:1998uv]) because the energy range of the smearing associated to the weight function $1/s^{n+1}$ in Eq. (\[momentdef2\]) decreases with $n$. We will therefore use $n=1$ for our final numerical analysis. Moreover, it is common practice to quote the $\overline{\rm MS}$ charm mass $\overline m_c(\overline m_c)$, i.e. for the scale choice $\mu_m=\overline m_c$. It is therefore useful to have a look at the quality of the perturbative behavior of the renormalization group evolution of the strong $\overline{\rm MS}$ coupling $\alpha_s$ and the $\overline{\rm MS}$ charm quark mass.

In Fig. \[fig:alphasevolutionexact\][^5] we have displayed $\alpha_s^{\rm N^3LL}(\mu)/\alpha_s^{\rm N^kLL}(\mu)$ using for $\alpha_s^{\rm N^kLL}$ the QCD and the respective exact numerical solution for $\alpha_s(3~\mbox{GeV})= 0.2535$ as the common reference point. We see that the convergence of the lower order results towards the 4-loop evolution is very good even down to scales of around $1$ GeV. The curves indicate that the remaining relative perturbative uncertainty in the 4-loop evolution might be substantially smaller than $1\%$ for scales down to $\overline m_c\sim 1.3$ GeV. It is also instructive to examine the evolution using a fixed-order expansion. In Fig. \[fig:alphasevolutionexpanded\] we display $\alpha_s^{\rm N^3LL}(\mu)/\alpha_s^{(m)}(\mu)$ where $\alpha_s^{(m)}(\mu)$ is the ${\cal O}(\alpha_s^{m+1})$ fixed-order expression for $\alpha_s(\mu)$ using the reference value $\alpha_s(3~\mbox{GeV})=0.2535$ as the expansion parameter. The convergence of the fixed-order expansion for $\alpha_s(\mu)$ towards the exact N${}^3$LL numerical solution $\alpha_s^{\rm N^3LL}(\mu)$ is somewhat worse compared to the renormalization group resummed results since the deviation of the ratio from one is in general larger. However, convergence is clearly visible. In particular there are not any signs of instabilities. It therefore seems to be safe to use renormalization scales down to the charm mass and associated renormalization scale variations as an instrument to estimate the perturbative uncertainties.

In Figs. \[fig:mcevolutionexact\] and \[fig:mcevolutionexpanded\] an analogous analysis has been carried out for the $\overline{\rm MS}$ charm quark mass. In Fig. \[fig:mcevolutionexact\] $\overline m_c^{\rm N^3LL}(\mu)/\overline m_c^{\rm N^kLL}(\mu)$ is plotted for $k=0,1,2,3$ using the exact numerical solutions of the renormalization group equations and $\overline m_c(\mu=3~\mbox{GeV})$ as the respective reference value.[^6] Compared to the Fig. \[fig:alphasevolution\] we observe a very similar convergence. In Fig. \[fig:mcevolutionexpanded\], finally, we show $\overline m_c^{\rm N^3LL}(\mu)/\overline m_c^{(m)}(\mu)$, where $\overline m_c^{(m)}(\mu)$ is the ${\cal O}(\alpha_s^m)$ fixed-order expression for $\overline m_c(\mu)$ using $\alpha_s(3~\mbox{GeV})=0.2535$ as the expansion parameter. Again, the convergence towards the exact N${}^3$LL evolved result is very similar to the corresponding results for the strong coupling, and we find again no evidence for perturbative instabilities. Of course the corrections are somewhat larger when the fixed-order expansion is employed. We therefore conclude that perturbative evolution and renormalization scale variations for the $\overline{\mbox{MS}}$ charm quark mass can be safely used down to scales above $\overline m_c\sim 1.3$ GeV. One should of course mention that the lines in Fig. \[fig:mcevolutionexpanded\] also give an indication about the expected size of scale variations depending on the range of the variations. Scales above $2$ GeV can lead to sub-MeV variations, while scales down to the charm mass will result at best in percent precision (i.e. ${\cal O}(10~\mbox{MeV})$).

Perturbative Uncertainties in the $\overline{\rm MS}$ Charm Mass {#subsectionmcerror}
----------------------------------------------------------------

In this section we discuss in detail the perturbative series for the determination of the $\overline{\rm MS}$ charm mass $\overline m_c$ and how to set up an adequate scale variation to estimate the perturbative uncertainty. In the previous subsections we have discussed four different ways to carry out the perturbative expansion and we presented the corresponding order-by-order analytic expressions. As described there, we can determine at each order of the perturbative expansion for the moments $M_n^{\rm th}$ a value for $\overline
m_c(\mu_m)$ which also has a residual dependence on $\mu_\alpha$, the renormalization scale used for $\alpha_s$. To compare the different mass determinations we then evolve $\overline m_c(\mu_m)$ to obtain $\overline
m_c(\overline m_c)$ using the 4-loop renormalization group equations for the mass and the strong coupling [@Tarasov:1980au; @Larin:1993tp; @vanRitbergen:1997va; @Chetyrkin:1997dh; @Vermaseren:1997fq].[^7] The obtained value of $\overline m_c(\overline m_c)$ thus has a residual dependence on the scales $\mu_m$ and $\mu_\alpha$, on the order of perturbation theory and on the expansion method.[^8] For the results we can therefore use the notation $$\begin{aligned}
\label{mcnotation1} &&
\overline m_c(\overline m_c)[\mu_m,\mu_\alpha]^{i,n}
\,,\end{aligned}$$ where $n=0,1,2,3$ indicates perturbation theory at ${\cal O}(\alpha_s^n)$ and $$\begin{aligned}
\label{methodnotation1}
i & = & \left\{
\begin{array}[c]{ll}
a \qquad & \mbox{(fixed-order expansion),} \\
b & \mbox{(linearized expansion),} \\
c & \mbox{(iterative expansion),} \\
d & \mbox{(contour improved expansion).}
\end{array}\right.\end{aligned}$$

To initiate the discussion, we show in Fig. \[fig:trumpet1\] results for $\overline m_c(\overline m_c)$ at ${\cal O}(\alpha_s^n)$ for expansions a–d using $\mu_m=\mu_\alpha$ (upper four graphs) and using $\mu_m=\overline
m_c(\overline m_c)$ (lower four graphs). For each method and order we have displayed the range of $\overline m_c(\overline m_c)$ values for a variation of $2~\mbox{GeV}\le \mu_\alpha \le 4~\mbox{GeV}$, which corresponds to the scale variation employed in Refs. [@Chetyrkin:2006xg; @Boughezal:2006px; @Chetyrkin:2009fv; @Kuhn:2007vp]. Their analysis used the fixed-order expansion with the setting $\mu_m=\mu_\alpha$ and is represented by graph 1. We make several observations:

-   Choosing $\mu_m$ and $\mu_\alpha$ both larger than $2$ GeV makes the $\overline m_c(\overline m_c)$ value decrease with the order of perturbation theory.

-   Choosing $\mu_m$ smaller than $1.5$ GeV and $\mu_\alpha$ larger than $2$ GeV makes the $\overline m_c(\overline m_c)$ value increase with the order of perturbation theory.

-   For most choices of the scale setting and the expansion method the spread of the $\overline m_c(\overline m_c)$ values from the variation of $\mu_\alpha$ does not decrease in any substantial way with the order. However, viewing all methods and scale setting choices collectively a very good convergence is observed.

We have checked that these statements apply also in general beyond the specific cases displayed in Fig. \[fig:trumpet1\].

Quite conspicuous results are obtained for the scale choice $\mu_m=\mu_\alpha$ for the fixed-order (graph 1) and linearized expansions (graph 2). Here, extremely small variations in $\overline m_c(\overline m_c)$ are obtained. They amount to $1.8$ MeV ($4$ MeV) and $0.6$ MeV ($1.4$ MeV) at order $\alpha_s^2$ and $\alpha_s^3$, respectively, for the fixed-order expansions (linearized expansions). We note that our scale variation for the fixed-order expansion at ${\cal O}(\alpha_s^3)$ is consistent with the corresponding numbers quoted in Ref. [@Chetyrkin:2006xg; @Boughezal:2006px], where the ${\cal O}(\alpha_s^3)$ corrections to the first moment were computed,[^9] but differs from the scale variations given in Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp] which also quoted numerical results from the fixed-order expansion. Interestingly the ${\cal O}(\alpha_s^2)$ and ${\cal O}(\alpha_s^3)$ variations we find do not overlap and the ${\cal O}(\alpha_s^3)$ ranges appear to be highly inconsistent with the ${\cal O}(\alpha_s^3)$ results from the iterative method (graph 3). A visual display of scale variations obtained from the four expansion methods at ${\cal O}(\alpha_s^3)$ with different types of variation methods is given in Fig. \[fig:trumpet2\].

An illustrative way to demonstrate how a small scale variation can arise is given in Fig. \[fig:mccontour1fixed\] - \[fig:mccontour1improved\]. For all four expansion methods contour curves of constant $\overline m_c(\overline m_c)$ are displayed as a function of (the residual dependence on) $\mu_m$ and $\mu_\alpha$. For the fixed-order (a) and the linearized expansions (b) we see that there are contour lines closely along the diagonal $\mu_m=\mu_\alpha$. For the fixed-order expansion (a) this feature is almost exact and thus explains the extremely small scale-dependence seen in graph 1 of Fig. \[fig:trumpet1\]. For the linearized expansion (b) this feature is somewhat less exact and reflected in the slightly larger scale variation seen in graph 2 of Fig. \[fig:trumpet1\]. On the other hand, the contour lines of the iterative (c) and contour improved (d) expansions have a large angle with respect to the diagonal $\mu_m=\mu_\alpha$ leading to the much larger scale variations of $9$ and $5$ MeV, respectively, at ${\cal O}(\alpha_s^3)$ visible in graphs 3 and 4 of Fig. \[fig:trumpet1\]. The contour plots shown in Fig. \[fig:mccontour1\] also show the variation along the line $\mu_m$ or $\mu_\alpha=\overline
m_c(\overline m_c)$ (border of gray shaded areas). Here, the contour lines are relatively dense leading to scale variation of around $15$ MeV at ${\cal O}(\alpha_s^3)$.

Overall, we draw the following conclusions:

-   The small scale variations observed for the fixed-order and linearized expansions for $\mu_m=\mu_\alpha$ result from strong cancellations of the individual $\mu_m$ and $\mu_\alpha$ dependences that arise for this correlation.

-   Other correlations between $\mu_m$ and $\mu_\alpha$ that do not generate large logarithms do not lead to such cancellations. One therefore has to consider the small scale variations observed for $\mu_m=\mu_\alpha$ in the fixed-order and linearized expansions as accidental.

-   For an adequate estimate of perturbative uncertainties specific correlations between $\mu_m$ and $\mu_\alpha$ that are along contour lines of constant $\overline m_c(\overline m_c)$ have to be avoided. Moreover, adequate independent variations of $\mu_m$ and $\mu_\alpha$ should not induce large logarithms.

As the outcome of this discussion we adopt for our charm mass analysis an independent and uncorrelated variation of $\mu_m$ and $\mu_\alpha$ in the range $$\begin{aligned}
\label{muammvariation}
\overline m_c(\overline m_c) & \le & \mu_m,\mu_\alpha \, \le \, 4~\mbox{GeV}\,,\end{aligned}$$ in order to estimate perturbative uncertainties. The excluded region $\mu_m,\mu_\alpha<\overline{m}_c(\overline{m}_c)$ in the plane in the contour plots of Fig. \[fig:mccontour1\] is indicated by the gray shaded areas. This two-dimensional variation avoids accidental cancellations from correlated variations, large logarithms involving ratios of the scales $\overline m_c(\mu_m)$, $\mu_m$, $\mu_\alpha$ and remains well in the validity ranges for the perturbative renormalization group evolution of the $\overline{\rm MS}$ charm quark mass and $\alpha_s$ (see Sec. \[subsectionrunning\]). The range of Eq. (\[muammvariation\]) is also consistent with a scale variation $\mu\sim 2\,\overline m_c\times (1/2,1,2)$ one might consider as the standard choice with respect to the to particle threshold located at $\sqrt{s}=2\,{\overline m}_c$. As a comparison the range $\mu=(3\pm 1)$ GeV corresponds to $\mu=2\,{\overline m}_c(0.8,1.2,1.6)$. We emphasize that the lower boundary $\overline m_c(\overline m_c)\sim
1.25$ GeV is also reasonable as it represents the common flavor matching scale where gauge coupling evolution remains smooth up to NLL order. We do not see any evidence for perturbation theory for the first moment $M_1$ being unstable at the charm mass scale, in contrast to claims made in Ref. [@Kuhn:2007vp].

In Fig. \[fig:mmcerror1\] we show the ranges of $\overline m_c(\overline m_c)$ at ${\cal O}(\alpha_s^{1,2,3})$ for the four expansion methods employing the scale variations of Eq. (\[muammvariation\]). We see that all four expansion methods now lead to equivalent results. The variations are also compatible with the overall variations shown in Fig. \[fig:trumpet1\]. At ${\cal
  O}(\alpha_s^3)$ we obtain a scale variation for $\overline m_c(\overline m_c)$ of around $20$ MeV. This is an order of magnitude larger than the perturbative uncertainties quoted in Refs. [@Chetyrkin:2006xg; @Boughezal:2006px; @Chetyrkin:2009fv; @Kuhn:2007vp]. In Fig. \[fig:m3cerror1\] we have displayed the corresponding results for $\overline m_c(3~\mbox{GeV})$. At ${\cal
  O}(\alpha_s^3)$ they also exhibit a scale variation of around $20$ MeV.

Experimental Data {#sectiondata}
=================

Data Collections {#subsectioncollections}
----------------

[**Narrow resonances**]{}\
Below the open charm threshold there are the $J/\psi$ and $\psi^\prime$ narrow charmonium resonances. Their masses, and electronic widths are taken from the PDG [@Nakamura:2010zzi] and are collected in Tab. \[tabpsidata\] together with the value of the $\overline{\rm MS}$ QED coupling at their masses. The total widths are not relevant since we use the narrow width approximation for their contributions to the moments. The uncertainty for the contribution to the moments coming from the masses can be neglected.

                                 $J/\psi$      $\psi^\prime$
  -------------------------- ---------------- ----------------
          $M$ (GeV)           $3.096916(11)$   $3.686093(34)$
     $\Gamma_{ee}$ (keV)        $5.55(14)$       $2.48(6)$
   $(\alpha/\alpha(M))^{2}$     $0.957785$       $0.95554$

  : Masses and electronic widths [@Nakamura:2010zzi] of the narrow charmonium resonances and effective electromagnetic coupling [@Kuhn:2007vp]. $\alpha=1/137.035999084(51)$ stands for the fine structure constant, and $\alpha(M)$ stands for the $\overline{\mbox{MS}}$ electromagnetic coupling at the scale $M$.\[tabpsidata\]

[**Threshold and data continuum region**]{}\
The open charm threshold is located at $\sqrt{s}=3.73$ GeV. We call the energies from just below the threshold and up to $5$ GeV the threshold region, and the region between $5$ GeV and $10.538$ GeV, where the production rate is dominated by multiparticle final states the data continuum region. In these regions quite a variety of measurements of the total hadronic cross section exist from BES [@Bai:1999pk; @Bai:2001ct; @Ablikim:2004ck; @Ablikim:2006aj; @Ablikim:2006mb; @:2009jsa], CrystalBall [@Osterheld:1986hw; @Edwards:1990pc], CLEO [@Ammar:1997sk; @Besson:1984bd; @:2007qwa; @CroninHennessy:2008yi], MD1 [@Blinov:1993fw], PLUTO [@Criegee:1981qx], and MARKI and II [@Siegrist:1976br; @Rapidis:1977cv; @Abrams:1979cx; @Siegrist:1981zp]. Taken together, the entire energy region up to $10.538$ GeV is densely covered with total cross section measurements from these 19 data sets.[^10] The measurements from BES and CLEO have the smallest uncertainties. They do, however, not cover the region between $5$ and $7$ GeV. Here CrystalBall and MARKI and II have contributed measurements albeit with somewhat larger uncertainties. The statistical and total systematical uncertainties of the measurements can be extracted from the respective publications. For some data sets the amount of uncorrelated and correlated systematical uncertainties is given separately (BES [@Bai:2001ct; @Ablikim:2006aj; @Ablikim:2006mb], CrystalBall [@Osterheld:1986hw; @Edwards:1990pc], CLEO [@:2007qwa], MARKI and II [@Siegrist:1981zp; @Abrams:1979cx], MD1 [@Blinov:1993fw]) while for all the other data sets only combined systematical uncertainties are quoted. All these data sets are shown in Figs. \[figdatacompilation\], where the displayed error bars represent the (quadratically) combined statistical and systematical uncertainties.

Interestingly, none of the previous charm mass analyses, to the best of our knowledge, ever used the complete set of available data. As examples, Bodenstein et al. [@Bodenstein:2010qx] used data sets [@Bai:1999pk; @Bai:2001ct; @Ablikim:2006mb; @CroninHennessy:2008yi] from BES and CLEO. Jamin and Hoang [@Hoang:2004xm] used the data sets of Refs. [@Bai:2001ct], [@Blinov:1993fw] and [@Ammar:1997sk] from BES, MD1 and CLEO, covering the regions $2~\mbox{GeV} \leq E \leq 4.8~\mbox{GeV}$ and $6.964~\mbox{GeV} \leq E \leq 10.538~\mbox{GeV}$. Boughezal et al. [@Boughezal:2006px], Kuhn et al. [@Kuhn:2001dm], and Narison [@Narison:2010cg] use only one data set from BES [@Bai:2001ct]. Kuhn et al. [@Chetyrkin:2009fv; @Kuhn:2007vp] used the data sets of Refs. [@Bai:2001ct; @Ablikim:2006mb] from BES covering the energy region $2.6~\mbox{GeV} \leq E \leq 4.8~\mbox{GeV}$.

We consider three different selections of data sets to study the dependence of the experimental moments on this choice:

-   The [*minimal selection*]{} contains all data sets necessary to cover the whole energy region between $2$ and $10.538$ GeV without any gaps and keeping only the most accurate ones. These 8 data sets are from BES [@Bai:1999pk; @Bai:2001ct; @Ablikim:2006aj; @:2009jsa], CrystalBall [@Edwards:1990pc], CLEO [@CroninHennessy:2008yi; @:2007qwa] and MD1 [@Blinov:1993fw] corresponding to the data sets 1, 2, 5, 6, 9, 12, 13, and 14 (see Tab. \[tab:datasets\] for references).

-   The [*default selection*]{} contains all data sets except for the three ones with the largest uncertainties. It contains 16 data sets and fully includes the minimal selection. It contains all data sets except for Mark I and II data sets 16, 17 and 19 from Refs. [@Siegrist:1981zp; @Rapidis:1977cv; @Siegrist:1976br].

-   The [*maximal selection*]{} contains all 19 data sets.

We use the default selection as our standard choice for the charm mass analysis, but we will also quote results for the other data selections.

[**Perturbative QCD region**]{}\
Above $10.538$ GeV there are no experimental measurements of the total hadronic that might be useful for the experimental moments. In this energy region we will therefore use perturbative QCD to provide estimates for the charm production . As a penalty for not using experimental data we assign a $10$% total relative uncertainty to the contribution of the experimental moments coming from this region, which we then treat like an uncorrelated experimental uncertainty for the combination with the moment contribution from lower energies. As we see in Sec. \[subsectionmoments\] the energy region above $10.538$ GeV contributes only $(6, 0.4, 0.03, 0.002)\%$ to $M^{\rm
  exp}_{1,2,3,4}$. In the first moment $M_1^{\rm exp}$ the total contribution is about three times larger than the combined statistical and systematical (true) experimental uncertainties from the other energy regions. So the 10% penalty we assign to this approach represents a subleading component of the final quoted uncertainty.[^11] In the second and higher moments $M_{n\ge 2}^{\rm exp}$ the contributions from above $10.538$ GeV and the corresponding uncertainty are negligible compared to the uncertainties from the lower energy regions.

As the theoretical formula to determine the moment contribution from the perturbative QCD region we use the ${\cal O}(\alpha_s^3)$ non-singlet massless quark cross section including charm mass corrections up to ${\cal O}(\overline m_c^4/s^2)$: $$\begin{aligned}
R_{cc}^{\rm th}(s) & = &
N_c \,Q_c^2\,R^{\rm ns}(s,\overline m_c^2(\sqrt{s}), n_f=4,\alpha_s^{n_f=4}(\sqrt{s}))\,,
\label{Rcchighdef1}\end{aligned}$$ where $$\begin{aligned}
\label{Rnsdef1}
\lefteqn{R^{\rm ns}(s,\overline m_c^2(\mu), n_f=4,\alpha_s^{n_f=4}(\mu),\mu)}\\
& = &1+\dfrac{\alpha_{s}}{\pi}
+\left(\dfrac{\alpha_{s}}{\pi}\right)^{2}\left[
1.52453-2.08333 L_s\right]
+\left(\dfrac{\alpha_{s}}{\pi}\right)^{3}
\left[-11.52034+34028 L_s^2 -9.56052 L_s\right]\nonumber\\
 & + & \dfrac{\overline m_c^{2}(\mu)}{s}
\Big\{ 12\,\dfrac{\alpha_{s}}{\pi}+
\left(\dfrac{\alpha_{s}}{\pi}\right)^{2}\left[109.167-49 L_s\right]
 +\left(\dfrac{\alpha_{s}}{\pi}\right)^{3}\left[634.957-799.361 L_s+
 151.083 L_s^2\right]\!\!\Big\} \nonumber\\
 & + & \dfrac{\overline m_c^4(\mu)}{s^2}\Big\{ \!\!
-6+\dfrac{\alpha_{s}}{\pi}\left[-22+24 L_s\right]
  + \left(\dfrac{\alpha_{s}}{\pi}\right)^{2}
\left[140.855+221.5 L_s-5.16667 L_m-73 L_s^2\right]\nonumber\\
  && \qquad +  \left(\dfrac{\alpha_{s}}{\pi}\right)^{3}\left[
3776.94-509.226 L_s -174.831 L_m-1118.29 L_s^2\right.\nonumber\\
  && \qquad\qquad + 10.3333 L_m^2+\left.42.1944 L_m L_s 
+198.722 L_s^3\right]\Big\}
\,,\nonumber\end{aligned}$$ with $$L_s \, \equiv \, \ln\Big(\frac{s}{\mu^{2}}\Big)\,,\qquad
L_m \, \equiv \, \ln\Big(\dfrac{\overline m_c^{2}(\mu)}{s}\Big)\,.$$ For the computation of the contribution to the experimental moments we determine $\overline m_c(\sqrt{s})$ and $\alpha_s(\sqrt{s})$ appearing in Eq. (\[Rcchighdef1\]) using $\overline m_c(\overline m_c)=1.3$ GeV and $\alpha_s(m_Z)=0.118$ as initial conditions.

It is instructive to examine for the moment contributions from $\sqrt{s}>10.538$ GeV terms related to charm production that we do not account for in Eq. (\[Rnsdef1\]). In Tab. \[tabneglected\] the relative size with respect to the full first four moments (in percent) of the most important neglected contributions are given. In the second column the size of the mass corrections up to order $\overline m_c^4$, which we have included in $R^{\rm th}_{cc}$, are shown as a reference. The third column shows the contributions coming from secondary $c\bar c$ radiation through gluon splitting. The fourth column depicts the contributions from the ${\cal O}(\alpha_s^3)$ singlet corrections (including the mass corrections up to order $\overline m_c^4$), which one can take as an rough estimate for the actual contributions from the charm cut. Finally in the last column we show the size of the Z-boson exchange terms integrated from threshold to $10.538$ GeV. This contribution represents the Z-exchange contribution that is contained in the data, but - by definition - not accounted for in the theory moments. We see that at least for the first two moments, the contributions neglected are much smaller than the charm mass corrections we have accounted for in the nonsinglet production rate, which are already constituting a very small effect. Overall the numerical effect on the charm mass of all these contributions is tiny considering the scaling $\overline m_c\sim M_n^{1/2n}$. Since we assign a $10$% error on the moments’ contribution from the energy region $\sqrt{s}>10.538$ GeV where we use theory input, our approach to neglect subleading effects is justified.

[**Non-charm background**]{}\
Experimentally only the total hadronic cross section is available. Although charm-tagged rate measurements are in principle possible [@CroninHennessy:2008yi] they have not been provided in publications. On the other hand, they would also exhibit sizable additional uncertainties related to the dependence on simulations of the decay of charmed mesons into light quark final states. So to obtain the charm production cross section from the data we have to subtract the non-charm background using a model based on perturbative QCD related to the production of $u$, $d$ and $s$ quarks. A subtle point is related to the secondary radiation of $c\bar c$ pairs off the $u$, $d$ and $s$ quarks from gluon splitting and to which extent one has to account theoretically for the interplay between real and virtual secondary $c\bar c$ radiation which involves infrared sensitive terms [@Hoang:1997ca]. Since in this work we define the moments from primary $c\bar c$ production (see Eq. (\[momentdef1\])), secondary $c\bar c$ production is formally counted as non-charm background. Thus for the model for the non-charm background for $\sqrt{s}>2\, \overline m_c$ we employ the expression $$\begin{aligned}
\label{Rudsdef}
R_{uds}(s) & = &
N_c(Q_u^2 + 2 Q_d^2)\,R^{\rm ns}(s,0,n_f=3,\alpha_s^{n_f=4}(s))\\
 &&+\,
N_c\, Q_c^2\left(\dfrac{\alpha_s^{n_f=4}(s)}{\pi}\right)^{2}\left(\rho^{V}+\rho^{R}+\dfrac{1}{4}\,\log\dfrac{\overline m_{c}^{2}(\overline m_{c})}{s}\right)
\,.\nonumber\end{aligned}$$ The second term on the RHS describes the contributions from real and virtual secondary $c\bar c$ radiation. The analytic expressions for $\rho^R$ and $\rho^V$ can be found in Eqs. (2) and (6) of Ref. [@Hoang:1994it]. We have checked that the numerical impact of real ($\rho^R$) and virtual ($\rho^V$) secondary radiation individually as well as the complete second term on the RHS of Eq. (\[Rudsdef\]) on the moments is negligible, see Tab. \[tabneglected\]. We use Eq. (\[Rudsdef\]) and fit the non-charm background including also data in the region $2~\mbox{GeV}\le E\le 3.73~\mbox{GeV}$ via the ansatz $R_{{\rm non}-c\bar c}(s) = n_{\rm ns}\,R_{uds}(s)$, where the constant $n_{\rm
  ns}$ represents an additional fit parameter.

Data Combination {#subsectioncombination}
----------------

Combining different experimental measurements of the hadronic cross section one has to face several issues: (a) the measurements are given at individual separated energy points, (b) the set of measurements from different publications are not equally spaced, cover different, partly overlapping energy regions and have different statistical and systematical uncertainties, (c) the correlations of systematical errors are only known (or provided) for the data sets within each publication, (d) there are a number of very precise measurements at widely separated energies.

In this section we discuss the combination of the experimental data from the threshold and the data continuum regions between $2$ and $10.538$ GeV using a method based on a fitting procedure used before for determining the hadronic vacuum polarization effects for $g-2$ [@Hagiwara:2003da]. In this work we extend this approach and also account for the subtraction of the non-charm background.

[**Combination method**]{}\
The method uses the combination of data in energy bins (clusters) assuming that the within each cluster changes only very little and can thus be well approximated by a constant. Thus clusters for energies where $R$ varies rapidly need to be small (in this case the experimental measurements are also denser). The in each cluster is then obtained by a $\chi^2$ fitting procedure. Since each experimental data set from any publication covers an energy range overlapping with at least one other data set, the clusters are chosen such that clusters in overlapping regions contain measurements from different data sets. Through the fitting procedure correlations are then being communicated among different data sets and very accurate individual measurements can inherit their precision into neighbouring clusters. Both issues are desirable since the hadronic is a smooth function with respect to the sequence of clusters.

To describe the method we have to set up some notation:

-   All measurements $R(E)$ are distinguished according to the energy $E$ at which they have been carried out.

-   Each such energy point having a measurement is written as $E_i^{k,m}$, where $k=1,\ldots,N_{\rm exp}$ refers to the $N_{\rm exp}$ data sets, $m=1,\ldots,N_{\rm cluster}$ runs over the $N_{\rm cluster}$ clusters and $i=1,\ldots,N^{k,m}$ assigns the of the $N^{k,m}$ measurements.

-   Each individual measurement of the is then written as $$\begin{aligned}
    \label{Rexp}
    R(E^{k,m}_i) & = & R_i^{k,m} \, \pm \, \sigma_i^{k,m} \, 
    \pm \, \Delta_i^{k,m}
    \,,\end{aligned}$$ where $R_i^{k,m}$ is the central value, $\sigma_i^{k,m}$ the combined statistical and uncorrelated systematical uncertainty and $\Delta_i^{k,m}$ the correlated systematical experimental uncertainty.

-   For convenience we define $\Delta
      f_i^{k,m}=\Delta_i^{k,m}/R_i^{k,m}$ to be the relative systematical correlated uncertainty.

As our standard choice concerning the clusters we use 5 different regions each having equidistant cluster sizes $\Delta E$. The regions are as follows:

-   [**non-charm region**]{}: has 1 cluster for $2~\mbox{GeV}\leq E \leq
      3.73~\mbox{GeV}$ ($\Delta E=1.73$ GeV).

-   [**low charm region**]{}: has 2 clusters for $3.73~\mbox{GeV}< E \leq
      3.75~\mbox{GeV}$ ($\Delta E=10$ MeV).

-   [**$\psi(3S)$ region/threshold region 1**]{}: has 20 cluster for $3.75~\mbox{GeV}< E \leq
      3.79~\mbox{GeV}$ ($\Delta E=2$ MeV).

-   [**resonance region 2**]{}: has 20 cluster for $3.79~\mbox{GeV}< E \leq
      4.55~\mbox{GeV}$ ($\Delta E=38$ MeV).

-   [**continuum region**]{}: has 10 cluster for $4.55~\mbox{GeV}< E \leq
      10.538~\mbox{GeV}$ ($\Delta E=598.8$ MeV).

We assign to this choice of $52+1$ clusters the notation (2,20,20,10) and later also examine alternative cluster choices demonstrating that the outcome for the moments does within errors not depend on them. The cluster in the non-charm region is used to fit for the normalization constant $n_{\rm ns}$ of the non-charm background contribution, see Eq. (\[Rudsdef\]).

Our standard procedure to determine the central energy $E_m$ associated to each cluster is just the weighted average of the energies of all measurements falling into cluster $m$, $$\begin{aligned}
\label{Eclusterdef}
E_m & = & 
\dfrac{\sum_{k,i}\dfrac{E_{i}^{k,m}}
{(\sigma_{i}^{k,m})^2+(\Delta_{i}^{k,m})^2}}
{\sum_{k,i}\dfrac{1}{(\sigma_{i}^{k,m})^2+(\Delta_{i}^{k,m})^2}}
\,.\end{aligned}$$ The corresponding for the charm cross section that we determine through the fit procedure described below is called $$\begin{aligned}
\label{Rclusterdef}
R_m & \equiv & R_{c\bar c}(E_m)
\,. \end{aligned}$$ We note that using instead the unweighted average or simply the center of the cluster has a negligible effect on the outcome for the moments since the clusters we are employing are sufficiently small.

[**Fit procedure and** ]{}\
We determine the charm cross section $R_{c\bar c}$ from a that accounts for the positive correlation among the systematical uncertainties $\Delta_i^{k,m}$ within each experiment $k$ and, at the same time, also for the non-charm background. To implement the correlation we introduce the auxiliary parameters $d_k$ ($k=1,\ldots, N_{\rm exp}$) that parametrize the correlated deviation from the experimental central values $R_i^{k,m}$ in units of the correlated systematical uncertainty $\Delta_i^{k,m}$, see Eq. (\[Rexp\]). In this way we carry out fits to $R_i^{k,m}+d_k \,\Delta_i^{k,m}$ and treat the $d_k$ as additional auxiliary fit parameters that are constraint to be of order one (one standard deviation) by adding the term $$\begin{aligned}
\label{chi2totald}
\chi^2_{\rm corr}(\{d_k\}) & = &
\sum_{k=1}^{N_{\rm exp}} \, d_k^2\,,\end{aligned}$$ to the . To implement the non-charm background we assume that the relative energy dependence of the non-charm cross section related to primary production of $u$, $d$ and $s$ quarks is described properly by $R_{uds}$ given in Eq. (\[Rudsdef\]). We then parametrize the non-charm background cross section by the relation $$\begin{aligned}
\label{Rncdef}
R_{{\rm non}-c\bar c}(E) &= & n_{\rm ns}\,R_{uds}(E)
\,,\end{aligned}$$ where the fit parameter $n_{\rm ns}$ is determined mainly from the experimental data in the first clusters below $3.73$ GeV by adding the term $$\begin{aligned}
\label{chi2noncharm}
\chi^2_{\rm nc}(n_{\rm ns},\{d_k\}) & = &
\sum_{k=1}^{N_{\mathrm{exp}}}
\sum_{i=1}^{N^{k,1}}\left(\dfrac{
R_i^{k,1}-(1+\Delta f_i^{k,1}\, d_k)\, n_{ns}\, 
R_{uds}(E_i^{k,1})}{\sigma_i^{k,1}}\right)^2\,,\end{aligned}$$ to the . The complete then has the form $$\begin{aligned}
\label{chi2total}
\chi^2(\{R_m\}, n_{\rm ns},\{d_k\}) & = &
 \chi^2_{\rm corr}(\{d_k\}) \, + \,
\chi^2_{\rm nc}(n_{\rm ns},\{d_k\})\,+\,\chi^2_{\rm c}(\{R_m\}, n_{\rm ns},\{d_k\})
\,,\end{aligned}$$ where[^12] $$\begin{aligned}
\label{chi2charm}
&&\chi^2_{\rm c}(\{R_m\}, n_{\rm ns},\{d_k\})  = \\[2mm]
&&\sum_{k=1}^{N_{\mathrm{exp}}}
\sum_{m=2}^{N_{\mathrm{clusters}}}\sum_{i=1}^{N^{k,m}}
\left(\dfrac{R_{i}^{k,m}-(1+\Delta f_i^{k,m}\, d_k)\,
(R_m+n_{ns}\, R_{uds}(E_i^{k,m}))}{\sigma_i^{k,m}}\right)^{2}
\,.\nonumber\end{aligned}$$ Note that in our approach the non-charm normalization constant $n_{\rm ns}$ is obtained from a combined fit together with the cluster values $R_m$.

This form of the is an extended and adapted version of the ones used in Refs. [@D'Agostini:1993uj; @Takeuchi:1995xe]. A special characteristic of the in Eqs. (\[chi2noncharm\]) and (\[chi2charm\]) is that the relative correlated experimental uncertainties $\Delta f_i^{k,m}$ enter the by multiplying the fit value $R_m$ rather than the experimental values $R_i^{k,m}$. This leads to a non-bilinear dependence of the on the $d_m$ and the $R_m$ fit parameters and avoids spurious solutions where the best fit values for the $R_m$ are located systematically below the measurements. Such spurious solutions can arise for data points with substantial positive correlation when with strictly bilinear dependences are employed  [@D'Agostini:1993uj; @Takeuchi:1995xe].[^13]

We also note that the implementation of the non-charm background subtraction given in Eq. (\[chi2charm\]) leads to a partial cancellation of systematical uncertainties for the $R_m$ best fit values for the charm cross section. Moreover, it is interesting to mention that in the limit where each cluster contains exactly one measurement (except below threshold, in which we always keep one cluster) the decouples, after performing the change of variables $R'_m=R_{m}+n_{ns}\, R_{uds}(E^{k,m})$, into the sum of two independent , one containing data below threshold and depending only on $n_s$ and $d_k$, and another one containing data above threshold $R^{k,m}$ (we drop the label $i$ because having only one data per cluster it can take only the value $1$) and depending only on $R'_m$. After minimizing the first one can obtain the best fit values for $n_s$ and $d_k$, denoted by $n_s^{(0)}$ and $d_k^{(0)}$, respectively. The second $\chi^2$ has a minimal value of $0$ and the best fit parameters read $R_m^{(0)}=R^{k,m}/(1+d_k^{(0)}\Delta^{k,m})-n_s^{(0)}R_{uds}(E^{k,m})$.

Close to the minimum the of Eq. (\[chi2total\]) can be written in the Gaussian approximation $$\begin{aligned}
\label{chi2Gauss}
\chi^2(\{p_i\}) & = & \chi^2_{\rm min}\,+\,  \sum_{i,j} (p_i-p_i^{(0)})V_{i,j}^{-1}(p_j-p_j^{(0)})
\, + \, {\cal O}\Big((p-p^{(0)})^3\Big)
\,,\end{aligned}$$ where $p_i=(
\{R_m\}, n_{\rm ns}, \{d_k\})$ and the superscript $(0)$ indicates the respective best fit value. After determination of the correlation matrix $V_{ij}$ by numerically inverting $V^{-1}_{i,j}$ we can drop the dependence on the auxiliary variables $n_{\rm nc}$ and $d_k$ and obtain the correlation matrix of the $R_m$ from the of $V_{ij}$ which we call $V_{m m^\prime}^R$. In order to separate uncorrelated statistical and systematical uncertainties from correlated systematical ones we compute the complete $V_{m m^\prime}^R$ accounting for all uncertainties and a simpler version of the correlation matrix, $V_{m m^\prime}^{R,u}$ accounting only for uncorrelated uncertainties. The latter is obtained from dropping all correlated errors $\Delta_i^{k,m}$ from the  (\[chi2total\]).[^14]

The outcome of the fit for the sum of the charm and the non-charm cross section in the threshold and the data continuum region using the standard data set explained above is shown in Figs. \[fig:fitcompilation\](a)-(f) together with the input data sets. The red line segments connect the best fit values and the brown band represents the total uncertainty. The clusters are indicated by dashed vertical lines. For completeness we have also given all numerical results for the $R_m$ values in Appendix A. There we also give results for the minimal and maximal data set selections.

Experimental Moments {#subsectionmoments}
--------------------

[**Narrow resonances**]{}\
For the $J/\psi$ and $\psi^\prime$ charmonium contributions to the experimental moments we use the narrow width approximation, $$\begin{aligned}
\label{Mnresonances}
M_n^{\rm res} & = & \dfrac{9\,\pi\,\Gamma_{ee}}{\alpha(M)^{2}M^{2n+1}}\,,\end{aligned}$$ with the input numbers given in Tab. \[tabpsidata\]. We neglect the tiny uncertainties in the charmonium masses as their effects are negligible.

[**Threshold and data continuum region**]{}\
For the determination of the moment contributions from the threshold and the continuum region between $3.73$ and $10.538$ GeV we use the results for the clustered $c\bar c$ cross section values $R_m$ determined in Sec. \[subsectioncombination\] and the trapezoidal rule. Using the relations ${\rm d}s={\rm d}E^2=2E{\rm d}E$ we thus obtain $$\begin{aligned}
\label{Mndata}
M_n^{\rm thr+cont} & = &\sum_{i=1}^{N_{\mathrm{clusters}}}\dfrac{R_{i}}{E_{i}^{2n+1}}(E_{i+1}-E_{i-1})+\dfrac{R_{0}}{E_{0}^{2n+1}}(E_{1}-E_{0})\\
&+&\dfrac{R_{N_{\mathrm{clusters}}+1}}{E_{N_{\mathrm{clusters}}+1}^{2n+1}}(E_{N_{\mathrm{clusters}}+1}-E_{N_{\mathrm{clusters}}})\,,\nonumber\end{aligned}$$ where $R_0$ and $E_0$ are the $R$ and energy values at the lower boundary of the smallest energy cluster, and $R_{N_{cl}+1}$ and $E_{M_{cl}+1}$ are the corresponding values of the upper boundary of the highest energy cluster. The values for $R_0$ and $R_{N_{cl}+1}$ are obtained from linear extrapolation using the $R_m$ values of the two closest lying clusters[^15] $m^\prime$ and $m^\prime+1$ with the formula $$\begin{aligned}
\label{REextrapol}
R(E) & = & 
\frac{R_{m^\prime+1}-R_{m^\prime}}{E_{m^\prime+1}-E_{m^\prime}}\, (E-E_{m^\prime})+ R_{m^\prime}
\,.\end{aligned}$$ For the computation of the moment contributions from subintervals within the range between $3.73$ and $10.538$ GeV we also use Eq. (\[Mndata\]) using corresponding adaptations for the boundary values at $m=0$ and $m=N_{cl}+1$.

[**Perturbative QCD region**]{}\
For the region above $10.538$ GeV where we use the perturbative QCD input described in Eqs. (\[Rcchighdef1\]) and (\[Rnsdef1\]) for the charm the contribution to the experimental moment is obtained from the defining equation (\[momentdef2\]) with a lower integration limit of $10.538$ GeV: $$\begin{aligned}
\label{MnQCD}
M_n^{\rm QCD} & = & \gamma \times\int_{(10.538\,\rm GeV)^2}^\infty {\rm d}s \,\dfrac{R_{cc}^{\rm th}(s)}{s^{n+1}}
\,.\end{aligned}$$ The variable $\gamma$ is an auxiliary variable used to parametrize the $10\%$ uncertainty we assign to the perturbative QCD contribution, $$\begin{aligned}
\label{gammadef}
\gamma & = & 1.0 \,\pm \, 0.1
\,.\end{aligned}$$

[**Correlations**]{}\
The experimental moments are obtained from the sum of the resonance, threshold plus continuum and perturbative QCD contributions described just above, $$\begin{aligned}
\label{Mnsum}
M_n^{\rm exp} & = &
M_n^{\rm res} \, + \,
M_n^{\rm thr+cont} \, + \,
M_n^{\rm QCD}
\,.\end{aligned}$$ To determine the uncertainties we account for the errors in the $e^+e^-$ widths of $J/\psi$ and $\psi^\prime$ and in the cluster values $R_m$, and for the $10\%$ assigned uncertainty in $M_n^{\rm QCD}$. For the evaluation we use the usual error propagation based on a $\bar m\times\bar m$ correlation matrix with $\bar
m=N_{cl}+3$. The correlation matrix of the experimental moments thus has the form $$\begin{aligned}
\label{CMcorr}
C_{nn^\prime}^{\rm exp} & = & \sum\limits_{i,j=1}^{N_{cl}+3} 
\Big(\frac{\partial M_n^{\rm exp}}{\partial \bar R_i}\Big)\,
\Big(\frac{\partial M_{n^\prime}^{\rm exp}}{\partial \bar R_j}\Big) 
\, V_{ij}^{\bar R}
\,,\end{aligned}$$ where we have $\bar R_i =
(\{R_m\},\Gamma_{e^+e^-}(J/\psi),\Gamma_{e^+e^-}(\psi^\prime),\gamma)$. The entries of $V^{\bar R}$ in the $R_m$ subspace are just the entries of the correlation matrix $V^R$ obtained from the cluster fit. The diagonal entries in the $\Gamma_{e^+e^-}$ subspace are the combined statistical and systematical uncertainties of the $e^+e^-$ widths and the $\delta\gamma=0.1$ for $M_n^{\rm QCD}$, respectively. We treat the latter uncertainty as uncorrelated with all other uncertainties. So all non-diagonal entries of $V_{ij}^{\bar R}$ for $i$ or $j=N_{cl}+3$ are zero. For the uncertainty of the $e^+e^-$ widths we adopt a model where the (quadratic) half of the error is uncorrelated and the other (quadratic) half is positively correlated among the two narrow resonances and to the $R_n$ cluster values as well. Thus for the corresponding non-diagonal entries of $V_{i,j}^{\bar R}$ with $i\in\{1,N_{cl}\}$ and we have the entries $V_{ij}^{\bar R}=\delta\Gamma^{1,2}_{e^+e^-}\Delta R_m/\sqrt{2}$, and for $i=N_{cl}+1$ and $j=N_{cl}+2$ we have $V_{ij}^{\bar R}=\delta\Gamma^1_{e^+e^-}\delta\Gamma^2_{e^+e^-}/2$, where $\delta\Gamma^{1,2}_{e^+e^-}$ are the respective $e^+e^-$ width total uncertainties for $J/\psi$ and $\psi^\prime$, respectively, and $\Delta R_m$ the correlated uncertainty of the respective cluster $R_m$ value.

The results for the moments showing separately the contributions from the resonances, various energy subintervals and their total sum using the defaults data set collection (see Sec. \[subsectioncollections\]) are given in Tab. \[tab:momres1\]. Using Eq. (\[CMcorr\]) it is straightforward to compute the correlation matrix of the moments, and we obtain $$\begin{aligned}
\label{CMtotal}
C^{\rm exp} & = &\left(
\begin{array}{cccc}
0.250 & 0.167 & 0.147 & 0.142\\
0.167 & 0.120 & 0.107 & 0.103\\
0.147 & 0.107 & 0.095 & 0.092\\
0.142 & 0.103 & 0.092 & 0.090
\end{array}
\right)\,,\end{aligned}$$ for the total correlation matrix accounting for all correlated and uncorrelated uncertainties. We remind the reader that all numbers related to the moment $M_n^{\rm exp}$ are given in units of $10^{-(n+1)}\,\mbox{GeV}^{-2n}$. To quote correlated and uncorrelated uncertainties separately it is also useful to show the correlation matrix that is obtained when only uncorrelated uncertainties are accounted for. The result it $$\begin{aligned}
\label{CMuncorr}
C^{\rm exp}_{\rm uc} & = &
\left(
\begin{array}{cccc}
0.041 & 0.035 & 0.034 & 0.034\\
0.035 & 0.034 & 0.034 & 0.035\\
0.034 & 0.034 & 0.035 & 0.036\\
0.034 & 0.035 & 0.036 & 0.037
\end{array}\right)
\,.\end{aligned}$$ These results can be used to carry out combined simultaneous fits to several of the moments. This is, however, not attempted in this work.

Examination {#subsectionexaminations}
-----------

We conclude this section with an examination of some of the choices and assumptions we have implemented for the treatment of the experimental uncertainties. Our defaults choices include

-   treating one (quadratic) half of $J/\psi$ and $\psi^\prime$ $e^+e^-$ partial width uncertainties as uncorrelated and the other half as positively correlated among themselves and with the correlated uncertainties of the $R_m$ cluster values;

-   treating the entire systematical uncertainties of the measurements as correlated for the data sets where only total systematical uncertainties were quoted;

-   defining the cluster energies $E_m$ through the weighed average of measurement energies falling into the clusters, see Eq. (\[Eclusterdef\]);

-   using $N_{cl}=52+1$ clusters distributed in groups of (2,20,20,10) clusters in the energy ranges bounded by $(3.73, 3.75, 3.79, 4.55, 10.538)$ GeV (see Sec. \[subsectioncombination\]) and

-   using the default data set collection consisting of all data sets discussed in Sec. \[subsectioncollections\] except for sets 16, 17 and 19 as defined in Tab. \[tab:datasets\].

In Tab. \[tab:compcorr1\] alternative correlation models are being studied. The second column shows, as a reference, the first four moments with our default settings. The third and fourth columns show the moments treating the uncertainties of the $J/\psi$ and $\psi^\prime$ partial widths being minimally correlated with the $R_m$ values[^16] and completely uncorrelated, respectively. In the fifth column we display the moments treating, for data sets 1, 6, 13, 15, 16 and 17, one (quadratic) half of systematical uncertainties for the uncorrelated and the other half correlated. In the sixth column all $R_m$ values of all data sets are treated as completely uncorrelated. We see that the central values depend only weakly on the correlation model for those data were the corresponding information is unknown. In particular, for the determination of the uncertainties the ignorance about the composition of the systematical uncertainties in the from data sets 1, 6, 13, 15, 16 and 17 is not essential. However, for quoting the final uncertainties it is important to account for all (known) correlations since treating all errors uncorrelated underestimates the uncertainties considerably.

In Tab. \[tab:compcorr2\] we examine the impact of modifying the definition of the cluster energy $E_m$ and of changes to the default cluster distribution (2,20,20,10). In the second column we display the resulting moments of the default setting. In the third and fourth columns we have shown the moments using for $E_m$ simply the mean of the energies and the center of the cluster, respectively. The resulting differences to the default definition is an order of magnitude smaller than the uncertainties and thus negligible. The fifth, sixth and seventh columns display the moments using some alternative cluster distributions. The deviations for the default choice illustrated in the table are much smaller than the uncertainties and typical for all modifications that satisfy the guidelines for viable cluster definitions we have formulated in Sec. \[subsectioncombination\]. This demonstrates that the choice of the cluster distribution does not result in a bias for the resulting experimental moments.

Finally, we also examine the dependence of the moments on the data set collections as described in Sec. \[subsectioncollections\]. In Tabs. \[tab:moments-results-minimal\] and \[tab:moments-results-maximal\] the results for the moments are displayed using the minimal and the maximal collections (with default choices for all other settings). We see that the differences in the central values to the default collection are about half the size of the systematical correlated uncertainties for the first moment $M_1^{\rm exp}$. For the higher moments the differences are much smaller than the uncertainties. Using, instead of the default, the minimal and maximal collections affects the uncertainties of $M_1^{\rm exp}$ by only about $10\%$. For the higher moments the differences decrease strongly and basically disappear for the fourth moment. Again, the results show that having an increased or decreased redundancy in the data set collection only has a minor impact on the final numbers for the experimental moments.

To summarize, we find that modifications to the choices and assumptions that go into the combined treatment of the experimental data from different publications and experiments lead to changes that are well within the experimental uncertainties we obtain from our combination method. We therefore consider these uncertainties as appropriate. An instructive comparison of the moments obtained in our analysis to those obtained in some previous publications is given in Tab. \[tab:mom-comparison\]. A graphical illustration of the results is shown in Fig. \[fig:Moment-comparison\].

Charm Quark Mass Analysis {#sectionanalysis}
=========================

   $\alpha_s(m_Z)$   ${\overline m_c}({\overline m_c})$   $\Delta_{\rm pert}$   $\Delta_{\rm stat}$   $\Delta_{\rm syst}$   ${\overline m_c}(3\,{\rm GeV})$   $\Delta_{\rm pert}$   $\Delta_{\rm stat}$   $\Delta_{\rm syst}$
  ----------------- ------------------------------------ --------------------- --------------------- --------------------- --------------------------------- --------------------- --------------------- ---------------------
       $0.113$                    $1.255$                       $0.013$               $0.006$               $0.013$                     $1.010$                     $0.013$               $0.006$               $0.014$
       $0.114$                    $1.259$                       $0.014$               $0.006$               $0.013$                     $1.006$                     $0.015$               $0.006$               $0.014$
       $0.115$                    $1.263$                       $0.015$               $0.006$               $0.013$                     $1.002$                     $0.016$               $0.006$               $0.014$
       $0.116$                    $1.267$                       $0.016$               $0.006$               $0.013$                     $0.999$                     $0.017$               $0.006$               $0.014$
       $0.117$                    $1.271$                       $0.017$               $0.006$               $0.013$                     $0.995$                     $0.019$               $0.006$               $0.014$
       $0.118$                    $1.275$                       $0.019$               $0.006$               $0.013$                     $0.990$                     $0.020$               $0.006$               $0.014$
       $0.119$                    $1.280$                       $0.020$               $0.006$               $0.013$                     $0.986$                     $0.022$               $0.006$               $0.014$

  : Results for the central values of ${\overline m_c}({\overline m_c})$ and ${\overline m_c}(3\,{\rm GeV})$ (second and sixth columns) and their perturbative (third and seventh columns), experimental statistical (fourth and eight columns) and experimental systematical (fifth and ninth columns) errors.[]{data-label="tabmcalphas"}

Since it is theoretically most reliable, we use the first moment $M_1$ for our final numerical charm quark mass analysis. As ingredients for the analysis we use

-   the iterative expansion method for the perturbative contribution of the theoretical moment at ${\cal O}(\alpha_s^3)$, see Eq. (\[eq:iterative-general\]),

-   the gluon condensate correction with its Wilson coefficient determined at ${\cal O}(\alpha_s)$ as described in Sec. \[subsectioncondensate\],

-   the first experimental moment $$\begin{aligned}
    \label{M1exp}
    M_1^{\rm exp} & = & 0.2138\pm 0.0020_{\rm stat}\pm 0.0046_{\rm syst}\end{aligned}$$ using our default settings as discussed in Sec. \[sectiondata\].

One important source of uncertainty we have not yet discussed is the value of the strong $\overline{\mbox{MS}}$ coupling $\alpha_s$. Since in the recent literature [@Abbate:2010xh; @Frederix:2010ne; @Gehrmann:2009eh; @Blumlein:2006be; @Blumlein:2010rn] $\alpha_s$ determinations with a spread larger than the current world average [@Bethke:2009jm] have been obtained, we carry out our numerical analysis for values of $\alpha_s(m_Z)$ between $0.113$ and $0.119$.[^17] The outcome of our analysis is shown in Tab. \[tabmcalphas\]. In Fig. \[figmcalphas\] the central values, perturbative, statistical and systematical uncertainties are displayed graphically. For the central value and the perturbative uncertainty, which show a significant dependence on $\alpha_s$, we can present a linear fit. For the statistical and systematical uncertainties the variation with $\alpha_s$ is smaller than $1$ MeV and we only quote constant values. We thus obtain $$\begin{aligned}
\label{mcfinalalphas}
\overline m_c(\overline m_c)& = &(0.788 + 4.13\times\alpha_s(m_Z)) 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.013)_{\rm syst}\\ &&
\, \pm \, (-0.127+1.23\times\alpha_s(m_Z))_{\rm pert} 
\, \pm \, (0.002)_{\langle GG\rangle}
\,,\nonumber
\\[2mm]
\overline m_c(3~\mbox{GeV})& = &(1.462 - 4.00\times\alpha_s(m_Z)) 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.014)_{\rm syst}\\ &&
\, \pm \, (-0.148+1.43\times\alpha_s(m_Z))_{\rm pert} 
\, \pm \, (0.002)_{\langle GG\rangle}\,.\nonumber\end{aligned}$$ Taking as an input $$\begin{aligned}
\label{alphaswa}
\alpha_s(m_Z) & = & 0.1184  \pm 0.0021\,,\end{aligned}$$ which is the current world Bethke average [@Bethke:2009jm] with a tripled uncertainty we obtain $$\begin{aligned}
\label{mcfinalalphaswa}
\overline m_c(\overline m_c) & = & 1.277 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.013)_{\rm syst} 
\, \pm \, (0.019)_{\rm pert}
\, \pm \, (0.009)_{\alpha_s} 
\, \pm \, (0.002)_{\langle GG\rangle}
\,,\nonumber
\\[2mm]
\overline m_c(3~\mbox{GeV}) & = & 0.998 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.014)_{\rm syst} 
\, \pm \, (0.022)_{\rm pert}
\, \pm \, (0.010)_{\alpha_s} 
\, \pm \, (0.002)_{\langle GG\rangle}
\,,\nonumber\\\end{aligned}$$ which represents, together with Eq. (\[mcfinalalphas\]), our final numerical result for the $\overline{\mbox{MS}}$ charm mass. Our result is in good agreement with other recent determinations of $\overline{m}_{c}(\overline{m}_{c})$. A summary of the numerical results is shown in Tab. \[tab:comparison-alpha\] and in Fig. \[fig:comparison-alpha\]. Compared to the analysis carried out in Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp] our experimental uncertainty is larger by $5$ MeV and our perturbative uncertainty is a factor of $10$ larger. Compared to Refs. [@Chetyrkin:2006xg; @Boughezal:2006px] the discrepancy in the perturbative error estimate is even larger. We note that the almost head on agreement in the central value of Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp] is accidental, since the central value we obtain is based on a different perturbative expansion (iterative expansion in our work versus fixed-order expansion in theirs), on the way how the central value is defined (average of the lowest and highest values in our work versus value for $\mu=3$ GeV) and on a different central value for the experimental moment ($M_1^{\rm exp}=0.2138$ in our work versus $M_1^{\rm exp}=0.2166$ in theirs). Whereas our definition of the central value makes $\overline m_c$ smaller, our smaller value for $M_1^{\rm exp}$ makes it bigger by almost the same amount.

![ Comparison of recent determinations of $\overline{m}_{c}(\overline{m}_{c})$. Red corresponds to our result, black and gray correspond to ${\mathcal O}\alpha_s^3)$ and ${\mathcal O}(\alpha_s^2)$ charmonimum sum rules analyses, respectively, green labels other kind of sum rule analyses (weighted finite energy sum rules [@Bodenstein:2010qx] and ratios of moments [@Narison:2010cg]), and blue stands for lattice simulations. \[fig:comparison-alpha\] ](alpha-comparison.eps){width="90.00000%"}

Conclusions and Final Thoughts {#sectionconclusions}
==============================

In this work we have used state-of-the-art ${\cal O}(\alpha_s^3)$ input from perturbative QCD to determine the $\overline{\mbox{MS}}$ charm quark mass from relativistic ($n=1$) charmonium sum rules using experimental data on the total hadronic cross section in $e^+e^-$ annihilation. The main aims were (i) to carefully reexamine perturbative uncertainties in the charm mass extractions from the moments of the charm vector current correlator and (i) to fully exploit the available experimental data on the hadronic cross section.

We carried out this work having in mind recent ${\mathcal O}(\alpha_s^3)$ sum rule analyses [@Chetyrkin:2006xg; @Boughezal:2006px; @Chetyrkin:2009fv; @Kuhn:2007vp], where perturbative errors of $2$ MeV or smaller were quoted. Moreover in Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp] an experimental uncertainty of $9$ MeV was quoted. Given these numbers we found it appropriate to reexamine their analysis. In their work the perturbative uncertainty estimate was achieved using a specific choice to arrange the perturbative expansion and by setting the renormalization scales in $\alpha_s$ ($\mu_\alpha$) and of the $\overline{\mbox{MS}}$ charm quark mass ($\mu_m$) equal. We demonstrated that this results in an accidental cancellation of $\mu_\alpha$ and $\mu_m$ scale variations that is not observed in other alternative ways to treat the perturbative expansion. Moreover, concerning the experimental input their work relied on perturbative QCD predictions instead of available data for energies $\sqrt{s}>4.8$ GeV which results in a rather strong dependence of the final experimental uncertainty on the model-uncertainty assigned to the theory input.

Concerning a proper assessment of perturbative and experimental uncertainties we established in our work the following improvements:

-   We demonstrated that for achieving an estimate of perturbative uncertainties based on scale variations that is independent of the perturbative expansion method one needs to vary $\mu_\alpha$ and $\mu_m$ independently, albeit with ranges that avoid large logs. As a result the perturbative uncertainty estimates using different ways to carry out the expansion in $\alpha_s$ become equivalent, which is not the case for $\mu_\alpha=\mu_m$.

-   Using a data clustering method similar to Refs. [@D'Agostini:1993uj; @Takeuchi:1995xe; @Hagiwara:2003da] we combined available data on the total $e^+e^-$ hadronic cross section from many different experiments covering energies up to $\sqrt{s}=10.538$ GeV to fully exploit the existing experimental information for the experimental moments. This avoids a strong dependence of the experimental moments on ad-hoc assumptions on the “experimental” model uncertainty being associated to the QCD theory input used for the experimental moments. As a result we were also able to quantify the correlation between different experimental moments.

Using $\alpha_s(m_Z)$ as an unspecified variable and the theoretically more reliable first moment $M_1$ for the fits we have obtained $$\begin{aligned}
\overline m_c(\overline m_c)& = &(0.788 + 4.13\times\alpha_s(m_Z)) 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.013)_{\rm syst}\\ &&
\, \pm \, (-0.127+1.23\times\alpha_s(m_Z))_{\rm pert}
\, \pm \, (0.009)_{\alpha_s}  
\, \pm \, (0.002)_{\langle GG\rangle}
\,,\nonumber
\\[2mm]
\overline m_c(3~\mbox{GeV})& = &(1.462 - 4.00\times\alpha_s(m_Z)) 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.014)_{\rm syst}\\ &&
\, \pm \, (-0.148+1.43\times\alpha_s(m_Z))_{\rm pert} 
\, \pm \, (0.002)_{\langle GG\rangle}
\,,\nonumber\end{aligned}$$ for the $\overline{\mbox{MS}}$ charm mass. At the level of uncertainties obtained in our work excellent convergence of perturbation theory was observed. Adopting $\alpha_s(m_Z) = 0.1184 \pm 0.0021$ we then obtain $$\begin{aligned}
\overline m_c(\overline m_c) & = & 1.277 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.013)_{\rm syst} 
\, \pm \, (0.019)_{\rm pert}
\, \pm \, (0.009)_{\alpha_s} 
\, \pm \, (0.002)_{\langle GG\rangle}
\,,\nonumber
\\[2mm]
\overline m_c(3~\mbox{GeV}) & = & 0.998 
\, \pm \, (0.006)_{\rm stat} 
\, \pm \, (0.014)_{\rm syst} 
\, \pm \, (0.022)_{\rm pert}
\, \pm \, (0.010)_{\alpha_s} 
\, \pm \, (0.002)_{\langle GG\rangle}
\,.\nonumber\\\end{aligned}$$ Our perturbative error of $20$ MeV is a factor of ten larger, and the experimental uncertainty of $14$ MeV is by $5$ MeV larger than in the most recent analysis of Ref. [@Kuhn:2007vp]. For estimating the perturbative error a range of scale variations between $\overline m_c(\overline m_c)$ and $4$ GeV was employed. Adding all uncertainties quadratically we obtain $$\begin{aligned}
\label{conclusionmc}
\overline m_c(\overline m_c) & = & 1.277 
\, \pm \, 0.026~\mbox{GeV}
\,,\nonumber
\\[2mm]
\overline m_c(3~\mbox{GeV}) & = & 0.998\,\pm \, 0.029~\mbox{GeV}
\,,\end{aligned}$$ giving an uncertainty that is twice the size of the one obtained in Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp]. This difference arises mainly from the more appropriate estimate of perturbative uncertainties obtained in this work. The agreement in the central value with Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp] is a coincidence, arising from an interplay of differences in the experimental input and theoretical treatment.

As a final thought one might ask which further improvements might be possible in the future. As can be seen from Tab. \[tab:momres1\], from the experimental side the biggest improvement could be made from more accurate measurements of the $J/\psi$ and $\psi^\prime$ electronic partial widths. The current relative uncertainties are $2.5$% and $2.4$%, respectively. Here some improvement might be conceivable with dedicated measurements. On theory side, viewing the uncertainties, it is not unreasonable to assume that the computation of ${\cal
  O}(\alpha_s^4)$ moments of the vector current correlator might further reduce the perturbative error below the level of $20$ MeV. Using the $\overline{\mbox{MS}}$ scheme for the charm mass the OPE states that the remaining perturbative infrared renormalon ambiguity is of order $\Lambda_{\rm
  QCD}^4/\overline m_c^3\sim {\cal O}(5 \mbox{\,-}15~\mbox{MeV})$. This expectation has also been confirmed by explicit bubble chain (large-$\beta_0$) computations [@Grozin:2004ez] and indicates that a further reduction of the perturbative uncertainty is not excluded.

However, to throw in some words of caution, at the level of the present perturbative uncertainties one should also remind oneself about possible loopholes still left in the charmonium sum rule method. An issue we would like to mention concerns the separation of charm and non-charm hadronic production rates needed to carry out the charmonium sum rule. On the theory side the issue is conceptually subtle due to the singlet and secondary radiation contributions which arise at ${\cal O}(\alpha_s^3)$ and ${\cal O}(\alpha_s^2)$, respectively. In this work (as well as in Ref. [@Chetyrkin:2009fv]) both contributions have been considered as non-charm although they contain terms belonging to the $c\bar c$ final state. This treatment might be justified since the size of the corresponding terms is quite small (see Tab. \[tabneglected\]) and since it is the common approach to determine the experimental charm production rate in the continuum region by subtracting theoretical results (or models) for the non-charm rate from the measured total hadronic . In our method to determine the experimental moments this subtraction involves a normalization constant multiplying the theoretical non-charm that is fitted within our clustering method as well accounting for data below and above the charm threshold. The result (see Eq. (\[eq:ns-results\])) reveals a disparity of $4$% between the theoretical non-charm and the data. Setting, in contrast, the normalization constant to unity results in a shift in the charm mass by $-15$ MeV.[^18] Since this shift is compatible with the overall systematical uncertainty in the experimental data we have not treated it as an additional source of uncertainty. On the other hand, the size of the shift could also be considered as an inherent conceptual uncertainty related to separating the charm from the non-charm , which is based on theory considerations rather than on experimental methods and which apparently cannot be improved simply by higher order perturbative computations. As an alternative, one might avoid the separation of the charm and the non-charm contributions altogether and use the total hadronic cross section for the charm mass fits. Apart from the shift mentioned above such an approach would, however, also lead to a substantially larger dependence on the uncertainties in $\alpha_s$. Given these considerations we believe that a substantial reduction of the uncertainties also relies on a resolution of the disparity mentioned above. This might certainly involve more precise measurements in the charm threshold and below-threshold regions, but also some deeper conceptual insight. Until then a substantial reduction of the uncertainties shown in Eq. (\[conclusionmc\]) appears hard to achieve without imposing ad-hoc assumptions.

This work was supported in part by the European Community’s Marie-Curie Research Networks under contract MRTN-CT-2006-035482 (FLAVIAnet), MRTN-CT-2006-035505 (HEPTOOLS) and PITN-GA-2010-264564 (LHCphenOnet), and by the U.S. Department of Energy under the grant FG02-94ER40818. V. Mateu has been partially supported by a DFG “Eigenen Stelle” under contract MA 4882/1-1 and by a Marie Curie Fellowship under contract PIOF-GA-2009-251174. S. Zebarjad thanks the MPI for hospitality while part of this work was accomplished. S. Zebarjad and V. Mateu are grateful to the MPI guest program for partial support. We thank S. Schutzmeier for confirmation of our numerical ${\cal O}(\alpha_s^3)$ fixed-order results. A. Hoang acknowledges discussion with H. Kühn and C. Sturm.

Appendix A: Results of the Fit Procedure {#appendix-a-results-of-the-fit-procedure .unnumbered}
========================================

\[ap:FitOutcome\]

In this appendix we present in some more detail the numerical results of our fit procedure. In Tabs. \[tab:FitFunction-Standard\], \[tab:FitFunction-Minimal\] and \[tab:FitFunction-Maximal\], the results for the cluster energies and the cluster charmed are shown for the standard, minimal and maximal selection of data sets, respectively, using our default setting for the correlations. We use the results for the standard data set selection for our final charm mass analysis. The numbers in the parentheses correspond to the statistical and systematical errors. The correlation matrices for the is available, but cannot be displayed due to lack of space. They can be obtained by the authors on request. For the three data selections, the fit gives the following minimal $\chi^2$ per degree of freedom, $$\dfrac{\chi_{\rm standard}^2}{\rm dof}=1.89\,, \qquad
\dfrac{\chi_{\rm minimal}^2}{\rm dof}=1.86\,,
\qquad \dfrac{\chi_{\rm maximal}^2}{\rm dof}=1.81\,,
\label{eq:chi2-results}$$ and the following normalization constants for the non-charm background $$\begin{aligned}
\!\!\!\!\!\!n^{\rm standard}_s&=&1.039 \pm 0.003_{\rm stat} \pm 0.012_{\rm syst},\;  n^{\rm minimal}_s=1.029 \pm 0.003_{\rm stat} \pm 
0.015_{\rm syst},
\label{eq:ns-results} \\
\!\!\!\!\!\!n^{\rm maximal}_s&=&1.023 \pm 0.003_{\rm stat} \pm 0.011_{\rm syst}.\nonumber\end{aligned}$$ Since the minimal $\chi^2/{\rm dof}$ values are not close to unity, one has to conclude that fit quality is not really very good. This is not at all visible from the agreement of the fit and the data for the total cross section (see Figs. \[fig:fitcompilation\]) and thus might be related to the disparity between the fits of charm versus non-charm production rates described in Sec. \[sectionconclusions\].

In Eq. (\[correlmatrices\]) we show for the correlation matrices of the first four experimental moments for the minimal and the maximal data set selection. The results for our standard selection are given in Eqs. (\[CMtotal\]) and (\[CMuncorr\]). All numbers are related to moments $M_n^{\rm exp}$ normalized to units of $10^{-(n+1)}\,\mbox{GeV}^{-2n}$. We show the results accounting for the full set of correlated and uncorrelated uncertainties and the correlation matrices accounting only for uncorrelated systematical and statistical uncertainties (subscript uc).

$$\begin{aligned}
\label{correlmatrices}
C^{\rm exp}_{\rm min}&=&\left(
\begin{array}{cccc}
0.304 & 0.194 & 0.165 & 0.157\\
0.194 & 0.133 & 0.115 & 0.109\\
0.165 & 0.115 & 0.100 & 0.096\\
0.157 & 0.109 & 0.096 & 0.092\\ 
\end{array}
\right)\,,\quad C^{\rm exp}_{\rm uc}\,=\,\left(
\begin{array}{cccc}
0.047 & 0.038 & 0.035 & 0.035\\
0.038 & 0.035 & 0.035 & 0.035\\
0.035 & 0.035 & 0.035 & 0.036\\
0.035 & 0.035 & 0.036 & 0.037
\end{array}\right)\,,\\\nonumber\\
C^{\rm exp}_{\rm max}&=&\left(
\begin{array}{cccc}
0.193 & 0.138 & 0.125 & 0.122\\
0.138 & 0.108 & 0.099 & 0.097\\
0.125 & 0.099 & 0.091 & 0.090\\
0.122 & 0.097 & 0.090 & 0.088\\
\end{array}
\right)\,,\quad 
C^{\rm exp}_{\rm uc}\,=\,\left(
\begin{array}{cccc}
0.037 & 0.034 & 0.033 & 0.034\\
0.034 & 0.033 & 0.034 & 0.035\\
0.033 & 0.034 & 0.035 & 0.036\\
0.034 & 0.035 & 0.036 & 0.037
\end{array}\right)\,.\nonumber\end{aligned}$$

Appendix B: On the Equivalence of $\chi^2$-Functions {#appendix-b-on-the-equivalence-of-chi2-functions .unnumbered}
====================================================

In this appendix we demonstrate that a in which the auxiliary fit parameters $d_k$, which describe the correlated deviation off the experimental central value within experiment $k$, multiplies only the experimental systematical uncertainties,[^19] $$\chi^2 = \sum_k\left[d_k^2+\sum_{i,m}\left(\dfrac{R_i^{k,m}+d_k\,
\Delta_i^{k,m}-R_m}{\sigma_i^{k,m}}\right)^2\right]\,, \label{additive}$$ is mathematically equivalent to the well known written solely in terms of the fit parameters $R_m$ and a correlation matrix, $$\bar \chi^2=\sum_{i,j,k,m,n}(R_i^{k,m}-R_m)(V_k^{-1})_{ij}^{mn}(R_j^{k,n}-R_n)\,,
\label{traditional}$$ where there is no correlation among the different experiments $k$ and $k^\prime$, and the correlation matrix within one experiment $k$ has the form $$(V_k)_{ij}^{mn}=\sigma_i^{k,m\,2}\delta_{ij}\delta^{mn}+\Delta_i^{k,m}\Delta_j^{k,n}
\,.
\label{tradcorr}$$ The matrix Eq. (\[tradcorr\]) can be inverted analytically giving $$(V_k^{-1})_{ij}^{mn} = \dfrac{\delta_{ij}\delta^{mn}}{(\sigma_i^{k,m})^2}-
\dfrac{1}{1+\dfrac{\left\langle \Delta^2\right\rangle_k}{\left\langle \sigma\right\rangle_k^2}}
\dfrac{\Delta_i^{k,m}}{(\sigma_i^{k,m})^2}\dfrac{\Delta_j^{k,n}}{(\sigma_j^{k,n})^2}\,,
\label{invtradcorr}$$ where we have defined the mean statistical error and the statistical average of the systematical error within one experiment as follows: $$\left\langle \sigma\right\rangle_k^2 \equiv 
\left(\sum_{i,m}\dfrac{1}{(\sigma_i^{k,m})^2}\right)^{-1},\qquad
\left\langle \Delta^2
\right\rangle_k\equiv\left\langle \sigma\right\rangle_k^2\sum_{i,m}
\dfrac{(\Delta_i^{k,m})^2}{(\sigma_i^{k,m})^2}\,.
\label{averages}$$

In case of a sizable positive correlation between measurements in the same experiment (such that $\Delta_i^{k,m}/R_i^{k,m}$ is sizable and approximately constant) is it known that the form of the in Eq. (\[additive\]) leads to best fit values $R_m$ that are systematically below the measurements [@D'Agostini:1993uj]. Our proof demonstrates that the standard of Eq. (\[traditional\]) has the same property for correlation matrices with the form of Eq. (\[tradcorr\]). This motivates to use the so-called minimal-overlap correlation model where the second term on the RHS of Eq. (\[traditional\]) is replaced by $\mbox{min}^2(\Delta_i^{k,m},\Delta_j^{k,n})$. In general, within the minimal-overlap model, the correlations are sufficiently reduced such that the unphysical effect described above does not arise.

We proceed by showing that one can “integrate out" auxiliary parameters $d_k$ in Eq. (\[additive\]) obtaining a new function $\tilde \chi^2(R^i)$ which yields the same results for the best fit for the $R_m$, as long as one works in the Gaussian approximation. The minimum of $\chi^2(R_i,d_k)$ is located at the best fit values (indicated by superscripts $(0)$) defined by the conditions $$\left.\dfrac{\partial\chi^{2}}{\partial R_{i}}\right|_{R_i^{(0)},d_k^{(0)}}=\left.\dfrac{\partial\chi^{2}}{\partial d_{j}}\right|_{R_i^{(0)},d_k^{(0)}}=0\,.\label{Mincond}$$ To invert the matrix of second derivatives we proceed in blocks $$\left(\begin{array}{cc}
\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}} & \dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial d_{k}}\\
\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial R_{j}} & \dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial b_{k}}\end{array}\right)_{\tilde{R},\tilde{b}}\left(\begin{array}{cc}
c_{jm} & b_{jr}\\
b_{mk} & a_{kr}\end{array}\right)=\left(\begin{array}{cc}
\delta_{im} & 0\\
0 & \delta_{lr}\end{array}\right).$$ where $c_{ij}$ and $a_{kr}$ are $N_{\rm cluster}\times N_{\rm cluster}$ and $N_{\rm exp}\times N_{\rm exp}$ square matrices, respectively, and $b_{jr}$ is a $N_{\rm cluster}\times N_{\rm exp}$ rectangular matrix. We find the following four matrix relations $$\begin{aligned}
&&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
\sum_{j}\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}}\, c_{jm}+\sum_{k}\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial d_{k}}\, b_{mk}=\delta_{im}\,, \quad 
\sum_{j}\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}}\,
b_{jr}+\sum_{k}\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial d_{k}}\,
a_{kr}=0\,,
\label{eq:general_block_1}\\
&&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
\sum_{j}\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial R_{j}}\,
c_{jm}+\sum_{k}\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial d_{k}}\,
b_{mk}=0\,, 
\qquad\;\, \sum_{j}\dfrac{\partial^{2}\chi^{2}}{\partial 
d_{l}\partial R_{j}}\, b_{jr}+\sum_{k}\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial d_{k}} a_{kr}=\delta_{lr}\,.\label{eq:general_block_2}\end{aligned}$$ Combining Eqs. (\[eq:general\_block\_1\]a) and (\[eq:general\_block\_2\]a) we find the inverse of the upper left , $$(c^{-1})_{ij}=2(V^{-1})_{ij}=\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}}-\sum_{kl}\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial R_{i}}\left[\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial 
d_{k}}\right]^{-1}\dfrac{\partial^{2}\chi^{2}}{\partial d_{k}\partial R_{j}}\,,
\label{eq:noEOM_Rj_bk}$$ where $\left[\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial
    d_{k}}\right]^{-1}$ stands for the $(l,k)$ element of the inverse matrix of $\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial d_{k}}$ (and not the inverse of the element). Combining Eqs. (\[eq:general\_block\_1\]b) and (\[eq:general\_block\_2\]b) one can obtain relations for the $a$ and $b$ submatrices: $$\begin{aligned}
b_{jr} & = & -\sum_{i,k}\left[\dfrac{\partial^{2}\chi^{2}}{\partial R_{j}\partial R_{i}}\right]^{-1}\dfrac{\partial^{2}\chi^{2}}{\partial d_{k}\partial R_{i}}\, a_{kr}\,,\\
\delta_{lr}&=&\sum_{k}\left(\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial d_{k}}-\sum_{i,j}\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial R_{j}}\left[\dfrac{\partial^{2}\chi^{2}}{\partial R_{j}\partial 
R_{i}}\right]^{-1}\dfrac{\partial^{2}\chi^{2}}{\partial d_{i}\partial R_{k}}\right)a_{kr}\,,\\
\delta_{im} & = & \sum_{j}\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}}\, c_{jm}-\sum_{k,l,s}\dfrac{\partial^{2}\chi^{2}}{\partial d_{k}\partial R_{i}}\,\left[\dfrac{\partial^{2}\chi^{2}}{\partial R_{m}\partial 
R_{s}}\right]^{-1}\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial R_{s}}\, a_{lk}\,,\end{aligned}$$ where again $\left[\dfrac{\partial^{2}\chi^{2}}{\partial R_{j}\partial R_{i}}\right]^{-1}$ stands for the $(j,i)$ element of the inverse matrix of $\dfrac{\partial^{2}\chi^{2}}{\partial R_{j}\partial R_{i}}$. We “integrate out" the auxiliary parameters $d_k$ by substituting their minimum conditions $d^{(0)}_k=\tilde d_k(R^{(0)}_i)$ (which is analogous to using the equation of motion when integrating out heavy particles): $$\left.\dfrac{\partial\chi^{2}}{\partial d_{j}}\right|_{d_k=\tilde d_k(R_i)}=0\,,\qquad \tilde \chi^2(R_i)=\chi^2(R_i,\tilde d_k(R_i))\,.\label{EOM}$$ Their first derivatives with respect to $R_{j}$ read $$\begin{aligned}
\dfrac{\partial}{\partial R_{j}}\left.\dfrac{\partial\chi^{2}}{\partial d_{i}}\right|_{\tilde d_m(R_{k})} & = & \left.\dfrac{\partial^{2}\chi^{2}}{\partial d_{i}\partial R_{j}}\right|_{\tilde 
d_m(R_{k})}+\sum_{l}\left.\dfrac{\partial^{2}\chi^{2}}{\partial d_{i}\partial d_{l}}\right|_{{\tilde d}_{m}(R_{k})}\dfrac{\partial {\tilde d}_{l}(R_{n})}{\partial R_{j}}=0\,,\nonumber \\
\dfrac{\partial \tilde d_i(R_{n})}{\partial R_{j}} & = & -\sum_l\left[\dfrac{\partial^{2}\chi^{2}}{\partial d_i\partial d_l}\right]_{\tilde d_m(R_{k})}^{-1}\,\left.\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial 
R_{j}}\right|_{\tilde d_m(R_{k})}
\,.\end{aligned}$$ The minimum of $\tilde \chi^2(R_i)$ is indeed located at $R^{(0)}_{i}$ because $$\dfrac{\partial\tilde\chi^2}{\partial R_{i}}=\sum_k\left.\dfrac{\partial\chi^{2}}{\partial d_k}\right|_{d_k=\tilde d_k(R_i)}\dfrac{{\rm d} \tilde d_k}{{\rm d} R_i}+
\left.\dfrac{\partial\chi^{2}}{\partial R_i}\right|_{d_k=\tilde
  d_k(R_i)}=\left.\dfrac{\partial\chi^{2}}{\partial R_i}\right|_{d_k=\tilde
  d_k(R_i)}
\,,\label{first_der}$$ and because the first term vanishes by the condition in Eq. (\[EOM\]). When evaluating Eq. (\[first\_der\]) for $R_i=R_i^{(0)}$ it vanishes because of Eq. (\[Mincond\]). Finally, let us calculate the inverse correlation matrix: $$\begin{aligned}
&&\left.\dfrac{\partial^{2}\tilde{\chi}^{2}}{\partial R_{i}\partial R_{j}}\right|_{\tilde{R}_{k}} = 
\left.\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}}\right|_{d^{(0)}_{m},R^{(0)}_{k}}+
\sum_{k}\left.\dfrac{\partial^2\chi^{2}}{\partial d_k\partial R_j}\right|_{d^{(0)}_m,R^{(0)}_k}
\left.\dfrac{\partial \tilde d_{k}(R_{i})}{\partial R_{i}}\right|_{\tilde{R}_{k}}\nonumber \\
 & = & \left.\dfrac{\partial^{2}\chi^{2}}{\partial R_{i}\partial R_{j}}\right|_{d^{(0)},R^{(0)}_{k}}-\sum_{k,l}\left[\dfrac{\partial^{2}\chi^{2}}{\partial d_{k}\partial 
d_{l}}\right]_{d^{(0)}_{m},R^{(0)}_{k}}^{-1}\,\left.\dfrac{\partial^{2}\chi^{2}}{\partial d_{l}\partial R_{i}}\right|_{d^{(0)}_{m},R^{(0)}_{k}}\left.\dfrac{\partial^{2}\chi^{2}}{\partial d_{k}\partial 
R_{j}}\right|_{d^{(0)}_{m},R^{(0)}_{k}},\label{eq:EOM_Rj_bk}\end{aligned}$$ which agrees with Eq. (\[eq:noEOM\_Rj\_bk\]).

We can now apply the previous results to Eq. (\[additive\]): $$\begin{aligned}
\chi^{2} & = & \sum_{k}\left[d_k^2+\sum_{i,m}\left(\dfrac{R_i^{k,m}+d_k\,
\Delta_i^{k,m}-R_m}{\sigma_i^{k,m}}\right)^2\right]\label{eq:with_bk_Rj}\\
 & = & \sum_k\left\{ d_k^2\left[1+
\dfrac{\left\langle \Delta^2\right\rangle_k}{\left\langle \sigma\right\rangle _k^2}\right]+2\, d_k\sum_{i,m}\dfrac{(R_i^{k,m}-R_m)\Delta_i^k}{\sigma_{i}^{k,m\,2}}+\sum_{i,m}
\left(\dfrac{R_i^{k,m}-R_{m}}{\sigma_i^{k,m}}\right)^2\right\} 
\,.\nonumber\end{aligned}$$ The equation of motion for $d_{k}$ reads $$d_k(R)=-\,\dfrac{1}{1+\dfrac{\left\langle \Delta^2\right\rangle_k}{\left\langle \sigma\right\rangle_k^2}}
\sum_{i,m}\dfrac{(R_i^{k,m}-R_m)\Delta_i^{k,m}}{\sigma_i^{k,m\,2}}
\,.$$ This renders for $\tilde\chi^2$ the form $$\begin{aligned}
\tilde\chi^2 & = & \sum_k\left[d_k(R)\sum_{i,m}\dfrac{(R_i^{k,m}-R_m)\Delta_i^{k,m}}{\sigma_i^{k,m\,2}}+
\sum_{i,m}\left(\dfrac{R_i^{k,m}-R_m}{\sigma_i^{k,m}}\right)^2\right]\\
 & = & \sum_k\left[\sum_{i,m}\left(\dfrac{R_i^{k,m}-R_m}{\sigma_i^{k,m}}\right)^2-
\dfrac{1}{1+\dfrac{\left\langle \Delta^2\right\rangle_k}{\left\langle \sigma\right\rangle_k^2}}\sum_{i,j,m,n}\dfrac{(R_i^{k,m}-R_m)\Delta_i^{k,m}}
{\sigma_i^{k,m\,2}}\dfrac{(R_j^{k,n}-R_n)\Delta_j^{k,n}}{\sigma_j^{k,n\,2}}\right]\,,\nonumber \end{aligned}$$ which reproduces Eq. (\[traditional\]), as we wanted to demonstrate.

We conclude this appendix presenting an alternative way to write Eq. (\[chi2total\]) after using the equations of motion for $d_k$. We again concentrate on the simpler case without subtraction of the non-charm contribution: $$\begin{aligned}
\chi^2 & = & \sum_k\left\{ d_k^{2}+\sum_{i,m}\left[\dfrac{R_i^{k,m}-(1+d_k\Delta f_k^{i,m})\, R_m}{\sigma_k^{i,m}}\right]^{2}\right\} \nonumber \\
 & = & \sum_k\left\{ \left[1+\sum_{i,m}
\left(\dfrac{\Delta f_k^{i,m}\, R_m}{\sigma_k^{i,m}}\right)^2\right]d_k^2-2d_k\sum_{i,m}\dfrac{\left(R_{k}^{i,m}-R_m\right)
\Delta f_k^{i,m}\, R_m}{(\sigma_k^{i,m})^2}\right.\nonumber \\
&+&\left.\sum_{i,m}\left[\dfrac{R_k^{i,m}-R_m}
{\sigma_k^{i,m}}\right]^2\right\} .\label{eq:chi2_auxiliary}\end{aligned}$$ The EOM for $d_k$ now reads $$\left[1+\sum_{i,m}\left(\dfrac{\Delta f_k^{i,m}\, R_m}{\sigma_k^{i,m}}\right)^2\right]
\tilde d_k(R_n)-\sum_{i,m}\dfrac{\left(R_k^{i,m}-R_m\right)\Delta f_k^{i,m}\,
  R_m}{(\sigma_k^{i,m})^2}=0
\,,$$ which upon insertion into Eq. (\[eq:chi2\_auxiliary\]) renders $$\begin{aligned}
\widetilde{\chi}^2 & = & \sum_k\left\{
  \sum_{i,m}\left[\dfrac{R_k^{i,m}-R_m}{\sigma_k^{i,m}}\right]^2-\tilde
  d_k(R_n)\sum_{i,m}\dfrac{\left(R_k^{i,m}-R_m\right)\Delta f_k^{i,m}\, R_m} 
{(\sigma_k^{i,m})^2}\right\} \nonumber\\
 & = & \sum_k\Biggr\{ \sum_{i,m}\left[\dfrac{R_k^{i,m}-R_m}{\sigma_k^{i,m}}\right]^2\nonumber\\
&-&
\dfrac{1}{\left[1+\sum\left(\dfrac{\Delta f_k^{i,m}\,
 R_m}{\sigma_k^{i,m}}\right)^2\right]}
\left[\sum_{i,m}\dfrac{\left(R_k^{i,m}-R_m\right)\Delta_k^{i,m}\, R_m}
{(\sigma_k^{i,m})^2}\right]^2\Biggr\}
\,.\label{eq:chi2_final} \end{aligned}$$ One can rewrite Eq. (\[eq:chi2\_final\]) in the matrix form $$\begin{aligned}
\widetilde{\chi}^2 & = & \sum_k\left\{ \sum_{i,j,m,n}\left[R_{k}^{im}-R_m\right]\left[V_{k}^{-1}\right]^{mn}_{ij}\left[R_{k}^{jn}-R_n\right]\right\} \,,\label{eq:Matrix-Form} \\
\left[V_{k}^{-1}\right]^{mn}_{ij} & = & \dfrac{\delta_{ij}\delta^{mn}}{(\sigma_{k}^{i,m})^{2}}-
\dfrac{\left(R_k^{i,m}-R_m\right)\Delta f_k^{i,m}\, R_m
\left(R_k^{jn}-R_n\right)\Delta f_k^{jn}\, R_n}
{\left[1+\sum\left(\dfrac{\Delta_k^{i,m}\, R_m}{\sigma_k^{i,m}}\right)^2\right]}\,.\nonumber\end{aligned}$$ In Eq. (\[eq:Matrix-Form\]) one can interpret the second term of the inverse correlation matrix, as a non-linear correlation among the measurements, where the correlation matrix itself depends on the value of the fit parameters. The total inverse correlation matrix is block diagonal, and the blocks correspond to $V_k^{-1}$.

[^1]: As a word of caution we mention that, except for the contributions of the $J/\psi$ and $\psi'$ resonances, no experimental separation of the charm and non-charm contributions in the hadronic cross section has been provided in available data, although charm-tagged measurements are possible, see e.g. Ref. [@CroninHennessy:2008yi] (CLEO collaboration). So the charm pair production rate from above the $J/\psi$ and $\psi^\prime$ that enters the charmonium sum rules in Eq. (\[momentdef2\]) is usually obtained partly from the measured total with theory motivated subtractions of the non-charm rate, and partly by using theory predictions for the charm production rate.

[^2]: Since the analyses of Refs. [@Chetyrkin:2006xg; @Boughezal:2006px] were based on outdated and less precise data for the $J/\psi$ and $\psi^\prime$ electronic partial widths [@Yao:2006px], we frequently only compare our numerical results with those of Refs. [@Chetyrkin:2009fv; @Kuhn:2007vp]. However the perturbative input of Refs. [@Chetyrkin:2006xg; @Boughezal:2006px; @Chetyrkin:2009fv; @Kuhn:2007vp] and ours is identical.

[^3]: This can to happen at any order for $n>1$ and for $\mu_m\sim 3$ GeV.

[^4]: This works in general as long as the path dependent $\mu_\alpha^c(q^2,\overline m_c^2)$ does not produce a spurious cut in $\alpha_s$ starting at $q^2=0$ and running towards $-\infty$. This condition is implemented into Eq. (\[mualphacontour\]) by $(\mu_\alpha^c)^2$ being negative along the physical cut of the vacuum polarization function above the charm pair threshold.

[^5]: In this examination and throughout our other analyses we use the $\overline{\rm
      MS}$ renormalization group equations with $n_f=4$ active running flavors.

[^6]: The numerical value of $\overline m_c(\mu=3~\mbox{GeV})$ is actually irrelevant, since the running involves $\overline m_c$ only in a linear way and exactly cancels in the ratio.

[^7]: For the discussions in this section we use $\alpha_s^{(n_f=5)}(m_Z)=0.118$ ($\alpha_s^{(n_f=4)}(4.2~\mbox{GeV})=0.2245$) as an input using five-to-four flavor matching at $4.2$ GeV. For the first moment employed for the charm mass fits we use $M_1=0.2138\,{\rm GeV}^{-1}$, corresponding to the standard selection of data sets, as described in Sec. \[sectiondata\].

[^8]: Of course the extracted mass depends on the moment considered as well. Since in our analysis we focus on the first moment only, for simplicity we drop that label.

[^9]: We are grateful to Thomas Schutzmeier for confirming agreement of the results of our ${\cal O}(\alpha_s^3)$ fixed-order code with theirs of Ref. [@Boughezal:2006px].

[^10]: There are 18 references quoted since Ref. [@Edwards:1990pc] provides results from two independent runs that we treat as two different data sets.

[^11]: The situation is quite different if perturbative QCD is used already from $4.8$ GeV. Here the contribution to the moments $M_{1,2,3,4}^{\rm exp}$ from energies above $4.8$ GeV amounts to $(30, 10, 3,
    1)\%$. For the first moment the $10$% penalty would then represent the largest source of uncertainty and correspond to an uncertainty in the charm mass of around $18$ MeV.

[^12]: We have checked that the effect of using $R_{uds}(E_m)$ instead of $R_{uds}(E_i^{k,m})$ is totally negligible.

[^13]: We prove in Appendix B the equivalence of a bilinear with fit auxiliary parameters to a bilinear without auxiliary fit parameters, but containing the standard correlation matrix.

[^14]: In a very good approximation, it can be also obtained by dropping in $V^{-1}_{i,j}$ the rows and columns corresponding to $d_k$ and inverting the resulting matrix. After that one also drops the raw and column corresponding to $n_{\rm nc}$. We adopt this simplified procedure for our numerical fits.

[^15]: For $R_0$ we have $m^\prime=1$ and for $R_{N_{cl}+1}$ we have $m^\prime=N_{cl}-1$.

[^16]: We use a modified version of the minimal correlation model. The non-diagonal entries of the correlation matrix are filled in with $\Gamma_i\, R_m\, {\rm Min}^2\left\lbrace
        {\Delta{\Gamma_i}/\Gamma_i, \Delta R_m/R_m}\right\rbrace
      $. Here $\Delta\Gamma_i$ and $\Delta R_m$ represent the systematical uncertainties of the width of the narrow resonance and the $R$ value of the cluster, respectively, and $i=1,2$ refer to $J/\psi$ and $\psi^\prime$.

[^17]: As our default we use $\alpha_s^{n_f=5}(m_Z)$ as the input, use the four-loop QCD beta-function for the renormalization group evolution and three-loop matching conditions to the $n_f=4$ theory at $\mu=4.2$ GeV.

[^18]: The same disparity was found in Ref. [@Kuhn:2007vp]. In their analysis the corresponding effect is $-5$ MeV since experimental data were used in the experimental moments only for energies $\sqrt{s}<4.8$ MeV. In this approach, however, the moments are strongly dependent on the uncertainty one assigns to the theory input used for $\sqrt{s}>4.8$ MeV.

[^19]: For the discussion in this appendix we do for simplicity of the presentation not account for fits of the non-charm contribution. It is, however, straightforward to generalize the presentation to this case.
