Determination of the charm-quark mass in the MS-bar scheme using charm production data from deep inelastic scattering at HERA

We determine the charm-quark mass mc(mc) in the MS-bar scheme using measurements of charm production in deep-inelastic ep scattering at HERA in the kinematic range of photon virtuality from 5 to 1000 GeV squared and Bjorken scaling variable from 0.0001 to 0.5. The extraction of charm quark mass from this process with space-like kinematics provides complementary information to results from hadronic processes. The QCD analysis of the HERA data yields a value of mc(mc) of 1.27 plus minus 0.05(exp) plus 0.06 minus 0.01(scale) GeV at next-to-leading order and mc(mc) of 1.36 plus minus 0.04(exp) plus 0.04 minus 0.00 (scale) plus minus 0.1(theory) GeV at approximate next-to-next-to-leading order. The results are consistent with and of comparable precision as the world average.


Introduction
Quark masses are fundamental parameters of the Standard Model. However, due to confinement no free quarks are observed in nature. Therefore, for the determination of heavy quark masses a careful theoretical description is needed for appropriate observables. This then enables the determination of heavy quark masses by comparing quark mass dependent theoretical predictions to experimental data. In such an analyses reference must be made to the specific scheme used to define the quark mass and it is mandatory to include radiative corrections in Quantum Chromodynamics (QCD) beyond the leading order (LO).
In QCD predictions for hard scattering processes employing the pole (or on-shell) scheme for the definition of the heavy-quark mass, m q is chosen to coincide with the pole of the heavy-quark propagator at each order in perturbative QCD. It is known since long, that the concept of pole masses in QCD has an inherent drawback, see, e.g., [1]. Since quarks are confined inside hadrons, there is no pole in the quark propagator of the full theory. As a consequence, the quark mass parameter in the pole mass scheme is limited to the perturbative domain with corrections of order O(Λ QCD /m q ) and its value depends strongly on the order of perturbation theory. Alternatively, QCD predictions for the charm production cross section in DIS can be considered as a function of the running mass m q (µ r ) at the renormalisation scale µ r in the MS scheme. The quark mass, m q (µ r ), is treated on the same footing as the strong coupling constant α s (µ r ) and represents a prominent example of the so-called short-distance mass definition, which requires the heavy-quark masses to be evaluated at the scale µ r much larger than the QCD scale Λ QCD , i.e., µ r ≫ Λ QCD .
To this date, the charm-quark mass has been measured very accurately from data in electronpositron annihilation, e.g., with the help of QCD sum rules, or from numerical simulations on the lattice, see [2]. These determinations are either entirely non-perturbative (lattice) or related to scattering processes in time-like kinematics only. It is therefore of particular interest to examine the possibility of alternative charm-quark mass determination from heavy-quark production in deepinelastic scattering (DIS). This reaction proceeds in the space-like kinematics with the distinct structure of the QCD amplitudes and in this way it provides an additional check of the QCD parameters' universality.
In this paper the theoretical framework developed in [3] is applied to determine the charmquark mass at next-to-leading order (NLO) and approximate next-to-next-to-leading order (NNLO). The recent results from the H1 experiment [4][5][6] on charm production in deep-inelastic ep scattering at HERA are used in addition to the data sets of [7]. By using the MS scheme the current analysis profits from improved stability of the perturbative expansion and the reduced theoretical uncertainty due to missing higher order contributions estimated by the variations of the renormalisation and factorisation scales. This paper is organised as follows: after a brief reminder of the theoretical basis for the treatment of heavy quarks in DIS we present the determination of the charm-quark mass in the MS scheme, based on recent charm production measurements from H1. Special emphasis is put on the analysis procedure to account for the charm-quark mass dependence of the D * ± production cross section in the fiducial kinematic range. Finally, the running charm-quark mass m c (µ r ) is determined at NLO and at NNLO in QCD. In the NNLO case, the theory predictions are available only to an approximation with a substantial uncertainty arising from missing information in the NNLO DIS Wilson coefficients. The results are compared to previous extractions of m c (m c ) from the DIS data in [3,8] and to the results of different determination methods, used for the world average [2].

Theory framework
Charm-quark production in neutral-current DIS proceeds by scattering of a charged lepton e off a proton P following the reaction in which a virtual boson of space-like momentum Q 2 = −q 2 = −(l − l ′ ) 2 is exchanged. Here l, p, l ′ , p 1 and p 2 denote the four-vectors of the incoming electron and proton and the outgoing electron, charm quark and anti-charm quark, respectively. The charm-quark pair, cc, in the final state is heavy, such that m 2 c ≫ Λ 2 QCD holds. We restrict the momentum transfer to be much smaller than the Z-boson mass Q 2 ≪ M 2 Z . The inclusive cross section is expressed in terms of the standard DIS proton structure functions F k with k = 2, L. The charm contributions to the inclusive structure functions F k are denoted by F cc k (x, Q 2 , m 2 c ). In referring to reaction (1) the structure function F cc k requires by definition a charm-quark pair in the final state, but it is otherwise a completely inclusive quantity, especially with respect to the proton initial state. Thus, the cross section for the reaction (1) reads Here α is the electromagnetic coupling constant and the DIS variables x and y are defined as x = Q 2 /(2p · q) and y = (p · q)/(p · l), respectively. For later use the reduced cross section is introduced as In the standard factorisation approach to perturbative QCD the structure functions F cc k can be written as a convolution of parton distribution functions (PDFs) and Wilson coefficients, see e.g., [9,10], where z max = 1/(1 + 4m 2 c /Q 2 ) and e c = 2/3 is the normalised charm-quark charge. The PDFs for the parton of flavor i are denoted as f i (x, µ 2 f ) and the sum in equation (4) runs over all flavor combinations, i.e., singlet and non-singlet, and the gluon. The Wilson coefficients C k,i depend on the kinematic variables z and ξ, with s denoting the partonic centre-of-mass energy. For the treatment of heavy quarks in DIS the so-called fixed flavour number scheme (FFNS) is chosen with the number of active quarks in the proton n f = 3 and massive charm quarks appearing exclusively in the final state. Moreover, the strong coupling constant is defined in the same scheme as α s (n f ) with n f = 3 for charm quark production. This description of QCD with one massive and n f light quarks can be related to QCD with (n f +1) light quarks by means of the standard matching conditions, cf. e.g., [11,12], and for a discussion of all-order resummations of logarithms in Q 2 /m 2 c in so-called variable flavor number schemes (VFNS) [13], see [14,15]. The C k,i are obtained in perturbative QCD as an expansion in α s = α s (µ r ), here the renormalisation and factorisation scales are set to µ = µ f = µ r . In order to estimate the uncertainties due to missing higher orders, the scale µ is varied by a factor two up and down around the nominal value 2 Conventionally, the Wilson coefficients are presented using the pole mass scheme for m c , see e.g., [9,16]. The necessary formulae for the conversion to the running mass m c (µ r ) in the MS scheme are well known [17][18][19] and the application of the MS scheme to the calculation of F cc 2 has been detailed in [3] (see also [20,21] for related work on heavy-quark hadro-production). For the inclusive cross sections at short distances an appropriate choice for the scale of the running mass m c (µ r ) is µ r = m c . The renormalisation group evolution for the scale dependence governed by the corresponding quark mass anomalous dimension [22,23] converges even for scales as low as the charm-quark mass, that is O(1) GeV.
The Wilson coefficients in equation (6) have been computed exactly to NLO and the functions c (1) k,i [9] are often used via the parameterisations of [10], see also [24]. At NNLO approximate results for the most important gluon and quark coefficient functions, c (2) 2,g and c (2) 2,q , are known [16] and denoted by NNLO approx in the following. These are based on recent partial NNLO improvements for the structure function F cc 2 which encompass various kinematic limits: The threshold approximation c (2) thr 2,g has been determined to the next-to-next-to-leading logarithm (NNLL) (see [25,26] for previous approximations at the level of the next-to-leading logarithm (NLL)). The function c (2) thr 2,q is consistent with zero to the accuracy considered. Likewise, fully analytic results for c (2) asm 2,i , i = g, q, in the asymptotic regime of Q 2 ≫ m 2 c have been obtained. The corresponding calculations make use of the formalism of [11,12] and a number of lowest even-integer Mellin moments for the necessary heavy-quark operator matrix elements at three loops [27][28][29], see also [30]. In the high-energy limit the expression for c (2) small−x 2,i is exact to leading-logarithmic (LL) accuracy at small-x due to [31] but is approximate only at NLL.
The combination of all available information leads to an expression for the NNLO Wilson coefficients, c (2) 2,i of the form: with a suitable matching function f (ξ) connecting the regions Q 2 ≫ m 2 c and s ≫ m 2 c , cf. equation (4.9) in [16].
The approach chosen has one caveat, though, which is a non-negligible theoretical uncertainty due to the poorly constrained NLL term at small-x in c (2) small−x 2,i and the limited number of known three-loop Mellin moments [27][28][29]. To account for these deficits two different scenarios c  [16]. An additional theoretical uncertainty on F cc 2 is estimated by the variation of the renormalisation and factorisation scales. This part of the theoretical uncertainty is obtained using exact results since the µ-dependence is fully known at NNLO [3,25] from the renormalisation group, see also [16].

QCD analysis of charm production measurements
To study the impact of the charm production measurements at HERA on the determination of the MS charm mass a variant of the ABM11 fit [7] is performed with the data of [4][5][6] included. The data were collected by the H1 experiment during the HERA II running period corresponding to an integrated luminosity of about 350 pb −1 and cover the kinematic range of photon virtualities 5 GeV 2 < Q 2 < 2000 GeV 2 . The charm mass m c (m c ) is a free parameter of the fit comparing equation (3) to the experimental data. Experimentally charm-quark production in DIS is tagged via fully reconstructed charmed mesons or by using secondary vertex information of tracks, exploiting the longevity and the large mass of charmed hadrons. These measurements are restricted (e.g. in the transverse momenta or angles of the produced particles) by the acceptance of the detector. The phase space, in which charmed hadrons can be fully reconstructed is usually referred to as the visible or fiducial phase space. Corrections for the non-measurable phase space and the fragmentation of charm quarks to charmed hadrons are applied in the fit. In the following, the measurements of charm production used as input and their treatment in the QCD analysis are described with particular emphasis on the corrections for non-measurable phase space.
The c-quark production in [4,5] is tagged via fully reconstructed D * ± -mesons in the decay mode D * ± → (D 0 → K ∓ π ± )π ± . These measurements are restricted to the visible kinematic range of the D * ± -meson's transverse momentum p T (D * ) > 1.25 GeV and pseudo-rapidities |η(D * )| < 1.8 at medium virtualities Q 2 < 100 GeV 2 . At high virtualities 100 GeV 2 < Q 2 < 1000 GeV 2 the visible range of p T (D * ) > 1.5 GeV and |η(D * )| < 1.5 is covered. Due to the phase space limitations equation (2) cannot be used directly in the analysis of these data. Therefore, for these measurements the invisible phase space region is accounted for in equation (2) through a factor where σ vis (D * ± ) and σ cc full are the QCD predictions for the D * ± meson cross section in the visible phase space and the charm production cross sections in the full phase space, respectively. These predictions are calculated at NLO in the FFNS with the fully exclusive program HVQDIS [24]. The contribution from b-quarks to the inclusive D * ± -meson production cross section, reported in [4,5], is estimated using HVQDIS and is subtracted. The D * ± cross sections are re-calculated using the recent branching ratio values [2].
The value of ε vis depends strongly on m c and its dependence on m c is taken into account in the fit. For this purpose ε vis is calculated in a first step for selected values of the charm quark mass, which correspond to the range of m c (m c ) scanned in our fit: obtained by mapping of the running-mass grid in equation (10) where a = α s (m c (m c ))/π. In the calculation of ε vis , the proton structure is described by special PDF sets, provided for this analysis. They correspond to the NLO variant of ABM11 fit [7] and are performed in the MS definition of charm quark mass. The PDFs depend on the value of m c through substantial correlations between α s , m c and the gluon PDF. In order to provide a fully consistent treatment of the charm-mass effects in ε vis one has to take into account this dependence. Therefore we use the PDFs, which exactly correspond to the current value of m c appearing in the fit. These PDFs are given by interpolation between those obtained in the variants of ABM fit with the m c settings of equation (10). The full dependence of ε vis on m c is provided by parabolic interpolation of the HVQDIS results obtained at the values of m c in equation (11) with the interpolation coefficients P, calculated independently for the LO and NLO terms in σ cc vis/full , Such a representation allows the determination of ε vis in terms of the MS mass. This is achieved by substituting m pole c in equations (13) and (14) with the matching condition of equation (12) similarly to the approach used earlier to derive the heavy-quark electro-production coefficient functions in the MS mass definition [3]. The terms of O(α 3 s ) and higher appearing after substitution exceed the NLO accuracy therefore they are dropped and the final expressions for σ cc vis/full employed in our analysis read for i = 0, 1, 2.
The cross section predictions σ vis (D * ± ) depend not only on the kinematics of the charm quark production mechanism but also on the fragmentation of the charm quark into D * ± mesons. The charm quark fragmentation function to D * ± mesons has been measured by H1 [32] using the production of D * -mesons with and without associated jets in DIS. In the calculation of σ vis (D * ± ) the longitudinal fragmentation is performed in the γ * − p centre-of-mass frame, using the fragmentation function of Kartvelishvili et al. [33] which is controlled by a single parameter, α K . This parameter has been determined for two different regions of the partonic centre-of-mass energy squared, s, depending on the jet requirements made in the different analyses. The fragmentation is observed to become softer with increasing s as expected from the perturbative evolution of the fragmentation function not implemented in HVQDIS. This is accounted for in the current analysis by using different values of α K corresponding to the measurements for two ranges of s, as listed in table 1. The limits on the ranges in s are determined with HVQDIS by applying the jet requirements of the individual analysis on parton level. The α K parameters and the s ranges are varied according to their uncertainties to evaluate the corresponding uncertainty on σ cc vis . The charmed hadrons also s range [GeV 2 ] α K measurement s < (70 ± 40) 6.1 ± 0.9 [32]  receive a transverse momentum,k T , with respect to the charm quark direction according to The transverse momentum averagek T is chosen as 0.35 ± 0.15 GeV 2 , in line with the experimental results on hadron production in e + e − collisions [34][35][36][37][38][39]. A fragmentation fraction of charm quarks into D * ± mesons of f (c → D * + ) = 0.2287 ± 0.0056, is used as determined by averaging the e + e − and ep results [40]. In order to evaluate the fragmentation model uncertainty on σ vis (D * ± )) the parametersk T , α K , and f (c → D * + ) are varied within the uncertainties quoted above and each variation is considered as a correlated uncertainty in ε vis . In a few cases, a symmetric variation of the model parameters results in an asymmetric uncertainty on the cross section. In such cases, the largest absolute value of the uncertainty is assigned. The dominant systematic uncertainty is arising from the variation of the fragmentation function parameter α K . The measurement [6] used in the fit is based on determination of the vertex displacement of the tracks and covers essentially the full phase space. However, the extrapolation to the cross sections of charm and beauty quark production is performed with a Monte-Carlo simulation and the determination of a correction at NLO as in case of D * ± analysis is not possible. We assume this correction to be small and no extrapolation factor is applied to these data.

Results
The NLO variant of our analysis applies O(α 2 s ) corrections to the heavy quark electro-production Wilson coefficients [9] and the NLO ABM11 PDFs. For the NNLO variant of our analysis we use the NNLO ABM11 PDFs and apply the corrections up to O(α 3 s ) as discussed in section 2. Since the NNLO Wilson coefficients c (2) 2,i in equation (8) are still approximate and affected by a residual uncertainty parameterised by c (2),A 2,i and c (2),B 2,i we define a particular shape for c (2) 2,i as a linear combination of the two envelopes with a parameter d N interpolating between these two options, The consistency of the prediction using results of the fit with the data of [4,5] on σ cc,VIS red = σ vis (D * ± )/ f (c → D * ± ) and of [6] on σ cc red is illustrated in figures 1 and 2, respectively. In the case of NLO, the fit quality is very good for each data set with the total value of χ 2 = 60 for the number of data points (NDP) equal to 60. The consistency of the NNLO prediction using the fit results with the data is equally good as for the NLO variant of the fit (cf. figures 1 and 2), with the best fit value of χ 2 /NDP = 63/60 being achieved for d N = −0.6. The variant of the NNLO Wilson coefficient given by c The experimental uncertainty is obtained from the propagation of all uncertainties in the data with account of their correlations. The fragmentation model uncertainties relevant for the data [4,5] are also included on the same footing. To estimate the influence of the PDF uncertainty on the precision of the m c (m c ) determination, the NLO analysis is repeated for each of the ABM11 PDF set members, representing the 1σ uncertainty in the fitted PDF parameters 3 . The differences between the values m c (m c ) obtained in this way and the ones obtained with the central PDF member are added in quadratures. The resulting uncertainty on m c (m c ) is about 10 MeV which is very small in comparison to the other uncertainties. As the fit is performed at fixed order in perturbation theory, one needs to account for missing contributions beyond the order considered. This uncertainty is calculated in the standard manner from the variation of the renormalisation and factorisation scales within the range given in equation (7). The scale variation is performed simultaneously in equation (3) and in the calculation of the extrapolation factors ε vis . For the NNLO case the parameter d N is fixed at the value of −0.6, which is preferred by the fit. Due to the sensitivity of ε vis on the scale variation, the resulting uncertainty on m c (m c ) does not improve significantly from the NLO prediction to the NNLO one. In case of the NNLO result, the theory error on m c (m c ) of 100 MeV due to lacking knowledge on the NNLO Wilson coefficient is significantly larger than the uncertainty due to scale variation. The necessary theory computations to remedy this unsatisfactory situation are discussed in [16].
The NNLO predictions for charm electro-production obtained for the two variants of the fit are shown in figure 3 for the kinematics of the HERA collider experiments. The two variants of the fit are not very different at large Q 2 which illustrates the moderate potential of the data [4][5][6] to separate between different shapes of the NNLO term in the heavy-quark Wilson coefficients. However, at small Q 2 the predictions diverge therefore the combined H1 and ZEUS data may provide an additional constraint on this term and can improve the accuracy of the NNLO determination of m c .
In the previous analysis [3,8], values of the charm quark MS mass of m c (m c ) = 1.26 ± 0.09(exp) ± 0.11(th) at NLO and m c (m c ) = 1.01 ± 0.09(exp) ± 0.03(th) at NNLO approx was obtained. While the NLO results of the present and the previous analyses are consistent, the present result at NNLO approx is significantly larger. This shift in the central value is due to several sources. First, in the present study we use a significantly improved theory of heavy-quark production in DIS. As mentioned above, we reconstruct the three-loop Wilson coefficient c (2) 2,i based on all available information in various kinematic limits [16], in contrast to earlier studies [3,8] which relied only on the threshold approximation c (2) thr 2,g . Also, the analysis of [3,8] was performed as a variant of the ABKM09 global analysis [14] of world DIS data, to which the HERA experiments, at that time, only contributed the inclusive data of [41,42]. The present determination of m c (m c ) is based on the ABM11 fit [7], which incorporates the new combination of the HERA run I inclusive DIS data from the H1 and ZEUS experiments [43] and to which we add the new precise data sets on DIS charm production [4][5][6] [47] inclusive spectra in semileptonic B-decays 1.28 ± 0.04 BLOSSIER [48] lattice QCD (n f = 2); hadron spectrum 1.273 ± 0.006 MCNEILE [49] lattice QCD (n f = 2+1); pseudo-scalar current correlator 1.279 ± 0.013 CHETYRKIN [50] e + e − → cc cross-section and QCD sum rules.
1.224 ± 0.017 ± 0.054 HOANG [54] global fit to inclusive B-decay data. QCD predictions have been applied, the extractions of m c (m c ) listed in table 2 are limited to processes with time-like kinematics, i.e., cross-section data from e + e − -collisions exposed to QCD sum rule analyses or data on B-decays. The latter method carries also an intrinsic dependence on the uncertainty of the b-quark mass, while the former, i.e., the QCD sum rule analyses fix the value of α s (M Z ) to the world average including the associated very small error 5 . It has been shown [55], that the systematic shift of m c (m c ) due to the value of α s (M Z ) in QCD sum rule analyses is quite sizable.

Conclusions
The high precision measurements of charm quark production in DIS at HERA offer an attractive possibility to extract the charm-quark mass in the theoretically well-founded MS scheme. The resulting experimental precision of the present determination based on DIS data is substantially improved with respect to previous analyses. It is now compatible with the theoretical uncertainty when comparing to the QCD predictions at NLO and it is significantly smaller than in case of the approximate NNLO QCD predictions. The latter suffer from missing information on the threeloop Wilson coefficients at small-x and small values of Q 2 and imply an additional theoretical uncertainty on m c of 100 MeV. In comparison to the measurements used to define the world average up to now, our results add complementary information from scattering processes with space-like kinematics and provide an important test of the QCD dynamics. The kinematic range covered by the DIS data allows for the extraction of m c in the MS scheme well within the regime of validity of perturbative QCD. This analysis accounts for the full correlation of the dependence on m c with other non-perturbative parameters, most prominently the value of strong coupling constant α s (M Z ) and the gluon PDF. Future improvements on the accuracy of the m c (m c ) extractions from DIS data rely mostly on the theoretical progress for the three-loop Wilson coefficients especially at small x. Also the experimental uncertainty can still be reduced by combining the charm production data from the H1 and the ZEUS collaborations.
In summary, the charm-quark mass extractions from DIS data may challenge the accuracy of QCD sum rules analyses once the uncertainties on all non-perturbative parameters are accounted for on equal footing.