Revisiting the D-meson hadroproduction in general-mass variable flavour number scheme

We introduce a novel realization of the open heavy-flavour hadroproduction in general-mass variable flavour number scheme at next-to-leading order in perturbative QCD. The principal novelty with respect to the earlier works is in the treatment of small-transverse-momentum limit, which has been a particularly challenging kinematic region in the past. We show that by a suitable choice of scheme, it is possible to obtain a well-behaved description of the open heavy-flavour hadroproduction cross sections from zero up to asymptotically high transverse momentum. We contrast our calculation with the available D$^0$-meson data as measured by the LHCb and ALICE collaborations at the LHC, finding a very good agreement within the theoretical and experimental uncertainties. We also compare our framework with other theoretical approaches.


Introduction
The hadroproduction of heavy-flavoured mesons at the LHC, in particular the D-and B-meson measurements at forward direction [1][2][3][4], has recently attracted a growing interest for its potential to provide information on partonic dynamics at low momentum fractions. Because of the finite heavy-quark mass m, the perturbative methods are applicable down to zero transverse momentum (p T ) of the observed meson, and the measurements provide opportunities e.g. to constrain the collinearly factorized gluon distributions at small momentum fractions in proton [5][6][7] or nucleus [8,9], or to test other scenarios like saturation physics [10,11], or k T factorization [12]. The D-meson production is also of great interest from the viewpoint of neutrino astrophysics as the secondary neutrinos from D mesons produced in scatterings of cosmic rays in the atmosphere form a significant background for the extraterrestrial neutrinos. Given that the D-meson measurements at LHCb [1][2][3] are kinematically close to the cosmic-ray-on-air scattering, the rates for secondary neutrinos can be constrained by the LHC data [13][14][15][16]. In heavy-ion collisions the measured open heavy-flavour data [17,18] provides opportunities e.g. to test the so-called dead-cone effect [19,20] in QCD medium [21,22].
Theoretically, there are several collinear-factorization-based ways to calculate cross sections for heavy-flavoured mesons in proton-proton (p-p) collisions, see e.g. Refs. [22][23][24] for reviews. On one hand, parton-level heavy-quark cross sections at fixed-flavour-number scheme (FFNS) [25][26][27] can be folded with phenomenological, scale-independent partonto-meson fragmentation functions (FFs), or the parton-level calculation is matched to a parton-shower [28,29] from a general-purpose Monte-Carlo event generator, such as Pythia 8 [30] or Herwig [31], and the showered event is then hadronized according to the hadronization model of the generator. Alternatively, one can work fully within the collinear factorization where the fragmentation is described with universal, scale-dependent FFs [32]. In this paper, we will focus on this latter approach.
The general framework in QCD to treat the heavy-quark production is the so-called general-mass variable flavour number scheme (GM-VFNS) [33][34][35][36][37][38][39]. In this framework, at low interaction scales Q 2 m 2 the heavy quarks are not treated as partons in PDFs but are considered only as massive objects in the final state. The full mass dependence is retained in the production cross sections, but the initial-state partons are restricted only to the light ones. These cross sections contain mass-dependent logarithmic terms which, towards higher interaction scales, will eventually dominate and diverge. In GM-VFNS these large logarithms are subtracted at a certain transition scale Q 2 t -typically the heavy-quark mass threshold -and resummed into the PDFs and scale-dependent FFs. At asymptotically high interaction scales Q 2 m 2 the result reduces (up to finite terms) to the calculation where the quark mass has been put to zero from the outset, the so-called zero-mass variable flavour number scheme (ZM-VFNS).
To obtain a well-behaved description for the heavy-flavoured mesons within GM-VFNS approach from zero to asymptotically large p T has, however, been a bit challenging. The difficulty is related to the intrinsic freedom in GM-VFNS to use the zero-mass formalism for the processes with heavy-quarks in the initial state or where the fragmenting parton is a light one. The massless coefficient functions display a divergent behaviour towards low p T and with a typical scale choice Q 2 ∼ m 2 + p 2 T their contribution dominates the cross sections immediately above p T = 0. Thus, the production cross sections diverge towards p T → 0.
A solution was proposed in Ref. [40]. In essence, the idea was to exclude the aforementioned divergent contributions at low p T by retaining the factorization and fragmentation scales at the threshold Q 2 = m 2 until large-enough p T . Formally, the difference with respect to a more natural choice Q 2 = m 2 + p 2 T is one order higher in QCD coupling than what one works at, but numerically the effect is large and the cross sections are rendered finite down to p T = 0. A relatively good description of the LHCb data can be obtained by tuning the scales in this manner [15,40], but the price to pay is that there will be a certain unphysical wiggle in the production cross section near the region where one decides to turn on the heavy-quark PDFs and light-parton FFs, see e.g. Figure 6 in Ref. [40]. An alternative strategy along this line would be to take the transition scale Q 2 t to be much higher than the heavy-quark mass [41]. However, this would lead to a discontinuity in the cross sections at the arbitrary point where one decides to make the transition. Higher-order calculations should decrease the transient effects in both cases, but will not cure them completely. Clearly, a different solution would be beneficial.
The option we propose here is to make use of the scheme dependence inherent to GM-VFNS. Physically, our choice of scheme is rooted in the observation that -in the absence of intrinsic charm component -the contributions from heavy-quark PDFs and light-parton FFs are simply an efficient way to resum diagrams where a heavy quark-antiquark (QQ) pair is dynamically produced. Being of the same origin, it is natural to require that these contributions respect the same kinematical constraints as the channels where the pair is explicit produced. These are formally O(m 2 ) effects and can be included in the definition of a scheme. However, the contributions from heavy-quark PDFs and light-parton FFs will no longer diverge in the p T → 0 limit, but are regulated by the heavy-quark mass. The production cross sections thus remain finite in the p T → 0 limit with arbitrary factorization and fragmentation scales.

Formalism
In this section we will describe our theoretical construction and its numerical implementation. As the GM-VFNS framework in hadroproduction of heavy quarks has been detailedly discussed in Refs. [38,42,43], we will here focus only on the most important features of our approach. However, enough details are still given so that our results can be reproduced.

General structure and kinematics
The process we study is an inclusive production of a hadron h 3 with momentum P 3 in collision of two hadrons h 1 and h 2 with momenta P 1 and P 2 , In the approximation where the masses of partons and produced hadron are neglected, the cross section differentiated with respect to the produced hadron's transverse momentum P T and rapidity Y can be written in the well-known factorized form [44], dσ ij→k (τ 1 , τ 2 , µ 2 ren , µ 2 fact , µ 2 frag ) dp T dy .
where the fragmenting parton's transverse momentum and rapidity are p T = P T /z and y = Y . Here f h 1,2 i (x 1,2 , µ 2 fact ) are the PDFs for parton species i in hadron h 1,2 and D l→h 3 (z, µ 2 frag ) is the parton-to-h 3 FF. The invariants τ i are defined as where p 1 and p 2 are the momenta of the incoming partons, p 3 is the momentum of the produced, outgoing parton and √ s is center-of-mass (c.m.) energy of the collision. The integration limits are given by

Partonic kinematics in the presence of mass
When a QQ pair is produced from light partons, the zero-mass partonic kinematics above should be adjusted to account for the heavy-quark mass m. In practice, this amounts to replacing the partonic transverse momentum p T in the x 1,2 integration limits and scaling variables τ 1,2 by the partonic transverse mass m T ≡ p 2 T + m 2 , These kinematics correspond to the inclusive heavy-quark production. When the produced parton is a heavy quark, the above replacements follow directly from the momentum conservation. However, in the case that the fragmenting parton is a light one or when there is a heavy quark in the initial state, these replacements are strictly speaking not necessary, but are part of our choice of scheme (SACOT-m T , explained in more detail later).
In the picture where the heavy quarks are generated perturbatively, the heavy-flavour PDFs and light-flavour FFs are merely an efficient way to resum diagrams where a heavy quark-antiquark pair is created. That is, the production of heavy-flavour pair is implicit in these contributions and motivates the usage of heavy-flavour kinematics.

Massive fragmentation variable
The zero-mass version of the fragmentation scaling variable z = P 3 /p 3 is ill-defined in the presence of massive quarks/hadrons, and the zero-mass relations Y = y and P T = zp T are no longer true. Here, we choose to define the scaling variable z in a Lorentz-invariant way as As indicated, in the c.m. frame of the colliding hadrons z can be interpreted as the fraction of partonic energy carried by the outgoing hadron [45]. Alternatively, the scaling variable could be defined e.g. in terms of light-cone momentum fractions [40,46]. From the above definition and considering the fragmentation to be collinear in the c.m. frame, we have two equations, where the hadronic transverse mass is defined as M T ≡ M 2 + P 2 T , M being the mass of the produced hadron. We can solve these equations for the hadronic transverse momentum and rapidity, The cross section corresponding to the above definition of z can be obtained as by integrating over p T and y. Using the relation we find again Eq. (2.1) where the partonic transverse momentum and rapidity are now given by 12) and the hadron mass corrects the lower limit of the z integration as Otherwise the cross-section formula is formally identical to the case of zero-mass partons and hadrons.

Partonic cross sections in SACOT and SACOT-m T schemes
The starting point in our GM-VFNS construction, is the next-to-leading-order (NLO) one-particle inclusive heavy-quark cross section in FFNS [25][26][27] where heavy flavour can be produced in three different partonic processes, g + g → Q + X, q + q → Q + X, q + g → Q + X . (2.14) In FFNS (and also GM-VFNS at low interaction scales), these are the only ways to produce heavy flavour. The heavy-quark mass m is kept finite in these processes and in the high-p T limit, the partonic cross sections develop logarithmic divergences ∼ log(p 2 T /m 2 ) coming from kinematic regions where the heavy quarks become collinear with other partons. These are the first terms in the whole series of large collinear logarithms which, in GM-VFNS framework, are resummed to heavy-quark PDFs and parton-to-hadron FFs when the interaction scale exceeds a chosen transition scale Q t . From now on identify Q t as the heavy quark mass, Q t = m. To avoid double counting when including the contributions also from heavy-quark PDFs and using the scale-dependent parton-to-hadron FFs, one has then to subtract these logarithmic pieces from the coefficient functions. In what follows, we will explain what are P q g P q g P q g P q q (a) (b) (c) (d) Figure 1. Diagrams illustrating the origin of four types of collinear logarithms in g + g → Q + X → h 3 + X channel. The relevant partonic splitting functions are indicated. Made with JaxoDraw 2.1 [47].
the added and subtracted terms in our case, using the g + g → (Q → h 3 ) + X channel as an explicit example. The gluon-fusion process g + g → (Q → h 3 ) + X at NLO entails four different sources of collinear divergences in the m → 0 (or equivalently p T → ∞) limit, illustrated in Figure 1: -one of the two initial-state gluons splits into a collinear heavy quark-antiquark pair, -an outgoing gluon splits into collinear heavy quark-antiquark pair, -an outgoing heavy quark emits a collinear gluon. A simple way to specify the GM-VFNS subtraction terms at NLO is to take as the starting point the leading-order (LO) contributions from channels where there are heavy quarks in the initial state or the fragmenting parton is a light one. Let us begin with the former case. Using Eq. (2.1), we write the leading-order contribution for process g This now uniquely determines the subtraction term which cancels the logarithmic term from diagrams like (b) in Figure 1. The expression for perturbative heavy-quark PDF, to the first order in strong coupling α s , reads where P qg (z) = T f z 2 + (1 − z) 2 with T f = 1/2, is the leading-order gluon-to-quark splitting function. Using this expression for f h 2 Q in Eq. (2.15) gives our definition of the subtraction term, When adding the leading-order contribution of Eq. (2.15), one must then compensate by subtracting Eq. (2.17). The difference contributes at O(α 4 s ) and is not considered at an NLO-level O(α 3 s ) calculation. Here, we also plainly see the origin of the scheme dependence in GM-VFNS: The exact form of |M gQ→Qg | 2 appearing in Eq. (2.15) and Eq. (2.17) is subject to a certain amount of arbitrariness. Indeed, the only requirement is that in the m → 0 limit |M gQ→Qg | 2 must tend to the zero-mass expression |M gq→qg | 2 so as to ensure that the corresponding collinear logarithm from g + g → (Q → h 3 ) + X process cancels. Otherwise we can choose it at will. Similarly, the exact expressions for the integration limits are irrelevant as far as the zero-mass expressions given in Eq. (2.3) are found in the m → 0 limit. The simplest option is to use the zero-mass matrix elements and kinematics from the beginning -this choice of scheme is usually dubbed as simplified ACOT, or SACOT scheme [48]. Here, we shall adopt a prescription where we use the zero-mass matrix elements but still retain the kinematic mass dependence. In other words, the integration limits and the invariants τ 1 and τ 2 are as in Eq. (2.4), and for the squared matrix element in Eqs. (2.15) and (2.17) we take, where |M gq→qg | 2 (τ 1 , τ 2 ) is obtained from the zero-mass expression [49], with C A = 3, C f = 4/3, D A = 8, and the "massive" Mandelstam variables being noŵ In practice, our prescription amounts to replacing the partonic transverse momentum p T in the zero-mass expressions by the transverse mass m T -hence we shall name the present implementation as SACOT-m T scheme. 1 The leading-order contribution and subtraction term for the Q + g → (Q → h 3 ) + X channel are defined in a similar manner as above, so let us then discuss the contributions from light-parton fragmentation and the corresponding subtraction terms. Proceeding as in the case of initial state, we define the leading-order contribution from g + g → (g → h 3 ) + X channel, originating from diagrams like (c) in Figure 1, by Together with the the perturbative expression for the gluon fragmentation function (considering that the only non-zero FF at the mass threshold is D Q→h 3 ), this defines the subtraction term where now The subtractions required to cancel the large logarithm originating from diagram (d) in Figure 1 goes slightly different than the above cases. The reason is that the contributions from g + g → (Q → h 3 ) + X channel (part of the inclusive heavy-quark cross sections) that we here use to determine the subtraction terms, are included using the full mass-dependence. Therefore, the subtraction term required to cancel the large logarithm that occurs when final-state heavy quark emits a collinear gluon is is the quark-to-quark splitting function, and the matrix element [25], now carries the full mass dependence. In order to recover the standard MS zero-mass results at high P T we must still compensate for the fact that the m → 0 limit in the massive calculation does not exactly match that of usual massless MS, but some finite differences remain as a relic of a different regularization procedure. As explained in Ref. [43], this can be effectively achieved by replacing Eq. (2.24) by where d QQ is the partonic fragmentation function [50,51], Also the renormalization procedure applied in FFNS calculations is slightly different than in the purely zero-mass case. Indeed, the FFNS results of Ref. [25] are obtained in a so-called decoupling scheme, where α s runs with only light partons (gluons + 3 light-flavour quarks). Above the transition scale Q t , also the heavy-quark is considered as being "active" in the running of α s , and the matching between the two schemes induces an additional contribution. Specifically, we must add a term as explained in Ref. [39]. The same line of reasoning is applied when defining the subtraction terms for q + q → Q + X and q + g → Q + X channels and the emerging leading-order contributions from q + q → (g → h 3 ) + X, q + Q → (Q → h 3 ) + X, and q + g → (g → h 3 ) + X channels. As in Eq. (2.27), the definition of the S qq→(Q→h 3 )+X subtraction term involves the partonic fragmentation function d QQ , and a term is added to recover the MS renormalization scheme [39]. In addition, our full results include the contributions from all other partonic subprocesses whose inclusion does not require a preparation of subtraction terms at the perturbative order we work at. The NLO O(α 3 s ) contributions, taken from Ref. [44] are included as well. We stress that when including these terms in our SACOT-m T scheme, we consistently retain the kinematics which they inherit from QQ pair-creation process as explained earlier. In practice this is done by trading the massless variables v and w used in Ref. [44] by their massive counterparts v → 1 − τ 1 and w → τ 2 /(1 − τ 1 ), see e.g. Sect. 2 of Ref. [43], and imposing the proper integration limits explained in Section 2.1. In this way, we already implicitly define the subtraction terms that would be required at a next-to-NLO (NNLO) -level calculation.
In comparison to the earlier works [42,43], the most notable advantage of the SACOTm T scheme is that the cross sections remain finite in the P T → 0 limit. Indeed, in Refs. [42,43] at least part of the contributions not coming directly from flavour-creation processes are included using purely zero-mass formalism, and give rise to a divergent P −n T behaviour at P T → 0 limit. The difficulty will not be completely resolved at NNLO either, though the divergences may be a bit "softer". In Ref. [40] these divergent contributions were excluded at small P T = 0 by maintaining the factorization and fragmentation scale at (or below) the heavy-quark mass threshold until large-enough P T . This procedure leads to finite cross sections in the P T → 0 limit, but causes certain unphysical slope change near the P T value where the factorization and fragmentation scales go above the mass threshold -we will come back to this in Section 3.2 (see also see Fig. 6 in Ref. [40]). In our case -and this applies also for the fixed-order calculations -the divergent behaviour is regulated by the heavy-quark mass and leads to finite cross sections even at P T = 0 (at any perturbative order) without a need to fine tune the scale choices. Technically, this happens because the lower limits for the scaling variables τ 1 and τ 2 appearing in the squared matrix elements are not zero but limited by the heavy-quark mass.

Numerical implementation
Our numerical realization of the GM-VFNS scheme described above is crafted around the public INCNLO [44,52] and Mangano-Nason-Ridolfi (MNR) [53,54] codes. The former provides the zero-mass matrix elements, and the latter one the one-particle inclusive heavy-quark cross section of Ref. [25]. As already noted in Ref. [55], in order to obtain reliable numerical results from INCNLO at high √ s away from the midrapidity |y| 0, the numerical stability of the original code has had to be improved, see p.30-32 in Ref. [56] for a detailed explanation. Schematically, we compute  where the inclusion of charge-conjugate contributions and shuffling between the initial-state partons is implicit. That is, from the full zero-mass result we subtract the zero-mass contributions of g + g → Q + X, q + q → Q + X, and q + g → Q + X channels which we add back using the full mass dependence. The subtraction terms provide the proper matching. Towards high P T , only the first sum term in Eq. (2.31) survives -others add up to zero. In the numerical evaluation we have used NNPDF31 nlo pch as 0118 variable-flavour-number PDFs and the corresponding running strong coupling α s [57]. This is the latest NNPDF fit assuming no intrinsic charm content in the proton. The PDFs are interfaced by using LHAPDF 6 library [58]. The introduced framework is applicable to production of any hadrons involving heavy quarks but in this work we consider only D-meson production due to good availability and precision of the experimental data from LHC experiments. In particular, we will focus on D 0 -meson production, with the data from LHCb [1][2][3], ALICE [62], and CMS [18] (though not yet available) extending to small P T region, which is where our SACOT-m T scheme mostly differs from other GM-VFNS implementations (the LHCb collaboration has also measured D ± at small P T ). For D * mesons there would be more recent FF analyses available [59,60] but for D 0 we use KKKS08 [61] FFs, which is the only available FF set for D 0 's. For the charm-quark mass we use m charm = 1.51 GeV in accordance with the used PDF set. The input charm mass in KKKS08 analysis was 1.50 GeV so the pairing with NNPDF3.1 is consistent. Our default scale choice will be µ ren = µ fact = µ frag = P 2 T + m 2 charm , and for the D-meson mass we use M = 1.87 GeV. The small contribution from b-quark fragmentation is retained in the calculation neglecting the finite b-quark mass. The D mesons from B-meson decays have been excluded from the LHCb and ALICE data we discuss later on, but as the KKKS08 FFs include these feed-down D mesons as well, there is no fully consistent way to exclude them without explicitly evaluating the D 0 meson spectra from B-meson decays and subtracting it from the fully inclusive cross section. However, the contributions from B decays are very small, less than 1% in the integrated inclusive D 0 -meson cross section of ALICE [62].
In order to compare with another popular approach, we have used here the Powheg method [63] in which the QQ production at FFNS is matched with the Pythia parton shower providing NLO accuracy for the matrix element generation and leading-log resummation from the parton shower. In practice, we have first generated cc events with the hvq part [29] of Powheg-Box generator [64]. The generated events are then fed into Pythia (version 8.230) [30] which generates the p T -ordered parton shower and hadronizes the events using the implemented Lund string model with parameter values from the default Monash tune [65]. The D 0 mesons (and its charge conjugate) are then picked up from the hadronized final state and binned in Y and P T . The same NNPDF3.1 PDFs as for the GM-VFNS calculations have been used for the event generation in Powheg and also in showering within Pythia. The sensitivity of the Pythia shower to PDFs is very mild as they affect only the initial-state emission probabilities and there only ratios of PDFs are involved.
In Powheg generation the default scale choice is µ ren = µ fact = p 2 T + m 2 charm with m charm = 1.5 GeV. We have not explicitly introduced the matching terms, Eqs. (2.29) and (2.30), at the heavy-quark mass thresholds as their effect has been found small in the P T range of LHCb data [66]. Indeed, with µ ren = µ fact the first matching term in Eq. (2.29) is zero, and the second term in Eq. (2.29) is small as the LO contribution of qq channels is small. As discussed in Ref. [66], Powheg+Pythia yields very similar results as e.g. FONLL [39,67] or Madgraph5 aMC@NLO [28,68] approaches in the kinematic domain of LHCb.

Results
In this section, we will first illustrate some features of our calculation that we have studied numerically and then compare with the available experimental LHC data.

Consistency checks and other trivia
We begin to fold out the numerical results by showing in the left-hand panel of Figure 2 contributions from the channels where the cc pair is explicitly produced. Here, we have   The results are normalized by the full GM-VFNS calculation including all the partonic channels. At high P T the solid and dashed curves merge which provides a non-trivial, strong check on the consistency of our implementation. Towards P T → 0 the two sets of curves, however, behave completely differently: Whereas all channels of the "massive" calculation yield a positive contribution at P T → 0 limit, even the overall result with zero-mass matrix elements remains negative. As can be seen from the left-hand panel of Figure 2, the overall contribution from the channels where the cc pair is explicitly produced, is only a few percents from P T ∼ 5 GeV onwards. In fact, almost the entire cross sections in this region accumulates from the partonic subprocesses with heavy quarks in the initial state or gluon fragmentation, around 50% coming from each of these two sources. This is demonstrated in the right-hand panel of Figure 2 where we plot the contributions from these channels, normalized to the full GM-VFNS result. The balance between the contributions shown in the left-and right-hand panels of Figure 2 depends rather strongly on the scale choices at low P T , and the pace at which the contributions in the right-hand panel begin to dominate can be controlled by adjusting the scales. Indeed, using a lower scale than our default choice, the contributions shown in the right-hand panel would begin to dominate at higher P T than now shown in Figure 2. As already mentioned in Section 1, it was exactly this property that was taken advantage of in Refs. [15,40] to suppress the divergent contributions at low P T . σ(y = Y, p In the left-hand panel of Figure 3 we estimate the effects of charm-quark and D-meson masses in our cross sections at Y = 0. The green curve corresponds to a ZM-VFNS calculation (but still using the aforementioned default scale choice) normalized with the full GM-VFNS result. In accord with what was seen in Figure 2, we observe that neglecting the charm mass leads to a lower cross section at low P T due to increasingly negative contributions from g + g → (c → h 3 ) + X and g + q → (c → h 3 ) + X channels in ZM-VFNS. The blue dashed curve corresponds to putting m, M = 0 in Eq. (2.11) and Eq. (2.12), that is, ignoring the mass dependence in the fragmentation variable z. We observe that this manoeuvre leads to increased cross sections. The origin of the effect can be understood relatively easily on the grounds of Eqs. (2.11) and (2.12) from which it follows that when m < M . That is, for fixed Y and P T the partonic cross sections are probed at larger y and larger p T in comparison to the massless kinematics. Since the partonic cross sections decrease steeply, particularly with increasing p T , also the hadronic cross sections are consequently lower. In our framework, this explains why the hadronic cross sections are suppressed in the presence of non-zero masses. This is in contrast to what has been found in Ref. [40] in the case of B mesons, though there a different version of the fragmentation variable z was used. Moreover, in Ref. [46] a very similar definition of z as in Ref. [40] was adopted and there, in turn, the mass effects led to suppressed cross sections (as in our case). To clear up the systematics of different definitions of the fragmentation variable warrants a separate study which is beyond our scope here. Nevertheless, the effects of finite hadron and quark masses can be non-negligible up to P T ∼ 20 GeV which signifies a possibly considerable source of theoretical uncertainty, given that the definition of fragmentation variable z is ambiguous. In the context of the present definition of z these effects will, however, get milder towards larger Y . This can be easily understood from Eq. (2.11) from which it follows that p T ∼ P T /z when Y 0. Thus, with the present definition of the fragmentation variable z, part of the significant effect found at Y = 0 can be expected to melt away. This is demonstrated in the right-hand panel of Figure 3 where we show the impact of massive fragmentation variable also in the forward direction. The differences between massless and massive fragmentation variable get clearly suppressed when moving to larger Y . At very large Y the effect starts to rise again as the partonic y spectrum gets steeper (near y ∼ 0 it is quite flat) and the condition |y| ≥ |Y | of Eq. (3.1) begins to matter increasingly. The validity of our calculation towards low P T → 0 could be potentially compromised by the unstable fixed-order NLO scale evolution of the fragmentation functions below z ∼ 0.1, stemming from singular α 2 s (log 2 z)/z terms in the time-like NLO quark-to-gluon and gluon-to-gluon splitting functions, see e.g. Ref. [70]. The proper treatment of this region requires resummation in both, splitting functions and Wilson coefficients [71]. To exclude contributions from the unstable region we have imposed a condition z > 0.1 when computing the cross sections. When doing so we must then make sure that this cut is not overly strict, i.e. that the contribution outside of the introduced cut is negligible. The reason why the P T → 0 limit could pose a problem, can be easily understood: As discussed e.g. in Refs. [45,55], approximating the convolution of partonic cross sections and PDF by a power law,  where n > 0 and C k does not depend on p T , one gets in the zero-mass approximation. If the partonic spectrum drops sufficiently strongly in p T (i.e. the exponent n is large enough), the factor z n−1 efficiently eliminates the contributions from the problematic low-z domain. However, in the low-P T region the LHC data [1-3, 62, 72] show that the hadronic D-meson cross sections tend to level off towards P T → 0, see Figure 7 ahead. That is, the exponent n in Eq. (3.2) decreases and the mechanism above is not as effective in suppressing the small-z contributions. In Figure 4 we show z distributions obtained directly from the full calculation for a few fixed values of P T . Unlike could have been expected on the basis of the above discussion, the cross sections are found practically inert to the small z region even at very small P T . Here, the explanation seems to be in the form of the D-meson fragmentation functions which at low µ frag are clearly suppressed in the small-z region as shown in Figure 5. Towards higher µ frag the small-z tails go up but then also the probed p T is larger (larger exponent n) and the contributions are suppressed by virtue of Eq. (3.3). The behaviour of the D-meson fragmentation functions are quite different in comparison to e.g. typical pion fragmentation function as demonstrated in Figure 5 as well. All in all, the cross sections get hardly any contributions from the small z region. At NNLO and beyond, also the PDF evolution becomes similarly unstable at small x, and resummation [73] appears to be required in order to optimally reproduce the small-x HERA data at low Q 2 [74,75]. However, at the NLO level we work at, these issues are not yet that pressing. Indeed, based on Ref. [74], the effects of resummation in PDFs are only modest at NLO, and thus we expect that the small-x resummation would lead to only subleading effects in comparison to the very large scale uncertainty in low P T D-meson production, see e.g. Figure 7 ahead.  Before comparing with the data, we wish to shortly discuss the predominant x ranges sampled by D meson production, in particular the small-x sensitivity within the LHCb acceptance. To this end, the upper left-hand panel of Figure 6 presents examples of x 2 distributions as obtained from our GM-VFNS calculation. The distributions are presented for √ s = 13 TeV with typical LHCb kinematics. While the lower limits for the x 2 distributions are always deep in the small-x domain, the distributions carry a long tail towards large xin all of the considered cases there is a clearly non-negligible contribution coming even from the x 2 10 −2 region. In part, this tail originates from the NLO contributions in processes where the cc pair is explicitly produced, but mostly it comes from the new partonic channels that "open" as a result of resumming the collinear logarithms into non-zero heavy-quark PDFs and light-parton FFs. To corroborate this point, the upper right-hand panel of Figure 6 shows various x 2 distributions for the bin Y ∈ [2, 2.5]. The LO calculation with the heavy-quark PDFs and light-parton FFs turned off (only "direct" cc) gets almost no contribution from large x, but when the NLO contributions are switched on, a smallish large-x tail develops. 2 The normalized x 2 spectra from the full LO and NLO calculations are mutually very similar, both getting a significant contribution from x 2 10 −3 unlike the contributions from "direct" cc production processes. This shows that the large-x tail mostly originates from other than explicit cc processes. To illustrate this point even further from a different viewpoint we show, in the lower panel of Figure 6, also predictions from pure LO Pythia simulations. On one hand, when the origin of the D meson is restricted to underlying cc events, the x 2 distributions are clearly suppressed at large x. On the other hand, when the D mesons are allowed to be created from all partonic QCD processes (here we omitted b quarks), the large-x tail emerges. All in all, the D-meson production at forward rapidity is sensitive to the PDFs at small-x, but the role of large-x contributions is still clearly non negligible. The importance of the large-x part is something that has maybe been a bit underrated in many recent articles [5,6,66], in part because their importance does not show up in fixed-order-based calculations. We note that a very similar large-x behaviour is present also in isolated photon production [76]. In inclusive pion production the large-x tail is even more pronounced [76] due to different behaviour of parton-to-pion FFs, see Figure 5.

Comparison with LHCb and ALICE data
In this section we will present comparisons of our calculation with the experimental LHCb and ALICE D 0 data taken in p-p collisions at √ s = 5 GeV [3], √ s = 7 GeV [1, 62] and √ s = 13 GeV [2]. We will consider two types of theoretical uncertainties, namely those related to the choices of the QCD scales and those related to PDF errors. In principle, there are also other sources of theory uncertainty in the input variables, like the value picked for the charm-quark mass or the value of α s at the Z-boson pole. However, to consistently vary these quantities they should be accompanied by PDFs and FFs extracted with the same variation. The are no error sets (unlike for most modern PDFs) available for the D 0 meson fragmentation functions, so the uncertainties in FFs are not considered either. The definitions of fragmentation variable and heavy-quark scheme are taken here as inherent to the presented calculation.
The PDF uncertainties are evaluated in the standard NNPDF way by computing the variance where dσ i is the cross section computed with the ith member out of the collection of N rep PDF replicas, and where        is the central prediction. The stability of the results against scale variations is quantified as e.g. in Ref. [55] by varying the three scales independently as where the parameters α i=fact,ren,frag vary between 0.5 and 2. The total scale uncertainty is taken as the maximum and minimum of the 16 cross section found in this way. With the above choice, the scales remain always above the charm mass-threshold and the potentially large contributions from log(µ ren /µ fact ) and log(µ ren /µ frag ) terms are suppressed by limiting the maximal difference of the respective scales to a factor of two. The results for absolute cross sections are presented in Figure 7 in the case of √ s = 13 TeV LHCb p-p data, Figure 8 in the case of √ s = 7 TeV LHCb p-p data, Figure 9 in the case of √ s = 5 TeV LHCb p-p data, and Figure 10 in the case of √ s = 7 TeV ALICE p-p data. In addition to our GM-VFNS results, we also show the central prediction from the Powheg+Pythia framework in all panels, and separately compare the scale uncertainties of these two different approaches in two rapidity bins at √ s = 13 TeV in Figure 11.
In all cases, the data are reproduced very well by the GM-VFNS calculations within the considered theory uncertainties, whereas the central values of the Powheg+Pythia calculations systematically fall short of the data. Essentially the same hierarchy has been observed e.g. in LHCb/ALICE papers [1][2][3]62] and elsewhere [23], though the GM-VFNS calculations do not extend to zero P T in these references. We believe that most, or at least D g h Figure 12. A contribution to heavy-quark production resummed by DGLAP equations in GM-VFNS calculation. The diagram has been drawn with JaxoDraw 2.1 [47]. a significant part, of the systematic difference between GM-VFNS and Powheg+Pythia setups can be explained by contributions from gluon fragmentation that are resummed in GM-VFNS but that are not accounted for in the Powheg+Pythia calculation. To better illustrate the point we show, in Figure 12, a diagram where the upper gluon line radiates one or more gluons before splitting into a QQ pair. Contributions of this type are implicitly included in the GM-VFNS calculation, resummed into the gluon FF D g→h 3 in collinear configuration. They correspond to roughly 50% of the D 0 -meson cross sections in P T and Y region investigated here (ALICE and LHCb acceptance). The Powheg+Pythia framework, however, does not allow for these contributions as the starting point are events in which the QQ pair has been produced in the first place. In other words, the possibility that the QQ pair is produced only later on by the Pythia parton shower is possible only for the hard processes where one QQ pair has already been produced in matrix-element level. On the other hand, the standard Pythia simulation with all hard-QCD processes subsumed does include also contributions like those in Figure 12 and, as was shown in the lower panel of Figure 6, enabling this possibility changes the x 2 distributions quite significantly. In the FONLL approach [39], the contributions like those in Figure 12 are resummed by partonic FFs, but the contribution of the resummed part is shrouded by a multiplicative factor G(m, p T ) = p 2 T /(p 2 T + 25m 2 ) engineered so as to suppress these contributions at low p T . As expected, the scale variation in GM-VFNS calculations, best visible in Figure 11, is maximal in the low-P T region and diminishes toward larger values. At the LHCb acceptance,  the scale uncertainty is clearly larger than the variation we found using massless/massive fragmentation variable, see back Figure 3. At low P T , the upper limit for the scale variation is set by the calculation with α fact,frag = 2 and α ren = 1, and the lower border traces the configuration α fact,frag = 0.5 and α ren = 1. Neither of these "extreme" cases reproduces the shape of the data particularly well at low P T : either the spectra rise too steeply (α fact,frag = 2) 3 , or they are too flat (α fact,frag = 0.5). The change of slope of the lower border between P T ∼ 3 GeV and P T ∼ 5 GeV is a consequence of the heavy-quark PDFs and light-parton FFs becoming 'active' as the scale goes above the mass threshold. That is, the contributions shown in the right-hand panel of Figure 2 become suddenly dominant. Below P T ∼ 3 GeV the set condition min(µ i ) > m charm halts the growth of scale uncertainty downwards. The "natural" choice α i = 1 matches much better with the shape of data. The central prediction goes somewhat above the LHCb data at low P T , which could be improved by using a bit lower scale. Indeed, we have checked that setting the scales as (α P T ) 2 + m 2 charm , with the parameter α < 1 [77], would improve the description with the central prediction. However, here our intention is not to fine tune the predictions but rather to present the calculation as it is "out of the box" with default settings. In the Powheg+Pythia case the scales µ ren and µ fact are varied independently by a factor of two wrt. central scale within 0.5 < µ fact /µ ren < 2 in the generation of Powheg cc events. As can be seen in Figure 11, the Powheg+Pythia approach is clearly more sensitive to the scale variations than GM-VFNS, especially at P T > 4 GeV, and the LHCb data still remain within these uncertainties, though at the very upper part of the band. The scale-uncertainty estimates of Powheg+Pythia method we find here are well in line with the error bands shown in the original data papers of LHCb and ALICE, and e.g. in Ref. [66]. Figure 13 presents ratios of cross sections measured by LHCb at different √ s. Here, a significant part of the theory uncertainties cancel and, sure enough, the central predictions of the both considered methods describe the data rather well even at low P T . However, the Powheg results are systematically below the GM-VFNS predictions in most of the cases, best visible in the ( √ s = 13 TeV)/( √ s = 5 TeV) panel. That is, the √ s dependence is stronger in the GM-VFNS calculation. We believe that this hierarchy follows mainly from the presence (absence) of gluon fragmentation in GM-VFNS (Powheg+Pyhtia) approach, similarly as in the case of absolute cross sections. As the Powheg generates events where the QQ pair is produced in the hard process or from the hardest emission, part of the increased phase-space for parton shower due to increased √ s is not in use for heavy-quark production. This is consistent also with the observation that the difference between the GM-VFNS and Powheg+Pythia cross sections decreases, though slowly, towards lower √ s. In Figure 14 we show the ( √ s = 13 TeV)/( √ s = 5 TeV) case without PDF errors and including also the predictions obtained with the zero-mass version of the fragmentation variable (z = P T /p T , y = Y ) to estimate the theoretical bias of this origin. We observe that the differences between our default definition and the zero-mass version are still clearly smaller than the scale uncertainties.
Another observable that has been discussed in the recent literature [6,7] is the normalized cross section dσ/(dP T dY ) dσ/(dP T dY ref ) , where Y ref is a fixed reference rapidity. Also here, large part of the theory uncertainties cancel upon taking the ratio. Our calculations for this quantity within the LHCb acceptance at √ s = 13 TeV in GM-VFNS framework are presented in Figure 15 taking 2 < Y ref < 2.5.
We also compare with the Powheg+Pythia approach and with the GM-VFNS predictions using the zero-mass fragmentation variable. The shown LHCb data points have been formed from the absolute cross sections adding all the uncertainties in quadrature (as if they were uncorrelated). Again, the scale variation appears more significant than the effect induced by using the zero-mass fragmentation variable (z = P T /p T , y = Y ). The Powheg+Pythia results are systematically above the GM-VFNS ones. As earlier, this seems to follow from the absence of gluon fragmentation in Powheg+Pythia: The "missing" contributions from gluon fragmentation are relatively larger for 2 < Y ref < 2.5 than for the more forward bins (since the phase space is larger for 2 < Y ref < 2.5), and thus the ratio is higher than in GM-VFNS.

Summary
In the present article, we have introduced a novel implementation of the GM-VFNS for hadroproduction of heavy-flavoured mesons. Here, the novelty amounts to a specific definition of scheme, SACOT-m T , which retains the kinematics of heavy quark-antiquark pair production also in contributions where the heavy-flavoured meson formally comes from light-flavour fragmentation or initial-state heavy quarks. As we have explained, this is physically a natural choice as the origin of these contributions is in the heavy quark-antiquark pair production, and as such it is analogous to the SACOT-χ scheme in deeply inelastic scattering. Within the SACOT-m T scheme, it is possible to compute the heavy-flavoured meson spectra down to P T = 0 with arbitrary choices for renormalization, fragmentation, and renormalization scales. In earlier works presented in the literature, a finite P T → 0 limit could only be achieved by setting the scales in a particular way. Comparisons with the available D 0 -meson data from LHCb and ALICE collaborations indicate that our calculation in SACOT-m T scheme performs well, though it must be admitted that the theoretical uncertainties at low P T are significant. Here, it is not only the scale, PDF, and FF uncertainties that matter, but also the scheme dependence and ambiguities in defining the fragmentation variable z in the presence of finite heavy-quark and heavy-flavoured meson masses. In particular, as we have shown, the latter can have a significant impact on the shape of the absolute spectra at low P T . Within our definition of the fragmentation variable we have found, however, that the mass effects are suppressed in the forward direction, especially within the LHCb acceptance. Nevertheless, a particular choice of fragmentation variable may still bias the usage of low-P T D-meson production e.g. as a constraint for PDFs. We note that this type of uncertainty is not inherent only to the GM-VFNS approach but the very same ambiguity is there also in FFNS calculations when scale-independent FFs are used for the c → D transition. We have also shown that in all considered rapidities there is a sizeable contribution from large-x 2 region which, as we have argued, appears to originate from gluon fragmentation. Thus, estimates based on fixed-order cc pair production overstate the small-x sensitivity of inclusive D-meson production. Though formally higher order in strong coupling, the effect of gluon fragmentation is numerically large and seems to explain why FFNS-based calculations are typically a factor of two below the LHC data.
In addition, we have observed that an approach (Powheg+Pythia) neglecting large part of the gluon fragmentation deviates systematically from the GM-VFNS predictions also in the case of normalized cross sections and ratios across different √ s. The size of these systematic differences competes or exceeds the scale uncertainty of GM-VFNS, though the scale uncertainties in the Powheg-based setup are presumably somewhat larger also for these observables. In the future, we plan to extend the SACOT-m T scheme also to the case of intrinsic/fitted charm, B-meson production, as well as nuclear collisions.