MSbar Charm Mass from Charmonium Sum Rules with Contour Improvement

A detailed error analysis is carried out for the determination of the MSbar charm quark mass $\bar m_c(\bar m_c)$ from moments at order alpha_s^2 of the charm cross section in e^+e^- annihilation. To estimate the theoretical uncertainties the renormalization scale is implemented in various ways including energy-dependent functions, which lead to ``contour-improved'' predictions. We obtain $\bar m_c(\bar m_c)=1.29\pm 0.07$ GeV which contains a substantial theoretical uncertainty.


Introduction
At the ongoing and future B-physics experiments a realistic estimate of the uncertainties in the values of the bottom and charm quark masses will become increasingly important for the measurements of the CKM parameters and the search for new physics, particularly from inclusive B-decay rates. [1,2] Moments of the respective e + e − R-ratios [3,4], where R qq = σ(e + e − → qq + X)/σ(e + e − → µ + µ − ), are the most important instrument at present to determine the bottom and charm quark mass parameters. (See also Refs. [5,6] for recent bottom quark mass extractions from spectral moments in semileptonic B-decays.) In general, one can distinguish between two regions in n, which require a different theoretical treatment. For low values of n the moments are dominated by relativistic dynamics and scales of order of the heavy quark mass m q . This allows the use of the usual expansion in the number of loops for the theoretical computations, and the MS scheme is an appropriate choice for the heavy quark mass parameter. However, the lack of data for R qq in the continuum regions above the quarkonium resonances introduces model-dependent errors that have not been quantified in most of the previous analyses. (See Refs. [7] for recent reviews.) On the other hand, for large values of n the continuum regions are suppressed and the moments become dominated by the quarkonium resonance region where good sets of data have been obtained in the past. However, the theoretical predictions of moments for large values of n is more complicated since the usual loop expansion breaks down and the size of non-perturbative effects increases. Here, summations of higher order contributions proportional to powers of (α s √ n) need to be carried out in order to capture the relevant non-relativistic perturbative information [8,9], and so-called threshold masses [10,11,12] are appropriate schemes for the heavy quark mass parameter. Moreover, there is an upper "duality" bound for the possible choices of n since the energy range contributing to the moments, which is of order m q /n, needs to be larger than the typical hadronization scale Λ QCD . [9,13] In the case of R bb and the determination of the bottom mass, this bound is around n = 10, and the two ranges in n are believed to be well separated with their boundary being approximately at n = 4. A good number of analyses exists for small and large values of n and respecting the cancellation of the dominant renormalon contributions associated to the choice of the quark mass definition. [14,15,16,17,18] (See also Ref. [19] for an early analysis using the MS bottom quark mass for larger values of n.) For R cc and the determination of the charm mass, the distinction between largeand low-n moments is more delicate because m c is not much larger than Λ QCD . Here, the upper duality bound for n is around 3 or 4. This leaves basically no space at all to carry out the non-relativistic summations that can be applied in the bottom quark case because the corresponding techniques are only valid for large n and as long as m q /n is larger than Λ QCD . On the other hand, even for n ≤ 4 the non-relativistic region close to the cc threshold can have a considerable contribution to the moments, while the model-dependences from the experimentally unknown continuum region can still be significant. For a reliable (error) analysis these issues need to be taken into account.
In this paper we determine the MS charm quark mass m c (m c ) from moments P n with n ≤ 4 using the usual loop expansion in powers of α s . In a previous work [14] an error of less than 30 MeV was obtained using results at order α 2 s and a fixed renormalization scale µ between m c and 4m c . This analysis assumed that the moments are dominated by scales of order m c . In our analysis we focus on the theoretical uncertainties coming from the low-energy region close to the cc threshold. It is well known that close to threshold, where the c.m. velocity β of the quarks is small, the dominant perturbative contributions to R cc are in fact governed by renormalization scales of the order of the relative quark momentum, µ ≃ m c β. [20,21] We believe that this effect should not be ignored in determinations of the charm quark mass since it accounts for important higher-order information. It is the main aim of this work to study the impact of these higher-order contributions. In particular, we find that they affect the estimate of the theoretical error. We also carry out a careful study of all sources of uncertainties including the continuum regions where no data is available.

Theoretical Moments
For the QCD parameters used in this work we adopt the MS renormalization scheme and the convention that the charm quark participates in the running (n f = 4). The theoretical moments are directly parametrized in terms of m c ≡ m c (m c ), and the masses of the up, down and strange quarks are set to zero. For the running of the strong coupling we use four-loop renormalization group equations, and three-loop matching conditions at the scale m b (m b ) = 4.2 GeV where we switch to n f = 4 flavors. We use three different methods to implement the renormalization scale for the theoretical predictions of the moments and apply order α 2 s results obtained earlier in Refs. [20,22,23,24].
The first method is simply based on a fixed energy-independent renormalization scale, µ 2 ≃ m 2 c , in analogy to previous analyses. In the OPE including the first nonperturbative contribution from the dimension-four gluon condensate, the moments have the form P n = P pert where P pert Here, β 0 = 11 − 2/3n f . The order α s terms were given in Refs. [3,4] and the coefficients for the condensate contribution were taken from Refs. [25,26] adopting the renormalization-group-invariant normalization of the gluon condensate. The numerical results for the coefficients in Eq. (2) for n = 1, 2, 3, 4 are collected in Tab. 1.
The numbers for f 0 n and f 1 n agree with results given in Ref. [14]. The numbers for f 2 n differ because we also include the contributions from secondary cc production, where the charm pair is produced through gluon splitting off primary massless quarks. [22] For the determination of m c (m c ) described in subsequent sections we use µ 2 = ξ 2 M 2 with M = 1.3 GeV as the renormalization scale and vary the parameter ξ in order to estimate the perturbative uncertainties. The second method is based on an energy-dependent renormalization scale of order of the c.m. energy, µ 2 ≃ s. We note that the same energy-dependent scale also needs to be used for the implementation of the MS charm mass m c (m c ) in order to avoid an incomplete cancellation of the large higher order corrections associated to the O(Λ QCD ) renormalon of the pole mass definition. A technical issue is related to the fact that in the MS scheme for the charm mass parameter, R cc is more singular at the threshold s = 4m c than in the pole scheme. This makes the computation of the perturbative contributions of the moments based on the dispersion integral in Eq. (1) somewhat cumbersome. It is therefore advantageous to carry out the computation using a contour integration in the positive complex s-plane (Re(s) > 0) around the cut of the vacuum polarization function Π cc associated to R cc , The same procedure can of course also be applied for the first method. The path C of the integration is illustrated in Fig. 1. For practical purposes we uses ≡ s/(4m c )  Table 2: Coefficients of the theoretical expressions for the moments P n at order α 2 s for massless light quarks based on energy-dependent renormalization scales as described in the text for ξ = 1 (ξ = 2). as the integration variable and µ 2 = ξ 2 M 2s as the renormalization scale. This choice allows to explicitly factor out the dependence of the moments on m c . The perturbative part of the moments can then be written in the form where the superscripts 0, 1, 2 refer to O(1), O(α s ) and O(α 2 s ), respectively. Exemplarily, for ξ = 1 (ξ = 2) and α n f =5 s (M z ) = 0.118 (i.e. α n f =4 s (4.2 GeV) = 0.225), the coefficients are given numerically in Tab. 2 The third method is based on an energy-dependent renormalization scale of order of the relative three-momentum of the charm quarks, µ 2 ≃ −sβ 2 , where β ≡ 1 − 4m 2 c /s. Again, the perturbative part of the moments is computed using the contour integration in Eq. (5). In the complex µ 2 -plane, α s (µ 2 ) has a cut along the negative real axis. It is therefore necessary to introduce the negative sign. With this choice the cut in α s (µ 2 ≃ −sβ 2 ) agrees with the one of the vacuum polarization function, and contours can be found such that β is never small along the path, and the perturbative description is appropriate. 1 In fact, this is in close analogy to the so-called "contour-improved" approach to compute the hadronic τ -decay width (see e.g. Ref. [27,28]). For our numerical computations we uses = s/(4m 2 c ) as the integration variable and µ 2 = ξ 2 M 2 (1 −s) as the renormalization scale. The perturbative part of the moments can then be written in the form Exemplarily, for ξ = 1 (ξ = 2) and α

Experimental Moments
For the contributions to the experimental moments from the J/ψ and the ψ ′ we use the most recent averages for masses and e + e − widths given by the PDG [29] (M J/ψ = 3.0969 GeV, M ψ ′ = 3.6860 GeV, Γ e + e − J/ψ = 5.14 ± 0.31 keV, Γ e + e − ψ ′ = 2.12 ± 0.12 keV). For the energy region between 3.73 and 4.8 GeV we employ the data for the total R ratio, R tot , obtained by the BES collaboration. [30] To obtain experimental data for R cc we use the last four data points below 3.73 GeV to fit the non-charm R ratio, R nc , assuming it to be energy-independent, and subtract it from the numbers for R tot . The effects from the Z exchange and the small logarithmic energy-dependence expected from theory are below 1% and can be neglected. The determination of the statistical uncertainties of R cc is standard. The determination of the systematical uncertainties in R cc is more subtle because only systematical errors for R tot (with all individual contributions being added quadratically) but no specific numbers for the D-meson final states above 3.7 GeV were given in Ref. [30]. Thus, naively subtracting the systematical errors of R nc from the ones in R tot underestimates the actual systematical errors in R cc . To obtain numbers for the systematical errors in R cc we make the assumption that the systematical errors for the non-charm and the charm final states are equally large, i.e. the systematical errors from R cc are obtained by multiplying those of R tot by a factor 1/ √ 2. This procedure is supported by the observation that the relative systematical errors given in Ref. [30] for energies 1 For the first and the third method one can also choose a closed path around the origin in the counter-clockwise direction.
, where for the electromagnetic coupling [α(3.1 GeV)] −1 = 134.3 has been adopted. above 3.7 GeV are in average approximately a factor 1.5 larger than expected from extrapolating the relative systematical errors for energies below 3.7 GeV.
In the region above 4.8 GeV the only other data available for R tot is from the MD1 experiment [31] and from the CLEO collaboration [32] in the energy region between 7.25 and 10.5 GeV. For the determination of R cc and the uncertainties we use the same method as described in the previous paragraph. The resulting relative systematical errors amount to about 8%.
There is no continuous experimental data for R cc in the regions between 4.8 and 7.25 GeV and above 10.5 GeV. As shown below, these contribution turn out to be not very important. For our analysis we use a model for these contributions based on the theory prediction of R cc at order α 2 s . This is reasonable since it is known that perturbative QCD agrees with data for cc production to several percent at M Z and at the level of a few 10% percent in the LEP2 region. [33]. The resulting modeldependence in the determination of m c is negligible for n > 1, as shown below. For n = 1 the contribution from the region between 4.8 and 7.25 GeV leads to a model-dependent error of about 10 MeV.
In Tab. 3 a collection of all contributions to the first four moments is given. For contributions coming from experimental data the first error is statistical and the second systematical. For the contributions from the J/ψ and the ψ ′ we used the averaged PDG [29] errors. Since detailed numbers for the various error sources are difficult to obtain, we make the simplified assumption that statistical errors (being later combined in quadrature) and systematical errors (being later combined linearly) are equal. This is consistent with information given in the original literature. For the continuum regions that have to be estimated from our theory-model the contributions from below and above M Z are distinguished in order to visualize the impact of these two regions. The respective errors given in Tab. 3

Uncertainties in m c (m c )
For the determination of m c (m c ) we fit single moments using the three different methods for the theoretical predictions described above. The contribution from the gluon condensate is taken from the first method for all cases. We use α  )  11  5  3  2  11  10  9  9  δξ  24  15  10  7  24  20  18  14  combined error  92  56  45  43  94  65  55  52   Table 5: Central values and individual errors for m c (m c ) in units of MeV based on the methods described in the text.
methods yield very similar results, particularly for n > 1. This is expected since the high-energy range, s ≫ 4m 2 c , where the renormalization scales of the two methods differ, is suppressed. The uncertainties from experimental data and from the treatment of the non-perturbative corrections are similar for the three theoretical methods. As mentioned before, the model-dependent error coming from the lack of experimental data in the continuum region is negligible for n > 1 and amounts to only about 10 MeV for n = 1. This is in contrast to determinations of m b (m b ) where the model-dependent error is significantly larger, particularly for n = 1 and n = 2. [15] The uncertainties coming from the error in α s (M Z ) for n > 1 are larger for the third method than for the first and the second one, since variations in the value of α s (M Z ) have a larger impact at lower scales. To obtain combined errors we have added all statistical uncertainties in quadrature, and all other uncertainties linearly. Since the statistical errors are smaller than the other sources of uncertainties, the reader should note that our total errors cannot be interpreted statistically.
It is instructive to have a closer look onto the perturbative uncertainties obtained by the different theoretical methods. For all methods the uncertainties obtained by varying ξ either in the range 3/4 ≤ ξ ≤ 3 or 1 ≤ ξ ≤ 4 become smaller for larger values of n. This behavior is, however, less pronounced for the third method. Thus, each of the methods yields a behavior in contrast to the expectation that the scaling uncertainties grow like the uncertainties from the non-perturbative corrections, since for larger n the moments become dominated by dynamics at lower energies. On the other hand, comparing the results for the central values obtained by the three theoretical methods it is found that they tend to decrease for larger values of n for the first and second method, while there is the opposite effect for the third method. In fact, the difference between the central values for n > 1 is larger than the uncertainties based on varying ξ. This feature is related to the different treatment of higher-order contributions and appears to be similar to the discrepancies found between the "contour-improved" and the "fixed-scale" approaches for predictions of the hadronic τ -decay rate. (See e.g. Refs. [34,35] for discussions on this issue). Interestingly, the convergence behavior of the moments for all three methods is quite good. For example, for n = 3, m c = 1.  35 for the perturbative expansion of the moment for the first, second and the third method, respectively. The discrepancy between the first two methods and the third one is considerably larger than an error estimate based on the size of the order α 2 s terms obtained from each of the methods. These discrepancies represent theoretical uncertainties that are not captured by conventional variations in the renormalization scale. However, they should be included in the estimate of the theoretical error. We include this type of theoretical uncertainty in our procedure to determine the combined results shown in Tab. 4.
As our final result we adopt m c (m c ) = 1.29 ± 0.07 GeV (8) obtained from our fits using the moments for n = 2, 3. They are sufficiently below the duality bound and, at the same time, have only small sensitivity to the continuum region above the ψ ′ resonance. Our result is compatible with earlier analyses based on moments of R cc [14,36,37,38,39], and with recent lattice determinations [40,41,42] as well as with the finite energy sum rule analysis of Ref. [43]. However, our error is much larger than in Refs. [14,37,38] due to our more pertinent treatment of theoretical and experimental uncertainties. On the other hand, the size of our error is comparable to Refs. [36,39] where summations of terms proportional to powers of (α s √ n) [8,9] were partially included to account for non-relativistic effects in the low-energy region s ≃ 4m 2 c . However, the resulting central values are systematically lower than ours by up to 100 MeV. This effect apparently arises from the summation of the higher order terms just mentioned and from the use of moments for n > 4. Both issues are problematic from the conceptual point of view. The size of our error is also comparable to the analysis of Ref. [43], but our central value is about 100 MeV lower. This discrepancy seems to indicate that there is a systematic difference between the sum rules based on Eq. (5) and the finite energy moments employed in Ref. [43]. We note, however, that the strong influence of systematic errors in the derivation of the non-charm R ratio R nc and possible ambiguities in the choice of the upper energy cut-off represent sources of uncertainties for the finite energy sum rule approach that are not easy to quantify. Based on the numbers given in Tabs. 4 and 5 it is straightforward to discuss how the uncertainty in m c (m c ) might be further reduced in the future. One improvement could be achieved by the computation of the order α 3 s corrections to the momentsif they indeed lead to a reduction of the discrepancies between "contour-improved" and "fixed-scale" predictions. At present the theoretical uncertainty coming from the different ways to implement the renormalization scale amounts to about 30 to 35 MeV. The uncertainties coming from the J/ψ e + e − -partial width, from the error in α s and from the treatment of non-perturbative effects are each at the level of 10 MeV. For the combination of these uncertainties a reduction of up to 15 MeV appears possible for an improved measurement of the e + e − -partial width and if smaller uncertainties in the determinations of the strong coupling are achieved. On the other hand, improvements in the data for ψ ′ and for the cc continuum will not improve the situation significantly.
We have also applied our method to estimate the theoretical uncertainties for determinations of m b (m b ) from low-n moments of R bb . Using the approach described above we find m b (m b ) = 4.22 ± 0.11 GeV. The error is only 20 MeV larger than the error obtained from an earlier analysis by one of the authors which was based only the first method with an energy-independent renormalization scale. This shows that the scale choice µ 2 ≃ −sβ 2 has a smaller impact in the bottom quark case due to its larger mass.