Heavy-quark mass dependence in global PDF analyses and 3- and 4-flavour parton distributions

We study the sensitivity of our recent MSTW 2008 NLO and NNLO PDF analyses to the values of the charm- and bottom-quark masses, and we provide additional public PDF sets for a wide range of these heavy-quark masses. We quantify the impact of varying m_c and m_b on the cross sections for W, Z and Higgs production at the Tevatron and the LHC. We generate 3- and 4-flavour versions of the (5-flavour) MSTW 2008 PDFs by evolving the input PDFs and alpha_S determined from fits in the 5-flavour scheme, including the eigenvector PDF sets necessary for calculation of PDF uncertainties. As an example of their use, we study the difference in the Z total cross sections at the Tevatron and LHC in the 4- and 5-flavour schemes. Significant differences are found, illustrating the need to resum large logarithms in Q^2/m_b^2 by using the 5-flavour scheme. The 4-flavour scheme is still necessary, however, if cuts are imposed on associated (massive) b-quarks, as is the case for the experimental measurement of Z b bbar production and similar processes.


Introduction
Parton distribution functions (PDFs) determined from global analysis [1][2][3][4][5] of fixed-target and collider data, mainly from deep-inelastic scattering (DIS), are an essential input to theory calculations for hard processes at hadron colliders such as the Tevatron and LHC. In addition to the fitted input PDFs, several other parameters enter into the global fits, which affect both the PDFs obtained and predictions for hadronic processes. One important example is the value of the strong coupling α S and its uncertainty [6][7][8]. Others, which are the focus of the present paper, are the values of the heavy-quark masses and the scheme used to calculate heavy-quark contributions to observables. In particular, while the precise values taken for m c and m b may not be crucial for the processes included in the global PDF analyses, calculations of processes with explicit c-and b-quarks might be more sensitive to these values, and it is therefore desirable to have PDFs available which have been consistently fitted and evolved with corresponding values of m c and m b .
We will study two topics, both of which concern the treatment of the heavy charm and bottom quarks in global parton analyses. The first topic, which is the subject of Section 3, concerns the sensitivity of the MSTW 2008 global parton analysis [1] to the values of the heavyquark masses m h , with h = c, b. In Ref. [1] these masses were taken to be m c = 1.40 GeV and m b = 4.75 GeV. Here we perform global fits for a range of m c and m b about these values, with α S (M 2 Z ) allowed to be a free parameter. In this way, we determine the values of m c and m b preferred by the data, and the correlations between these values of m h and α S (M 2 Z ). Due to a significant correlation between m c and α S (M 2 Z ), we also perform fits with varying m c but with α S (M 2 Z ) held fixed. We study how the cross sections for W , Z and Higgs boson production at the Tevatron and the LHC depend on the choice of m c and m b , and we provide a recommendation on how to include the uncertainty arising from m h in a general cross-section calculation.
The second topic, described in Section 4, is the generation of sets of 3-and 4-flavour PDFs corresponding to the MSTW 2008 [1] 5-flavour sets of PDFs. We follow the same procedure previously used to generate the 3-and 4-flavour versions [9] of the 5-flavour MRST 2004 PDFs [10], i.e. the input PDFs (and α S ) at Q 2 0 = 1 GeV 2 determined from fits in the 5-flavour scheme are used as the initial condition for 3-or 4-flavour evolution. However, going beyond Ref. [9], we also provide 3-and 4-flavour eigenvector PDF sets to enable calculation of PDF uncertainties, and we provide 3-and 4-flavour PDF sets for a wide range of m c and m b values, respectively. As an example of the use of the central 4-flavour PDFs for the default quark masses, we compare the total cross sections for Z production at the Tevatron and LHC in the 4-and 5-flavour schemes. However, first we begin with a short résumé of the alternative "schemes" to treat heavy flavours 1 , and in particular, an explanation of what precisely we mean by 3-, 4-, and 5-flavour PDFs.

Schemes for the treatment of heavy flavours
It is appropriate to briefly recall the various schemes for the treatment of heavy flavours in global parton analyses. In PDF analyses it is common to start evolving upwards from Q 2 = Q 2 0 ∼ 1 GeV 2 with the distributions of the three light quarks (u, d, s), assuming that they are massless. As we evolve upwards, we have the choice to turn on the heavy-quark (c, b, t) distributions as we pass through their respective transition point, for which a convenient choice is µ 2 F = m 2 h . As we pass through a transition point, the number of active quarks is increased from n to n + 1, and the parton densities are related to each other perturbatively, i.e.
where the perturbative matrix elements A jk (µ 2 F /m 2 h ) contain ln(µ 2 F /m 2 h ) terms which are known to O(α 2 S ) [11] and O(α 3 S ) [12]. The "x" arguments have been suppressed in Eq. (1), and the symbol ⊗ is shorthand for the convolution Eq. (1) relates the f n i (µ 2 F ) set of partons to the f n+1 i (µ 2 F ) set, guaranteeing the correct evolution for both the n and n + 1 regimes. We make the simplest choice, µ 2 F = Q 2 , for the factorisation scale.
Hence, we have to decide whether or not to keep a heavy quark as just a final-state particle, and not as a parton within the proton. We may choose to keep just the three light flavours as parton distributions. We will call this the 3-flavour scheme (3FS), though it is often referred to as the fixed flavour number scheme (FFNS). Alternatively, we may include the c-quark in the evolution above Q 2 = m 2 c and generate 4-flavour PDFs in a 4-flavour scheme (4FS). Actually, in the global MRST/MSTW parton analyses we also include the b-quark distribution in the evolution above Q 2 = m 2 b , but not the t-quark above Q 2 = m 2 t , so we generate 5-flavour sets of parton distributions in a 5-flavour scheme (5FS). So to be precise, in our definition of n f -flavour parton sets, n f refers to the maximum number of quark flavours in the evolution.
In each n f -flavour scheme (n f FS) the structure functions are given by the usual convolution of coefficient functions and parton distributions: where the sum j is over the gluon and the (variable, depending on Q 2 ) number of active quark flavours, n ≤ n f . We have a choice in how to choose n f and define the coefficient functions. One simple choice is to fix n = n f = 3. For the heavy flavours, all the m h dependence then occurs in the coefficient functions, and these are called the FFNS coefficient functions. The structure functions may be written in the form F (x, Q 2 ) = j=u,d,s,g However, the sum of the α k S ln m (Q 2 /m 2 h ) terms, with m ≤ k, is not included in the perturbative expansion. Thus the accuracy of the expansion becomes increasingly unreliable as Q 2 increases above m 2 h . In addition, there is the problem that the full mass dependence of the coefficient functions is known up to NLO [13], but is not completely defined at NNLO, i.e. the α 3 S coefficient, C FF,3,(3) 2,hg , for F 2 is not fully known, see Ref. [12]. (Here, the outer subscript "hg" denotes the g → h contribution to the heavy-flavour structure function F h 2 .) As an aside, we note that it would be possible to treat the charm quark as light, and the bottom quark as heavy. Then it would be possible to express the structure functions and cross sections in terms of four light quarks with all the mass dependence of the bottom quark contained in the coefficient functions. This is sometimes called the 4-flavour FFNS. However if one needs to use scales Q 2 < m 2 b and in the vicinity of m 2 c , as for example in the global PDF analysis, then the charm-quark mass dependence, and the charm transition point (if Q 2 ≤ m 2 c ), should be included. Thus, this 4-flavour FFNS is only applicable in a restricted range of Q 2 , and it will not be considered further here.
The alternative, and better, approach is to use the 4-flavour (and 5-flavour) PDFs of the variable flavour number schemes (4FS, 5FS). The simplest variable flavour number evolution procedure is to evolve treating some (or all) of the heavy quarks as massless, but to turn on the distributions at the appropriate transition points, Q 2 = m 2 h . That is, to assume that these heavy-quark distributions evolve according to the splitting functions for massless quarks. Thus, the resummation of the large logarithms in Q 2 /m 2 h is achieved by the introduction of heavyflavour parton distributions and the solution of the evolution equations. It is motivated by the observation that at high scales the massive quarks behave like massless partons, and the coefficient functions are simply those in the massless limit, e.g. for structure functions where, as in Eq. (3), the sum j is over the gluon and the (variable) number n ≤ n f of quarks that are active during the evolution. This is the so-called zero-mass variable flavour number scheme (ZM-VFNS). However, the ZM-VFNS has the failing that it simply ignores O(m 2 h /Q 2 ) corrections to the coefficient functions, and hence it is inaccurate in the region where Q 2 is not so much greater than m 2 h . Thus, strictly speaking, the ZM-VFNS does not define a scheme, but rather an approximation.
So we have two approaches where the treatment of heavy flavours is relatively simple. The 3-flavour scheme (or FFNS), appropriate to the region Q 2 m 2 c , and the ZM-VFNS, appropriate to the region Q 2 ≫ m 2 h . Clearly, for precision parton analyses, we must use a so-called general-mass variable flavour number scheme (GM-VFNS) which smoothly connects these two well-defined regions, so as to reduce to the (3-flavour) FFNS in the low Q 2 limit and to the ZM-VFNS in the high Q 2 limit (up to possible terms of higher order in α S ). In particular, in Eq. (3) the value of n increases by one each time we reach Q 2 = m 2 h in the evolution, and the coefficient functions C nGMVF for the maximum value of n = n f . There is some freedom in how one does this, and in the MSTW 2008 analysis [1] we use the definition of the GM-VFNS as described in Ref. [14]. 2 We also note that actually a 5-flavour GM-VFNS is used, since we do not include the top quark in the evolution. Note that this means that processes such as associated production of Higgs bosons with top quarks should be calculated using the 5-flavour PDFs with the full 2 → 3 subprocess matrix elements, i.e. gg, qq → Htt. In a 6FS, one would calculate this via tt → H, introducing a top-quark PDF. Although in principle this method would include resummed higher-order corrections ∼ [α S ln(M 2 H /m 2 t )] n , we would expect that for all practical purposes 3 at the LHC, the 5FS approach, gg, qq → Htt, supplemented by NLO corrections, would be adequate, see also Section 4.2.
As just noted, the use of a GM-VFNS is nowadays essential for the determination of a set of precision PDFs from the data. 4 However, there are processes where the full mass dependence of the coefficient functions for c and b production processes are only known in the case where the heavy quark is treated as a final-state particle, and not as a parton in the proton. For example, FFNS parton distributions are needed for use with the hvqdis [19,20] and mc@nlo [21] programs. Thus in this paper we make available 3-and 4-flavour PDFs. We also give an example of their use. This is the subject of Section 4.

Dependence of PDFs on the heavy-quark masses
In the MSTW 2008 global parton analysis [1] we used m c = 1.40 GeV and m b = 4.75 GeV, where these are the pole mass values. Some limited justification for choosing these values was briefly given at the end of Section 3 of Ref. [1]. To summarise, there was little data constraint on m b , so it was simply fixed at a value of 4.75 GeV, close to the calculated MS mass transformed to the pole mass. The fixed value of m c = 1.40 GeV was close to the best-fit value at NLO if treated as a free parameter, but a little higher than the best-fit value at NNLO, and lower than the calculated MS mass transformed to the pole mass. We discuss limitations in constraints on the pole masses by transforming from the calculated MS masses below in Section 3.3.
Here, using exactly the same data sets, we repeat the global analysis for different fixed values of m c and m b , but with α S (M 2 Z ) left as a free parameter. The MSTW 2008 global fit [1] used a wide variety of data from both fixed-target experiments and the HERA ep and Tevatron pp colliders. Neutral-current structure functions (F 2 and F L ) were included from fixedtarget lepton-nucleon scattering experiments (BCDMS [22,23], NMC [24,25], E665 [26] and 2 At NNLO this involves some modelling of the O(α 3 S ) FFNS coefficient functions which are not fully calculated. In the GM-VFNS we only need these in the low-Q 2 regime and we approximate them using the known results of the small-x [15] and threshold limits [16] (see Refs. [17,18] for recent refinements). 3 Note that for M H 1 TeV, α S ln(M 2 H /m 2 t ) 0. 35. 4 In fact we only use the GM-VFNS for structure functions in DIS, whereas hadron collider cross sections are calculated in the ZM-VFNS under the assumption that mass effects are negligible at the typically large scales. (Low-mass fixed-target Drell-Yan production is also calculated in the ZM-VFNS, but here heavy flavours contribute ≪ 1% of the total, so inaccuracies induced by using the ZM-VFNS are in practice irrelevant.) (a) NLO:  SLAC [27][28][29]), low-mass Drell-Yan cross sections from the E866/NuSea experiment [30,31], and charged-current structure functions (F 2 and xF 3 ) and dimuon cross sections from neutrinonucleon scattering experiments (CCFR/NuTeV [32,33] and CHORUS [34]). From the HERA experiments, H1 and ZEUS, data were included on neutral-and charged-current reduced cross sections (σ NC r and σ CC r ) [35][36][37][38][39][40][41][42][43], the charm structure function (F c 2 ) [44][45][46][47][48][49][50], and inclusive jet production in DIS [51][52][53]. From the Tevatron experiments, CDF and DØ, Run II data were included on inclusive jet production [54,55], the lepton charge asymmetry from W decays [56,57] and the Z rapidity distribution [58,59]. A more detailed description of the treatment of each of these data sets can be found in Ref. [1]. Note that more precise H1 and ZEUS data on F c 2 (and the beauty structure function, F b 2 , not included originally) are now available, but we stick to fitting exactly the same data sets as in the MSTW 2008 analysis.

Dependence on charm-quark mass m c
The sensitivity to the charm mass of the data in the global PDF analysis, and the subset of F c 2 data, in terms of the goodness-of-fit measure χ 2 , is shown in Table 1  Here, the 1σ uncertainty is obtained in exactly the same way as described for the uncertainty on α S (M 2 Z ) in Ref. [6], i.e. by examining the χ 2 profile of each data set included in the global fit. The distinguishing power of the various data sets is shown in Fig. 2. Here, the points (•) indicate the values of m c for which the χ 2 for each data set is minimised (within the context of the global fit), while the inner error bars extend across the 68% confidence-level (C.L.) region and the outer error bars extend across the 90% C.L. region (see Section 4 of Ref. [6] for the precise definitions of the 68% and 90% C.L. limits). As may be expected, the F c 2 (x, Q 2 ) data of the H1 and ZEUS collaborations are the most discriminating. Indeed, they dominate the determination of m c , as can be seen from Fig. 3, which shows χ 2 versus m c for the F c 2 data alone. For the HERA inclusive data, changes in the value of m c are partially compensated by changes in the gluon and in α S . The NMC data prefer a lower value of m c giving a quicker evolution near threshold. Similarly the BCDMS data mainly prefer a lower m c , but only due to the correlation with a lower value of α S . Supplementary plots of the χ 2 profiles versus m c for all data sets in the global fit are available from Ref.
[60]. The experimental uncertainty on m c , given in Fig. 1 and indicated by the horizontal dashed lines in Fig. 2, is chosen to ensure that all data sets are described within their 68% or 90% C.L. limits.
The dependence of the PDFs on the value taken for m c is shown in Fig. 4. The plots show the ratio of the NLO PDFs obtained with m c = 1.3 GeV and m c = 1.5 GeV to those of the MSTW 2008 analysis (which assumed m c = 1.4 GeV) at scales of Q 2 = 4 GeV 2 and Q 2 = 10 4 GeV 2 . We see that the ratios lie well within the 90% C.L. PDF uncertainty. The exception is the charm distribution at Q 2 = 4 GeV 2 where, as expected, the different charm transition points have an appreciable effect.

Dependence on bottom-quark mass m b
We repeated the whole exercise for m b . That is, we performed a series of global fits for different values of m b . The results are shown in Table 2. Now there is less sensitivity to the value of the quark mass, and essentially no correlation between m b and α S (M 2 Z ). Indeed, the global χ 2 stays fairly flat all the way down to m b = 3 GeV. For the lower values of m b , there is a slightly better description of the HERA data, including F c 2 (x, Q 2 ). A similar conclusion holds at NNLO, but with about half the change in global χ 2 , that is, even less sensitivity to m b .
Note that the HERA F b 2 (x, Q 2 ) data [46,47,[61][62][63] are not included in the global fits, neither in MSTW 2008 [1], nor in the fits described above. In Fig. 5 we compare the predictions of the NNLO fits with varying m b values to the F b 2 data from H1 [61]. We also show a few of the (less precise) ZEUS data points [62]. In Table 2 we give the χ 2 for the H1 F b 2 data accounting 5 for all 24 sources of correlated systematic uncertainty, and we also give the simple addition with 5 We use Eqs. (38)(39)(40) of Ref. [1], noting that there is a typo in Eq. (40) of that paper and σ uncorr.      The horizontal dashed lines in the plots indicate the 68% and 90% C.L. limits, determined according to a "hypothesis-testing" criterion, see Section 4 of Ref. [6].      Table 2. We conclude that the global data do not meaningfully constrain m b , and that, in view of Table 2

Pole masses of heavy quarks
Note that in our analyses m c refers to the pole mass of the charm quark. Both the transition matrix elements defining the boundary conditions for heavy-quark evolution [11] and the heavyquark coefficient functions used to define the GM-VFNS [11] are defined using the on-shell renormalisation scheme for the heavy quark. Unlike α S (µ 2 ), the pole mass m c is, in principle, a physically-defined (spoilt by confinement which gives power corrections, see below) quantity which is stable to the order of perturbation theory. So the values that we obtain in the NLO and NNLO analyses should, in principle, be the same. However, we see from Fig. 1 that the values are only marginally compatible. This is due to a variety of reasons. The perturbative expansion of heavy-flavour coefficient functions is not as convergent as one might hope, with large corrections near threshold and in the small-x limit. So even though the pole mass might, in principle, be the same quantity at two different perturbative orders, the calculation at lower 13 orders is missing significant corrections at higher orders which will affect the value of the pole mass obtained in an extraction from data. Perhaps most important in this respect is the large negative boundary condition (at Q 2 = m 2 h ) for the heavy-flavour distribution at NNLO, which has the effect of flattening the slope ∂F h 2 /∂ ln Q 2 (see Ref. [14] for details) and implies a lower mass to correct for this. There is also more potential variation in the GM-VFNS definition at NLO than NNLO, and there would be slightly more stability in preferred pole mass in, for example, the "optimal" scheme introduced in Ref. [64] (although this does not remove the variation and is the subject for a future study). Finally, as briefly noted above and discussed further below, the definition of the pole mass is contaminated by non-perturbative corrections. These are reflected in the divergence of the perturbation theory and power corrections are, in practice, difficult to disentangle from higher-order corrections, so extractions at different orders will have a variation from this source. We now examine this issue in more detail.
Let us consider the values we extract and use as default. In the Review of Particle Physics [65], the authors quote values for the "running" masses, m c (µ = m c ) and m b (m b ), in the MS scheme: which they evaluate from a whole series of independent determinations. This implies a precise determination of masses. 6 The best precision of the pole mass should be obtained by using this MS value (the best determinations coming from NNNLO calculations) with the most accurate conversion factor, also at NNNLO. However, the conversion of the "running" mass to the pole mass is problematic, since it relies on a very weakly convergent perturbative expansion. Using our NNLO value of α S , the term-by-term conversion factor (obtained from the formula in Ref. [65]), starting at zeroth order and going to the highest known order of α 3 S [68,69], for the bottom-quark mass is (1 + 0.095 + 0.045 + 0.036) and for charm is (1 + 0.16 + 0.14 + 0.18). 7 The former implies that convergence is ceasing, while the latter has no real convergence at all, since it relies on a much higher value of α S (and larger coefficients due to less cancellation between individual terms as the number of light quark flavours decreases [68,69]). Indeed, the contamination by renormalons of the perturbative series for the transformation of a running quark mass in the MS scheme to the pole mass has been known for a long time [71,72]. These papers estimate the intrinsic ambiguity of the perturbative series due to infrared contributions, translating into an uncertainty on the pole mass of at least 0.05 GeV, but more likely of order 0.1-0.2 GeV. The best estimate one can obtain from the perturbative series 6 There are numerous individual determinations of both charm-and bottom-quark masses in the MS scheme which are actually even more precise than the above-for recent examples see Refs. [66,67]. However, these can vary from each other by considerably more than the quoted uncertainties, so in the Review of Particle Physics [65] account is taken of the spread, as well as the individual uncertainties. 7 Use of an NNNLO value of α S would make little difference. As shown in Ref. [70], NNNLO extractions of α S (M 2 Z ) are very similar to the corresponding NNLO values, and the slightly quicker running at NNNLO would marginally increase α S at the scale of the charm-and bottom-quark masses, i.e. the convergence of the perturbative series would likely be very slightly worse. is usually defined by including all terms in the series until they cease to fall in magnitude from one order to the next. For the bottom quark, under the plausible assumption that the unknown O(α 4 S ) term is similar in size to the calculated O(α 3 S ) term, this would give a pole mass of m b = 4.9 GeV. The uncertainty on this value from the conversion is either the estimate from Ref. [72] or roughly the size of the last term included in the series, i.e. 0.15 GeV. The good correspondence of the two estimates of the uncertainty is expected, and suggests that the series is indeed truncated at the correct point. If this uncertainty is added (approximately) in quadrature with the uncertainty on the MS scheme mass, Eq. (6), we obtain m b = 4.9±0.2 GeV.
The determination of the charm-quark pole mass by conversion from the MS scheme mass is more problematic since, as noted above, the perturbative series displays no clear convergence at all. However, as discussed in Ref. [73], the leading renormalon contribution is the same for different masses, i.e. δm 2 ∝ Λ 2 QCD , independent of quark flavour. This means the difference between m b and m c can be well determined, and turns out to be 3.4 GeV with a very small uncertainty [73]. Using this result we obtain a best determination of the charm-quark pole mass of m c = 1.5 ± 0.
are slightly high compared with our default values (1.4 GeV and 4.75 GeV) and our best-fit values of m c , particularly at NNLO, but even in the most discrepant case the error bars overlap. Combining our best-fit m c values and the pole mass determination given in Eq. (7), it seems that a range for m c of about 1.25-1.55 GeV at 68% C.L. and 1.15-1.65 GeV at 90% C.L. would seem reasonable at present. With improved HERA data on F c 2 to come [74] this range can hopefully be narrowed in the near future. Turning now to the bottom mass, the limited information that we obtain from comparison to data, and from the determination of the pole mass given in Eq. (7), suggest that a range of m b from about 4.65-5.05 GeV at 68% C.L. would seem sensible. As we will see in the next section, LHC cross sections that do not involve explicit b-quarks are not very sensitive to the PDF variation with m b , so using our grids at 4.5 and 5 GeV should give a good estimate of the uncertainty due to m b at 68% C.L., and the grids at 4.25 and 5.25 GeV should give a conservative estimate of this source of uncertainty at 90% C.L. As with charm, improved HERA data on F b 2 are likely to limit the allowed spread in future. If we reach the stage where we become confident that the data are starting to constrain the pole masses to an accuracy clearly better than the renormalon ambiguity, it may become preferable to transform to a scheme where MS definitions are used instead, even though the mass is less directly related to a physical variable in this case. This would be appropriate if we become strongly constrained by data with Q 2 ∼ m 2 h and data close to the kinematic threshold, W 2 = Q 2 (1/x − 1) = 4m 2 h (for neutral-current DIS), where non-perturbative effects 15 and the interplay between the leading-twist power-series and the power corrections become very important, and the ambiguities may be reduced in a different renormalisation scheme. Indeed, we have already noted in Section 9.2 of Ref. [1] that the lowest-Q 2 EMC data on F c 2 [75] imply a non-perturbative correction to the cross section.
3.4 Impact on W , Z and Higgs production at the Tevatron and LHC Tables 3 and 4 show how the predictions for the "standard candle" NNLO W and Z cross sections at the Tevatron and LHC change using the PDFs obtained when the charm-and bottom-quark masses are varied in the global fit. The cross sections are calculated in the 5flavour ZM-VFNS as described in Section 15 of Ref. [1], e.g. using the PDG 2008 [65] electroweak parameters. These changes are not due primarily to the changes in the heavy-quark PDFs themselves, which would in any case be in the opposite direction. Rather they are due to the changes in the light quarks, which evolve slightly more rapidly at small x to compensate for the slower turn-on of the heavy-flavour contribution to the structure functions when the heavy-quark masses are increased. In the free coupling case this is achieved mainly by an increase in α S , while for fixed coupling it occurs from an increase in the small-x input gluon distribution. Additionally, the input sea quarks, mainly the less well-constrained strange quark distributions, also increase with increasing m h to compensate for the suppressed heavy-flavour contribution to small-x structure functions. At the Tevatron, varying the charm mass by ±0. In Tables 3 and 4 we also give the changes to the cross section for the production via gluongluon fusion through a top-quark loop of a Standard Model Higgs boson of mass 120 GeV, calculated in the 5-flavour ZM-VFNS as described in Section 6.2 of Ref. [6]. From Table 3 we see that the variation of the Higgs cross sections with m c is different depending on whether α S is (a) allowed to vary or (b) kept fixed. This can be understood from the behaviour of the gluon distribution as m c is varied, see Fig. 7(b). When α S is fixed, the change in the Higgs cross section reflects the change in the gluon distribution at the relevant x values (indicated). Thus the {14 TeV LHC, 7 TeV LHC, Tevatron} Higgs cross sections are {correlated, uncorrelated, anticorrelated} with m c . On the other hand, for variable α S , the additional correlation of α S with m c (see Table 1) results in {14 TeV LHC, 7 TeV LHC, Tevatron} Higgs cross sections that are {correlated, correlated, almost uncorrelated} with m c , i.e. at the Tevatron the anticorrelation of the gluon distribution largely cancels the correlation of α S .
From Table 4 we see that the W , Z and Higgs production cross sections are much less dependent on the value of m b . As seen from Table 2, the correlation between m b and α S (    18  Table 2. is negligible, so we do not explicitly show results with fixed α S (M 2 Z ) since they would be almost the same as in Table 4. In all cases, for both varying-m c and varying-m b , the cross-section ratios σ W /σ Z (and also σ W + /σ W − at the LHC) are almost completely unaffected.
In order to assess the total uncertainty on a hadron collider cross section arising from variations in PDFs, α S and the heavy-quark masses, we need to define a prescription for combining the latter with the former. Our recommendation is to vary m c in the range m c = 1.40±0.15 GeV at 68% C.L. and ±0. 25 GeV at 90% C.L., centred on the default value for fixed α S , and to vary m b in the range m b = 4.75 ± 0.25 GeV at 68% C.L. and ±0.50 GeV at 90% C.L., and then to separately add the resulting cross-section variations in quadrature with the usual "PDF+α S " uncertainty. Of course, this prescription does not account for all possible correlations between PDFs, α S and m c,b , but it should be a sufficiently good approximation. The range of m c values is based on the slightly contrasting pulls of the pole mass determination from the MS conversion and our fit. As an example, in Table 5 we show the increase in uncertainty on the W , Z and Higgs total cross sections due to m c and m b variations, compared to the "PDF only" [1] and "PDF+α S " [6] uncertainties.
In Table 6 we compare the values of the heavy-quark masses, m c and m b , used in the MRST/MSTW analyses, with the values taken in the most recent public PDF fits of other groups. The other two global PDF fitting groups, CTEQ [3,4] and NNPDF [5], do not attempt to quantify the uncertainty coming from m c,b , neither do the fits of "dynamical" PDFs by GJR/JR [77,78]. The HERAPDF1.0 fit [79] uses the same central values as MSTW08 [1], but additional PDF fits are also provided with m c = {1.35, 1.65} GeV and m b = {4.3, 5.0} GeV, and this additional model uncertainty is recommended to be added in quadrature with the other uncertainties. The ABKM09 analysis [80] uses fixed values of m c = 1.5 GeV and m b = 4.5 GeV to determine the central fit. However, additional pseudo-measurements of m c,b are then added, with values given in the last line of Table 6, and m c and m b are taken as free parameters to calculate the covariance matrix used for the final error propagation. This means that each of Tevatron,   Table 5: NNLO predictions for W , Z and Higgs (M H = 120 GeV) total cross sections at the Tevatron, 7 TeV LHC and 14 TeV LHC, with PDF uncertainties only [1], with the combined "PDF+α S " uncertainty [6], then finally also including the uncertainty due to m c and m b . The 68% C.L. uncertainties are given in all cases. We take µ R = µ F = M W,Z,H . 1.50 ± 0.10 4.50 ± 0.50 Table 6: Values of m c and m b used in various PDF fits. The 1σ uncertainties are given only for the three PDF groups who attempt to account for uncertainties on m c,b in their analyses.
the public eigenvector PDF sets will be associated with different values of m c and m b , but these values are not readily accessible.

3-flavour and 4-flavour scheme parton distributions
As well as looking at the variation in the parton distributions as a function of the heavy-quark masses, it is also interesting to consider the PDF sets obtained in the framework of a different maximum number of active quark flavours. Hence, in this section we will consider our PDFs when charm becomes an active parton but bottom does not-the 4-flavour scheme-and when both the charm and bottom quarks only appear in the final state-the 3-flavour scheme. We have argued on various occasions (e.g. in Section 4 of Ref. [1]) that the use of a GM-VFNS, i.e. with up to five active quarks (or even six if we include top), is preferable. However, there are cases where the cross section has only been calculated with finite mass dependence for the fixed flavour number case; see, for example, the hqvdis program [19,20] for details of final states in heavy-quark production in DIS. For this reason we make available sets of 3FS and 4FS distributions with a variety of both charm and bottom masses.

Obtaining 3-flavour and 4-flavour scheme parton distributions
It might be thought preferable to obtain these lower active quark number PDF sets from a fit performed using FFNS scheme coefficient functions. However, as argued in Refs. [9,81], it is not actually so obvious that this is the case. This is largely because rather few of the data sets included in a truly global PDF fit can be kept, even at NLO, in a fit using the FFNS, due to lack of the full coefficient functions (even charged-current DIS coefficients are not known with full mass dependence at order α 2 S except for Q 2 ≫ m 2 h [82]). Hence, the central values of the PDFs are likely to be influenced as much by the lack of data as by the change of scheme. However, it is also a consideration that the lack of resummation of the large logarithms in Q 2 /m 2 h potentially affects the stability of the fit compared to the presumably more stable GM-VFNS fit. Ultimately a correct GM-VFNS will provide results very similar to the FFNS near to the transition points Q 2 = m 2 h anyway, so we deem it best to simply obtain the 3FS and 4FS PDFs from the inputs for the full fits performed using the GM-VFNS.
When obtaining the 3FS and 4FS PDFs it is vital to make a self-consistent definition of the strong coupling constant α S . It is generally the case that coefficient functions calculated in the 3-or 4-flavour schemes are made in a renormalisation scheme where the contribution of the heavy quark decouples and the coupling itself does not include the heavy quark as an active flavour. On this basis we define the PDFs using this definition of the coupling in the splitting functions. It is certainly possible to use a different definition of the coupling, but this must be applied universally, in both the PDFs and the coefficient functions. As illustrated in Ref. [9], the error, made from not doing so, can be a few percent. Indeed, the change in the  10809. This is illustrated in Fig. 8, and more clearly in Fig. 9(a), which shows the ratio of the 3-and 4-flavour α S to the 5-flavour one.
We make available both the 3FS and 4FS NLO PDF sets for the full variety of fits with varying charm and bottom masses discussed in the previous section, including the case of varying m c with fixed α S (M 2 Z ). We do not vary the quark mass at LO since the quality of the fit is already poor at LO. However, we do provide the 3FS and 4FS PDFs for the default masses at LO since these may be useful in some Monte Carlo event generators. We also make available the full range of 3FS and 4FS sets at NNLO, although the application of these is distinctly limited at present due to the lack of cross sections with full mass dependence calculated at this order. We will comment more on this in the following section. The ratios of the 3FS and 4FS NLO PDFs to those in the 5FS are shown in Fig. 10 at Q 2 = 10 4 GeV 2 . In both cases, the gluon is larger in the 4FS and, particularly, in the 3FS, due to the fact that it splits into a smaller number of quark-antiquark pairs; see also Fig. 9(b). The effect of the increased growth of the gluon distribution is countered exactly for the leading term in α S (Q 2 ) ln(Q 2 /m 2 h ) by the quicker decrease of the coupling. This is illustrated in Fig. 9, where the change in the gluon distribution comparing different flavour numbers is, to a good approximation, the inverse of  that for α S . This behaviour is only violated by higher-order corrections, so the gluon-driven change in the light quarks at small x is rather minimal. At higher x, the onset of the smaller coupling in the lower flavour number schemes means a slower evolution and consequently larger light-quark distributions, as seen in Fig. 10.
Finally, we also make available new eigenvector PDF sets for the 3-and 4-flavour PDFs using the default quark masses at LO, NLO and NNLO. These are evolved from the saved PDF parameters at Q 2 0 = 1 GeV 2 for the default MSTW fits, i.e. the "dynamical" tolerance values, T = (∆χ 2 global ) 1/2 , are those determined from the 5-flavour fit [1]. For the first time, this will allow PDF uncertainties to be consistently included in a 3FS or 4FS calculation. We have not yet generated the additional eigenvector PDF sets with varying α S values needed for the "PDF+α S " uncertainty calculation [6], although these could be provided at a future point if necessary. Note, however, that 3FS or 4FS calculations generally have large higher-order corrections, and an associated large (renormalisation and factorisation) scale dependence, which is common in processes with multiple hard scales. The theory uncertainty due to neglected higher-order corrections is therefore likely to dominate over any PDF or "PDF+α S " uncertainty. Precise calculations of the latter quantities are therefore relatively less important than in typical 5FS calculations where higher-order corrections are more readily available.

Comparison of Z total cross-section predictions in 4FS and 5FS
We have already noted in Section 2 that while the parton distributions of a nFS, with appropriate coefficient functions, can give an adequate description of the production of an (n + 1)th parton near threshold, the accuracy of the perturbative expansion becomes increasingly unreliable as the scale of the hard scattering process, Q 2 , increases above m 2 h . In this section we illustrate this by considering the total Z cross section at hadron colliders. In particular, we compare the cross-section predictions in the 5FS (≡ MSTW 2008 [1]) and 4FS defined previously. Since in this case Q 2 = M 2 Z ≫ m 2 b , we would expect to see quantitative differences between the predictions due to the resummed [α S log(M 2 Z /m 2 b )] n terms which are implicit in the 5FS, via the evolution of the b-quark PDF and α (5) S , but absent from the 4FS. For reasons which will become apparent below, we consider the Z cross section at both NLO and NNLO.
In the 4FS the b quarks are not considered as partons in the initial state, but contribute to the Z cross section via real and virtual contributions which first appear at O(α 2 S ) in perturbation theory. Therefore the only difference in the predicted cross sections at NLO are (i) explicit bparton contributions b+b → Z(+g) and g +b(b) → Z +b(b) which only contribute in the n f = 5 scheme (5FS), (ii) small differences in the light-quark and gluon distributions arising from the slightly different evolution in the two schemes, see S (M 2 Z ), see Figs. 8 and 9(a), which affect the size of the NLO K-factor. We would expect the Tevatron Z cross sections to be more similar than those at the LHC, since the b-quark contributions to the 5FS cross section are much smaller at the lower collider energy. We do not explicitly include PDF and/or α S uncertainties on the calculated cross sections, as Tevatron,     these are presumed to be more or less the same for the 4FS and 5FS calculations. Contributions from top quarks are also not included, as these have been shown to be very small in Ref. [83], and in any case should be added to both the 4FS and 5FS cross sections.
The results at NLO are shown in Table 7. 8 The 4FS cross section is smaller by 1.0%, 2.8% and 3.9% at the Tevatron, LHC(7 TeV) and LHC(14 TeV) respectively. The O(α 1 S ) contributions differ by more than the O(α 0 S ) contributions, due to the ∼ 5% differences in the α S (M 2 Z ) values in the two schemes. Also shown in Table 7 (final column) are the contributions to the 5FS cross sections from the b +b → Z(+g) and g + b(b) → Z + b(b) processes, which evidently account for 0.4%, 2.0% and 3.0% of the total at the Tevatron, LHC(7 TeV) and LHC(14 TeV) respectively. Thus these contributions account for the bulk of the differences in the 4FS and 5FS cross sections at the LHC energies.
As noted above, when making a comparison between the 4FS and 5FS calculations at NNLO, we must include in the former the explicit real and virtual b-quark contributions which first appear at O(α 2 S ) in perturbation theory. When calculated with a non-zero b-quark mass 9 these contributions are finite, and simply add to the 4-flavour contributions. Sample Feynman   ) contributions to the total Z 4FS NNLO cross section in nb (multiplied by leptonic branching ratio) at the Tevatron and LHC arising from real and virtual b-quark processes.
diagrams for the various real and virtual b-quark contributions are shown in Fig. 11, and analytic expressions are given in Refs. [83,84]. We summarise these in Appendix A, and derive small-mass (i.e. m 2 b /M 2 Z ≪ 1) expansions which are useful in practice. The results are presented in Table 8 Table 9 we add these real and virtual b-quark contributions to the bulk of the NNLO cross section coming from 4FS light quarks and gluons, and compare the total with the benchmark 5FS NNLO results presented in Ref. [1], in which the b quark is treated as a massless parton in the subprocess cross sections, i.e. both in the initial and final states and in loop contributions. Also shown in Table 9 are the (bb, bg, . . .) contributions in which the Z couples directly to a b quark and in which there is at least one b orb quark in the initial state. 10 These represent 0.03%, 1.5% and 2.4% of the total 5FS NNLO cross sections at the Tevatron, LHC(7 TeV) and LHC(14 TeV) respectively.
From Table 9, we see that at the Tevatron the 5FS and 4FS+b-quark total cross sections are the same to within 1%. At the LHC, however, the additional b-quark contributions to the 4FS cross section are at the +1% level and so do not completely compensate the "missing" bb → Z contributions to the 5FS cross section. To be precise, at √ s = 14 TeV (7 TeV), which results in the "full" 4FS total cross section being 2.3% (1.7%) smaller than the the 5FS cross section. We interpret this as meaning that the DGLAPresummed [α S ln(M 2 Z /m 2 b )] n contributions that are absorbed into the b parton distribution are numerically important.
Our results are in broad agreement with the study of Ref. [85], in which the different calculational methods for heavy particles produced in association with b quarks, gg → Xbb (4FS) versus bb → X (5FS), were studied in detail for X = H, Z. In particular, it was shown Tevatron,   that with the canonical choice of scale µ = M Z , the LO gg → Zbb cross section at the LHC was a factor of two smaller than the NNLO bb → Z cross section, consistent with our results in Eq. (8) above. It was also shown in Ref. [85] that the agreement between the 4FS and 5FS calculations improves if the scale µ is reduced to around M Z /3 (found by choosing µ ∼ √ −t near the end of the collinear plateau in the quantity −t dσ/dt for the process gb → Zb). The NNLO bb → Z cross section is approximately scale independent, while the LO gg → Zbb cross section increases with decreasing scale, primarily because of the overall α 2 S factor. Of course the explicit gg → Zbb contribution corresponds only to the n = 1 term in the resummed [α S ln(M 2 Z /m 2 b )] n perturbation series implicit in the bb → Z 5FS calculation. Complete agreement between the two schemes would be obtained in a fully all-orders perturbative calculation. Note that since at all collider energies ∆ b σ Z is dominated by the contributions involving the Zbb final state, we would expect that higher-order corrections to the Zbb production process, which is here calculated only at leading order, will generate the [α S ln(M 2 Z /m 2 b )] n terms implicit in the b-quark PDF. The NLO (i.e. O(α 3 S )) corrections to the Zbb total cross sections with m b = 0 have recently been calculated [86][87][88][89]. Although the results presented in Ref. [89] impose a minimum p b T (> 5 GeV), the K-factor 11 is evidently rather independent of p b T at small p b T , suggesting that the NLO/LO K-factor for the fully inclusive Zbb cross section is approximately 1.5 for the LHC at √ s = 14 TeV. It is therefore plausible that even higher-order perturbative corrections can account for the factor of 2 difference in the 4FS and 5FS cross sections, Eq. (8). This conclusion is supported by the fact that the scale dependence of the 4FS calculation for Zbb production at NLO is only mildly weaker than at LO [86][87][88][89].
In Ref. [91] it was shown that for x in the region of 0.01-0.05, the relevant region for Z + bb production at the 14 TeV LHC, the ratio of the GM-VFNS structure function F h 2 to the FFNS structure function was ∼ 1.5 at LO at high scales. This represents the effect of either resumming the [α S ln(Q 2 /m 2 h )] n contributions, or keeping only the contribution of fixed order in α S for one parton. For hadron-hadron processes we would expect the difference to be about 1.5 2 > 2, exactly as observed. At NLO for structure functions at this x the ratio is reduced to ∼ 1.1, so inclusion of the extra ln(Q 2 /m 2 h ) removes much of the discrepancy present at LO. However, for hadron-hadron processes, NLO in the fixed flavour scheme only contains the extra large logarithm for one of the two incoming partons, so the ratio between the 5FS and 4FS would be roughly 1.5 × 1.1 ≈ 1.6, again as we expect to see in practice. It would only be at NNLO in the 4FS, when the double logarithm for both incoming PDFs is included, that we would expect to see the reduction to roughly 1.1 2 ≈ 1.2 in the ratio of the 5FS to 4FS cross sections. This is a general feature, i.e. the 4FS (or 3FS) will converge to the resummed 5FS results more slowly for hadronic processes than for those in DIS.
In summary, the 5FS PDFs are clearly the most appropriate to use for inclusive quantities such as σ Z (or σ dijet , etc.) at high Q 2 scales, where resummation of the [α S ln(Q 2 /m 2 b )] n contributions is evidently important. However, for more exclusive quantities like σ Zbb , where 11 The K-factor is here defined as the ratio of cross sections calculated using the dynamical scale µ 2 = M 2 Z + (p b,1 T ) 2 + (p b,2 T ) 2 , and with the standard 5-flavour CTEQ6M/CTEQ6L1 PDFs [90] at NLO/LO. the b quarks are measured in the final state, the 4FS parton distributions are more appropriate since the 5FS calculation gives no information on the spectator b andb quarks which must be present in the final state for the b-andb-initiated processes. Note that if only the total cross section is required, without cuts imposed on the b-quarks, then a 5FS is still better, e.g. for Zbb a 5FS calculation can be used for bb → Z at NNLO, where the b-quarks couple directly to the Z, and so there are implicitly also two b-quarks in the final state [85]. However, if cuts must be applied to the b-quarks, as is the case in the experimental measurement, then a 4FS calculation is more appropriate. Similar remarks apply to the calculation of Hbb production [92][93][94] and other processes where b-quarks are detected in the final state. In a recent study [95], the production of a charged Higgs boson in association with a t quark was considered in both the 4FS (gg → tbH − etc.) and 5FS (gb → tH − etc.) to NLO in pQCD, using the appropriate MRST 2004 PDF sets [9,10]. The central predictions in the 5FS were shown to be approximately 40% larger than those in the 4FS. Even taking the scale uncertainty into account the 4FS and 5FS NLO cross sections are barely consistent.
An ideal calculation would combine the best features of the 4FS and 5FS so as to resum [α S ln(Q 2 /m 2 b )] n terms while also retaining a finite m b dependence in the partonic cross section (rather than setting m b to zero as done in the 5FS). This matching has, of course, been done for structure functions in DIS using different variants of the GM-VFNS, but applications to hadron collider cross sections are more involved and have so far been limited (see Ref. [96] for an application of the GM-VFNS to the p T spectrum in heavy-flavour hadroproduction). However, for processes where the hard scale is, for example, Q 2 ∼ M 2 Z , then the GM-VFNS calculation will differ from the ZM-VFNS (5FS) only by terms O(m 2 b /M 2 Z ∼ 0.3%). We would therefore expect the complete GM-VFNS calculation to give results very close to the pure ZM-VFNS (5FS) for the total cross section.
Note that rather than producing separate 4-flavour PDFs for use in a 4FS calculation, an alternative approach (e.g. used at NLO in Refs. [97,98]) is to use the conventional 5-flavour PDFs, then pass to the 4-flavour scheme using counterterms given in Ref. [96]. However, these counterterms [96] are equivalent to using the inverse of transition matrix elements, but only out to order α S . One could indeed use the transition matrix elements themselves to go from 4-flavour to 5-flavour PDFs, except this would not sum the logarithmic terms, [α S ln(Q 2 /m 2 b )] n , in the PDF evolution. 12 Hence, the use of counterterms [96] is a less complete way of going from a 5FS to a 4FS, and instead, we recommend that dedicated 4-flavour PDFs be used in 4FS calculations. 13 Previously, a major advantage of using 5-flavour PDFs with counterterms was that eigenvector PDF sets to calculate PDF uncertainties were not made available for existing 4-flavour PDFs [9]. However, we have now provided eigenvector PDF sets also for the 4-flavour PDFs, therefore this advantage no longer holds.

Conclusions
We have repeated the NLO and NNLO MSTW 2008 global PDF analyses [1] for a range of heavy-quark masses about their "default" values m c = 1.40 GeV and m b = 4.75 GeV. For the charm quark, we found that the global data prefer the values m c = 1.45 (1.26) GeV at NLO (NNLO). The most discriminating data are, as anticipated, the HERA data for F c 2 [44][45][46][47][48][49][50]. On the other hand, for the bottom quark, the data included in the global fit (excluding F b 2 ) do not put a meaningful constraint on the value of m b , while the HERA F b 2 data slightly favour m b ≈ 4.75-5 GeV. We pointed out that precise determinations of the heavy-quark masses in the MS scheme are affected by poorly convergent perturbative series in the conversion to the pole masses, particularly for the case of the charm quark. Recent precise combined HERA data on σ NC r [79] and F c 2 [74] will in future be able to narrow the favoured range of the charm-quark pole mass m c . Note, however, that uncertainties from the choice of GM-VFNS [64] mean that the favoured value of m c will be correlated to some extent with the particular choice of GM-VFNS, although this correlation will be much smaller at NNLO than at NLO, as will other uncertainties arising from the choice of GM-VFNS [64].
We explored the effect of the values of the heavy-quark masses on W , Z and Higgs production at the Tevatron and LHC. Varying the charm mass by ±0.15 GeV changes the cross sections by ±1% or less at Tevatron energies and by ±2% at the LHC energy of √ s = 14 TeV.
The various weak boson cross-section ratios are essentially unchanged. The predictions for W , Z and Higgs cross sections are much less dependent on the value taken for m b . We provided a recommendation on how to include the uncertainty arising from the choices of m c and m b in a generic cross-section calculation.
We also presented PDF sets obtained in a framework with different active numbers of quark flavours, as done previously in the context of the MRST 2004 analysis [9]. Explicitly, we determined 4-flavour PDF sets in which charm becomes an active parton, but bottom does not, and 3-flavour PDF sets where charm and bottom are not partons, but only appear in the final state. The analogous 5-flavour parton sets are simply those of MSTW 2008 [1]. Of course, the latter, which in the absence of top corresponds to PDFs of a (general-mass) variable flavour number scheme, are generally to be preferred, particularly for inclusive quantities at high Q 2 scales where the resummation of the [α S ln(Q 2 /m 2 b )] n contributions is essential. However, for more exclusive processes, such as Zbb, Hbb, . . ., where b-quarks are observed in the final state, the 4-flavour parton distributions are more appropriate. For illustration, we computed the various components of the Z production cross section to O(α 2 S ) at the Tevatron and LHC, and compared the predictions obtained using the 4FS and 5FS (≡ MSTW 2008 [1]) parton sets.
The additional grids for all PDF sets discussed in this paper are made publicly available [60], for use either with the standalone MSTW interpolation code or via the lhapdf interface [99]. To be precise, grids for the following PDF sets are made available: • For the default quark masses, m c = 1.40 GeV and m b = 4.75 GeV, we provide LO, NLO and NNLO grids for 3-and 4-flavour PDFs (central set and 40 eigenvector sets at both 68% and 90% C.L.). These grids complement the existing grids for the 5-flavour PDFs.
• For m c in the range 1.05 GeV to 1.75 GeV (in steps of 0.05 GeV), we provide NLO and NNLO grids for 3-and 5-flavour PDFs (central set only) for both free and fixed α S (M 2 Z ).
• For m b in the range 4.00 GeV to 5.50 GeV (in steps of 0.25 GeV), we provide NLO and NNLO grids for 4-and 5-flavour PDFs (central set only) for free α S (M 2 Z ) only.
These additional grids should prove to be useful in future for detailed studies of a variety of collider processes involving heavy quarks.

A Appendix
We consider the b-quark contributions to the 4FS calculation, closely following the discussion in Ref. [83]. First, we have the two-loop corrections to the q +q → Z Born process, examples of which are shown in Fig. 11(a) and (b). In the notation of Ref. [83], with ρ b = m 2 b /M 2 Z ≈ 2.7 × 10 −3 . The functions F and G 1 are defined in Refs. [83] and [84] respectively. Since in practice ρ b ≪ 1, we can use the following small-mass expansions: with L = ln(1/ρ) ≈ 5.9 for ρ = ρ b . The contributions to the total Z cross section from these loop diagrams are given in the first row of Table 8. They are numerically very small and negative (positive) at the LHC (Tevatron). 14 Next we have the one-loop corrections to the 2 → 2 processes q +q → Z + g and g + q(q) → Z + q(q), examples of which are shown in Fig. 11(c) and (d): whereŝ =τ −1 M 2 Z and the functions H 1 , J 1 and J 2 are defined in Ref. [84]. The following small-mass approximations are again useful: Finally, we have the 2 → 3 processes q +q → Z + b +b, g + g → Z + b +b, (A.7) see Fig. 11(e) and (f). The calculation of these is in principle straightforward, with matrix elements squared integrated over three-body phase space. The non-zero b-quark mass makes these contributions infra-red and collinear finite, i.e. the singularities present for massless quarks are replaced by terms with logarithmic (ln(M 2 Z /m 2 b )) behaviour in the small m b limit. We use the mcfm implementation [100] of the matrix elements and phase space, with exactly the same parameter choice as in the previous MSTW calculations. The results are given in Table 8. Note that the gg contribution is larger than the qq contribution at the LHC, while the converse is true at the Tevatron.