Heavy-flavor parton distributions without heavy-flavor matching prescriptions

We show that the well-known obstacle for working with the zero-mass variable flavor number scheme, namely, the omission of O(1) mass power corrections close to the conventional heavy flavor matching point (HFMP) mu_b=m, can be easily overcome. For this it is sufficient to take advantage of the freedom in choosing the position of the HFMP. We demonstrate that by choosing a sufficiently large HFMP, which could be as large as 10 times the mass of the heavy quark, one can achieve the following improvements: 1) above the HFMP the size of missing power corrections O(m) is restricted by the value of mu_b and, therefore, the error associated with their omission can be made negligible; 2) additional prescriptions for the definition of cross-sections are not required; 3) the resummation accuracy is maintained and 4) contrary to the common lore we find that the discontinuity of alpha_s and pdfs across thresholds leads to improved continuity in predictions for observables. We have considered a large set of proton-proton and electron-proton collider processes, many through NNLO QCD, that demonstrate the broad applicability of our proposal.


Introduction
It is well established that a proton p undergoing an inelastic collision at a scale Q > m p will exhibit non-trivial heavy quark content. Heavy quarks are the ones whose mass m is large enough that the strong coupling at that scale is perturbative, i.e. α s (m) 1. In practice, only the charm and bottom quarks are considered massive while the top quark is too massive to be relevant for currently accessible collider energies (which might change at future high energy colliders). In this work we consider the simplified situation of a single massive flavor which for convenience we choose to be the bottom. We will return to the issue of including charm and top in the Conclusions.
This work is based on the collinear factorization approach [1][2][3] within which one identifies two types of contributions to the heavy quark parton distribution function (pdf): a perturbative and an intrinsic one. In this study we will not consider the intrinsic component since it is very small and, in fact, all recent attempts [4][5][6][7][8] to derive the intrinsic charm component of the proton from experimental data have been consistent with vanishing intrinsic charm (intrinsic bottom would be even smaller). We note that not including intrinsic heavy flavor component is mostly for convenience, however, and this is a problem which is orthogonal to the scope of the current work. Therefore, for the rest of this work we will simply speak of heavy flavor pdfs and will always have in mind the perturbative component.
The heavy flavor decoupling theorem [9][10][11] provides the natural framework for discussing heavy flavor pdfs. When a quark of mass m is probed at a scale Q one has two unambiguous limiting behaviors: if Q m then the quark is effectively massless; all O(m) contributions can be neglected, large quasi-collinear logs ∼ ln n (Q/m) are resummed by the DGLAP evolution and the quark behaves exactly as the light quarks. In the opposite limit, when Q m the quark is very heavy and can be integrated out from the theory. Corrections behaving like O(1/m) are neglected and, effectively, the heavy quark disappears from the theory.
There exists a substantial interval of scales Q not very different from the mass m, Q ∼ m, which is not covered by the two above-mentioned asymptotic regimes and where power corrections of m are important. The theory provides us with only a minimal guidance about this intermediate region and it is this inherent ambiguity which is the subject of this work as well as the vast majority of past works on including heavy flavors in pdfs.
In the interval of energies from Q ∼ m down to Q m the dependence on the mass m can be reconstructed, at least at the conceptual level, in the following way. In any observable, there are three basic ingredients where the dependence on the mass m appears: the pdfs, the strong coupling α s and the partonic cross-section dσ. In the limit Q m the pdfs and α s evolve with n f = n l active flavors. As emphasized extensively in ref. [3], since they both are renormalized in the MS scheme -which is mass independent -one can extend their evolution all the way up to scales Q ∼ m. On the other hand the mass dependence in dσ can also be included by computing all diagrams that include the heavy quark in the final state. This way, one can have a prediction with full dependence on the mass m from scales Q ∼ m down to Q m. This is the usual fixed-flavor (FF) picture for describing charm, bottom and top production.
The difficulties arise when one tries to extend this result towards high energies Q m. Following the decoupling approach discussed above, at some scale Q ≡ µ b ∼ m one switches from the n f = n l description to a n f = n l + 1 one, where the massive flavor is now considered as massless. In the following we will refer to the point µ b as heavy flavor matching point (HFMP). 1 We will avoid calling it threshold to avoid confusion with physical thresholds. For all scales Q > µ b then one evolves the pdfs and α s with n l + 1 flavors and introduces pdfs for the heavy flavor and anti-flavor. As usual, matching relations for the coupling and pdfs at the scale µ b have to be imposed. Those are fully known through NNLO [12][13][14] and partly known beyond [15][16][17][18] for the pdfs, and through four loops [19,20] for α s . The above construction represents what is known as a variable flavor number scheme (VFNS) [21] or, specifically, a zero-mass one (ZM-VFNS). The ZM-VFNS described above is consistent. However, it has an important shortcoming: the intermediate region Q ∼ m in the n l + 1 scheme misses power corrections O(m) that are important numerically (indeed, possibly dominant) in that region. A number of approaches have been proposed in the past, starting with the ACOT proposal [21], which allows for the power corrections present in the FF result to be incorporated in the ZM-VFNS prediction. Many other such proposals exist [22][23][24][25][26][27][28][29][30] and approaches of this type are known as general mass VFNS (GM-VFNS). These proposals have been refined in a number of phenomenological applications [31][32][33][34][35][36][37] and in studies on the rationale behind the choice of either scheme [38][39][40][41].
Schematically, the GM-VFN schemes work in the following way (see ref. [28] for details): dσ (GM−VFNS) = dσ (n l +1) (Q, m) + dσ (n l ) (Q, m) − dσ (n l ,0) (Q, m) , (1.1) where dσ (n l +1) (Q, m) is the ZM-VFNS result, dσ (n l ) (Q, m) is the one in the n f = n l scheme with full mass dependence in the perturbative coefficient function and dσ (n l ,0) (Q, m) is its massless limit (which contains ln(m) and m-independent terms but no O(m) ones). Such a construction naturally converges to the ZM-VFNS one for Q m and has the added benefit that at scales Q ∼ m also O(m) power corrections that are part of the massive n f = n l calculation are included. Note that since in most GM-VFNS pdfs and α s are transformed into the same scheme, either n f = n l or n f = n l + 1, the notation in eq. (1.1) emphasizes the number of flavors in the partonic cross-sections and not the ones in α s and pdfs.
While the above construction is undoubtedly an improvement over the ZM-VFNS for scales Q ∼ m, GM-VFN schemes suffer from certain ambiguities. First, O(m) terms in reactions initiated by heavy flavors are usually not introduced. Presumably those are suppressed numerically by the smallness of the heavy-flavor pdf; we revisit these power corrections in sec. 2.2. Second, dσ (GM−VFNS) , as defined in eq. (1.1), does not behave as desired for Q ∼ m because the difference dσ (n l +1) (Q, m)−dσ (n l ,0) (Q, m) does not automatically vanish there. To restore the expected dσ (GM−VFNS) → dσ (n l ) (Q, m) behavior in this limit, one typically suppresses this difference by hand 2 . In the existing literature this suppression is implemented in two ways: either by multiplying that difference by an arbitrary function which vanishes below Q = m or by introducing a specially designed DIS-inspired rescaling (known as χ-rescaling) of the partonic variable x. See ref. [28] for a detailed discussion on this point.
In this work we demonstrate that a numerically accurate prediction can be achieved within the simplest possible variable flavor scheme: ZM-VFNS. To that end only one modification needs to be made compared to the conventional formulation: move the matching point µ b that separates the n f = n l and n f = n l + 1 schemes towards higher values, i.e. from the conventional value µ b = m to µ b = κ · m with the parameter κ as high as 10. The detailed justification for, and implications of, such a choice are given in the following section. The remainder of the paper is then devoted to the phenomenological study of our proposal.

The proposed idea
The single most relevant piece of information about our work is the position of the matching point µ b . In order to explain why and how we change it we need to first revisit its status. 2 While in the calculation of DIS structure functions all schemes implement a suppression, such a suppression is not necessarily introduced in calculations of Tevatron or LHC processes where the typical scales are much larger than the heavy quark mass m. See for example refs. [32,33,36,37,42].
All existing pdf sets and virtually all papers on the subject (with handful of exceptions we discuss below) set µ b = m. We have been unable to trace this choice to any specific proposal in the literature. It seems to us that this choice stems from a combination of two results. First, as discussed in sec. 1, decoupling implies that the matching point µ b has to be of the order of the mass m. This, of course, does not imply that µ b has to be set equal to m and, in fact, there is a considerable range of values around m that satisfy this requirement. It is precisely this range that we will be exploring in this work.
Second, the requirement for continuity of the pdfs across the HFMP was historically very influential. The matching conditions for α s and the full set of pdfs f i read where the sum over j goes over all flavors but the heavy one and ⊗ stands for the usual integral convolution in the variable x. It is obvious from eqs. (2.1) that at the leading order (LO) both the coupling and pdfs are continuous for any value of µ b . At higher orders the coefficients c (n) and K (n) ij are polynomials in log(µ 2 b /m 2 ). A peculiar feature of the NLO result is that the "constant" (i.e. log(µ b /m)independent) terms in both c (1) and K (1) ij are zero. This implies that the coupling and pdfs are continuous also at NLO if the matching scale is chosen to be the mass, i.e. if µ b = m. Naively, the requirement for continuity seems desirable and much attention to it has been devoted in the past [22,23]. However, this is not so as we demonstrate in the following.
First, it turns out that the continuity of α s and pdfs across the matching point µ b for µ b = m is accidental. It does not persist at higher perturbative orders in the space-like region [13,14] or for related quantities in the time-like case [43][44][45][46]. Therefore, one has to accept that the argument for continuity by a special choice of the matching point µ b is invalid.
Second, pdfs are not observables and their continuity is not a formal theoretical requirement. Perhaps somewhat surprisingly, in this work we demonstrate precisely the opposite: the presence of discontinuities across µ b in α s and pdfs is, in fact, beneficial for the continuity of predictions for observables at higher perturbative orders. Technically this happens because the discontinuities in α s and various pdfs conspire in such a way that they tend to largely compensate each other. This happens in a variety of processes and observables (both inclusive and differential). Unsurprisingly, this observation reminds us that the theory we work with is incredibly self-consistent! At this point we can introduce our proposal in the following way: Proposal: • For the observable of interest (which could be inclusive or fully differential) choose the value (or functional form) of the factorization scale; • Decide on the value of µ b . In this work we explore values for µ b as large as 10m; • Events with kinematics for which µ F ≤ µ b are computed in the n f = n l scheme and full mass dependence is to be retained; • Events with kinematics for which µ F > µ b are to be evaluated in a scheme with n f = n l + 1 active flavors where the mass m is set to zero.
In effect this is a ZM-VFNS scheme with a HFMP that is set to a (much) higher value than in all existing pdf sets. For short, we will sometimes call it "variable µ b approach". With the exception of appendix A, in the rest of this work we focus on the bottom pdf. For added clarity, we make our proposal explicit with the following b-pdf specific equation:

Addressing some obvious questions
In the following we address a number of questions about our construction.
1. What about the mass power corrections O(m) neglected above the HFMP? Our proposal is not exact and indeed for scales above the point µ b it misses power corrections O(m). However, unlike existing VFNS realizations, we have a parameter that controls the size of the error we make by neglecting these power corrections. For a given threshold the actual error that is made is O(m 2 /µ 2 b ) 3 and could be as small as 1% for µ b = 10m. This is a negligible effect at the current level of experimental and theoretical precision. Even for µ b = 5m the error is around 4%, i.e. very small. We would like to stress that none of the existing schemes have a way of parametrically estimating the error they make.
2. Can this proposal be complemented with existing, or future, GM-VFNS schemes like ACOT [21] or FONLL [28]? Yes. Our proposal does not preclude the inclusion of power corrections O(m) above the HFMP. However, by construction, such corrections will be small and therefore typically there is no need for them at the current accuracy level (we verify this explicitly in appendix A). One of the main goals behind our proposal is to have a scheme which is both accurate and very straightforward to implement.
3. How high should the value of µ b be? Clearly it should not be very high because if µ b m then collinear resummation will be spoiled. On the other hand, the power corrections missed in the massless calculation above the HFMP are suppressed by O(m 2 /µ 2 b ) and the motivation for choosing large values for µ b is that they can be made tiny. In fact, as we demonstrate at length in the rest of this paper, a value of µ b that is as high as µ b = 10m is allowed and, therefore, for such a choice any power corrections are about 1% effect and thus rendered phenomenologically irrelevant. 4. Is resummation accuracy maintained? Yes, if the HFMP is not chosen too high. We have checked that even for µ b = 10m the asymptotic behavior of cross-sections for Q m is maintained. 5. What is the distinction between choosing µ b and the factorization scale µ F ? There is a fundamental difference between these two scales. While µ b is totally process independent -in effect it only knows about the proton and its choice should reflect that -the factorization scale µ F is process and observable dependent. This scale implements the distinction between short-and long-distance physics specific to the way the proton is probed in a given measurement. For this reason the choice for µ F should be made prior to, and independently of, that for µ b .
6. What about intrinsic heavy flavor? In the context of our work a possible intrinsic heavy flavor contribution will only affect the boundary condition at the HFMP used subsequently by the DGLAP evolution. If needed it could be implemented on top of our considerations.
7. What is the relation to fragmentation functions? Although we have not explored the implications of our proposal for the case of fragmentation functions we see no reason why it would not apply directly there, too. One should note that our situation corresponds to the case of heavy flavor contribution to the fragmentation of light hadrons. The fragmentation of heavy flavored mesons is different and may have to be analyzed with more care.
2.2 Is the value of µ b a matter of choice or part of the theoretical uncertainty?
At this point one may wonder if the value of µ b is a matter of making a suitable choice or if it should be considered apriori undetermined and the ambiguity due to its variation around the point µ b = m considered as part of the theoretical systematics. The latter possibility has been explored in three recent papers [36,37,47] which, to the best of our knowledge, are the only ones where values for µ b different from the standard one µ b = m have been considered. For physical predictions, refs. [36,37] pick the canonical value µ b = m as the central choice and use variations about this point to estimate the associated theoretical uncertainties. Our motivation for moving away from the canonical choice µ b = m is different from the ones in the above-mentioned references. We explain it next.
As far as α s and pdfs are concerned any choice µ b ∼ m is equally correct and therefore the choice of µ b is unimportant 4 . This "translational invariance" of µ b is, however, not respected by the perturbative differential cross-sections (i.e. coefficient functions). While in the n f = n l scheme full mass dependence is retained, this is not the case for the n f = n l + 1 scheme where power corrections O(m) are missed. It is these missing power corrections that make the theory predictions depend on the position of the HFMP.
In this work we take the viewpoint that the correct position of the HFMP is a matter of choice: it is the position where these corrections are minimized. Such a position may not be unique; it may belong to a finite range whose size is determined by the uncertainty tolerance level. Values of κ in the interval κ = 5 − 10 seem to satisfy this requirement for b-production. Once κ is chosen within such a range, its variation within that range would be indicative of missing higher-order corrections related to perturbative matching. Indeed, in this work we see that the inclusion of higher-order corrections in the pdf matching conditions (and to a lesser degree in the pdf evolution) leads to a systematic reduction of the sensitivity of observables to the position of the HFMP µ b for sufficiently high values of µ b .
Finally, we would like to stress that there are both technical and conceptual reasons why O(m) corrections are always missing. In existing GM-VFNS implementations this is due to the introduction by hand of rescaling or damping functions. These are checked in known cases but their applicability for any process or to higher-orders is not obvious or automatic. O(m) corrections are also neglected in higher-order calculations with incoming massive quarks. The reason for this is twofold: on the technical side such calculations are uncommon and technically challenging [48]. On the conceptual side, even if one makes the effort to compute a cross-section with full mass dependence for initial-state massive quarks and then subtracts all quasi-collinear singularities as appropriate for a massless MS subtraction, it is known that collinear factorization is violated starting at NNLO with two initial state massive fermions [49][50][51]. More about this topic can be found in the textbook by Collins [52].

The effect of changing HFMP on α s and pdfs
In this work, we produce 6 sets with varying b−HFMP's, in which the value of µ b is set to κ · m with κ = 1, 2, 4, 6, 8, 10. For reasons of technical convenience we base our study on the NNPDF3.0 family of pdf sets, although it could be performed with any existing pdf set. The pdf set with κ = 1 coincides with the standard NNPDF3.0 pdf set with α s (m Z ) = 0.118. The strong coupling constant is fixed at the conventional reference scale, i.e. the Z-pole mass, with m Z = 91.2 GeV and for n f = 5. All pdfs are parametrized below the charm threshold (m c = 1.275 GeV) at the scale Q 0 = 1.0 GeV, with n f = 3.
Following the NNPDF setup, in our analysis the strong coupling constant is evolved from the initial scale m Z to the scale Q of the data included in the fit, starting with n f = 5 and crossing all HFMP between m Z and Q. In turn, the pdfs are evolved from the initial scale Q 0 crossing at least the charm HFMP and up to the scale at which data are available. 5 This is illustrated schematically in fig. 1.
Notice that in the NNPDF3.0 default set, the top mass is set to infinity, thus the maximum number of active flavors in both pdfs and α s evolution is n f = 5. The hadronic observables included in the NNPDF3.0 analysis are computed in the ZM-VFN scheme, i.e. no charm or bottom quark mass effects are included in the computation of Drell-Yan, jets or vector boson production cross sections. Instead, in the computation of DIS observables, the FONLL scheme is adopted, in which the massive calculation in the vicinity of the heavy quark HFMP is matched to the massless computation far above it.
In order to fully define the sets with µ b = κ · m, with κ > 1, one needs to specify their values at the initial scale Q 0 = 1 GeV. This could be done in two ways. One may wish to obtain the boundary condition for each set with κ > 1 by refitting the data as described above and evolving consistently with HFMP set at µ b = κ · m, as appropriate for that set. This way, in general, one will obtain initial conditions which are different for pdf sets with different κ. Alternatively, one may simply require that the initial condition for any value of κ be the same (and in particular be the same as the canonical case κ = 1). The rationale for the latter possibility is that one may imagine a lattice QCD-based prediction (see e.g. refs. [53,54] for latest developments) which presumably will be insensitive to the heavy quark content at high scales. 6 One may wonder if in practice choosing one or the other approach leads to significant differences in the pdf boundary conditions. We have checked this explicitly by producing pdf sets based on both approaches. We observe very small differences, with PDFs at the initial scale deviating from each other by less than a tenth of the pdf uncertainty bands for all values of x. As a result we have decided to use in this work only pdfs whose initial value is the same as the one of the standard NNPDF3.0 set. We next study the effect from changing the HFMP on the scale evolution of α s and all pdfs. We plot them as functions of the factorization scale Q at NLO and NNLO. The pdfs are shown for four values of the partonic momentum fraction x.
The evolution of α s is shown in fig. 2. In fig. 13 fig. 15 (fig. 16), respectively. In all cases we observe significant and progressively increasing discontinuities as the value of κ increases. As we demonstrate in the following sections, however, these discontinuities cancel each other in observables.
Although this may appear somewhat counterintuitive, we would like to point out that the discontinuity in the matching conditions has a beneficial effect on the b-quark pdf. As can be seen from fig. 14 going from LO towards NNLO the discontinuity in the matching condition increases and the bottom pdf at large scales becomes less sensitive to the position of the threshold. This happens because the large discontinuity at NNLO allows the b-quark pdf to DGLAP-evolve starting from a non-zero value immediately after the HFMP. As an example for how pdf continuity negatively affects pdfs one can take the LO case, not shown here, where the matching condition is continuous for any position of the HFMP. For this reason the b-quark pdf has to start from zero for any position of the HFMP, which means any two pdfs derived with different HFMPs will only very slow converge towards each other at large Q.
Finally, we also notice that in all cases, for very large Q all pdf sets with different HFMPs converge towards each other and this convergence is well within the pdf uncertainty. This is an important validation of our approach which needs to preserve the asymptotic behavior for large Q (i.e. to preserve the DGLAP evolution and the resummation of terms ∼ ln n (Q/m)).

Testing the approach with LHC processes
Our idea -to work with pdfs with high heavy flavor matching points -is motivated formally, not phenomenologically. Therefore, in this section we want to verify that it works in practice by quantifying to what extent the position of the HFMP affects a broad range of LHC processes and observables. We will be paying particular attention to the discontinuities in inclusive and differential cross-sections across the HFMP and the effect of changing HFMP on precision LHC observables.

Effect of varying HFMP on processes sensitive to b-quark pdf
We consider the following three processes that are sensitive to the b-content of the proton: the total t-channel single top cross-section at LHC 13 TeV with µ F = µ R = m t /2, the total cross-section for the process bbZ at LHC 13 TeV and the differential cross-section of Z + b-jet at LHC 13 TeV 7 as a function of p T (J b ) with The coefficient 1/4 in the equation above is chosen so that for small p T,j the scale µ F is small enough and crosses as many values of µ b as possible (see fig. 9). At leading order, the process bbZ contains bb → Z in the 5 flavor scheme (FS) while in the 4FS it is given by gg → Zbb. At higher orders it is defined in such a way that it involves all diagrams containing a bbZ vertex dressed with QCD radiation, as appropriate.
All calculations through NLO QCD are derived with the MadGraph5 aMC@NLO code [55]. All public and private pdf sets we use are incorporated through the LHAPDF library [56]. The 5FS NNLO corrections to the t-channel single top cross-section are not computed from first principles but are obtained by rescaling the NLO 5FS results with an NNLO Kfactor derived from the results of refs. [57,58]. Specifically, the proxy to the NNLO t-channel single-top cross-section in the 5FS is where the ratio σ NNLO /σ NLO is constructed from the 13 TeV inclusive cross-section numbers in table 1 of ref. [58] (i.e. for κ = 1) and σ NLO (NNLO pdf, κ) is calculated exactly, using the relevant pdf set with κ ≥ 1. Given the extremely flat scale dependence and small NNLO correction such an approach is more than adequate for our purpose. For this reason we have not extended the NNLO curves for scales outside the µ F range considered in refs. [57,58], which is sufficient for our goals. The 5FS cross-section for bbZ has been calculated at NNLO by using a private code [59], which has been cross-checked at LO and NLO against MadGraph5 aMC@NLO. The 4FS cross-section has been computed with MadGraph5 aMC@NLO. In both calculations the coupling of the Z boson to light quarks is set to zero, so that only bottom-initiated diagrams are considered.
We first consider the cross-sections with absolute normalization. They are plotted in figs.      this result should be expected, it is a nice check that the approach works as anticipated. Second and more important, all 5FS curves that correspond to different values of µ b become more and more tightly packed together as we go to higher perturbative orders. This means that at higher orders the predicted cross-sections become less sensitive to the position of µ b . These two observations apply universally for all three different processes.
In order to clarify the origin of the second feature mentioned above we perform further checks. In figs. 6, 7, 8 we plot the 5FS curves as the ratio: , κ = (1, 2, 4, 6, 8, 10) . Strikingly, despite the very different nature of the three processes, there is an extreme similarity between all of them. Clearly, we observe a process-independent feature: by increasing the order of the perturbative cross-section the slope of the curves becomes smaller (i.e. they become less scale-dependent) while by increasing the order of the pdfs the curves become closer to each other. The improvement from the inclusion of higher order pdfs is very significant and, at NNLO, the dependence on the position of the HFMP is dramatically reduced irrespectively of the order of the perturbative cross-section.

Discontinuities in cross-sections across HFMPs
The size of the discontinuity of hadronic cross-sections during the transition between 4FS and 5FS at the scale µ b is an important criterion for our work. We first consider the three processes already studied in sec. 4.1. In fig. 9 we show the discontinuities in the predicted     cross-sections due to the transition from 4FS to 5FS, as a function of the matching point µ b . The discontinuity across threshold, as a function of µ F , is defined as where M stands for the relevant mass parameter for each process (m Z or m t ) and is very small. In the discontinuity eq. (4.4) above (as well as in eq. (4.5) below), the 5FS cross-section is evaluated just to the right of the matching point while the 4FS one is computed just to the left of it. Ideally, the discontinuity when going across the HFMP should be zero. There are two important features in fig. 9. First, by going from LO to NLO the dis- LO 5F/4F disc. NLO 5F/4F disc. Figure 9. Discontinuity between 5FS cross-section above threshold and 4FS cross-section below threshold for LO (dotted) and NLO (dashed) for t-channel single top cross-section (left), bbZ production cross-section (center) and the p T (J b ) differential distribution in Z + b (right). All are for LHC at 13 TeV as a function of κ. The discontinuity is defined in eq. (4.4) and, in the ideal case, it should be zero.
continuity decreases drastically. Second, we note that the 4FS-to-5FS discontinuity tends to decrease when the HFMP is increased. The details are process dependent but the trend is significant and generic. We, therefore, conclude that both the inclusion of higher order corrections in observables and the increase in the value of the HFMP lead to a decrease in the discontinuities in observables. We stress that in the case of the p T (J b ) distribution only three b-quark matching points are crossed as p T (J b ) sweeps the interval between zero and infinity. This is due to the specific functional form (4.1) adopted for the factorization scale. The processes just considered have two limitations as far as our study is concerned. First, the 4FS predictions for these processes are not known at NNLO. Second, as the matching point is varied in the range κ = 1 to κ = 10 the factorization scale µ F , which is kept equal to µ b , can significantly deviate from the natural scale 8 in the respective process. This is not ideal since when the factorization and/or renormalization scales are taken to be very different from their natural value the convergence of the perturbation series gets affected.
In order to demonstrate that the above two limitations do not alter our conclusions we extend our study as follows. We construct two families of (unphysical) processes for which both the 4FS and 5FS predictions can be derived at NNLO. These families of processes are constructed in such a way that for any value of the matching point µ b the factorization scale is always close to the natural scale for that process. Specifically, we study the cross-sections for Z-like and tt-like production. By Z-like we mean Z production but with appropriately chosen Z mass, such that m Z = µ b as we change µ b from m to 10m. Same for tt-like production but with m t = 2µ b . The factorization and renormalization scales are taken to be µ F = µ R = µ b thus avoiding large ratios between the natural scales and µ F,R .
The 4FS predictions for these processes are not known at NNLO but we devise an approximation which is sufficiently accurate for our purpose: we take the partonic-cross section to be the 5F NNLO one for massless b-quark and we evaluate it with a 4F α s and convolute The discontinuity is computed in both the canonical approach ("Standard", in green) and the variable HFMP approach advocated here ("Variable", in blue).  Table 1. Discontinuity eq. (4.5) for Z-like production at LHC 13 TeV at LO, NLO and NNLO in the variable HFMP approach ("Variable") advocated here and in the canonical ("Standard") approach.
it with 4F pdfs. We expect this to be a good approximation because the extra diagrams present in the 5F NNLO partonic cross-sections are multiplied by vanishing pdfs and thus do not contribute. Terms proportional to the b-quark mass that originate from 4F diagrams with b-quark loops or b-quark pair emission are also missed in the massless-b 5F diagrams but these terms should contribute little, especially at large scales. For these reasons we expect the error we make with this approximation (i.e. the difference to the full 4F NNLO results) to be at the sub-percent level. With the help of explicit calculations we have checked that at NLO the inclusive tt and Z cross-sections are affected at the permille level. The definition of the discontinuity is adapted from eq. (4.4):  Table 3. Dependence of the total Z cross section at LHC 13 TeV on the threshold scale µ b = κ · m (recall that κ = 1 represents the standard choice in all publicly available pdf sets). Shown is the 5FS cross-section predicted at LO, NLO and NNLO for m Z = µ R = µ F = 91.1876 GeV.
with M = m Z or m t , as appropriate for the process.
The results for Z-like and tt-like production are given in table 1 and table 2, respectively. Both are plotted in fig. 10. The discontinuities are given at LO, NLO and NNLO in the variable HFMP approach advocated here ("Variable") as well as in the canonical approach ("Standard"). The numbers in the canonical approach are derived in the following way: we take standard 4FS and 5FS pdf sets (with µ b = m, as usual) and evaluate the ratio (4.5) at the points κ · m. In other words we extend the 4FS pdf set from scales ∼ m all the way up to κ · m which reflects the current practice within the canonical approach.
The way to read fig. 10 is to recall that in a given process the value of M is fixed (with M being m Z or m t ). Thus, one should choose between computing the cross-section for that process in either the standard or variable approaches and at LO, NLO or NNLO.
It is evident from fig. 10 that the variable µ b approach advocated here has smaller discontinuities than the standard approach for the whole range of masses M considered here and for any perturbative order (except the Z-like cross-section for small values of M at NLO where numerical effects play a role). The discontinuity at NLO is smaller than the NNLO one for tt-like production with small masses. However, as the mass m t increases the NNLO discontinuity becomes competitive with the NLO one. For Z-like production, the NNLO discontinuity in the variable approach is smallest for the whole range of m Z values.  1.649 Table 6. As in table 3 but for the t-channel single-top total cross-section and the ratio R t/t of single top versus single antitop cross-sections with m t = 173.3 GeV and µ F,R = m t /2.

Effect of changing the value of µ b on standard LHC candles
Another important test for our approach is how increasing the value of the HFMP affects standard precision LHC candles. We note that these processes have typical scales that are much larger than any of the b-quark matching points considered here. In  while in table 6 we show the t-channel single top cross-section and the ratio of single top versus single antitop cross-sections, R t/t . All results are for LHC at 13 TeV. We observe that for all processes the changes between the canonical approach κ = 1  and any of the variable ones with κ as large as 10 at NNLO are small compared to the theoretical and experimental errors. We note that the increased order of the prediction significantly decreases this difference. The most sizable effect is observed in single top, where for µ b = 10m we have a NNLO prediction that we estimate to be almost 5% lower than the current prediction. This change is significantly larger than the NNLO scale error (which is around 1%) but the two are compatible within the experimental error (which currently is below 10%). The long-term prospective for measuring the single top cross-section at the LHC points towards 5% precision for this observable. This means that it may not be easy to use single top LHC data to discriminate between these two approaches. An alternative may be the measurement of the ratio σ(tt)/σ(Z) which, as can be seen in fig. 11, is currently known with accuracy of about 2.8% [62]. Future improvements in this ratio may become the leading candidate for distinguishing the predictions based on these two approaches.

Conclusions
In this work we advocate a new approach to constructing heavy-flavor pdfs, namely, a standard ZM-VFNS but with heavy flavor matching point µ b which is taken to be significantly higher than the conventional value µ b = m. We extensively test our proposal on a range of NLO and NNLO precision observables at the LHC. We find that our approach is competitive with the current state-of-the-art GM-VFNS approaches. Its main advantage over existing GM-VFN schemes is its transparency and simplicity; it is straightforward to formulate and implement in practice for any process at hadron colliders and it avoids the need for adding by hand rescaling or damping factors. Our proposed approach typically leads to smaller discontinuities in observables across the heavy flavor matching point compared to conventional approaches. It is straightforward to implement in practice thanks to the existence of public tools like xFitter [63], APFEL [64] and LHAPDF [56].
We demonstrate that our proposal satisfies all requirements for constructing a good pdf set: first, it maintains collinear resummation provided the HFMP µ b is not chosen to be too large; we consider the range µ b = 5m − 10m to be optimal. Second, power corrections O(m) of the heavy quark mass, which constitute the main problem in constructing heavyflavor pdfs, are under control. Despite the fact that in the n f = n l + 1 scheme we only use purely massless coefficient functions (which greatly simplifies scheme's implementation) the neglected power corrections can be made negligible, by construction. We have verified this explicitly by a direct comparison with FONLL predictions for the case of DIS production of bottom (at NLO and NNLO) and charm (at NNLO); see appendix A for details.
We pay close attention to the issue of continuity of observables when one crosses the HFMP. As the HFMP is moved away from the point µ b = m and both α s and pdfs become more discontinuous across that point, one may naively expect that observables may also become more discontinuous (which would be bad). Our analysis shows precisely the opposite: typically, using pdfs of higher orders or with higher HFMPs, leads to smaller discontinuities in observables. We suspect that pdfs with full N 3 LO accuracy will further improve continuity, mainly as a result of higher-order corrections to the heavy-flavor matching conditions. The results of ref. [61] suggest that large shifts in predictions compared to N 2 LO should not be expected. The significant cancellation between various contributions we observe is in line with the expectation that observables should be continuous to all orders. This demonstration of the self-consistency of the theory is, in our opinion, not as apparent in other schemes.
For the sake of simplicity in our consideration we have mostly focused on the case of a single heavy flavor which we have taken to be the bottom. However, this need not be the case; all heavy flavors (charm, bottom and top) can be treated this way and their corresponding matching points may be set independently. Our ideas are particularly well suited for applications in top production. No calculations for top production exist in the 6F scheme, 9 despite the fact that top quarks with p T in the TeV range are already routinely measured at the LHC. This fixed-flavor-like treatment clearly contradicts the spirit of all existing GM-VFNS. On the contrary, our proposal indicates that, especially when working with the factorization scale proposed in ref. [65], top production should be described within the 5FS for p T 's as large as 3 TeV if κ t = 10 is chosen.
Similarly, the recently introduced pdfs with full SM content at very-high energies [73] represent another natural application for the ideas introduced in the present work.
Our proposal may be useful also in the context of fragmentation functions, especially when treating heavy quark contributions to the fragmentation of light hadrons (this problem was recently considered in ref. [66] within a FONLL-inspired GM-VFNS framework). While typically the uncertainties in fragmentation function analyses are larger than in pdf ones, recent advances in this subject have introduced NNLO QCD precision into global fragmentation function fits [67,68] which suggest that a more refined treatment of HFMPs in fragmentation functions may also be beneficial.
Finally, we would like to emphasize that in order to fully explore the phenomenological implications of our proposal, a more complete analysis along the lines of the global fits performed by the pdf fitting collaborations will be required. This will allow one to precisely determine the impact of our proposal on precision observables like Z and single top production where we see differences with respect to standard approaches that are significant yet not sufficient at present to differentiate between the two approaches. We hope that future work will ultimately help clarify the proper interpretation of such differences, i.e. if they should be considered as due to difference in approach or as a genuine uncertainty within existing approaches not fully exhibited until now.

A Example: charm and bottom fits from DIS data
In this appendix we explicitly demonstrate in the context of inclusive DIS that the power corrections O(m) indeed become negligible as the HFMP is increased. We perform fits to inclusive DIS data for bottom and charm production at NLO and NNLO. The size of the power corrections can be established by comparing pdf fits based on our approach and on the FONLL approach which is designed to contain the bulk of those power corrections.
Specifically, we perform a series of pdf fits by implementing the method proposed in this paper and using the open-source fitting code xFitter (former HERAFitter) [63] interfaced to the APFEL code [64]. The dataset included in these fits comprises the combined HERA1+2 H1 and ZEUS inclusive DIS cross-section data [69], the combined H1 and ZEUS charm production cross-section measurements [70] and the separate bottom production cross sections from H1 [71] and ZEUS [72]. This dataset provides sufficient information for a reliable determination of pdfs. Furthermore, thanks to the inclusion of heavy-quark tagged data, it is also sensitive to heavy-quark mass effects which allows us to assess the impact of different prescriptions for describing heavy flavors. In order to validate our method, we perform fits both at NLO and NNLO moving separately the charm and bottom HFMPs µ c,b by a factor κ = 1, 2, 5, 10 with respect to the masses m c,b . We then compare to the corresponding fits performed in the FONLL scheme. In the case of charm, given that m b /m c ∼ 3, in all fits we set µ b = 5m b in such a way that the bottom HFMP is always sufficiently above the charm one even for κ c = 10. The aim of these fits is to show that, as κ increases, the fits using our method and the FONLL scheme become closer. In order to assess the difference between our method and the FONLL scheme, in fig. 12 we display the quantity: as a function of the threshold rescaling parameter κ when the charm threshold is moved at NLO (red curve) and NNLO (orange curve) and when the bottom threshold is moved at NLO (blue curve) and NNLO (green curve).
First we observe that, as expected, the FONLL scheme at κ = 1 provides a much better description than the conventional ZM-VFNS. As the value of κ increases, however, the difference between our prescription and the FONLL scheme starts to decrease. Notably, in the case of bottom, a low value of κ = 2 is sufficient to bring δχ 2 close to zero both at NLO and NNLO. This is partly due to the fact that the value α s (µ b ) is sufficiently small.
The picture for charm is more complicated. At NLO the FONLL fit is systematically better than our prescription even for large values of κ. The large value of α s (µ c ) plays a role here since likely it enhances the discontinuity around the HFMP. As can be seen from Heavy quark Pert. order Scheme  Table 7. Values of the total χ 2 normalized to the number of degrees of freedom for the fits with different values of κ for the charm and bottom HFMPs, at NLO and NNLO, derived using the FONLL scheme and the proposal in this work. table 7, however, one should keep in mind that the quality of the FONLL fit is significantly degraded for NLO charm compared to the other cases we consider. This most likely indicates that the perturbative description of charm is simply problematic at low orders. The situation is greatly improved at NNLO where the convergence between the fits derived in our method and in the FONLL scheme is reached around κ 5, fully in line with our expectations.
In table 7 we give the values of the total χ 2 normalized to the number of degrees of freedom for each of the fits discussed above. In order to assess the significance of deviations, it is useful to keep in mind that all fits discussed here have been performed with the same dataset and using the same pdf parametrization. Therefore, the number of degrees of freedom is common to all fits and equals 1207. This means that variations of the order of a few permill at the level of the normalized χ 2 's are statistically significant. 10 For the case of charm at NLO it is clear that going from κ = 1 to κ = 10 the χ 2 deteriorates significantly, by about 7%, for FONLL and about 4% for our fits. As a results of this we believe that the χ 2 values in the first two lines of table 7 do not provide a basis for assessing the goodness of our prescription; see also the related discussion in ref. [47].
The variation of the χ 2 between κ = 1 and κ = 10 for charm at NNLO is much smaller, about 1.5% for FONLL. This reduced sensitivity to the charm threshold position allows one to reliably estimate the quality of the method proposed in this work. Indeed, as mentioned above, the fits based on our prescription approach the FONLL ones for κ = 5 which is fully in line with our expectations of κ being in the interval κ = 5 − 10.
As far as the bottom is concerned, the picture at NLO and NNLO is essentially the same as for the charm at NNLO. In particular, the fits in the FONLL scheme are very stable upon displacement of the HFMP. Correspondingly, the quality of the fit within our method quickly approaches that of the FONLL scheme as κ increases and already at κ = 2 the two are essentially equivalently good.    Figure 14. As in fig. 13 but for the bottom pdf.  Figure 15. As in fig. 13 but for Σ = u +ū + d +d + s +s + c +c + b +b.  Figure 16. As in fig. 13 but for Σ = u +ū + d +d + s +s + c +c.