A determination of mc(mc) from HERA data using a matched heavy-flavor scheme

The charm quark mass is one of the fundamental parameters of the Standard Model Lagrangian. In this work we present a determination of the MSbar charm mass from a fit to the inclusive and charm HERA deep-inelastic structure function data. The analysis is performed within the xFitter framework, with structure functions computed in the FONLL general-mass scheme as implemented in APFEL. In the case of the FONLL-C scheme, we obtain mc(mc) = 1.335 +- 0.043(exp) +0.019 -0.000(param) +0.011 -0.008(mod) +0.033 -0.008(th) GeV. We also perform an analogous determination in the fixed-flavor-number scheme at next-to-leading order, finding mc(mc) = 1.318 +- 0.054(exp) +0.011 -0.010(param) +0.015 -0.019(mod) +0.045 -0.004(th) GeV, compatible with the FONLL-C value. Our results are consistent with previous determinations from DIS data as well as with the PDG world average.


Introduction
The masses of the heavy quarks, charm, bottom and top, are fundamental parameters of the Standard Model [1]. A precise determination of their values is of utmost importance; as an example, the fate of the electroweak vacuum depends crucially on the exact value of m t [2]. In the case of the charm quark, since its mass m c is larger than the scale Λ QCD of Quantum Chromodynamics (QCD), its value is a direct input of many perturbative calculations involving charm quarks in the initial and/or in the final state.
Differences in the value of the charm quark mass and in the treatment of its effects in deep-inelastic-scattering structure functions can lead to differences in modern analyses of parton distribution functions (PDFs) [3][4][5][6][7], with implications for precision phenomenology at the Large Hadron Collider (LHC). As a consequence, a high-precision determination of the charm quark mass is of interest both in principle, as a fundamental test of the Standard Model and a measurement of one of its fundamental parameters, and in practice, as input for LHC calculations.
The current global-average value of the charm mass in the MS renormalization scheme is m c (µ R = m c ) = 1.275 ± 0.025 GeV [8], where the result is dominated by high-precision data from charm production in e + e − collisions. It is therefore interesting to provide alternative determinations of the charm mass from other processes, both to test the robustness of the global average and to attempt to further reduce the present uncertainty.
A process directly sensitive to the charm mass is open-charm production in leptonproton deep-inelastic scattering (DIS). This process has been measured with high accuracy at the HERA collider and the results of different measurements implying various charmtagging techniques are combined [9]. The charm contribution to the inclusive structure functions can be determined through the measurement of the charm-pair production cross section. In addition, the final combination of inclusive measurements from Runs I and II at HERA has been recently presented in [7].
DIS structure functions can be described using a variety of theoretical schemes, including the fixed-flavor number (FFN) scheme, where charm mass effects are included to a fixed perturbative order, the zero-mass variable-flavor number (ZM-VFN) scheme that neglects power-suppressed terms in the charm mass but resums to all orders large collinear logarithms, and the so-called matched general-mass variable-flavor-number (GM-VFN) schemes, which interpolate smoothly between the two regimes. A recent discussion and summary of the application of these schemes to heavy-flavor data at HERA can be found e.g. in [10].
The original formulation of the FONLL general-mass scheme for DIS structure functions was derived in the pole (on-shell) heavy quark scheme [11]. In Ref. [26] it was shown how DIS structure functions in the FFN scheme can be modified to include MS heavyquark masses. The same scheme conversion can be applied to any GM-VFN scheme, and in this work we provide the relevant expressions for FONLL structure functions with MS running masses. The main advantage of the use of MS masses is the possibility of direct connection with the precise determinations from low-energy experimental data [8].
In this work we will use the xFitter open-source framework [27] (previously known as HERAfitter) to extract the MS charm mass from a PDF fit to the most up-to-date inclusive and charm data from HERA. Structure functions are computed using the FONLL scheme as implemented in the APFEL [28] code. Our results have been obtained employing the most accurate perturbative calculations presently available and will include a detailed characterization of the different sources of uncertainties on m c (m c ) from data, theory and fitting methodology. As we will show, the results are consistent with the global PDG average as well as with previous determinations based on the FFN [9,[23][24][25] and in the S-ACOT [29] schemes. 1 The uncertainty in our results turns out to be competitive with that of previous determinations based on DIS structure functions.
The outline of this paper is the following. In Sect. 2 we discuss how FONLL can be formulated in terms of MS masses and present a benchmark of its implementation in APFEL. In Sect. 3 we describe the settings of the PDF fits and the treatment of the uncertainties. Results for the determination of m c (m c ) are presented in Sect. 4, where we also compare with previous determinations. We conclude and discuss possible next steps in Sect. 5.

FONLL with MS heavy-quark masses
In this section we discuss how the FONLL general-mass variable-flavor-number scheme for DIS structure functions can be expressed in terms of MS heavy-quark masses. We also describe the subsequent implementation in the public code APFEL, and present a number of benchmark comparisons with other public codes.
In general, higher-order calculations are affected by ambiguities in the prediction for the physical quantities due to the choice of the subtraction scheme used to remove divergences. In fact, different prescriptions imply different numerical values of the parameters of the underlying theory.
As far as the mass parameters are concerned, the pole mass definition is usually more common in the calculation of massive higher-order QCD corrections to heavy-quark production processes. The main reason for this is that the pole mass is, by its own definition, more closely connected to what is measured in the experiments. On the other hand, it is well known that observables expressed in terms of the pole mass present a slow perturbative convergence. This is caused by the fact that the pole mass definition suffers from non-perturbative effects which result in an intrinsic uncertainty of order Λ QCD [31]. The MS scheme, which stands for modified minimal subtraction scheme, is instead free of such ambiguities and as a matter of fact massive computations expressed in terms of heavyquark masses normalized in this scheme present a better perturbative convergence [26]. As a consequence, the results obtained in the MS scheme are more appropriate to achieve a reliable determination of the numerical value of the charm mass.
The FONLL scheme, as any other GM-VFN scheme, aims at improving the accuracy of fixed-order calculations at high scales by matching them to resummed computations. In DIS this results in the combination of massive (fixed-order) calculations, that are more reliable at scales closer to the heavy-quark masses, with resummed calculations that are instead more accurate at scales much larger than the heavy-quark masses. However, in the original derivation, the massive component of the FONLL scheme was expressed in terms of the pole masses [11].
It is then one of the goals of this paper to provide a full formulation of the FONLL scheme applied to DIS structure functions in terms of MS masses. A detailed discussion on such a formulation is given below in Sect. 2.1. Here, we limit ourselves to describing the main steps needed.
The generic form of the DIS structure functions in the FONLL approach applied to charm production is: where x, Q, and m c are the Bjorken variable, the virtuality of the photon, and the mass of the charm quark, respectively. In eq. (2.1) the three-flavor structure function F (3) is evaluated retaining the full charm-mass dependence and with no charm in the initial state. The four-flavor structure function F (4) is instead computed by setting m c to zero and allowing for charm in the initial state, and its associated PDF reabsorbs the mass (collinear) divergences which are in turn resummed by means of the DGLAP evolution. Finally, F (3,0) represents the massless limit of F (3) where all the massive power corrections are set to zero and only the logarithmically enhanced terms are retained. This last term is meant to subtract the double counting terms resulting from the sum of F (3) and F (4) . In fact, the role of F (3,0) is twofold: for Q m c , by definition F (3) and F (3,0) tend to the same value so that the FONLL structure function reduces to F (4) . By contrast, in the region where Q m c it can be shown that F (d) becomes subleading in α s reducing the FONLL structure function to F (3) up to terms beyond the nominal perturbative accuracy.
It should be noticed that, even though F (d) in eq. (2.1) becomes subleading in the low-energy region, it might become numerically relevant and it is advisable to suppress it. To this end, the term F (d) in eq. (2.1) is usually replaced by: where the function D(Q, m c ) is usually referred to as the damping factor and has the explicit form: 3) The role of the damping factor is clearly that of setting F (d ) to zero for Q < m c , suppressing it for Q m c , and reducing it to F (d) for Q m c . It should be pointed out that the particular functional form of the damping factor given in eq. (2.3) is somewhat arbitrary. In fact, any function D such that F (d ) and F (d) only differ by power-suppressed terms, namely: is a formally suitable choice. In the results section we will also consider the effect of varying the functional form of the damping factor in order to estimate the associated theoretical uncertainty on m c (m c ). Given the possible different perturbative structure of the elements that compose the FONLL structure function in eq. (2.1), two possibilities for the definition of the perturbative ordering are possible: the relative and the absolute definitions. In the relative definition F (4) and F (3) are combined using the same relative perturbative accuracy, that is LO with LO, NLO with NLO, and so on. The absolute definition, instead, is such that LO refers to O(α 0 s ) (parton model), NLO to O(α s ), and so forth. This issue is relevant in the neutralcurrent case where the lowest non-vanishing order is O(α 0 s ) for F (4) and O(α s ) for F (3)2 such that the relative and absolute orderings lead to different prescriptions.
Beyond LO, there are currently three possible variants of the FONLL scheme, all of them implemented in APFEL: • the FONLL-A variant adopts the absolute ordering at O(α s ) and thus only terms up to this accuracy are included. This variant is formally NLO and thus also PDFs should be evolved using the same accuracy in the DGLAP evolution.
• The FONLL-B variant is instead computed using the relative ordering at NLO. Therefore, F (4) is computed at O(α s ) and combined with F (3) at O(α 2 s ). F (3,0) is instead computed dropping the non-logarithmic O(α 2 s ) term to match the accuracy of F (4) in the low-energy region. PDFs are again evolved at NLO.
• Finally, the FONLL-C scheme adopts the absolute ordering at O(α 2 s ). This is formally a NNLO scheme thus PDFs should be evolved using the same accuracy.
Presently, no other variant beyond FONLL-C can be pursued because the O(α 3 s ) massive coefficient functions are not known yet. Approximate NNLO corrections valid near the partonic threshold, in the high-energy (small-x) limit, and at high scales Q 2 m 2 have been derived in Ref. [34] and they are currently employed by the ABM group to determine NNLO PDFs [6].
As clear from the description above, the computations for the three-flavor structure functions F (3) and F (3,0) depend explicitly on the charm mass, while F (4) does not. In addition, as already mentioned, the expressions needed to compute F (3) and F (3,0) are usually given in terms of the pole mass. As a consequence, one of the steps required to achieve a full formulation of the FONLL structure functions in terms of MS masses is the adaptation of the heavy-flavor contributions to the structure functions. A thorough explanation of the procedure adopted to perform such transformation can be found in Ref. [26] for both neutral-and charged-current structure functions. In Sect. 2.1 we re-derive the main formulae and report the full expressions for the relevant coefficient functions. It should be pointed out that the derivation presented in Ref. [26] is performed assuming µ R = m c (m c ), µ R being the renormalisation scale, and the renormalisation scale dependence of α s is restored only at the end using the expansion of the solution of the relative RG equation. Such a procedure implies that the heavy-quark mass is not subject to the relative RG equation: in other words, the mass running is not expressed explicitly. The reason is that in the running of the heavy-quark mass in MS one can resum logarithms of µ R /m c (m c ) and this is not required in a fixed-order calculation. On the contrary, when dealing with a GM-VFN scheme like FONLL, such a resummation is an important ingredient and thus should be consistently incorporated into the derivation. For this reason, the transition from pole to MS masses of the massive structure functions presented in Sect. 2.1 is done at the generic renormalisation scale µ R and the connection between m c (m c ) and m c (µ R ) is established solving the appropriate RG equation.
A further complication that arises in FONLL as a VFN scheme is the fact that the involved running quantities, that is PDFs, α s and the mass itself, have to be properly matched when crossing a heavy-quark threshold in their evolution. The matching conditions for PDFs and α s are presently known up to O(α 2 s ) [35] and O(α 3 s ) [36], respectively, but those for PDFs are given in terms of the pole mass. In the next section we will show how to express them in terms of the MS mass up to the relevant accuracy. As far as the matching of the mass is concerned, the expressions for the matching conditions are given in Ref. [37] up to O(α 3 s ) also in terms of MS mass.

Implementation
In this section we will describe in some detail the implementation of the FONLL scheme in terms of the MS heavy-quark masses in APFEL. Starting from the more usual definition of structure functions in terms of pole masses, our goal is to consistently replace them with the MS mass definition.

MS mass vs. pole mass
The (scale independent) pole mass M and the (scale dependent) MS mass m(µ) arise from two different renormalization procedures and, as already mentioned, in perturbation theory they can be expressed one in terms of the other. The relation connecting pole and MS mass definitions has been computed in Ref. [31] up to four loops. However, in the following we will only need to go up to one loop and thus we report here the corresponding relation: with: where C F = 4/3 is one of the usual QCD color factors. Moreover, we have defined: and: In the following we will use eq. (2.5) to replace the pole mass M with the MS mass m(µ).

Solution of the RGE for the running of the MS mass
In order to evaluate the running of m(µ) with the renormalization scale µ we have to solve the corresponding renormalization-group equation (RGE): whose first three coefficients can be taken from Ref. [38] 3 : where N f is the number of active flavors. In addition, the RGE for the running of α s reads: with: Combining eqs. (2.9) and (2.11) we obtain the following differential equation: whose solution is: In order to get an analytical expression out of eq. (2.14), one can expand the integrand in the r.h.s. using the perturbative expansions of γ m (a s ) and β(a s ) given in eqs. (2.9) and (2.11). This allows us to solve the integral analytically, obtaining: where we have defined: and a ≡ a s (µ) and a 0 ≡ a s (µ 0 ). Eq. (2.15) represents the NNLO solution of the RGE for the MS mass m(µ).
Of course, the NLO and the LO solutions can be easily extracted from eq. (2.15) just by disregarding the terms proportional to a 2 and a 2 0 for the NLO solution and also the terms proportional to a and a 0 for the LO solution 4 . 3 The following expressions have been adjusted taking into account our definition of as which differs by a factor of 4 with respect to that of Ref. [38]. 4 In order to be consistent, the evaluation of a and a0 eq. (2.15) must be performed at the same perturbative order of m(µ). So, for instance, if one wants to evaluate the NNLO running of m(µ) also the value of a and a0 must be computed using the NNLO running.

Matching conditions
When working in the context of a VFN scheme, all running quantities are often required to cross heavy-quark thresholds when evolving from one scale to another. Such a transition in turn requires the matching different factorization schemes whose content of active flavors differs by one unit. In other words, if the perturbative evolution leads from an energy region where (by definition) there are N f − 1 active flavors to another region where there are N f active flavors, the two regions must be consistently connected and such a connection can be evaluated perturbatively. This goes under the name of matching conditions. In general, matching conditions give rise to discontinuities of the running quantities at the matching scales and in the following we will report the matching conditions up to NNLO in terms of the MS heavy-quark thresholds for: α s (µ), m(µ) and PDFs.
Matching of α s (µ) The matching conditions for α s were evaluated in Ref. [36] to three loops. We report here the relation up to two loops (again taking into account the factor 4 coming from the different definitions of a): M being the pole mass of the n-th flavor. From eq. (2.5) we can easily infer that: (2.18) Therefore, it is straightforward to see that: so that: consistently with eq. (20) of Ref. [37].
In order to simplify this expression, it is a common procedure to perform the matching at the point where the logarithms vanish. In this particular case, choosing µ = m(µ) = m(m), we get: which can be easily inverted obtaining: It is interesting to observe that, in order to perform the matching as described above, one just needs to know the value of m(m). This is the so-called RG-invariant MS mass.

Matching of m(µ)
The running of the MS masses also needs to be matched at the heavy-quark thresholds. In particular, one needs to match the (N f − 1)-with (N f )-scheme for the mass m q (µ), with q = c, b, t, at the threshold m h (µ), where h = c, b, t. From Ref. [37] we read: where: .

(2.24)
Exactly as before, if we choose to match the two schemes at the scale µ = m h (µ) = m h (m h ), the logarithmic terms vanish and we are left with:

Matching of PDFs
To conclude the section on the matching conditions, we finally consider PDFs. One can write the singlet and the gluon in the (N f )-scheme in terms of singlet and gluon in the (N f − 1)-scheme at any scale µ as follows: where the form of the functions entering the transformation matrix above are given in Appendix B of Ref. [39] in terms of the pole mass. We omit the matching conditions for the non-singlet PDF combinations because they have no O(a s ) correction and the first correction appears at O(a 2 s ). This leaves the conversion from the pole to the MS mass scheme unaffected up to NNLO.
In order to replace the pole mass M with the MS mass m(µ), we just have to plug eq. (2.19) into eq. (2.27). In doing so, only the O(a s ) terms proportional to ln(µ 2 /M 2 ) play a role in the conversion up to NNLO. Since the functionsÃ S,(1) hg and A S, (1) gg,h can be written as:Ã where: replacing M with m in eq. (2.28) using eq. (2.19) leads to: (2.30) Therefore eq. (2.27) in terms of m becomes: (2.31) As usual, we choose to match the (N f )-scheme to the (N f − 1)-scheme at µ = m(µ) = m(m) so that all the logarithmic terms vanish, obtaining:

Renormalization scale variation
The scale µ that appears in a s and m q is the renormalization scale, which we will now denote as µ R . The scale that explicitly appears in the PDFs is instead the factorization scale, which we will now denote with µ F . In principle, renormalization and factorization scales are different but one usually takes them to be proportional to each other, as µ R = κµ F , where κ can be any real number 5 . The most common choice when matching the (N f − 1)-scheme to the (N f )-scheme is to set µ F equal to heavy-quark thresholds (M c , M b and M t in the pole-mass scheme and m c (m c ), m b (m b ) and m t (m t ) in the MS scheme). In doing so, the logarithmic terms in the PDF matching conditions are assured to vanish. However, if κ is different from one, the logarithmic terms in the matching conditions for a s (µ R ) and m q (µ R ) do not vanish anymore. In the following we will show how the matching conditions for a s and m q change for κ = 1.
Let us start with α s . Inverting eq. (2.20) we obtain: where: Setting µ F = κµ F , we have that: It should be noticed that in the case κ = 1 PDFs acquire an implicit dependence on µR that comes from a redefinition of the splitting functions that in turn derives from the expansion of αs(µR) around µR = µF .
As usual, the matching scale is chosen to be µ F = m(m), so that: But using eq. (2.14), it is easy to see that: so that: It should be noticed that in the eq. (2.38), since a Therefore, setting µ = µ R = κm(m) = κm in eq. (2.20) and using eq. (2.38), one gets: (2.39) whose inverse is: (2.40) Now let us turn to m q . In this case there is not much to do. In fact, for an arbitrary matching point the matching condition of the MS mass starts at O(α 2 s ) (cfr. eq. (2.23)), therefore writing L µm in terms of ln κ would give rise to subleading terms up to NNLO (see eq. (2.38)). As a consequence, we have that: whose inverse is:

Structure functions
We finally turn to discuss how the DIS massive structure functions change when expressing them in terms of the MS masses. We will first consider the neutral-current (NC) massive structure functions up to O(α 2 s ), which is the highest perturbative order at which corrections are known exactly, and then we will consider the charged-current (CC) massive structure functions again up to the highest perturbative order exactly known 6 , that is O(α s ). In order to shorten the notation, we will adopt the following definitions:

Neutral current
Dropping all the unnecessary dependences, the NC massive structure functions up to O(a 2 s ) have the form: The goal is to replace explicitly the pole mass M with the MS mass m using eq. (2.5). To this end, following the procedure adopted in Ref. [26], we expand F (0) (M ) and F (1) (M ) around M = m: Finally, we have that: We now need to evaluate explicitly the derivative in eq. (2.46). First of all we observe that: where g is the gluon distribution and we have used the following definitions: (2.48) Defining: the derivative of eq. (2.47) can be written as: where G(z, M ) is the primitive of G(z, M ) with respect to z (i.e. ∂ G/∂z = G). But: It can be shown that the boundary term in eq. (2.52) vanishes (see Ref. [26]), thus it can be omitted. Gathering all pieces and taking into account that: we have that: Finally, considering that: and using eqs. (2.46) and (2.54), one can explicitly write down the full structure of the massive structure functions (F 2 and F L ) in terms of MS masses up to O(α 2 s ) as follows: (2.56) In order to carry out the implementation, we need to evaluate explicitly the derivative of C (0) g in eq. (2.56) and this must be done separately for F 2 and F L . We consider F 2 first. The explicit expression of C (0) 2,g is the following: where: From the definitions in eq. (2.60), we obtain: (2.61) Therefore: To find the explicit expression, we just need to evaluate the derivative of I q and J q starting from eqs. (2.58) and (2.59) which is easily done: (2.63) In the end we get: (2.64) The implementation of the FONLL scheme given in eq. (2.1) requires the massless limit of the massive structure functions. In practice this means that one needs to compute the limit M → 0 of the massive coefficient functions retaining the logarithmic enhanced terms. In order to apply this recipe to eq. (2.64), we observe that: and that: We now turn to consider F L . In this case the the gluon coefficient function takes the simpler form: (2.68) Therefore, using eq. (2.61), we immediately get: It is finally easy to realize that:

Charged current
In this section we consider the CC massive structure functions. The treatment follows the exact same steps as the NC structure functions, with the only difference being that in the CC case the first non-vanishing term is O(a 0 s ). This means that, truncating the perturbative expansion at O(a s ), we have: with k = 2, 3, L. Therefore, expanding F (0) and F (1) around M = m and keeping only the terms up to O(a s ), one obtains: The leading-order contribution can be written as follows: where: where we have also defined: (2.75) Therefore: that can be conveniently rewritten as: so that, using eq. (2.74), we have that: (2.78) Finally, we notice that in the massless limit, where λ → 1, all expressions in eq. (2.78) vanish, with the consequence that the CC massive structure functions up to O(a s ) in terms of the pole mass M or the MS mass m are exactly the same.

Benchmark
In order to validate the implementation in APFEL, we have benchmarked it against public codes. To the best of our knowledge, there exist no public codes able to compute structure functions in the FONLL scheme with MS masses. For this reason the best we could do is to benchmark the various ingredients separately.
As a first step, we present the benchmark of the running of PDFs, α s and m c 7 in the VFN scheme with MS heavy-quark thresholds. The difference with respect to the more common pole-mass formulation arises from the fact that the matching of the evolutions at the heavy-quark thresholds needs to be adapted to take into account the different scheme used to renormalize the masses. The full set of such matching conditions for PDFs, α s and m c has been collected in Sect. 2.1.
We start with the DGLAP PDF evolution in the VFN scheme with MS heavy-quark thresholds. A careful benchmark was already presented in the original APFEL publication.  Figure 1: Comparison between APFEL v2.7.0 and HOPPET v1.1.5 for the VFNS DGLAP evolution at NNLO with MS heavy-quark thresholds. The evolution settings, i.e. initial scale PDFs, reference value of α s , and heavy-quark thresholds, are the same as used in the Les Houches PDF evolution benchmark [42]. The upper inset shows the gluon PDF xg, the valence up and down PDFs xu v ≡ xu − xu and xd v ≡ xd − xd, respectively, and the total strangeness xs + ≡ xs + xs at µ F = 100 GeV as functions of the Bjorken variable x as returned by APFEL. In the lower inset the ratio to HOPPET is displayed showing a relative difference of 10 −4 or better all over the considered range.
In particular, the APFEL evolution has been checked against the HOPPET code [41] v1.1.5, finding a very good agreement at the O 10 −4 level or better. Since then, APFEL has undergone several changes and improvements and thus we repeated the benchmark using the same settings and finding the same level of agreement with HOPPET, as shown in Fig. 1 for a representative set of combinations of PDFs 8 .
Although the benchmark of the DGLAP evolution already provides an indirect check of the evolution of α s , we have also performed a direct check of the VFNS evolution with MS heavy-quark thresholds of α s along with the evolution of the MS charm mass. To this end, we have used the CRunDec code [43], which is the C++ version of the Mathematica package RunDec [37]. In Fig. 2 we show the comparison between APFEL and CRunDec for the three-loop evolution (NNLO) of the strong coupling α s (left plot) and the charm mass m c (right plot). As is clear from the lower insets, the agreement between the two codes is excellent. Also the one-and two-loop evolutions have been checked finding the same level of agreement.
Finally, we benchmarked the implementation of massive DIS structure functions (i.e. F (3) in eq. 2.1) with MS masses against the public code OPENQCDRAD v1.6 [44]. OPENQCDRAD   implements DIS structure functions in terms of the MS heavy-quark masses following the formalism discussed in Ref. [26]. However, as already mentioned above, such a procedure does not directly correspond to what is needed for the implementation of the FONLL scheme. In order to make the comparison with OPENQCDRAD possible, we have implemented in APFEL a variant of the FONLL scheme with MS masses where, as done in OPENQCDRAD , the RG running of the heavy-quark masses is expanded and truncated to the appropriate order. In Fig. 3 we show the comparison between APFEL and OPENQCDRAD for the exclusive charm neutral-current structure functions F c 2 (left plot) and F c L (right plot) at O(α 2 s ) for three different values of Q 2 and over a wide range of x. As is clear from the lower ratio plots, the agreement is typically at the per-mil level except in the very large-x region where, due to the smallness of the predictions, the relative difference tends to increase but maintains a good level of absolute accuracy.
To conclude this section, we observe that, referring to eq. (2.1), the introduction of the MS masses does not affect the four-flavor structure function F (4) . The structure function F (3,0) is instead affected by the transition from pole to MS masses. Since we are not aware of any public code that computes such structure functions, a direct bechmark has not been possible. However, as a sanity check we have checked that F (3,0) and F (3) for large values of Q 2 tend to the same value, as the definition of F (3,0) requires.

QCD fit settings
The QCD fits were performed to the combined H1 and ZEUS charm production crosssection measurements [9] together with the combined HERA1+2 H1 and ZEUS inclusive  Figure 3: Comparison between APFEL v2.7.0 and OPENQCDRAD v1.6 for the neutral-currents massive charm structure functions with MS heavy-quark masses at O(α 2 s ). As an input PDF set we have used MSTW2008nlo68cl nf3 [45] from which also the numerical values of α s and m c are taken. The upper insets show F c 2 (left) and F c L (right) as functions of x for Q 2 = 10, 100, 1000 GeV 2 as returned by APFEL. In the lower insets the ratios to OPENQCDRAD are displayed showing a relative difference at the per-mil level except in the very large-x region where, due to the smallness of the predictions, the relative differences tend to increase but maintain a good level of absolute accuracy.
DIS cross-section data [7], accounting for all given sources of systematic uncertainties.
The kinematic region covered by HERA is constrained by the invariant mass of the hadronic system of W > 15 GeV and the Bjorken scaling variable of x < 0.65, therefore target mass corrections are expected to have negligible effects and are not discussed in this paper. The settings of the QCD fits in xFitter closely follow those used for the HERAPDF2.0 PDF extraction [7], with a few differences related to the specifics of the current analysis which are motivated in the following.
The nominal result is extracted using the FONLL-C variant of the FONLL scheme discussed in Sect. 2. It should be pointed out that, while being accurate at NNLO for the inclusive DIS cross sections, the sensitivity to mass corrections of the FONLL-C scheme is actually NLO. The reason is that at O(α 0 s ) the FONLL scheme reduces to the parton model which is insensitive to heavy-quark mass effects. Therefore, the first mass-sensitive term is O(α s ) which is the accuracy of the FONLL-A scheme which would thus provide a LO determination of the charm mass. Both the FONLL-B and the FONLL-C schemes, instead, include the O(α 2 s ) massive corrections and thus would both produce determinations of the mass of the charm accurate at NLO. The advantage of FONLL-C with respect to FONLL-B is that it is accurate at O(α 2 s ) also in the massless sector and thus it is supposed to provide a better description of the data. In other words, FONLL-C is the most accurate variant of the FONLL scheme presently available and as such it will be employed for our determination of m c (m c ).
The result obtained in the FONLL scheme is accompanied by an analogous determination of m c (m c ) obtained using the FFN scheme with MS masses [6] at NLO. Access to the structure functions calculated with the FFN scheme is possible via the xFitter interface to the OPENQCDRAD program [44] using the QCDNUM program for the PDF evolution [46].
The procedure to determine the MS charm mass follows closely the methodology described in Ref. [9]. It involves a series of fits in each of which a set of PDFs is determined corresponding to numerical values of charm mass ranging between m c (m c ) = 1.15 GeV and m c (m c ) = 1.60 GeV with steps of 0.05 GeV. For each value of m c (m c ) a value of global χ 2 is obtained. The best fit value of m c (m c ) is determined from the minimum of the parabolic fit to the resulting χ 2 distribution and the associated 1-σ uncertainty, which reflects the sensitivity of the data set to the charm mass, is determined as the ∆χ 2 = 1 variation around the minimum.
We now discuss the settings of the nominal fits and the variations that we performed to assess the different sources of uncertainty deriving from: the PDF parametrization, the model parameters, and the theoretical assumptions.
The assumption that heavy-quark PDFs are dynamically generated via gluon splitting at the respective thresholds requires that the starting scale Q 0 at which PDFs are parametrized is below the charm threshold, which in turn is identified with m c (m c ). Given the range in which the scan of m c (m c ) is done (from 1.1 to 1.6 GeV), we have chosen to set Q 0 = 1 GeV to allow all fits to be parametrized at the same starting scale. The combinations and the relative functional forms of the initial scale PDFs have been chosen following the parametrization scan procedure as performed for the HERAPDF2.0 determination [7], and the optimal configuration has been found to be: There are 14 free parameters, since additional constraints were applied as follows. The QCD sum rules are imposed at the starting scale and constrain the normalisation parameters A g , A uv , A dv . The light-sea quark parameters that affect the low-x kinematic region BŪ and BD, as well as the normalisation parameters AŪ and AD, are constrained by the requirement thatū →d as x → 0, leading to the following constraints: with f s being the strangeness fraction ofD assumed at the starting scale, i.e. f s =s/D, because HERA data alone are not able to provide a precise light-sea flavor separation. The strangeness fraction for the nominal fits is set to f s = 0.4, as in the HERAPDF2.0 analysis [7]. In order to estimate the uncertainty associated to the PDF parametrization, we have considered the following variations with respect to the nominal configuration: • we have moved up the initial scale Q 0 from 1 to √ 1.5 GeV. In the FONLL scheme, this restricted the m c (m c ) range in which we did the scan because we could not use values of the charm mass such that m c (m c ) < √ 1.5 GeV. We were however able to perform the parabolic fit in order to find the best fit value of m c (m c ). This complication does not arise in the FFN scheme in which there is no threshold crossing.
• In the xu v distribution we have included an additional linear term so that the last factor in second line of eq. (3.1) reads (1 + D uv x + E uv x 2 ). After trying different variations of the parametrization, we found that this particular choice leads to the largest differences. The uncertainty associated to model parameters will be estimated by considering the following variations: • the bottom mass has been moved up and down by 0.25 GeV, i.e. m b (m b ) = 3.93 GeV and m b (m b ) = 4.43 GeV. The magnitude of the variation is actually much larger than the present uncertainty on the bottom mass and thus our choice is meant to provide a conservative estimate of the associated uncertainty.
• The variation of the strong coupling follows the recent PDF4LHC prescription [47]. We finally turn to the theory assumptions and their variations. These mostly concern unknown higher-order corrections and the most common way to estimate them is by varying the renormalization and the factorization scales µ R and µ F . As nominal scales in our analysis we have chosen µ 2 R = µ 2 F = Q 2 for both the FONLL 9 and the FFN scheme analyses. Another possible source of theoretical uncertainty in the FONLL scheme is the presence of the damping factor discussed in Sect. 2 which is meant to suppress unwanted subleading terms and whose explicit form in the nominal fits is given in eq. (2.3).
The theoretical uncertainty associated to the missing higher-order corrections has been estimated as follows: 9 A scale choice involving the heavy-quark mass would lead to technical complications with the FONLL matching as implemented in APFEL. However, we have checked that the more commonly used scales µ 2 R = µ 2 F = Q 2 + 4mc(mc) 2 produce a very marginal difference in the determination of mc(mc) in the FFN scheme.
• the factorization and renormalization scales were varied by a factor 2 up and down with respect to the nominal values, that is choosing µ 2 R = µ 2 F = Q 2 /2 and µ 2 R = µ 2 F = 2Q 2 . Such variations have been applied only to the heavy-quark components of the structure functions, while the light part has been left unchanged. The reason for this is that, in order to estimate the theoretical uncertainty associated to the determination of m c (m c ), we want to perform scale variations only in the part of the calculation sensitive to this parameter, which is clearly the charm structure function (for consistency, the same variation was applied also to the bottom structure functions).
• As already mentioned, the FONLL damping factor represents a further source of uncertainty. It has the role of suppressing unwanted subleading terms but the particular way in which this suppression is implemented is somewhat arbitrary. To assess the impact of our particular choice on the determination of m c (m c ), we have changed the suppression power around the nominal one, considering the following functional form: with p = 1, 4.
In addition, to assure the applicability of perturbative QCD and to keep higher-twist corrections under control, a cut on Q 2 is imposed on the fitted data. Our nominal cut is Q 2 > Q 2 min = 3.5 GeV 2 . The choice of the value of Q 2 min requires some care; an extensive discussion on the impact of varying it on the determination of m c (m c ) is given in Sect. 4.3.
To conclude this section, we observe that the self-consistency of the input data set and the good control of the systematic uncertainties enable the determination of the experimental uncertainties in the PDF fits using the tolerance criterion of ∆χ 2 = 1.

Results
In this section we will present the result for our the determination of the value m c (m c ) in the MS renormalization scheme using the FONLL scheme with its associated set of uncertainties.
The parabolic fit to the global χ 2 as a function of m c (m c ) is shown in Fig. 4 and yields a best fit value and its 1-σ experimental uncertainty equal to m c (m c ) = 1.335 ± 0.043 GeV. An estimate of the parametric, model, and the theoretical uncertainties, performed following the procedure described in Sect. 3, is summarised in the second column of Tab. 1 and leads to our final result: An illustration of the deviations, again determined through parabolic fits, caused by the variations employed to determine the parametric, model, and theoretical uncertainties is given in Fig. 5.   Table 2: χ 2 's resulting from the fit in the FONLL-C scheme using the best fit value of the charm mass m c (m c ) = 1.335 GeV. The partial χ 2 's per data point along with the total correlated χ 2 , the logarithmic penalty, and the total χ 2 / d.o.f. are reported, as defined in Ref. [48].
After we have determined the best fit value of the charm mass in eq. (4.1), we have used the central value to perform a further fit in the FONLL-C scheme (nominal fit). In Tab. 2 we report the partial χ 2 's over the number of data points for each subset along with the total correlated χ 2 , the logarithmic penalty, and the total χ 2 per degree of freedom.
As an illustration, the singlet and the gluon PDFs extracted from the nominal fits are compared with other GM-VFNS PDF sets: CT14 [5], HERAPDF2.0 [7], MMHT14 [49], NNPDF3.0 [3]. They are shown in Fig. 6 at the scale Q 2 = 10 GeV 2 , where the the experimental uncertainties from the nominal fits on PDFs are estimated using Monte Carlo procedure with the root mean square estimated from 500 replica. An overall good agreement is observed.
The FONLL determination of m c (m c ) presented above is supported by an analogous determination in the FFN scheme at NLO. The corresponding parabolic fit with the associated experimental uncertainty is shown in Fig. 7. Also in this case a full characterization of the non-experimental uncertainty has beed achieved by carrying out the same parametric, model, and theory variations (except for the variation of the damping factor which is specific of the FONLL scheme). The results of the variation in the FFN scheme are reported in the third column of Tab which is in agreement with the FONLL determination given in eq. (4.1). It is interesting to notice that we observe a reduced scale dependence in the FONLL scheme as compared to the FFN scheme. We ascribe this effect to the fact that the leading contributions in the FONLL scheme involve both gluon-and quark-initiated processes; typically the contributions from gluon processes decrease with the scale, while the contributions from quark processes tend to increase. Conversely, the FFN scheme is mostly driven by gluon processes the contributions of which (along with α s ) tend to be monotonic in µ leading to larger scale variations 10 .
As discussed Sect. 2.1.3, the running of the MS heavy-quark masses in the VFN scheme, exactly like the running of α s and PDFs, is not univocally defined at the heavy-quark thresholds due to the presence of the so-called matching conditions. In particular, when giving the value of the mass at one of the heavy-quark thresholds, one should also specify whether this corresponds to the value immediately below or above the threshold itself.

Comparison to other results
It is interesting to compare our results with the past determinations of MS charm mass m c (m c ) using a similar methodology (also see Ref. [10,25,29] for previous comparisons).
The analysis of Ref. [24] was performed in the ABM11 framework [50] using the FFN scheme at NLO and at approximate NNLO and based on world data for DIS from HERA, and fixed-target DIS experiments and Tevatron Drell-Yan data. While the analysis in Ref. [24] was performed including the same exclusive charm cross-section data used in this study, it did not include the HERA1+2 combined inclusive cross-section data set which was not available at the time, but used instead the HERA combined data from run 1 only. An earlier analysis [23] used a partial charm dataset only, with correspondingly larger uncertainties, while a subsequent analysis [25] investigated the correlation between the measurement of m c (m c ) and the strong coupling constant.
The analysis of Ref. [29] is instead based on the CT10NNLO global analysis, and uses the S-ACOT-χ GM-VFN scheme discussed, e.g., in Ref. [17]. It is based on a slightly wider data set as it includes LHC jet production data and also a set of older F c 2 measurements at HERA [51] that are not included in the more recent combined charm data. The authors of Ref. [29] provide a set of four determinations deriving from different strategies to convert the pole-mass definition into MS. They also provide a separate estimate of the uncertainty due to the O(α 3 s ) corrections for one of the four strategies essentially by varying the parameter that governs a generalized version of the rescaling variable χ.
Finally, a determination of the charm mass m c (m c ) was produced by the H1 and ZEUS collaborations in the framework of the HERAPDF QCD analysis in the same publication in which the charm cross-section measurements employed in our study were presented [9]. That determination also used only the HERA combined inclusive data from run 1 [52].
In Tab. 3 we report the numerical values for the m c (m c ) determinations listed above along with our results and the world average value [53]. A short clarification about the nomenclature of the uncertainties reported in Tab. 3 is in order. In Sect. 3 we discussed extensively the meaning of the uncertainties associated to our determinations. In doing so, we tried to be consistent with the previous determinations, nevertheless some differences remain. As far as the determination in Ref. [9] is concerned, while their definition of "(exp)" and "(param)" essentially coincides with ours, their "(model)" uncertainty includes the variation of the cut in Q 2 (that we will discuss separately in Sect. 4.3) but does not include the α s variation, which is instead quoted separately. In addition, the authors do not quote any scale variation uncertainty. The nomenclature of Ref. [24] is also different from ours. Apart from the common "(exp)" uncertainty, for the NLO determination the authors only quote the "(scale)" uncertainty, which essentially coincides with our "(th)" (even though the FONLL "(th)" uncertainty also accounts for the variation of the damping factor), while for the approximate NNLO determination they also quote a "(th)" uncertainty which, differently from our nomenclature, accounts for the uncertainty on the approximated expressions used at O(α 3 s ). Finally, the determinations in Ref. [29] only quote the experimental uncertainty (the asymmetric uncertainties are due to the use of a generic second-degree polynomial to fit the χ 2 profiles). A graphical representation of Tab. 3 is shown in Fig. 8 where the inner error bars display the experimental uncertainty while the outer error bars (when present) are obtained as a sum in quadrature of all uncertainty sources. The blue vertical band represents the world average and provides a reference for all other determinations. It is clear that, while the spread of the current determinations of m c (m c ) from DIS data covers a pretty large range, they are generally in agreement with the world average. As far as our determinations in particular are concerned, we observe that, apart from being consistent with each other and with the world average, they also present competitive uncertainties. This is particularly relevant for the FONLL determination because this is the first time that this scheme is employed for a direct determination of the charm mass. Fig. 8 shows that our determinations tend to be larger than the world average while most of the previous determinations place themselves below it. Detailed investigations show that the largest contribution to this difference arises from the use of to the new combined HERA1+2 combined inclusive cross section measurements that are employed for the first time to determine the charm mass and that, as we will discuss in Sect. 4 [53] 1.275 ± 0.025 Table 3: List of the recent determinations of m c (m c ) from fits to DIS data along with the determinations extracted in this work. The PDG world average value is also reported for reference.

Cross-checks
It is worth mentioning that we have also employed the variants A and B of the FONLL scheme discussed in Sect. 2 to determine m c (m c ). While the FONLL-A scheme is accurate to LO in the massive sector and thus does not produce a reliable determination of the charm mass, the FONLL-B has the same formal accuracy in the massive sector as FONLL-C and indeed it leads to a determination comparable to that given in eq. (4.1) both for the central value and the uncertainties. It is interesting to notice that the FONLL-B scheme in the low-energy region resembles very closely the FFN scheme at NLO. In particular, both schemes are accurate to O(α 2 s ) in the massive sector and to O(α s ) in the light sector. As a matter of fact, we find that the experimental uncertainty associated to the FONLL-B determination is very close to the FFN one quoted in eq. (4.2), which in turn is around 20% larger than that associated to the FONLL-C determination. This suggests that the O(α 2 s ) corrections to the light sector that are present in the FONLL-C scheme, which depend on the heavy-quark mass by means of diagrams in which a gluon plits into a pair of heavy quarks, provide a further constraint on m c (m c ).
Finally, we have also attempted a determination in the FFN scheme using the approximate NNLO massive structure functions as implemented in OPENQCDRAD. However, we did not pursue a full characterization of the uncertainties because we believe that this determination, while giving a quantitative indication of the effect of the NNLO corrections, cannot claim an NNLO accuracy and thus does not add anything to our NLO determinations.
4.3 Discussion on the Q 2 min dependence of the mass determination Our determination of m c (m c ) given in eq. (4.1) was obtained cutting off all data with Q 2 < Q 2 min = 3.5 GeV 2 . The necessity of such a cut stems from the fact that low-energy data are hard to describe for two main reasons: the large value of α s with consequent large higher-order corrections, and sizable higher-twist corrections. In addition, as pointed  out in Ref. [54], the low-Q 2 region (low-x, in fact) might be affected by deviations from the fixed-order DGLAP evolution whose description might require small-x perturbative resummation. The dependence on Q 2 min of fits to HERA data has already been discussed in the context of the inclusive measurements only. In this section, we will address this issue considering also the HERA charm production data.
The particular value of Q 2 min used in our analysis (3.5 GeV 2 ) was determined by requiring a good fit quality but maintaining a good sensitivity to m c (m c ). This is illustrated in Fig. 9 where the global χ 2 per degree of freedom is plotted as a function of Q 2 min in the left panel while the best fit of m c (m c ) is plotted as a function of Q 2 min in the right panel. Looking at the left panel it is clear that, as expected, the global χ 2 improves as more and more low-energy data are excluded from the fit. On the other hand, the right plot shows that the experimental uncertainty associated to m c (m c ) gets larger and larger as Q 2 min increases indicating that, again as expected, the sensitivity to m c (m c ) deteriorates if lowenergy data are excluded. In the light of the plots in Fig. 9, we conclude that Q 2 min = 3.5 GeV 2 represents a good compromise between a good description of the full data set and a  good sensitivity to m c (m c ).
In this context, it is interesting to look at the behaviour of the partial χ 2 's as a function of Q 2 min of the charm and inclusive cross-section data separately to assess in a more specific way which nominal value of Q 2 min is more convenient. Since the meaning of "degrees of freedom" is unclear for a subset of the full data set, in order to quantify the degree of improvement in the partial χ 2 's, we consider the following quantity: which provides an estimate of the improvement of the χ 2 per data point with respect to our lowest cut Q 2 min = 2.5 GeV 2 . If for a given value of Q 2 min this quantity is larger than one, this means that that specific cut leads to an improvement of the χ 2 which is larger than the degrees of freedom subtracted by excluding a given number of data points and thus the excluded data points with respect of the reference cut (2.5 GeV 2 ) are poorly described. On the contrary, if the quantity in eq. (4.3) is smaller than one, this means that the excluded data points are better described than the fitted ones. In the left panel of Fig. 10 we show the behaviour of the contribution to the global ∆χ 2 /∆N points originating from the charm data points only. It is clear that any cut between 3.5 and 5 GeV 2 improves drastically the partial χ 2 while cuts above 5 GeV 2 either cause a much less significant improvement or even lead to a deterioration. This provides a further confirmation of the fact that our nominal cut (3.5 GeV 2 ) is a sensible choice.
It is also interesting to look at the best fit values of m c (m c ) and the relative uncertainty preferred by a given subset as a function of Q 2 min to quantify the sensitivity to m c (m c ) as more and more data are excluded from the fit. This is plotted in the right panel of Fig. 10 for the charm cross-section data. It is clear that this particular subset of data tends to prefer values of m c (m c ) around 1.23 GeV which is substantially lower than the global value is remarkable and, as expected, the experimental uncertainty tends to increase for larger value of Q 2 min indicating a loss of sensitivity. Finally, we have done the same exercise for the HERA1+2 inclusive cross-section data and in Fig. 11 we present the relative plots. In the left panel we observe that the χ 2 of this subset improves essentially monotonically as Q 2 min increases while from the right panel it is clear that the preferred value of m c (m c ) of the inclusive cross sections is substantially larger than that preferred by the charm cross sections with, again, uncertainties than become broader for larger values of Q 2 min . It is finally clear that our best value for m c (m c ) quoted in eq. (4.1) is a compromise between the lower value preferred by the exclusive charm data and the larger value preferred by the inclusive data. It is clear from the right panels of Figs. 10 and 11 that the exclusive charm and inclusive data subsets prefer somewhat different values of m c (m c ). However, the values shown in these figures are clearly correlated because they were obtained in a simultaneous fit to all data. In order to investigate a possible tension, we have performed a fit to the inclusive data only using both the FONLL-C and FFN schemes. The χ 2 profiles are shown in Fig. 12. In contrast to Figs. 4 and 7, in both schemes the scan in m c (m c ) of the fits to inclusive data only yielded a shallow χ 2 dependences with a minimum around 1.7 GeV. This demonstrates that the inclusive data alone cannot constrain m c (m c ) reasonably well, but also why this data exerts an upwards pull on the m c (m c ) value in the combined fit. Furthermore, since Figs. 9, 10, and 11 in Sect. 4.3 present an overall remarkable stability of the central value of m c (m c ) for different values of Q 2 min , the observed feature cannot be attributed to the low Q 2 part of the inclusive data.

Conclusions
In this work we have presented a new determination of the MS charm quark mass m c (m c ) obtained by fitting HERA charm and inclusive DIS data. In particular, we included in our fits the combined H1 and ZEUS charm production cross-section measurements [9] and the final combination of HERA1+2 H1 and ZEUS inclusive DIS cross-section data [7], the latter being used in this work for the first time for the extraction of the charm mass. Our determination is based on the FONLL general-mass variable-flavor-number scheme, and has required the generalization of the FONLL structure functions, originally constructed in the pole-mass scheme, in terms of MS heavy quark masses.
A detailed estimate of the various sources of uncertainty that affect our determination of m c (m c ) has been performed. In particular, we estimated the uncertainties due to the choice of the PDF parametrization, the model parameters used as input for the theoretical computations, and the missing higher-order corrections. We found that those sources of uncertainty are smaller than the experimental uncertainty, resulting in a competitive determination of the charm mass.
We complemented the FONLL extraction of the charm mass with an analogous determination based on the fixed-flavour number scheme at next-to-leading order, finding a good agreement between the two. In addition, we compared our results with previous determinations also based on fits to DIS data and with the PDG world average finding again a generally good agreement. We find that the values extracted in this work, although compatible within uncertainties, tend to be slightly higher than previous determinations from HERA data. This feature seems to be associated to the final HERA1+2 combined inclusive dataset, which tends to prefer larger values of m c (m c ) as compared to the charm structure function data, and thus increases the best-fit value.
In the future, it would be interesting to repeat the FONLL determination in the context of a global PDF analysis, since, in addition to the inclusive and charm HERA data, other experiments are expected to have some sensitivity to the value of the MS charm mass. In addition, the use of a wider dataset might lead to a reduction of the experimental uncertainties of the m c (m c ) determination. Moreover, our analysis is based on the standard assumption that the charm PDF is dynamically generated by collinear splitting from gluons and light quarks. In this respect, it would be useful to redo the determination of m c (m c ) in the presence of a possible non-perturbative charm PDF, for which the generalized FONLL structure functions accounting for a fitted heavy quark PDF are available [32].