Heavy-flavor hadro-production with heavy-quark masses renormalized in the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document}, MSR and on-shell schemes

We present predictions for heavy-quark production at the Large Hadron Collider making use of the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} and MSR renormalization schemes for the heavy-quark mass as alternatives to the widely used on-shell renormalization scheme. We compute single and double differential distributions including QCD corrections at next-to-leading order and investigate the renormalization and factorization scale dependence as well as the perturbative convergence in these mass renormalization schemes. The implementation is based on publicly available programs, MCFM and xFitter, extending their capabilities. Our results are applied to extract the top-quark mass using measurements of the total and differential tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t\overline{t} $$\end{document} production cross-sections and to investigate constraints on parton distribution functions, especially on the gluon distribution at low x values, from available LHC data on heavy-flavor hadro-production.


Introduction
Charm-anticharm and bottom-antibottom pair-production are among the most frequent inelastic processes occurring in hadronic collisions at the Large Hadron Collider (LHC), with cross-sections smaller than for the dijet case but well above those for the top-antitop case, as follows from the hierarchy of the heavy-quark masses, the available phase-space, and the structure of the Standard Model (SM) Lagrangian.
Heavy quarks are not observable as free asymptotic states. Charm-and bottom-quarks hadronize, due to confinement, whereas the top-quark decays before hadronizing, due to its large decay width. Collider experiments implement procedures for reconstructing topquarks from their decay products and are able to detect the products of the hadronization / fragmentation of intermediate-mass quarks, i.e. heavy mesons and baryons, as well as the jets associated to them, i.e. b-jets and c-jets. B-hadrons and D-hadrons are reconstructed by their decay products, with decay vertices displaced with respect to the primary vertex. This operation requires good tracking, vertexing and particle identification capabilities. The measurements are indeed easier to perform in case of B-hadrons than for D-hadrons, because the proper lifetimes of the first ones are longer than those of the latter.
On the theory side, analytical formulae for the hadro-production of massive quark pairs are known since many years in quantum chromodynamics (QCD) perturbation theory at the next-to-leading order (NLO) accuracy, both for the total cross-sections as well as for one-particle inclusive differential distributions [1][2][3][4][5]. More recently, next-to-next-to-leading order (NNLO) QCD predictions have been computed for the total cross-sections of heavyquark pair-production [6][7][8][9]. On the other hand, differential predictions at NNLO are available for top-quark pair production [10,11], but not yet for the case of charm-quark JHEP04(2021)043 pair. Very recently, a first computation of differential cross-sections for bottom-quark pair production including NNLO QCD corrections has appeared [12].
Beyond fixed-order perturbation theory, the resummation of logarithms in the ratio p T /m, relevant when the transverse momentum of the heavy-quark p T is much larger than its mass m, up to the next-to-leading-logarithmic accuracy, performed through the fragmentation function approach [13], has been combined with NLO predictions [14][15][16]. The resummation of recoil logarithms related to soft gluon radiation from initial state partons, as well as the one of threshold logarithms and high-energy logarithms, have also been worked out (up to various degrees of accuracy) and presented in a number of papers (see e.g. [17][18][19][20][21][22]). The transition from quarks to hadrons is described either by fragmentation functions (FFs) [23][24][25], or by matching NLO matrix elements to parton shower and hadronization approaches [26][27][28][29][30].
One of the inputs of all aforementioned calculations are the values of the heavy-quark masses. The SM by itself does not predict the values of the quark masses, which are thus fundamental parameters and subject to renormalization. The ultraviolet divergences appearing in the heavy-quark self-energies, to be eliminated by renormalization, require to fix a specific renormalization scheme for relating the bare masses to the renormalized ones. The most common choice is the mass in the on-shell scheme, defining the pole mass by the relation that the inverse heavy-quark propagator with momentum p vanishes on-shell, i.e. for p 2 = (m pole ) 2 , and it is known at four loops in QCD [31,32]. Another option, also known at four loops [33,34], is the MS prescription. In complete analogy to the renormalization of the strong coupling constant α S , only divergent terms are absorbed such that the quark propagator becomes finite after wave-function and mass renormalization. Finally, the MSR scheme [35] defines a scale-dependent short-distance mass, that interpolates smoothly between a valid mass definition at low scales much below the mass and the MS mass for scales larger than the mass, using infrared renormalization group flow. Thus, predictions for physical observables in perturbation theory carry scheme dependence through the choice of a particular mass renormalization scheme. For a given observable, the behavior of the truncated expansion in perturbative QCD, including the apparent convergence and the residual scale dependence, can vary significantly between different schemes employed.
In this paper we describe the phenomenological effects of the use of the MS and MSR schemes for the renormalization of the heavy-quark masses in charm, bottom and top production at hadron colliders. We investigate the perturbative convergence in these schemes, by providing comparisons between physical quantities calculated at various levels of accuracy, and we discuss applications concerning the extraction of mass values and parton distribution functions (PDFs) from collider data.
Various phenomenological and experimental investigations on top quark hadro-production with top mass renormalized in the MS scheme already exist in the literature (see e.g. [36][37][38][39][40][41][42]). On the other hand, the MSR mass is discussed in various theory papers, mainly focusing on its definition and properties (see e.g. [43][44][45]). As for the top quark case, the main novelties of our work are the phenomenological predictions for NLO differential cross-sections in pp collisions using the MSR mass, which we compare to those with the MS and pole masses in a consistent framework, the implementation of a dynamical mass renor-JHEP04(2021)043 malization scale in the computation of transverse momentum distributions of the heavyquarks, as an alternative to the static m(m) case, and the extraction of a top MSR mass value from the analysis of state-of-the-art LHC triple-differential cross-section data, previously used for extracting m(m) [46]. The extraction is done preserving correlations with the strong coupling constant α S and the PDFs, fitted simultaneously to the heavy-quark mass.
Total cross-sections for charm and bottom production at hadron colliders with the mass renormalized in the MS scheme are available since a while [47]. The MS scheme has also been used for heavy-quark production in deep-inelastic scattering (DIS), which has allowed extractions of the charm mass m c (m c ) and bottom mass m b (m b ) (see e.g. [48][49][50]). On the other hand, differential cross-sections for charm and bottom production at hadron colliders with MS and MSR masses have never been presented in a dedicated phenomenological paper before this one, to the best of our knowledge. We show examples of them here for the first time and discuss the relative uncertainties due to scale and mass variation, working at NLO accuracy. In case of MS cross-sections we also show the role of variations of the heavy-quark mass renormalization scales, comparing their effects to those induced by the variation of other scales appearing in fixed-order computations, i.e. the factorization and α S renormalization scales.
In recent years heavy-flavor hadro-production data at the LHC have been used to constrain gluon and sea quark distributions at small x in PDF fits at NLO, working in the on-shell mass renormalization scheme [51][52][53][54]. In this work we comment on the simultaneous extraction of proton PDFs and the charm mass in the MS scheme, as an alternative to the pole scheme, using LHCb and HERA data and a fitting procedure were everything else is the same, except the charm mass renormalization scheme. Our findings at NLO support the use of masses renormalized in the MS scheme for further fits of this kind (see e.g. [55,56]). We also show for the first time how LHC charm production data extrapolated to the full phase-space can be used to constrain NNLO PDF fits, suggesting a reduction of the uncertainties in the gluon distributions of various widely used global PDF fits [57][58][59].
The implementation of the MS and MSR schemes for renormalizing the heavy quark masses, as an alternative to the on-shell scheme, is described in section 2. The obtained theoretical predictions for differential distributions including NLO QCD corrections and mass renormalization in these three schemes, are presented in section 3, together with considerations on the convergence of the perturbative expansion in the strong coupling constant. The results are applied in section 4 to investigate the impact of LHC data on possible extractions of PDFs from collider data and on determinations of the heavyquark mass values in different mass renormalization schemes. In section 4.3, we also use predictions for total cross-sections of charm hadro-production up to NNLO accuracy in QCD. Finally, our conclusions are summarized in section 5.

Implementation of heavy-quark mass renormalization schemes
In this work light quarks are assumed to be massless. For the heavy-quark masses, on the other hand, different renormalization schemes can be adopted and we briefly recall the relevant relations for the above mentioned cases of the MS, MSR and on-shell schemes.

JHEP04(2021)043
Other choices for the mass renormalization are possible. For physical observables inherently connected to the heavy-quark production threshold, for instance, the potential subtracted mass [60] was suggested as a possibility to produce an improved perturbative convergence at energies slightly above the quark-pair production threshold and the 1S mass has been presented [61] as a way of stabilizing the position of the peak of the vector-current-induced total cross-section for tt production in e + e − collisions, as a function of the center-of-mass energy √ s, for √ s ∼ 2m. In boosted kinematics, limited to the case of top-quarks, the top-quark jet-mass [62], was introduced in the framework of Soft Collinear Effective Theory. Further mass renormalization schemes are described, e.g., in refs. [44,45] and references therein.
The on-shell mass coincides with the pole in the propagator of the renormalized quark field and is known up to four loops in QCD [31,32]. Thus, it is the same at all scales and infrared finite to all orders in perturbation theory. This definition of the pole mass m pole , although being gauge invariant, has its short-comings [63,64]. It does not extend beyond perturbation theory, i.e. to full QCD, since it is based on the (unphysical) concept of colored quarks as asymptotic states. Therefore, m pole must acquire non-perturbative corrections, which leads to an intrinsic uncertainty in its definition of the order O(Λ QCD ) related to the renormalon ambiguity [65]. The latter manifests itself as a linear sensitivity to infrared momenta in Feynman diagrams, leading to poorly convergent perturbative series for the observables expressed in terms of m pole .
On the other hand, short-distance mass definitions such as the MS or the MSR schemes are renormalon-free. In general, such short-distance masses m sd are related to the pole mass through the relation where the term δm pole-sd removes the renormalon and the dependence of the short-mass definition on long-distance aspects of QCD. Here, µ R denotes renormalization scale, connected with the ultraviolet divergences, whereas the scale R is associated with the infrared renormalization group equation (RGE) [35]. In many short-distance mass renormalization schemes, R coincides with the renormalized mass itself, as for instance in the MS scheme, where R = m(µ R ). However, the possibility to consider other choices of R through the associated RGE allows to improve the stability of the conversion between short-distance mass schemes characterized by different values of R.
In the MS scheme, the renormalized mass of the heavy quark evolves with the RGE in the renormalization scale µ R , governed by the mass anomalous dimensions γ(α S (µ R )), where the perturbative expansion of is known at four loops [33]. Precise determinations of the MS masses for charm-and bottom-quarks are summarized by the Particle Data Group (PDG) [23]. For the MS mass of the top-quark, see, e.g., refs. [23,40,41,49,66]. The conversion to the on-shell scheme proceeds in the standard manner 3)

JHEP04(2021)043
where the first numerical coefficients c i read [67][68][69], . (2.5) Here, ζ denotes the Riemann zeta-function, L ≡ ln(µ 2 R /m(µ R ) 2 ) and the function ∆ accounts for all quarks with masses m i smaller than the heavy-quark one. As the light quarks are taken massless here, i.e., m i = 0, the ∆ term vanishes. The strong coupling α S is evaluated at the scale µ R and renormalized in the MS scheme with the number of active flavors set to n f = n lf + 1 at and above the heavy-quark threshold scale, which is assumed to be equal to its running mass. The number of light flavors is n lf = 3, 4, 5 for charm, bottom and top production, respectively.
For the particular choice m(m) ≡ m(µ R = m(µ R )), i.e. the MS mass renormalized at the specific scale µ R = m(µ R ), the logarithmic terms L cancel and eq. (2.3) evaluates numerically (up to terms O(α 4 S )) as [70] m pole = m(m) 1 + 1.333 The infrared renormalon ambiguity in the conversion in eqs. (2.3), (2.6) manifests itself in practice as factorially growing terms in the perturbative expansion, that spoil convergence. The sizable coefficients in eq. (2.6) indicate the poor convergence of m pole for the case of charm and bottom, when α S at low scales is large. For top-quarks, the convergence is better due to the smaller value of α S at larger scales. Including the four-loop QCD results [33,34], the residual uncertainty in m pole for top-quarks, including renormalon contributions, is estimated of the order of a few hundred MeV [71], i.e., of the order of O(Λ QCD ). All available relations for scheme changes from m(µ R ) to m pole and vice versa have been summarized in the programs CRunDec [72] and RunDec [73]. While α S in eqs. (2.3), (2.6) is renormalized in the MS scheme, the matrix elements, as well as the PDFs and the associated α S evolution used in the fixed-order massive calculations presented in this paper are all defined with a fixed number of light flavors n lf = 3 for charm and bottom production 1 and n lf = 5 for top production, even at scales well above the heavy-quark mass value. For α S renormalization in the decoupling scheme [77], subtractions in graphs with light-quark loops are done at zero mass, as in the MS scheme, whereas those in graphs with heavy-quark loops are done at zero momentum. The heavy quarks therefore do not contribute as active flavors to the running of α S in the effective JHEP04(2021)043 theory. Factorization is performed in the same scheme, which implies the use of nonvanishing PDFs only for light partons, and a calculation of partonic cross-sections done consistently. The latter is achieved by including contributions of virtual amplitudes corresponding to Feynman diagrams with massive heavy-quark fermionic loops. For the case of bottom production in the 3-flavor scheme, we verify that including or not the relevant charm-loop diagrams in our computations, assuming m c = 0, produces differences well within 1% on the total cross-section of bb hadro-production. We argue that including a charm mass different from zero, as appropriate for a fully consistent calculation in the 3-flavor scheme, makes also a small effect, well below the uncertainties on cross-sections due to scale variations. This effect compensates only partially the differences due to the modified α s value, that, at renormalization scales of the order of the bottom quark mass, is a few percent lower in the scheme with n f = 3 active flavors, with respect to the n f = 4 case. Possible advantages of computing bottom production in the 3-flavor scheme have also been claimed in the framework of DIS calculations [78,79], where the effects of including a finite charm-mass in the fermionic loop corrections turn also out to be small. We observe that our predictions for bottom production typically differ by some percent from those of fully consistent calculations with matrix-elements, PDFs and α S in the 4-flavor scheme, depending on the input.
Using the standard decoupling relations, it is possible to relate the PDFs, α S and the partonic cross-sections in the MS and decoupling schemes, once the matching scale is fixed. In this way, also eqs.

vanishes.
In practice, although the perturbative expansion in eqs. (2.3), (2.6) is known up to four loops [33], we truncate it in this work to two loops (order α 2 S ) for computing the NNLO cross-sections and to one loop (order α S ) for computing NLO cross-sections, respectively, unless stated otherwise. In addition, the evolution of α S as a function of µ R and the corresponding α S values entering in eqs. (2.3), (2.6) and other parts of the fixed-order computation are evaluated retaining three loops for producing NNLO cross-sections and two loops for producing the NLO ones, respectively, unless stated otherwise.
The MSR mass is a specific realization of the short-distance mass introduced in eq. (2.1). It is obtained, e.g., by considering the difference between m pole and m(m), see eq. (2.3), and substituting m(µ R ) with R in the terms proportional to α S to determine the difference between m pole and m MSR (R) as where the numerical coefficients a i are given in ref. [35]. The evolution of the MSR mass with the R scale follows the RGE   In the following we use what has been called practical MSR mass in ref. [43], in contrast to the natural MSR mass. For our purposes the numerical differences between those definitions are mostly negligible.
The evolution of the MS heavy-quark masses with renormalization scale is shown in the left panel of figure 1. It is calculated at one loop using the CRunDec program [72,73] (n lf = 3 for charm and bottom, n lf = 5 for top). The evolution of the MSR heavy-quark masses with the R scale at one loop is shown in the right panel of figure 1. It is obtained by solving the RGE in eq. (2.8): where expressions for the first few coefficients of the anomalous dimension γ M SR can be found in refs. [35,43]. Here, eq. (2.9) is expanded up to the lowest non-vanishing order of α s . As is visible in figure 1

Predictions for differential cross-sections
Predictions for cross-sections of heavy-quark production with different mass renormalization schemes can be obtained from those in the widely used on-shell scheme by substituting eqs. (2.3) and (2.7) in the cross-sections and performing a subsequent perturbative expan-JHEP04(2021)043 sion in α s , see refs. [36,38]. In particular, the cross-sections converted to the MS or MSR mass schemes are determined to NLO accuracy as follows Here σ 0 is the Born contribution to the cross-section (proportional to O(α 2 s )), and the differences m(µ R ) − m pole and m MSR (R) − m pole are calculated up to the lowest nonvanishing order in α s , such that all terms of order O(α 4 s ) are dropped in eq. (3.1), as appropriate for a NLO calculation at order O(α 3 s ). Formulae for scheme conversions up to NNLO haven been given in refs. [36,38].
We are considering theory predictions for stable heavy quarks (in case of bottom and charm augmented by FFs to describe the final state B-and D-hadrons, as for the applications in section 4). The additional impact of parton showers and the dependence of the quark mass parameter on their cutoff [80] as well as the study of renormalon effects in observables with cuts leading to corrections of order Λ QCD in the extracted quark masses [81] are subject of ongoing theory research.
We have computed double-differential cross-sections as functions of the transverse momentum p T and rapidity y of the heavy quark Q, and single-differential cross-sections as a function of the invariant mass M QQ of the heavy-quark pair in the on-shell, MS and MSR mass renormalization schemes at NLO using both frameworks, MCFM [82,83] with modifications [38,84] and xFitter [85]. In both cases, the original NLO calculations are done in the pole mass scheme [3,4,86]. The modified MCFM program [38] is capable of converting the NLO calculations using a pole mass into those with the heavy-quark mass renormalized in the MS mass scheme, in case of single-differential distributions in p T of the heavy quark, y of the heavy quark and invariant mass M QQ of the heavy-quark pair. On the other hand, the developed xFitter framework implements the calculation of one-particle inclusive cross-sections (i.e., with the other particle integrated over), i.e., it is capable to compute the double-differential cross-sections as a function of p T and y of the heavy quark, but not as a function of the invariant mass M QQ of the heavy-quark pair. It converts the pole-mass NLO cross-sections into the MS and MSR mass schemes for fully differential distributions. The derivative of the Born contribution appearing in eq. (3.1) is calculated semi-analytically in MCFM (see [38]), whereas it is computed numerically in xFitter, which allows for cross-checks of both methods. The differential cross-sections in different schemes which can be computed by MCFM and xFitter are summarized in table 2. For all crosssections calculated with both programs, MCFM and xFitter, (i.e. all those in the pole mass scheme and the p T and y differential distributions in the MS mass scheme), agreement within one percent accuracy is found. 3 xFitter is also interfaced to other codes, like e.g.

JHEP04(2021)043
Cross-section pole mass scheme MS mass scheme MSR mass scheme Table 2. Summary of the capabilities of the MCFM (M) and xFitter (X) frameworks to compute differential cross-sections for heavy-quark hadro-production in different mass schemes.
aMC@NLO, and thus it can be used for computing cross-sections with heavy-quark masses renormalized in the MS scheme (instead of the pole scheme implemented in the standalone standard version of these codes) for a wider range of processes (e.g. tt+j hadro-production, see section 4.1 for the application of this interface to a phenomenological study).
In the calculations of the differential distributions presented in the following, we use the PDG values [23]

JHEP04(2021)043
the range of a few percent to ∼ 40%, when using the MS or MSR mass scheme instead of the pole mass scheme, and they are more evident in the bulk of the phase space. However, this is a small effect compared to the size of scale uncertainties at NLO. The latter amount to a factor of ∼ 2 in the bulk of the phase space, decreasing slightly towards large p T values. It turns out that the scale uncertainties are very similar in all mass schemes for variations around the chosen nominal scale µ R = µ F = 4m 2 c + p 2 T . Modifying the value of the charm-quark MS mass, which is set to the PDG value in figure 2, to the value extracted in the ABMP16 NLO fit, produce results qualitatively similar, shown in figure 3. Differences between predictions in different mass renormalization schemes in figure 3 are smaller than in figure 2, due to the fact that the ABMP16 MS masses are smaller than the PDG ones.
In figure 4 the same comparison of NLO differential cross-sections in the various mass renormalization schemes is presented for bottom-quark production. In this case, the impact of converting the pole mass calculations into the MS or MSR schemes vary from a few percents to 25%, which is still small compared to the NLO scale uncertainties of the order of 50%. With the choice for the nominal scale µ R = µ F = 4m 2 b + p 2 T , the scale uncertainties are similar in the pole and MSR mass schemes, whereas they are more asymmetric and slightly smaller at low p T in the MS mass scheme. Again, modifying the value of the bottom-quark MS mass, which is set to the PDG value in figure 4, to the value extracted in the ABMP16 NLO fit, leads to results qualitatively similar, shown in figure 5, with slightly smaller differences (up to ∼ 20% among predictions in different mass renormalization schemes.
Finally, figure 6 displays the same comparison for top-quark production. In this case, the impact of converting the pole mass calculations into the MS mass scheme is about 20% at low p T , which is no longer small compared to the NLO scale uncertainties. It decreases towards higher p T values. When converting the cross-sections from the pole to the MSR mass scheme, the impact is below 10% and is within the NLO scale uncertainties for variations around the nominal scale µ R = µ F = 4m 2 t + p 2 T . The scale uncertainties in the MS mass scheme are slightly smaller than in the pole mass scheme, as was already reported previously [38], while the scale uncertainties in the MSR and pole mass schemes are very similar. Again, modifying the value of the top-quark MS mass, which is set to the PDG value in figure 6, to the value extracted in the ABMP16 fit, leads to predictions qualitatively similar, shown in figure 7.
In general, the differences between predictions in different mass renormalization schemes slightly increase with the rapidity, as can be seen in all figures 2-7. In

JHEP04(2021)043
affected by an intrinsic renormalon ambiguity of the order of Λ QCD , as already mentioned in section 2. If we made a more conservative assumption for the uncertainties of charm and bottom quark pole masses, using a value twice as large, i.e. ± 0.5 GeV, the size of the pole mass variation uncertainty bands in figures 8 and 10 would approximately be doubled with respect to the one shown. This is due to the fact that the cross-sections for charm (bottom) hadro-production in the on-shell scheme scale approximately linearly with the charm (bottom) mass (decreasing with increasing masses), at least for values of the charm (bottom) mass close enough to the central one considered in this work. Therefore, calculations in the MS or MSR mass schemes afford substantially smaller uncertainties (in particular at low p T ) due to precise input quark masses. In figure 12 the single-differential cross-sections as a function of the invariant mass M QQ of the heavy-quark pair in the pole and MS mass scheme are shown, as calculated using MCFM (no implementation of the MSR scheme is available for this distribution). The impact of changing from the pole to the MS mass scheme is largest at the lowest values of M QQ close to the production threshold. At a technical level, this is due to the derivative term in eq. (3.1) becoming dominant in this kinematic region. However, this implies that the term δm pole-sd = m pole − m sd in eq. (2.1) for the conversion of m pole to a short distance mass grows parametrically as δm pole-sd ∼ m sd α S , hence is no longer small either. This situation is realized for the MS mass definition and it persists even when changing the MS mass value, as follows from the comparison of figure 12 with figure 13, where different m(m) values are employed. This excludes the MS scheme from being a suitable mass renormalization scheme for observables very close to threshold, cf. ref. [38] for a detailed analysis for the top-quark pair invariant mass distribution. Alternative mass renormalization schemes for observables dominated by the production threshold have been mentioned in section 2.
For comparison to current experimental data on pair-invariant mass M QQ distributions in hadro-production, however, this aspect is of minor relevance. For instance, for top production at the LHC [46,66] the size of the lowest M tt bin is large, extending to O(50) GeV above threshold, so that sensitivity to threshold dynamics is significantly reduced and the MS mass scheme is still applicable in analyses of those data.
In figure 14 we show the impact of using different PDF sets (together with their α S (M Z ) value) in the rapidity distributions for charm-, bottom-and top-quarks (see also ref. [84]). We fix the heavy-quark MS masses to the PDG values. Slight changes in the normalization of the distributions can be ascribed to the fact that different   normalization and in shapes are related to the different behaviour of different PDFs as a function of x and µ F . In particular, in case of charm production, the shape of the rapidity distribution obtained with the central set of the MMHT14 PDF fit [57] for pp collisions at √ s = 7 TeV is much wider with respect to that obtained with the central PDF sets from other widely used fits. This is particularly evident when using the MS heavy-quark mass, instead of the pole mass, in the computation, due to the lower value of the first one with respect to the second one, and is related to the peculiar and very flexible MMHT14 PDF parameterization and the particular behaviour of the gluon distribution at small x. JHEP04(2021)043 At the scales relevant for the calculation, the MMHT14 NLO central gluon distribution steeply rises for smaller x and displays large uncertainties, in absence of data capable of constraining it for x < 10 −4 in the fit, see also ref. [87]. On the other hand, in case of top and bottom production, the differences among predictions making use of different PDF sets are smaller than for the charm case, because, for fixed rapidity values, these processes probe larger (x, Q 2 ) values, where more data have been used to constrain the various PDFs. In particular, the predictions obtained by different PDF sets, turn out to be within the scale uncertainty band computed using the ABMP16 NLO PDF nominal set, at least for rapidities away from the far-forward region.

JHEP04(2021)043
Additionally, in this paper we explore the possibility of using a dynamical scale in the heavy-quark MS mass renormalization, as an alternative to the static value m Q (m Q ) and its variations used in the previous distributions and in ref. [38]. There, the p T distribution of the top-quark at NLO was computed for static central scales µ R = µ F = µ m = m t (m t ), varying them simultaneously by factors 1/2 and 2 around their central value and finding that the scale uncertainty band was reduced with respect to the case when µ R and µ F are varied and µ m is fixed to m t (m t ). In general, we expect that dynamical scales, catching the different kinematics of different events, provide a more accurate description of differential distributions. Thus, in the following we consider the case when the central values for the renormalization and mass renormalization scales are chosen dynamically and coincide, i.e., We fix the central factorization scale to the same value µ F = µ 0 . For this configuration, we compute scale uncertainties, by two different procedures. In the first procedure (i) we fix µ R = µ m even in the scale variation but we still vary independently µ R in the interval [µ R,1 , µ R,2 ], where µ R,1 = 0.5 p 2 T + 4m 2 Q (µ R,1 ) and µ R,2 = 2 p 2 T + 4m 2 Q (µ R,2 ), and µ F in the interval [1/2, 2] around the chosen (mass) renormalization scale, excluding the (µ R , µ F ) extreme combinations (2, 1/2) and (1/2, 2), but keeping all the others, as in the conventional 7-point scale-variation procedure. These variations implicitly also encode a heavy-quark mass variation, with the mass value spanning the interval [m(µ R,2 ), m(µ R,1 )]. In the second variation procedure (ii), which is more general than (i), we fix µ R = µ m = µ F = µ 0 = p 2 T + 4m 2 Q (µ m ) in the central predictions as before, but we vary µ R , µ F and µ m independently from each other, each by factors 1/2 and 2 around µ 0 , excluding the extreme scale combinations as in the conventional scalevariation procedure. In other words, we release the constraint µ R = µ m during the variation of these scales. This procedure leads to a 7-point (µ R , µ F ) scale variation band at fixed  In figure 16 we also add leading order (LO) predictions to all panels, obtained by using the same values of m(µ m ) and m(m) as in the NLO ones, and the same NLO PDFs and α S (M z ) value. The central LO predictions are not always included in the NLO uncertainty bands, however the LO uncertainty bands (not shown) are much larger than the NLO ones, cf. also ref. [38] for top production.
Another example of dynamical mass renormalization scale choice was shown in ref. [42], where the MS mass renormalization scheme was studied as an alternative to the pole mass scheme for producing predictions for top-quark related distributions at NNLO. There the tt-pair invariant mass distribution was studied at NNLO, using the MS mass at a scale µ m ∼ M tt /2, setting µ R = µ F = M tt /2 and making a 15-point scale variation of factors (1/2, 2) around the (µ R , µ F , µ m ) central value. Predictions were compared to the case when µ m = m t (m t ), seeing small differences. On the other hand, p T and rapidity distributions were computed using static µ m values. To the best of our knowledge, our paper is the first work where the use of a dynamical µ m scale in computing p T distributions for heavy-quark hadro-production is investigated.
We checked that our NLO predictions are consistent with those reported in ref. [42], when using their configuration. In figure 17 we present the p T distribution of the antitopquark for tt production in pp collisions at √ s = 13 TeV, using as input the NNPDF3.1 NLO PDF set with its α S (M z ) default value and α S evolution, i.e. one of the configurations already considered in ref. [42], and multiple choices for the (µ R , µ F , µ m ) scales. For fixed µ m = m t (m t ) = 163.7 GeV, we observe that central predictions using µ 0

JHEP04(2021)043
have larger (µ R , µ F ) uncertainty bands (especially in the peak region) and have smaller absolute values than those using central scales µ R = µ F = p 2 T + m 2 t (m t ) or µ R = µ F = m t (m t ), with differences between central values at the peak amounting to ∼ 10%, as shown in the lower panel. The latter two scale choices can be considered better scale choices (i.e. scale choices leading to a faster perturbative convergence) for tt production than the first one, as also proven by the fact that NNLO corrections, reported in ref. [42] for the (µ R , µ F ) case in comparison to the NLO ones, are quite small. On the other hand, the central predictions we obtained using µ R = µ F = p 2 T + 4m 2 t (µ R ) are in much better agreement with the previous ones than those with (µ 0 R , µ 0 F ), as shown in the upper panel, and have smaller scale uncertainty bands (not reported in the plot), which shows that the use of m t (µ R ) instead of m t (m t ) in the dynamical scale definition improves the perturbative convergence of the calculation, corresponding to smaller NNLO corrections. The predictions with (µ R , µ F ) are larger than those with (µ 0

Phenomenological applications
The use of the theory results can be illustrated with a number of applications in phenomenology, determining the strong coupling constant α S (M Z ) and values of the top-quark mass in the different renormalization schemes as well as constraints on PDFs by using available LHC data.

t + α S (M Z ) from differential tt cross-sections at NLO
The top-quark mass can be extracted using measurements of the total or differential tt production cross-sections. As an example, we use the recent CMS measurement [46] of normalized triple-differential tt cross-sections as a function of invariant mass and rapidity of the tt pair, and the number of additional jets. These observables provide decent sensitivity to the values of m t (m t ) and m MSR t in a simultaneous fit with α s (M Z ) and the PDFs, i.e. the complete set of input theoretical parameters of fixed-order calculations for stable top-quark JHEP04(2021)043 pair production. We compare the results with the ones obtained in the CMS analysis [46]. In particular, the distributions of the tt invariant mass and the additional jet multiplicity are sensitive to the top-quark mass through threshold and cone effects [37].
The QCD analysis setup follows the original CMS analysis [46], and the main settings are summarized in the following paragraph. 6 The QCD analysis is done using the xFitter framework [85]. Theoretical predictions for the tt data are obtained at NLO in the pole mass scheme using the MadGraph5_aMC@NLO program [91], interfaced with aMCfast [92] and ApplGrid [93] to store the calculated cross-sections bin-by-bin in the format which is suitable for PDF fits with xFitter. The dependence of the theoretical predictions on the top-quark mass is taken into account by generating several sets of predictions with different values of this parameter and smoothly interpolating them in the fit. The HERA combined inclusive DIS data [89] are included in the fit to provide constraints on the valence and sea quark distributions and to probe the gluon distribution and α s through scaling violations, while the CMS tt data provide direct constraints on the gluon PDF and α s , as well as on the top-quark mass as discussed in ref. [46].
In our analysis, we convert the NLO calculations for the tt production cross-sections from the pole mass scheme into the MS or MSR mass schemes according to eq. As the analysis of triple-differential tt cross-sections requires NLO predictions not only for inclusive tt production (N jet ≥ 0), but even for inclusive tt+1 jet production (N jet ≥ 1), MadGraph5_aMC@NLO is the only public code, among those providing such calculations, that is already interfaced to ApplGrid. In general, also other frameworks implementing NLO QCD corrections could be adopted, even beyond the fixed-order studies considered here, but they are not yet interfaced with ApplGrid.
The fit results obtained using different mass schemes are given in table 3. The values of χ 2 characterize the fit quality. These values are very similar in all fit variants and illustrate a general good description of the tt data. To estimate uncertainties, we follow the procedure from ref. [46] and determine fit, model, parametrization and scale variation uncertainties. As in the CMS analysis, the scales are varied coherently in all bins of the measured crosssections. As shown in table 3, in the pole mass scheme, switching from the dynamic scale H to the static scale m pole t modifies the extracted pole mass by about 0.6 GeV, a value still smaller than the fit uncertainties amounting to 0.7 GeV, but enlarges the scale uncertainties substantially. Therefore, the larger scale uncertainties obtained in this analysis using the MS or MSR mass schemes, as compared to ref. [46] Table 3. The values for α S (M Z ) and the top-quark mass in different mass schemes obtained in ref. [46] and in this work by fitting the CMS data on tt production and the HERA DIS data [89] to theoretical predictions. The fit, model (mod), parametrisation (par) and scale variation (µ) uncertainties are reported. Also the values of χ 2 are reported, as well as the partial χ 2 values per number of degrees of freedom (dof) for the tt data (χ 2 tt ) for 23 tt cross-section bins in the fit. The scale H is defined in the text. the MSR mass m MSR t (3 GeV) does not affect the scale uncertainties significantly. On the other hand, if in the future one would know the value of the MS masses very precisely (from some other measurement), one could use them to get accurate predictions for differential cross-sections with smaller heavy-quark mass uncertainties, while the pole mass would be affected by O(Λ QCD ) uncertainties.
In the light of these observations, it will be worth to implement the transition to the other mass schemes directly in the MadGraph5_aMC@NLO program 8 and in further Monte Carlo integrators/event generators, such that predictions for differential tt cross-sections in association with jets can be obtained in the format which is suitable for PDF fits in different mass schemes and with dynamical scales. The advantages of the latter for the running masses have been illustrated in the previous section 3.
Furthermore, in table 3 we do not observe any noticeably larger theoretical uncertainty when fitting the MS running mass instead of the pole mass, contrary to what was reported in refs. [40,41]. A direct comparison of our analysis to those ones is not possible because of the different data sets used. Our analysis is based on triple differential distributions in the invariant mass and rapidity of the tt-pair and in the number of light jets, using recent  [23,40,41,49]. The world average labelled as 'PDG2018, appr. NNLO' is based on a single determination of the D0 collaboration [90]. precise data obtained by the CMS collaboration during the LHC Run 2 at √ s = 13 TeV, whereas the studies in refs. [40,41] make use of a different differential distribution (dσ/dρ s , where ρ s is an observable related to the inverse of the M ttj invariant mass), as measured by the ATLAS collaboration, at different center-of-mass energies ( √ s = 7 and 8 TeV, respectively), during Run 1. Switching from the pole mass scheme to the MSR mass scheme with R = 3 GeV changes the extracted mass value by 0.3 GeV only, which is well within the current experimental and theoretical uncertainties. On the other hand, the value of α s (M Z ) extracted from the fit does not change significantly when using different schemes, as also shown in table 3. The obtained values are compatible with α s (M Z ) = 0.1191 ± 0.0011 extracted in the ABMP16 fit at NLO [56] within two standard deviations. 9 The extracted value of m t (m t ) is compared with several other determinations in figure 18. In the ABMP16 analysis, the running top-quark mass was determined from measurements of total top-quark pair and single-top production cross-sections in a global QCD fit at NNLO [49]. In ref. [41] ATLAS extracted a m t (m t ) value at NLO from their measurement of tt + 1 jet production cross-sections, while ref. [40] has obtained m t (m t ) at NLO using the ATLAS measurement of tt + 1 jet production [94] on the basis of LHC Run-1 data. Currently, the world average value of m t (m t ) by the PDG [23] is based on a single determination of this parameter by the D0 collaboration at approximate NNLO [90]. When comparing to the other determinations of m t (m t ) displayed in figure 18, it is worth to note that only the results of this work and of the ABMP16 analysis are obtained in a JHEP04(2021)043 simultaneous fit of m t (m t ), α s (M Z ) and PDFs, preserving correlations among these quantities, while the other determinations were done using a value of α s (M Z ) and a PDF set fixed a-priori.
In line of principle, the applied methodology can be extended to the extraction of the m c (m c ) and m b (m b ) mass values from measurements of charm and bottom production in association with jets at colliders. However, this is a great challenge from the experimental point of view, because measuring jets at low p T , where the sensitivity to the charm-and bottom-quark mass would be particularly large, is hard.

NLO PDF fits with differential charm hadro-production cross-sections
The application of differential distributions for charm hadro-production with the MS mass definition allows for an update of PDF fits which use heavy-flavor measurements from the LHC, to constrain the gluon distribution at low values of the longitudinal momentum fraction x [51][52][53]. In particular, constraints at the lowest x values explored nowadays (x 10 −6 ) can be obtained by considering the charm hadro-production process at high rapidities (|y| 4.5) at the LHC, whereas the bottom hadro-production process at similar rapidities at the LHC is sensitive to slightly larger x values (x 10 −5 ), with a region of sensitivity that partially overlaps with the one of charm data. Because of the large scale dependence of the NLO calculations for charm hadro-production, it is customary to include in such fits only ratios of cross-sections, which are constructed using measurements at different values of rapidity and/or transverse momentum, or at different center-of-mass energies.
As an example, in the PROSA analysis [51] charm and bottom hadro-production crosssections [74,95] as a function of rapidity were used in ratios to the respective cross-section in the rapidity interval 3 < y < 3.5 for each p T bin, together with the inclusive DIS data [96] and the heavy-flavor production DIS data [97,98] from HERA. These ratios feature a reduced scale dependence, but, at the same time, they have reduced sensitivity to the heavy-quark mass. We repeat this PROSA analysis using the MS heavy-quark masses in the calculations of both the DIS structure functions [48] and the charm and bottom hadro-production cross-sections, instead of pole masses, while all other settings are as in ref. [51]. As a result, we observe only a small impact on the χ 2 value and the fitted PDFs, with a new central PDF that is well within the previously found PDF uncertainties. These small differences are driven mainly by the change in the predictions for the HERA data, because the LHCb data used in the format of normalised cross-sections do not provide any notable sensitivity neither to the heavy-quark mass scheme, nor to the value of the heavy-quark mass. As a result, the fitted MS heavy-quark masses are determined as  [97][98][99][100], where the HERA data alone were analyzed to determine the heavy-quark MS masses.
The MS masses are also in better agreement with the world average values [23], than the pole masses of eqs. (4.3), (4.4), indicating that the latter carry a significant intrinsic theoretical uncertainty. Therefore in our most recent PDF analysis [55] we have solely adopted heavy-quark running masses.

NNLO PDF fits with total charm hadro-production cross-section
The NLO predictions for differential charm hadro-production at the LHC have very large scale uncertainties (> 100% in some phase space regions), as illustrated in section 3. The lack of theory predictions for differential cross-sections on charm and bottom hadroproduction at NNLO prevents including the corresponding existing data in the state-ofthe-art PDF fits, which nowadays are mostly provided at NNLO accuracy. In this context measurements of the total charm hadro-production cross-section would be beneficial, because they can be confronted in the PDF fits with the already available inclusive NNLO predictions [6][7][8][9] which have significantly reduced scale uncertainties. However, no such measurements have been performed to date. On the other hand, the ALICE [101,105], ATLAS [102] and LHCb [95,103,104] experiments have provided measurements of charm production in different kinematic regions which cover more than one half of the phase space. One can reliably determine the total cross-section by extrapolating these measurements to the full phase space. The extrapolation procedure is analogous to that adopted for extracting reduced cross-sections for charm production in ep collisions at HERA [100] from measurements in a fiducial phase space. These reduced cross-sections are then routinely used in global PDF fits. In the following, we perform such extrapolations and provide the inferred values of the total cc production cross-section at different center-of-mass energies and their ratios, together with experimental and theoretical uncertainties arising from the extrapolation procedure. We then compare the results to theoretical predictions at NNLO in QCD which are obtained using different PDF sets, and demonstrate how these data can help to reduce PDF uncertainties.
The existing most precise LHC measurements of open charm production are summarized in table 4. The ALICE measurements at √ s = 5 and 7 TeV cover the central region |y| < 0.5, the LHCb measurements at 5, 7 and 13 TeV provide coverage of the forward region 2 < y < 4.5, and the ATLAS measurement at 7 TeV essentially bridges the gap by providing data at |η| < 2.1. However, while both ALICE and LHCb provide measurements nearly in the full p T range starting from 0 GeV, ATLAS reports the cross-sections only for p T > 3.5 GeV, thus leaving the bulk of the corresponding p T kinematic range unmeasured. Furthermore, it turns out that the most precise data of ALICE and LHCb among all open Dmeson data are those for D 0 production, while this final state was not measured by ATLAS.
Given these arguments, we extrapolate ALICE and LHCb measurements of D 0 production at 5 and 7 TeV to the full phase space. In order to maintain the least dependence on the theoretical predictions, both ALICE and LHCb measurements are extrapolated to nearby regions of y, namely to 0 < |y| < 1.5 and |y| > 1.5, respectively: where Here σ ALICE and σ LHCb denote the ALICE and LHCb data on fiducial cross-sections, respectively, and σ NLO in different rapidity ranges are the theoretical predictions. The factor 2 in the second term takes into account that the LHCb data are provided for only 2 < y < 4.5 and need to be extrapolated to 2 < |y| < 4.5. We exploit the symmetry around y = 0 and assume that the cross-sections for 2 < y < 4.5 and −4.5 < y < −2 are equal, as reasonably expected in case of pp collisions. Also the measurements are extrapolated into the full range of p T (not shown in eqs. (4.5), (4.6) for brevity), which implies only a 1% correction for the LHCb data at 7 TeV provided for 0 < p T < 8 GeV, and even smaller corrections for the ALICE data sets. This procedure is used to obtain the total cross-section for D 0 production at collision energies √ s = 5 and 7 TeV, while at 13 TeV we extrapolate solely the LHCb measurement since no other data are available at this energy. 10 We calculate the total charm production cross-section from the D 0 production crosssection dividing the latter by the fragmentation fraction from ref. [107]: The factor 2 in eq. (4.7) accounts for the fact that both c andc fragment into charmed hadrons. We assume f (c The uncertainty on f (c → D 0 ), which amounts to 1%, is neglected. We also compute ratios of cross-sections at different center-of-mass energies R 7/5 = σ 7 TeV /σ 5 TeV and R 13/7 = σ 13 TeV /σ 7 TeV , which benefit from a partial cancellation of theoretical uncertainties [108]. 10 Preliminary predictions on D 0 production in pp collisions at √ s = 13 TeV were reported by the ALICE collaboration in a conference proceeding [106] in 2018, but they have neither been confirmed yet nor further refined in a regular article. In addition, the data are presented in plots, but no numerical tables are provided in ref. [106].

JHEP04(2021)043
The theoretical predictions σ NLO in eqs. (4.5), (4.6) are computed using the MS masses as described in the previous sections. The hard-scattering cross-sections for heavy-quark hadro-production are supplemented with phenomenological non-perturbative FFs to describe the c → D 0 transition. The factorization and renormalization scales are chosen to be µ R = µ F = µ 0 = 4m 2 c (m c ) + p 2 T and varied by a factor of two up and down (both simultaneously and independently for µ R and µ F ) to estimate the scale uncertainties with the conventional 7-point scale variation, leaving out the combinations (µ R , µ F ) = (0.5, 2)µ 0 and (2, 0.5)µ 0 . The MS charm-quark mass is set to m c (m c ) = 1.275 ± 0.030 GeV [23].
The proton is described by the PROSA PDF set [51], which is expected to have a reliable gluon distribution at low x thanks to the heavy-quark data used for its determination. Furthermore, to estimate the PDF uncertainties, the extrapolation is performed using the ABMP16 [56], CT14 [58], MMHT2014 [57], JR14 [109], NNPDF3.1 [88] and HERAPDF2.0 FF3A [89] NLO PDF sets. Then, the envelope covering the PROSA PDF uncertainties and the difference obtained using any of the additional PDF sets is constructed. This conservative procedure is essential, because the theoretical calculations for the highest y values involve the gluon PDF at the lowest x values (up to 4 × 10 −8 ), which are not directly covered by data in any of the PDF fits (not even in the PROSA fits which include the charm data up to y = 4.5 as measured by LHCb, for which PDFs at x < 10 −6 and their uncertainties are extrapolated from the results obtained up to x ∼ 10 −6 , using built-in procedures in the LHAPDF library [110]).
All theoretical uncertainties are assumed to be fully correlated for cross-sections in different kinematic regions and at different center-of-mass energies. The robustness of the extrapolation procedure is checked by varying the boundary y between the kinematic regions into which the ALICE and LHCb measurements are extrapolated by ±0.5 (at the same time, this variation tests consistency of the ALICE and LHCb data). As a further check of the method, we have computed predictions for the ALICE and LHCb data using NLO matrix elements matched, according to the Powheg method [111,112], to parton shower and hadronization implemented in PYTHIA8 [113], and found these predictions to be consistent with our NLO + FF predictions within theoretical uncertainties.
The results of the extrapolation are reported in table 5. The scale, mass, PDF and fragmentation uncertainties are added in quadrature to obtain the total theoretical uncertainty assigned to the extrapolated results. The experimental uncertainties of the input data are propagated to the extrapolated cross-sections and reported separately. The experimental uncertainties of the input data sets are assumed to be fully uncorrelated. 11 The experimental and theoretical extrapolation uncertainties are approximately of the same size. The total uncertainties are obtained by adding the experimental and theoretical uncertainties in quadrature, and amount to ≈ 10%. Our value for the total charm cross-section at 7 TeV  Table 5. Extrapolated total charm production cross-sections and their ratios at different centerof-mass energies together with uncertainties from parametric variations of the scales at NLO, the mass m c (m c ) ± 0.03 GeV, α K ± 1.7, PDFs and the rapidity y ALICE,LHCb ± 0.5. The correlation factor between R 7/5 and R 13/7 is −0.61. α S uncertainties are negligible compared to the PDF ones, computed using as a baseline the CT14 PDF set of eigenvectors at NLO.
is in agreement with the extrapolated cross-sections reported in refs. [101,102,114] within uncertainties.
While the central values for the extrapolation factors in eqs. (4.5), (4.6) were obtained at NLO, their uncertainties are calculated such that they should cover potential deviations from the unknown true QCD result. Therefore the resulting total cross-sections, with these uncertainties included, can be compared to calculations in any QCD scheme to any order. Furthermore, for determining these extrapolation factors, only the shape of the predictions for the p T and y differential cross-sections is relevant, while a large part of theory uncertainties related to normalization cancels.
The extrapolated cross-sections and their ratios are compared to NNLO predictions obtained using the NNLO PDF sets ABMP16 [49], CT18 [59], MMHT2014 [57], JR14 [109], NNPDF3.1 [88] and HERAPDF2.0 [89]. The cross-sections are computed using the Hathor program [47] interfaced in xFitter [85]. The factorization and renormalization scales are chosen to be µ R = µ F = 2m c (m c ) and µ R and µ F are varied by a factor of two up and down according to the 7-point scale variation procedure to estimate scale uncertainties. The MS charm-quark mass is set to m c (m c ) = 1.275 GeV [23].
In figure 19 we show the extrapolated cross-sections and their ratios compared to NNLO predictions. For the NNLO predictions, the theoretical uncertainty arising from scale variations and the PDF uncertainty are shown separately. All theoretical predictions agree with the data within uncertainties, but noticeably the MMHT2014, HERAPDF2.0 (and CT14, not plotted in the figure) PDF sets have uncertainties which are larger than both scale and data extrapolation uncertainties for some of the observables. In particular, the MMHT2014 and HERAPDF2.0 predictions for the cross-sections at √ s = 13 TeV are consistent with negative values within uncertainties (see also ref. [87]). Predictions based on the new CT18 PDFs (and unlike those using the previous PDF set CT14) do not show anymore a large positive uncertainty which greatly exceeds the extrapolated cross-section. These PDF sets could benefit from including in their fits data on charm production crosssections or on their ratios. Remarkably, also the scale uncertainties appear to be different when using different PDF sets. Among the considered observables, the most conclusive one is R 7/5 for which both data extrapolation and theoretical scale uncertainties are moderate (≈ 10%). Our extrapolated value for this observable can be used in future NNLO PDF fits to constrain the gluon PDF at low x. The other ratio R 13/7 has a larger extrapolation uncertainty suffering from a lack of experimental measurements of charm production in the central rapidity region at 13 TeV. We are confident that this lack will be solved by the data which will appear in forthcoming experimental studies at the LHC (see footnote 10).
As a demonstration that the provided observables can indeed constrain the PDFs, we employ a profiling technique [115] based on minimizing the χ 2 function built from data and theoretical predictions, taking into account both data and theoretical uncertainties arising from PDF variations. The analysis is performed using the xFitter program [85]. We consider the MMHT2014 PDF set at NNLO and the ratios R 7/5 , which exhibits the least scale uncertainties, and R 13/7 . The correlation of R 7/5 and R 13/7 due to the common input of 7 TeV data sets is taken into account. The PDF uncertainties are included in the χ 2 functional through nuisance parameters, and the values of these nuisance parameters at the minimum are interpreted as optimised (or profiled) PDFs, while their uncertainties determined using the tolerance criterion of ∆χ 2 = 1 correspond to the new PDF uncertainties.
The original and the profiled MMHT2014 gluon PDF are shown in figure 20 at the scales Q 2 = 10 and 100 GeV 2 . The profiled distribution exhibits greatly reduced uncertainties at low x, and in this region the distribution is shifted towards larger values of the gluon density. In case of the MMHT2014 set, the original gluon PDF is negative at low x values, while the profiled one remains positive down to at least x ∼ 5·10 −6 , thanks to the constraint realized by adding the ratios of charm data in the PDF fit. We emphasize that the strong impact at low x is obtained as well when working with other PDF sets. As an example,  figure 20, but for the CT18 PDF at Q 2 = 10 GeV 2 (up left) and Q 2 = 100 GeV 2 (up right) and for the CT14 PDF at Q 2 = 10 GeV 2 (down left) and Q 2 = 100 GeV 2 (down right), all at NNLO. in figure 21 the CT14 and CT18 gluon distributions are shown before and after profiling. For these sets the gluon PDF is always positive in the entire x range for all eigenvectors by construction. In case of CT14, adding the aforementioned data strongly reduces PDF uncertainties at low x, whereas the effect is milder for CT18, but still sizable at low Q 2 .

Conclusions
The hadro-production of heavy-quarks is an important class of processes at LHC. Not only for top, but also for bottom and charm, a wealth of very precise high-statistics data has been collected by the LHC collaborations, differential in the relevant kinematic variables, such as the transverse momentum p T , the rapidity y or the pair-invariant mass M QQ of the JHEP04(2021)043 heavy-quarks (or of the respective heavy hadrons). In comparison to theory predictions in perturbative QCD, these data can be directly used for the extraction of heavy-quark masses, which are typically correlated with the value of the strong coupling constant α S (M Z ). The data also have an impact on fits of fundamental non-perturbative QCD parameters such as PDFs, where they give unique kinematic constraints.
In order to provide meaningful determinations of heavy-quark masses, the value of α S (M Z ) and PDFs, QCD predictions with good accuracy are needed. Theory predictions are available at NNLO for top-quark hadro-production since some time, also for differential distributions. Very recently differential predictions have also appeared for bottom pair hadro-production, but not for charm pair. In the latter case, the predictions at NLO accuracy are generally not sufficiently precise enough considering the large theoretical uncertainties, which stem predominantly from scale variations. In view of the much smaller experimental uncertainties reached in modern analyses improvements in the theoretical descriptions are clearly required.
One such aspect, which has been studied in this paper, is the choice of a suitable renormalization scheme for the heavy quark masses. We have investigated different heavy-quark mass renormalization schemes with emphasis on the MS and MSR masses as representative short-distance mass definitions. The choice of a particular mass scheme as well as the values for the scales µ R and µ F have an impact on the rate of apparent convergence of the perturbative expansion of the cross-sections. We have investigated a range of dynamical scale choices for the cross-sections and, in case of the MS mass, also for the mass renormalization scale µ m . In particular, we have found that dynamical renormalization and factorization scales of the type (µ R , µ F ) p 2 T + κ m 2 Q (µ m ) for heavy-quark p T -distributions with the running of mass m Q (µ m ) included and µ m set equal to µ R , when m Q (µ m ) does not reduce to a value below about 1 GeV, can lead to reduced residual scale uncertainties, with respect to the use of the analogous functional form with m Q (µ m ) replaced by m Q (m). The amount of reduction depends on the input parameters of the computation and on the convention adopted for scale variation. The most general scale variation procedure that we have considered corresponds to independent variations of µ R , µ F and µ m in the range [1/2, 2] around their respective central values. The resulting 7-point envelope in (µ R , µ F ) keeping µ m = µ 0 m dominates the total uncertainty band, i.e. adding an independent µ m variation for a 15-point envelope in (µ R , µ F , µ m ) increases the size of the band only by a moderate amount. This confirms previous findings by other authors in case of top quark pairs. The maximum amount of scale uncertainty reduction that we have observed when adopting this scale variation procedure, when comparing p T distributions with m Q (µ m ) to those with m Q (m) with the input configuration of figure 16 (κ = 4), occurs in case of top-antitop production and amounts to a few tens percent in the region around the peak of the top-quark p T distribution. At NLO accuracy in QCD scale uncertainties are, however, in general still large for all mass schemes, but theory predictions using MS or MSR masses carry smaller parametric uncertainties in the mass values, being theoretically well-defined and free of renormalon ambiguities.
We have demonstrated these features in extractions of the top-quark MS and MSR masses at NLO from recent differential distributions measured by CMS, finding good con-JHEP04(2021)043 sistency with other determinations. These extractions were performed simultaneously to the one of α S and PDFs, preserving correlations between these quantities. Using differential charm hadro-production cross-sections we have also been able to improve available constraints on PDFs and, using the MS mass scheme, to decrease extrapolation uncertainties when determining total cross-section from open charm data measured in fiducial regions of phase space by the LHC collaborations. In the latter case, ratios of cross-sections are particularly useful observables to cancel residual theoretical uncertainties. In order to carry out theses studies, we have developed software frameworks using the MCFM and xFitter programs to determine differential distributions at NLO in QCD efficiently.
Avenues for theoretical improvements include the obviously needed QCD predictions for charm hadro-production at NNLO accuracy, possibly combined with the resummation of large logarithms in specific kinematics, but also an improved description of charm-and bottom-quark fragmentation to mesons, an issue which has been side-stepped in the present analysis. In addition, further systematic studies of different (µ R , µ F , µ m ) dynamical scale choices for different differential distributions of heavy-quark hadro-production are desirable.
The extended xFitter program, implementing the MSR and MS mass renormalization schemes, as an alternative to the on-shell scheme in heavy-quark hadro-production, is publicly available on the web, and further extensions of the MCFM and Hathor programs used to perform calculations in this paper are available from the authors upon request.