Heavy-flavor PDF evolution and variable-flavor number scheme uncertainties in deep-inelastic scattering

We consider a detailed account on the construction of the heavy-quark parton distribution functions for charm and bottom, starting from $n_f=3$ light flavors in the fixed-flavor number (FFN) scheme and by using the standard decoupling relations for heavy quarks in QCD. We also account for two-mass effects. Furthermore, different implementations of the variable-flavor-number (VFN) scheme in deep-inelastic scattering (DIS) are studied, with the particular focus on the resummation of large logarithms in $Q^2/m_h^2$, the ratio the virtuality of the exchanged gauge-boson $Q^2$ to the heavy-quark mass squared $m_h^2$. A little impact of resummation effects if found in the kinematic range of the existing data on the DIS charm-quark production so that they can be described very well within the FFN scheme. Finally, we study the theoretical uncertainties associated to the VFN scheme, which manifest predominantly at small $Q^2$.


I. INTRODUCTION
The process of deep-inelastic scattering (DIS) of leptons off a nucleon target provides important information on the nucleon structure and the parton content. Therefore, it plays a central role in the determination of the parton distribution functions (PDFs), especially for the proton PDFs [1]. At large values of Bjorken x the DIS data constrain the valence quark distributions, while at small x they are sensitive to the sea-quark and gluon distributions. In addition, the DIS cross sections at small x contain substantial contributions from charm and bottom quarks. The virtuality Q 2 of the exchanged gauge-boson is the other important kinematic variable in DIS. It offers a wide range of scales to probe, for instance, in electron-proton scattering the parton dynamics inside the proton. Depending on the value of Q 2 , different theoretical descriptions of DIS within Quantum Chromodynamics (QCD) may be applied. This concerns in particular the number n f of active quark flavors and the treatment of the heavy quarks, as charm and bottom.
At low scales, when Q 2 is of the order of the heavy quark mass squared m 2 h , one typically works with n f = 3 massless quark flavors. Then, the proton structure function is composed only out of light-quark PDFs for up, down and strange and of the gluon PDF. Massive quarks appear in the final state only or contribute as purely virtual corrections. At higher scales, for Q 2 m 2 c , m 2 b compared to the charm and bottom quark masses squared, additional dynamical degrees of freedom lead to theories with n f = 4 or 5 effectively light flavors, depending on whether charm is considered massless, or even both, charm and bottom. The massive renormalization group equations rule these dynamics and provide the corresponding scale evolution, linking to n f = 3 massless quarks at very low virtualities. In general, the transition for the flavor dependence of the strong coupling α s , i.e., α s (n f ) → α s (n f + 1), is achieved at some matching scale µ 0 with the decoupling relations [2] which can be implemented perturbatively in QCD. These decoupling relations introduce a logarithmic dependence on the heavy quark masses m c , m b . In a similar manner, this is realized for the PDF f i of a parton i with the help of suitable heavy-quark operator matrix elements (OMEs) [3,4], which, in the perturbative expansion, also depend logarithmically on the heavy-quark masses. The transition f i (n f ) → f i (n f +1), again at a matching scale µ 0 , implies also the introduction of new heavy-quark PDFs for charm or bottom when they become effectively light flavors, and can then be considered as effective dynamical degrees of freedom inside the proton.
For a given fixed value of n f , and having decoupled the heavy quarks in an appropriate manner, one may define the fixed-flavor number (FFN) scheme. In the FFN scheme used in the ABMP16 global fit of proton PDFs [5], only light quarks and gluons are considered in the initial state, while heavy quarks appear in the final state as a result of the hard scattering of the incoming massless partons. Existing data on the heavy-quark DIS production are well described by the FFN scheme with n f = 3, see Refs. [1,6]. However, many PDF fits, like those of CT18 [7], MMHT14 [8] and NNPDF3.1 [9] employ various different versions of the so-called variable-flavor-number (VFN) scheme. In the VFN scheme the quark flavors charm and bottom are considered also in the initial state from a certain mass scale onward and are dealt with as partonic components in the proton. As a consequence, the original distributions f i (n f ) are mapped into the distributions f i (n f + 1) at a chosen scale µ 0 , cf. [3]. In addition, the VFN scheme effectively performs a resummation of logarithms in the ratio Q 2 /m 2 c (or Q 2 /m 2 b ) through the parton evolution equations for the charm (or bottom) PDF [10], although the corresponding logarithms are not necessarily large for realistic kinematics.
The difference in modeling of the heavy quark contribution, i.e., the choice for the FFN or the VFN scheme, has an impact on the PDFs obtained in global fits [1,11]. Therefore, a detailed comparison of the two approaches is mandatory in view of the use of the respective PDF sets in QCD precision phenomenology.
A particular prescription for a VFN scheme has been proposed in [3], commonly referred to as BMSN-scheme. In PDF fits, the VFN scheme using this approach yields results which are not very different from the ones in the FFN scheme [12]. This happens due to a smooth transition between the n f -and (n f + 1)-flavor regimes at the matching scales µ 0 = m c and µ 0 = m b , respectively, which is imposed in the BMSN ansatz. However, the BMSN prescription is based on heavy-quark PDFs, i.e., charm and bottom, which are derived with the help of fixed-order matching conditions. Therefore, the results of our previous study [12] cannot be directly compared to the PDF fits in Refs. [7][8][9] which apply heavy-quark PDF evolution.
In the present article we study the phenomenology of a modified BMSN prescription, which also includes the scale evolution of heavy-quark PDFs in order to clarify the basic features of such VFN schemes. Our studies are limited to the case of DIS charm-quark production, since this process is most essential phenomenologically and, at the same time, a representative case.
The paper is organized as follows. Basic features of QCD factorization, the VFN scheme and the BMSN prescription are outlined in Sec. II. In Sec. III we describe the particularities introduced by the heavy-quark PDF evolution and Sec. IV contains the benchmarking of various factorization schemes based on existing data for DIS charm-quark production. We address implications of VFN schemes for predictions at hadron colliders Sec. V and conclude in Sec. VI. Technical details of the various implementations of heavy-quark schemes are summarized in App. A.

II. HEAVY-QUARK PDFS
The dynamics of massless partons in the proton are parameterized in terms of the PDFs f i with i = u, d, s, g for up, down, strange quarks and the gluon. These define the set of flavor-singlet quark and gluon PDFs, q s and g, where µ denotes the factorization scale and we suppress the dependence on the momentum fractions x here and below. Using standard QCD factorization 1 , the PDFs for the heavy quarks charm and bottom (h = c, b) at the scale µ in the MS scheme and using on-shell renormalization for the mass m h are then constructed from the quark-singlet and gluon PDFs in Eq. (1) and the heavy-quark OMEs A i j as follows [3,4] f h+h (n f + 1, µ 2 ) = A ps hq n f , where h = c, b and '⊗' denotes the Mellin convolution in the momentum fractions x. Typically, the matching conditions are imposed at the scale µ 0 = m h , and we further assume that f h+h = 0 at scales µ ≤ m h . In addition, the transition {q s (n f ), g(n f )} → {q s (n f + 1), g(n f + 1)} for the set of the light-quark singlet and the gluon distributions with the respective heavy-quark OMEs has to account for operator mixing in the singlet sector g(n f + 1, with h = c, b, see Refs. [3,4], also for matching relations for the non-singlet distributions. The perturbative expansion of the OMEs in powers of the strong coupling constant α s reads (using a s = α s /(4π) as a short-hand), where the expressions a (k,0) i j contain the information, which is genuinely new at the k-th order. The leading-order (LO) and next-to-leading order (NLO) contributions to the OMEs are given by the coefficients at order a s and a 2 s in Eq. (5), respectively. They have been determined analytically in closed form in Refs. [3,[13][14][15] 2 . At next-to-next-to-leading (NNLO) the heavy-quark OMEs are known either exactly or to a good approximation [4,5,[16][17][18][19]. This includes specifically the non-singlet and pure-singlet constant parts a (3,0) hq in Eq. (5) and the term a (3,0) hg at order a 3 s . In the latter case an approximation based on fixed Mellin moments [4] with a residual uncertainty in the small-x region has been given in Ref. [5,17].
It should be noted that, the decoupling relations in Eqs.
(2)-(4) assume the presence of a single heavy quark at each step only. Thus, the bottom-quark contributions are ignored in the transition from n f = 3 to 4 and in the construction of the charm-quark PDF. However, starting at two-loop order, the perturbative corrections to the heavy-quark OMEs contain graphs with both, charm-and bottom-quark lines. With the ratio of masses (m c /m b ) 2 ≈ 1/10, charm quarks generally cannot be taken as massless at the scale of the bottom-quark. Such two-mass contributions to the heavy-quark OMEs have been computed recently [20][21][22]. At three-loop order (and beyond), these corrections can neither be attributed to the charm-nor to the bottom-quark PDFs separately. Rather, one has to decouple charm and bottom quarks together at some large scale and the corresponding VFN scheme, i.e. the simultaneous transition with two massive quarks, f i (n f ) → f i (n f + 2), has been discussed recently in Ref. [23]. This proceeds in close analogy to the simultaneous decoupling of bottom and charm quarks in the strong coupling constant α s , see for instance Ref. [24]. We will elaborate on these aspects further below.
First, we will limit our studies to the case of the charm-quark PDF and apply Eqs. (2)-(4) to change from n f = 3 to 4. At LO only the heavy-quark OME A s hg contributes and the coefficients are i.e., the constant term of the unrenormalized massive OME A s, (1) hg vanishes and the logarithmic one with T f = 1/2 is proportional to the LO quark-gluon splitting function P (0) qg in the normalization of Ref. [25]. For four active flavors we abbreviate the charm PDF in Eq. (2) as c(x, µ 2 ) ≡ f c+c (4, x, µ 2 ) and consider its perturbative expansion where the LO term c (1) (x, µ 2 ) has a particularly simple form, Here, g denotes the gluon PDF in the 3-flavor scheme. This expression is used in the BMSN prescription [3] of the VFN scheme and determines the charm-quark distribution at all scales µ ≥ m c at fixed-order perturbation theory (FOPT). On the contrary, other VFN prescriptions, like ACOT [26][27][28], FONLL [29] or RT [30] use Eq. (8) as a boundary condition for c(x, µ 2 ) at µ = m c and derive the scale dependence with the help of the standard QCD evolution equations (DGLAP) for massless quarks. The evolution resums logarithmic terms to all orders, so that the charm-quark distribution acquires additional higher order contributions which are not present in the FOPT expression in Eq. (8). In order to illustrate the numerical difference between these two approaches, we consider the derivative of c(x, µ 2 ) whereġ(x, µ 2 ) ≡ dg(x, µ 2 )/d ln µ 2 . The first term in Eq. (9) corresponds to the right hand side of the standard DGLAP evolution equations, recall Eq. (6), i.e., a (1,1) hg is proportional to P (0) qg . The second and the third term, however, account for the difference between the FOPT distributions and the evolved ones. These terms vanish at the matching scale µ 0 = m c as they should by definition. For scales µ > m c the second term proportional to the QCD β function is negative, since da s /d ln µ 2 = β(a s )/(4π) < 0. However, the net effect of the difference between the FOPT and the DGLAP evolved distributions shown in Fig. 1 on the left is positive at small x and driven byġ in the third term. Only at large x, where the gluon PDF is negligible, the term proportional to β(a s ) dominates and the net difference between the FOPT and the DGLAP evolved distributions is negative.
The matching conditions for the charm-quark at NLO are more involved. The NLO term c (2) (x, µ 2 ) in Eq. (7) has the form It includes the NLO corrections to the massive OMEs A s, (2) hg and A ps, (2) hq for n f = 3, see Eq. (5), the gluon and the quark-singlet PDFs, g and q (s) , are taken in the 3-flavor scheme again, cf. Eq. (1). Since a (2,0) hg and a (2,0) hq in Eq. (5) are non-zero in the MS scheme, c(x, µ 2 ) at NLO does not vanish anymore at the matching scale µ 0 = m c , see the off-set in Fig. 1 on the right. The comparison of the charm-quark FOPT distributions at NLO based on Eqs. (8) and (10) and the evolved ones, using c(x, µ 2 ) only as the boundary condition at the matching scale, shows in Fig. 1 qualitatively the same pattern as at LO, although the numerical differences are smaller now. At small x, driven by the scale derivativeġ of the gluon PDF, the FOPT distributions are larger while at large x the terms proportional to β(a s ) dominate and the DGLAP evolved distributions are larger. These observations can be expressed in quantitative form through the scale derivative of the NLO term c (2) (x, µ 2 ), which reads where, again,ġ(x, Here, the first two terms in the right hand side contain the expressions used in the standard DGLAP equations to evolve the charm-quark PDF, since the NLO splitting functions P (1) qg and P (1) qq appear in the terms a (2,1) hg and a (2,1) hq , cf. [3,4]. However, there are also other contributions, since the heavy-quark OMEs enjoy their own (massive) renormalization group equation. In addition, the full expression dc(x, µ 2 )/d ln µ 2 at NLO contains, of course, also the terms fromċ (1) (x, µ 2 ) in Eq. (9) expanded to higher order in a s , for example the term proportional to β(a s ). In summary, these terms are responsible for decreasing the difference between the FOPT and the evolved distributions at NLO in Fig. 1.
As a further variant in the study of the DGLAP evolved charm-quark PDF, one can perform the evolution using the full NNLO splitting functions P (2) i j of Ref. [25] starting at the matching scale µ 0 = m c from the boundary condition for c(x, m 2 c ) at NLO in Eqs. (8) and (10). We denote this variant as NNLO * , since there is a mismatch in the orders of perturbation theory between the heavy-quark OMEs and the accuracy of the evolution equations. The difference with the NLO variant is due to terms which are formally of higher order, but nevertheless have significant numerical impact at small x as shown in Fig. 1. There, the FOPT distributions at NLO and the evolved ones at NNLO * accuracy are very similar in the entire µ-range. Only at large x, the increased order in the DGLAP evolution is negligible.
In Fig. 2 we display the scale derivatives of the charm-quark PDFċ(x, µ 2 ) ≡ dc(x, µ 2 )/d ln µ 2 calculated using Eqs. (9) and (11). We consider the difference ofċ(x, µ 2 ) determined in FOPT, c FOPT , and the one evolved with the standard DGLAP evolution,ċ evol , choosing n f = 4 and starting from the expressions in Eqs. (8) and (10) at the matching scale µ 0 = m c = 1.4 GeV. Evidently, at LO the differenceċ FOPT −ċ evol has to vanish at the matching scale, while at NLO or in the NNLO * variant some finite off-set at µ 0 = m c = 1.4 GeV remains. Remarkably, the results at NLO and at NNLO * , i.e., using NLO boundary conditions from Eqs. (8) and (10) and NNLO splitting functions in the evolution ofċ evol are very different at low factorization scales and only converge above µ 2 10 2 . . . 10 3 GeV 2 , depending on the value of x. These large scales, however, at which the NLO and the NNLO * variants become of similar size, are typically well outside the kinematic range of the HERA collider, whose upper limit is indicated by the vertical arrow. These findings indicate that there is a substantial numerical uncertainty in the VFN prescriptions ACOT [26][27][28], FONLL [29] or RT [30] due to the order of the QCD evolution applied. In particular the additional higher order terms in the NNLO * variant do have a sizable effect within the µ-range covered by experimental data on DIS charm-quark production and, hence, on the quality of the description of those data in a fit using such VFN prescriptions.
An additional source of uncertainty in the VFN scheme concerns the choice of the matching scale µ 0 . Conventionally it is set to the corresponding heavy-quark mass, m c and m b for the 4and 5-flavor PDFs, respectively. A variation of µ 0 leads to a modification of the shape of the evolved heavy-quark PDFs, whereas, in contrast, the FOPT ones remain unchanged by construction. Therefore, for µ 0 > m h a difference between the FOPT and the evolved heavy-quark PDFs is generally becoming smaller, in particular within the phase-space region covered by existing data, cf. Fig. 3. Such a variation of µ 0 also implies the use of the FFN scheme to describe data in a wider kinematic range, e.g., up to µ 0 = 2m c instead of µ 0 = m c for the illustration in Fig. 3. Therefore, the uncertainty due to a variation of matching scale is not completely independent from those related to the choice of heavy-quark PDFs employed in the VFN scheme.

III. BMSN PRESCRIPTION OF THE VFN SCHEME
The heavy-quark distribution derived using the matching conditions Eqs. (2), (3) enter the zeromass VFN scheme (ZMVFN) expression for F 2,h where C (k) 2,i are the massless DIS Wilson coefficients at the k-th order, which are known to nextto-next-to-next-to-leading order (N 3 LO) [31]. This expression is valid at asymptotically large  momentum transfer Q 2 m 2 h , while it is unsuitable for scales Q 2 m 2 h since the heavy-quark decoupling is not applicable. Therefore, a realistic implementation of the VFN scheme commonly includes a combination of the ZMVFN expression in Eq. (12) with the FFN one where H (k) 2,i are the Wilson coefficients for the DIS heavy-quark production, all known exactly at NLO [32] and H (3) 2,g to a good approximation at NNLO [5,17]. Furthermore, in order to avoid double counting, a subtraction has to be carried out when combining Eqs. (12) and (13). For the BMSN prescription of the VFN scheme [3] this subtraction arises from the asymptotic FFN expression as follows where H (k),asy 2,i is derived from H (k) 2,i taken in the limit of Q 2 m 2 h . In summary BMSN prescription then reads where a factorization scale µ F = m h is used throughout. The asymptotic Wilson coefficients H asy 2,i can be expanded into a linear combination of the massless Wilson coefficients C 2,i and the massive OMEs [3,4,13]. For this reason, the asymptotic expression of Eq. (14) coincides with the ZMVFN one of Eq. (12), when the FOPT matching conditions Eqs. (2)-(4) are employed, up to the subleading non-singlet terms and the difference between a k s (n f + 1) and a k s (n f ) [12]. The latter exhibits a small discontinuity at Q 2 m 2 h beyond  . In summary, the BMSN prescription Eq. (15) provides a smooth transition between FFN scheme at small momentum transfer to the ZMVFN scheme at large scales, cf. Fig. 4.
A version of the BMSN prescription based on the NLO evolution of the (n f + 1)-flavor PDFs also allows for a smooth matching with the FFN scheme at Q 2 = m 2 h , because the NLO-evolved PDFs do not have a discontinuity with respect to the FOPT ones at the matching point, cf. Figs. 1 and 5. For the variant denoted NNLO * which uses NNLO-evolved PDFs the trend is different: the slope of F 2,h at Q 2 = m 2 h predicted by the BMSN prescription is substantially larger than the one obtained with the FFN scheme. This is in line with the difference between the NNLO and FOPT PDFs. Obviously, this difference is not explained by the impact of the resummation of large logarithms, but rather by the mismatch in the perturbative order of the matching conditions and evolution kernels employed to obtain the NNLO heavy-quark PDFs. Therefore, the difference between the NLO and NNLO * variants of the VFN scheme should essentially quantify its uncertainty due to the missing NNLO corrections to the massive OMEs. A choice of the matching scale µ 0 = m h , i.e., at the mass of the heavy quark, is a matter of convention rather than appearing as a consequence of solid theoretical arguments. Also note, that for DIS charm production, the matching scale µ 0 cannot be significantly shifted to scales much lower than m c , because in this case the matching would be performed at scales well below 1 GeV, where QCD perturbative expansions are not converging anymore. When µ 0 is shifted upwards, e.g., µ 0 = 2m h , the difference between the NLO and NNLO * variants of the VFN scheme is becoming less significant. This is particularly due to the fact, that then essential parts of the problematic small-Q 2 region are left for a theoretical description within the FFN scheme, cf. Fig. 6.
The impact of scheme variations and the choice of the matching scale are qualitatively similar for the c-and b-quark production. Nonetheless, the effects are less pronounced for the b-quark case [35], mainly because of the smaller numerical value of strong coupling at the scale m b . For this reason, and also due to more representative kinematics of data, all our phenomenological comparisons are focused on the c-quark contribution.

IV. BENCHMARKING OF THE FFN AND VFN SCHEMES WITH THE HERA DATA
To study the phenomenological impact of the VFN scheme uncertainties we consider several variants of ABMP16 PDF fit [5], which include the recent HERA data on heavy-flavor DIS production [6]. Furthermore, the inclusive neutral-current DIS HERA data used in the ABMP16 fit are excluded in order to illuminate the impact of the scheme variation on the PDFs extracted from the fit. For the same reason we exclude the collider data on W ± -and Z-boson production, which provide an additional constraint at small-x on the PDFs in the ABMP16 fit. However, in order to keep the different species of quark flavors disentangled, we add data on DIS off a deuteron target, analogous to an earlier study in Ref. [36]. For all variants we employ the NLO massive Wilson coefficients [32] and the pole-mass definition for the heavy-quark masses, so that a consistent comparison of the FFN scheme with the original formulation of the BMSN prescription and its modifications is possible. For the same purpose we take the factorization scale µ With these settings, the FFN scheme provides a good description of the c-quark production data, cf. Fig. 7. The agreement of the fit with the data is equally good, both at small and at large Q 2 , underpinning the fact, that any additional large logarithms cannot improve the theoretical data treatment within the range of kinematics covered by the HERA data. This observation is indeed long known [38].
In order to check this aspect in greater detail we also compare predictions of various versions of the VFN scheme with the data. Let us consider the VFN predictions for the heavy-quark production cross sections which are computed by using the BMSN prescription of Eq. (15) for F 2 , while still keeping the FFN scheme for F L . The justification of this approach derives from the small numerical contribution of F L as compared to F 2 . In addition, the modeling of F L within the VFN framework is conceptually problematic [12], because the effects of power corrections in m 2 h /Q 2 cannot be disregarded for this observable [13]. The PDFs used in this comparison are the ones obtained in the FFN version of the fit. Therefore the obtained pulls display the impact of the scheme variation only.
As expected, predictions of the VFN scheme based on the BMSN prescription and the FOPT heavy-flavor PDFs are close to the FFN ones. The same applies to the case of NLO-evolved PDFs, which are smoothly matched with the FFN ones at small scales, cf. Fig. 5. In contrast, an excess with respect to the small-Q 2 data appears for the variant of the fit with the NNLO * PDFs employed. This excess is clearly related to the mismatch between the FFN scheme and this variant of the VFN one. At large Q 2 the impact of the resummation of large logarithms is marginal, in particular given the size of the data uncertainties. The latter is true also for the case of NLO-evolved PDFs.
The HERA data on c-quark production used in the present analysis are accurate enough to provide a sensible constraint on the small-x gluon distribution. Moreover, the latter demonstrate sensitivity to the choice of the factorization scheme, cf. Fig. 8. The FFN scheme and the BMSN scheme with the FOPT and the NLO-evolved PDFs are in qualitative agreement, while a much lower small-x gluon distribution is preferred in the variant based on the NNLO * PDFs. This is in line with the trends observed for the pull comparison, cf. Fig. 7.
The difference between gluon and quark distributions obtained in the NLO-and NNLO * -based fits is pronounced at small x due to kinematic correlations with the small-Q 2 region, where the difference between these two approaches is localized, and reaches ∼ 30% at x = 10 −4 . The description of the small-x inclusive DIS data is also sensitive to the scheme choice due to a substantial contribution of the heavy-quark production. In order to check this quantitatively, we consider variants of the fits with various VFN scheme prescriptions and the HERA inclusive data [39] added. In line with the recent update of the ABMP16 fit [40], we impose strong cuts on the momentum transfer  Q 2 > 10 GeV 2 and on the hadronic mass W 2 > 12.5 GeV 2 , which allow to avoid any impact of higher twist corrections, cf. [41]. The PDF uncertainties are improved due to the additional data included. However, the sensitivity of the resulting gluon distribution to the choice of heavy-quark PDF evolution still reaches ∼ 30% at x = 10 −4 , cf. Fig. 9. Such a spread induces quite essential uncertainties in the small-x VFN predictions, in particular in the c-and b-quark input distributions for scattering processes at hadron colliders. The gluon distribution obtained using the BMSN prescription with the NLO-evolved PDFs is increased with respect to the FOPT one at x ∼ 0.01, which gives a hint on the impact of the resummation of large logarithms at these kinematics. No further substantial change in the gluon distribution at x 0.01 is observed, when the NNLO corrections to the evolution are taken into account. Therefore, one should expect a minor impact of the logarithmic terms at higher order (higher powers) on the description of the existing DIS data, although the comparison is somewhat deteriorated by the uncertainty from the mismatch in perturbative orders in the NNLO * fits appearing at x 0.01. In this context it is also instructive to consider the results of the FFN fit performed with account of the NNLO corrections, which include the terms up to O(ln 2 (µ 2 /m 2 h )) [5] and MS masses m c (m c ) = 1.27 GeV, m b (m b ) = 4.18 GeV [42]. The gluon distribution obtained with these settings is similar to the VFN ones at x 0.01, located between the NLO and NNLO * fit results at x ∼ 10 −4 and lower by ∼ 5% than both of these variants at x ∼ 0.01, where they agree with each other, cf. Fig. 9. This plot also yields an upper limit on the estimate of the impact of missing large logarithms in the NNLO FFN fit. On the other hand, a comparison with the NNLO VFN fit at small x is inconclusive due to the large uncertainties in the VFN scheme appearing at these kinematics. A more accurate estimate requires the NNLO VFN fit with a consistent boundary condition based on OMEs at NNLO accuracy [4,5,[16][17][18][19]. Nonetheless, at the present level of data accuracy this upper limit is comparable with the experimental uncertainties in the gluon distribution obtained from the fit.
Finally, considering a variation of the matching scale for the 4-flavor PDFs from µ 0 = m c to µ 0 = 2m c , leads to VFN heavy-flavor predictions being closer to the FFN ones, cf. Fig. 6. The phenomenological effect of such a variation is more substantial at small Q 2 and x due to kinematic characteristics of the existing DIS experiments. Therefore, the corresponding change of the gluon distribution due to a matching point variation is significant mostly at x 10 −3 , cf. Fig. 10. It is comparable in size with the VFN scheme uncertainty related to the boundary conditions for the evolution. However, strictly speaking, these two uncertainty sources should not be considered independently since the impact of the matching scale variation also manifests itself through the scheme change.

V. IMPLICATIONS OF VFN SCHEMES FOR PREDICTIONS AT HADRON COLLIDERS
The contribution of heavy flavors to the hadro-production of massive states, like W ± -, Z-and Higgs-bosons, t-quarks, etc., are commonly taken into account within the 4-or 5-flavor scheme. This allows for great simplifications of the computations, since the VFN PDFs employed in this case contain resummation effects, which are generally rising with the factorization scale, cf. Fig. 1. Therefore, the VFN scheme provides a relevant framework for the phenomenology of heavy particle hadro-production.
The NNLO 4-and 5-flavor PDFs still suffer from the uncertainty due to the yet unknown exact NNLO corrections to the massive OMEs. Moreover, for the NNLO PDFs derived from the VFN fit including the small-x DIS data this uncertainty is enhanced, since the part of those DIS data, which provides an essential constraint on the PDFs, also populates the matching region. The observed spread in the 5-flavor gluon distributions, which are obtained from the VFN fits with varying treatments of the matching ambiguity is somewhat reduced with increasing scales due to the general properties of the QCD evolution. However, it is still comparable to experimental uncertainties at x ∼ 0.01 and substantially larger at x ∼ 10 −4 , cf. Fig. 11.
Altogether, this implies an uncertainty in predictions of the production rates of the Higgsboson and t-quark pairs at the Large Hadron Collider (LHC) within a margin of few percent and somewhat larger at the higher collision energies discussed for future hadron colliders. Note that in the ABMP16 fit [5], which is based on a combination of both, DIS and hadron collider data, the FFN and the 5-flavor VFN schemes are used for the theoretical description of these samples, respectively. This allows to keep the advantages of the VFN scheme at large scales, while avoiding its problems concerning the DIS data. Nevertheless, the NNLO massive OMEs [4,5,[16][17][18][19] are still necessary to generate NNLO PDFs free from the matching ambiguity.
In closing the studies of VFN schemes we wish to address a conceptual problem of the 5flavor scheme definition due to the fact, that the b-quark mass m b is not too much larger than m c . This relates to the inherent limitations of the VFN schemes due to the successive decoupling of one heavy quark at a given time. As discussed above, starting from the two-loop order, the DIS structure functions also receive contributions which contain two different massive quarks [23]. At two loops, they are given by one-particle reducible Feynman diagrams, while one-particle irreducible graphs appear at the three-loop order for the first time, cf. [20][21][22].
Here we will consider the two-loop effects, which arise from virtual corrections with both, charm and bottom quarks. Thus, no production threshold is involved. For the structure function F 2 one obtains which is to be added to Eqs. (12) or (13), with e h denoting the fractional heavy-quark charge and using T F = 1/2. The effect of the 2-mass contributions rises at small x and large Q 2 , being more pronounced for the case of b-quark production, cf. Fig. 12. For the kinematics of the proposed lepton-proton LHeC collider it reaches up to ∼ 3%, which has impact on the phenomenology of heavy-quark production.
As demonstrated in Ref. [23], the two-mass diagrams at the two-loop order have the largest effects for the b-and c-quark distribution at large Q 2 . The respective PDFs can be obtained by adding the two-mass contributions to the OMEs in Eq. (2). Comparing the heavy-quark PDFs with and without the two-mass effects included, one finds that the relative size of the effect is negative: b-quark distributions with the two-mass contributions included are decreased by -2% to -6% in the range for Q 2 from 30 to 10000 GeV 2 at small x, x = 10 −4 ; c-quark distribution the relative variations are smaller, amounting to -1% to -4% for Q 2 = 100 GeV 2 to 10000 GeV 2 and x = 10 −4 . In precision fits these two-mass effects have consequences for all PDFs and require the use of a different VFN scheme compared to those with the decoupling of a single heavy quark at the time, cf. [23]. At this point, however, we leave detailed studies of VFN schemes with two massive quarks, i.e., the simultaneous transition f i (n f ) → f i (n f + 2) for PDFs for future studies.

VI. CONCLUSIONS
The precise description of the parton content in the proton across a large range of scales is a an important ingredient in precision phenomenology. The treatment of heavy quarks with a mass m h requires adapting the number of light flavors in QCD to the kinematics under consideration, set by the factorization scale µ, which is typically associated with the hard scale of the scattering process. Within the ABMP16 global PDF fit, the FFN scheme with n f = 3 light flavors provides a good description of the existing world DIS data, while the LHC processes are typically described with n f = 5 massless flavors by implementing decoupling of heavy quarks and a transition from 3to 4-or 5-flavor PDFs, including the possibility for the resummation of large logarithms in Q 2 /m 2 h . To check the effects of such a resummation on the analysis of existing DIS data we have studied the c-quark PDF, constructed with the help of massive OMEs in QCD, and we have quantified differences between the use of perturbation theory at fixed order and subsequent evolution. We have found that the impact of the PDF evolution as used in the BMSN prescription of VFN scheme is sizable and rather x-dependent than Q 2 -dependent, showing little impact on the large-log resummation on the heavy-quark production at realistic kinematics. Moreover, these differences must be considered an inherent theoretical uncertainty of VFN schemes since using NLO or NNLO accuracy for the evolution leads to significantly different results due to mismatch in the orders of perturbation theory between the heavy-quark OMEs and the accuracy of the evolution equations. Likewise, and related, the choice of the matching point position employed in the VFN schemes has the impact on heavy-quark PDFs and therefore brings additional uncertainty.
With the help of variants of the ABMP16 PDF fit, we have confronted the FFN scheme and different realizations of VFN schemes (FOPT, evolved at NLO, evolved at NNLO) in the BMSN approach with the combined HERA data and DIS c-quark production. The FFN scheme delivers a very good description of those data and we have found little need for the additional resummation of large logarithms in the kinematic range covered by HERA. From the fit variants, we have also determined the gluon and the total light-flavor sea quark distributions, illustrating again the sizable numerical differences, obtained by adopting the respective VFN scheme variants. Depending on the value of x, the observed differences for the gluon PDF are well outside the experimental uncertainties at low factorization scales and persist as well as at high scales of O(100) GeV. The VFN scheme choices are, therefore, highly relevant for LHC phenomenology and affect the predictions for the hadro-production of massive particles within a margin of few percent.
In summary, despite being applicable in a limited kinematic range, the FFN scheme works very well for the modern PDF fits and contains much smaller theoretical uncertainty than the VFN schemes currently available. As an avenue of future development, the latter will benefit from improving the perturbative accuracy of the massive OMEs used, including their NNLO corrections, which are known exactly or to a good approximation. Other features of VFN schemes to be improved concern the simultaneous decoupling of bottom and charm quarks, which is advisable due to the close proximity of the mass scales m b and m c . We leave these issues for future studies. the convolutions to x 1 + 4m 2 h Q 2 ≤ z ≤ 1 with the Bjorken variable x. The slow rescaling is motivated by its properties to model energy conservation in the DIS production of heavy final states. Ref. [44] also explores a wider family of rescaling choices, which interpolate smoothly between z and χ.

MSTW
MMHT14 [8] uses the RT VFN scheme [30], specifically the TR' prescription from Ref. [45] for PDF fits at NNLO. The RT scheme requires as a constraint the continuity of physical observables in the threshold region, i.e., for the expression for F FFN 2,h in Eq. (13) below and F ZMV FN 2,h in Eqs. (12) above threshold. To that end, the derivative of the structure function, dF 2 /d ln Q 2 is supposed to be continuous at the matching point Q 2 = m 2 h in the gluon sector. To achieve this modeling constraint, a Q 2 -independent term is added above the matching point to the expression for F Z MV FN 2,h to maintain continuity of the structure function. The TR' prescription specifies this procedure up to NNLO [45]. NNPDF NNPDF3.1 [9] uses the FONLL VFN scheme [29], which has been devised to combine the heavy-quark DIS structure functions and the ZMVFN expressions in analogy to Eq. (15). FONLL suppresses the difference of F ZMV FN 2,h in Eq. (12) and the necessary subtraction term, i.e., the expression analogous to F asy 2,h in Eq. (14), which is needed to avoid double counting, with a kinematical damping factor 1 − Q 2 m 2 2 . In this manner, it is guaranteed, that only F FFN 2,h of Eq. (13) remains for virtualities Q 2 m 2 h near threshold. The variant FONLL-C is used to determine the PDFs at NNLO [29].