Bottom-quark production at hadron colliders: fully differential predictions in NNLO QCD

We report on the first fully differential calculation of the next-to-next-to-leading-order (NNLO) QCD radiative corrections to the production of bottom-quark pairs at hadron colliders. The calculation is performed by using the $q_T$ subtraction formalism to handle and cancel infrared singularities in real and virtual contributions. The computation is implemented in the Matrix framework, thereby allowing us to efficiently compute arbitrary infrared-safe observables in the four-flavour scheme. We present selected predictions for bottom-quark production at the Tevatron and at the LHC at different collider energies, and we perform some comparisons with available experimental results. We find that the NNLO corrections are sizeable, typically of the order of $25$-$35\%$, and they lead to a significant reduction of the perturbative uncertainties. Therefore, their inclusion is crucial for an accurate theoretical description of this process.

At the theoretical level, heavy-quark production at hadron colliders is one of the most classic tests of perturbative QCD. The cross section to produce a pair of heavy quarks with mass m Q is computable as a power series expansion in the QCD coupling α S (µ R ), where the renormalisation scale µ R has to be chosen of the order of m Q . In the case of the bottom (b) quark the relatively low mass, m b ∼ 4 − 5 GeV, leads to a slow convergence of the perturbative expansion and, therefore, to large theoretical uncertainties.
Theoretical predictions for bb production at next-to-leading order (NLO) in QCD have been available for a long time [19][20][21], including calculations [22] of bb correlations and generic infrared safe (IR) observables. NLO studies based on different schemes for the renormalisation of the bottom-quark mass are presented in Ref. [23]. At high transverse momenta p T of the bottom quark, large logarithmic terms of the form ln(p T /m b ) need to be resummed to all perturbative orders [24]. In the case of single-particle (b quark or antiquark) inclusive cross sections the resummation can be performed by introducing the perturbative fragmentation function [25] of the bottom quark (which can also be supplemented with all-order soft-gluon effects [26]). Predictions obtained by matching such resummed computations to the NLO calculation (the so called "FONLL" prediction) [27][28][29][30] have become the standard reference for the comparison with experimental data. Perturbative predictions for the inclusive production of bare bottom quarks can then be folded with non-perturbative functions [31,32] describing the fragmentation into the triggered b-hadrons. The parameters that control such fragmentation functions are typically extracted from LEP data (see, e.g., Ref. [33]). The variable-flavournumber scheme [34,35] is another procedure that is used to combine the NLO calculation with high-p T resummation effects for single-inclusive b-hadron production, and ensuing data-theory comparisons have been presented in the literature (see, e.g., Refs. [36,37]). Phenomenological studies on the impact of bottom production measurements on parton distribution functions have been presented in Refs. [38][39][40].
At the next-to-next-to-leading order (NNLO) in QCD, theoretical predictions for bb production are available only for the total cross section. Indeed, the bb results can be directly derived by exploiting the corresponding NNLO theoretical calculation [41][42][43][44] of the total cross section for top-quark pair production. The calculation of the bb total cross section at NNLO is implemented in the numerical program Hathor [45,46].
In this paper we report on the first fully differential NNLO QCD calculation for bb production at hadron colliders, and we also present comparisons with some inclusive b-hadron data obtained at the Tevatron and at the LHC. Our computation, which follows the analogous calculation carried out for top-quark pair production [47][48][49], is implemented within the Matrix The NNLO differential cross section for bottom-pair production, dσ bb NNLO , is obtained within the q T -subtraction method according to the following main formula, where dσ bb+jet NLO represents the bb+jet cross section at NLO accuracy, which we evaluate by using the dipole subtraction method [60][61][62]. As is customary in the Matrix framework, the NLO cross section dσ bb NLO is also computed by using dipole subtraction, while q T subtraction is actually used only to evaluate the NNLO correction dσ bb NNLO − dσ bb NLO . The expression in Eq. (1) is completely analogous to the tt case [48]. We remind the reader that the term in the square bracket of Eq. (1) is formally finite in the limit q T → 0 (q T is the transverse momentum of the bb pair), but each of the two contributions in the square bracket is individually divergent. Therefore a technical cut is introduced in the dimensionless quantity r = q T /m bb (m bb is the invariant mass of the bb pair). The r cut → 0 extrapolation is performed following the procedure of Ref. [50], and it is also applied in the computation of differential distributions on a bin-by-bin basis.
The core of Matrix is the Monte Carlo program Munich 1 , which contains a fully automated implementation of the dipole subtraction method for massless [60,61] and massive [62] partons and a general implementation of an efficient phase space integration. The required spin-and colour-correlated tree-level and one-loop amplitudes are obtained by using Open-Loops [51][52][53], with the exception of the four-parton tree-level colour correlators for which we rely on an analytic implementation. The two-loop amplitudes are obtained via an interpolation routine, based on the numerical results presented in Refs. [54,63].

Results
In the following we present our results for total cross sections and differential distributions for bb production at the Tevatron ( √ s = 1.96 TeV) and at the LHC ( √ s = 7 and 13 TeV). We use the NNPDF31 parton distribution functions (PDFs) [64] with n f = 4 massless-quark flavours and the value of the QCD coupling α S (m Z ) = 0.118. The pole mass of the bottom quark is fixed to m b = 4.92 GeV, consistently with the value that is used in the chosen PDF set. Predictions at the n-th perturbative order are obtained by using the corresponding N n LO PDF set and the evolution of α S at (n + 1)-loops. Since the NNPDF31 set is not available at LO with n f = 4, we instead use the corresponding NNPDF30 set [65] for our LO predictions. Perturbative uncertainties are estimated with the customary 7-point variation of the renormalisation (µ R ) and factorisation (µ F ) scales by a factor of two around a common central value µ 0 , with the constraint 0.5 ≤ µ R /µ F ≤ 2. The value of µ 0 is chosen of the order of the characteristic hardscattering scale of the process. The total cross section is controlled by scales of the order of m b , while each differential distribution is characterised by a different hard-scattering scale that has to be specified accordingly. The fully differential nature of our calculation also allows us to use dynamic scales.
We point out that starting from NNLO the inclusive production of a bb pair receives contributions from tree-level diagrams with an additional bb pair in the final state. Since the bottom-quark mass is kept non vanishing, these four-bottom contributions are separately finite and may be included or not in the perturbative calculation, according to the actual (theoretical and experimental) definition of the inclusive cross section. We find that these contributions are generally rather small. They are completely negligible at the Tevatron while at the LHC they typically contribute at the per mille level, reaching about 0.5% at large transverse momenta of the bottom quarks. Since the theoretical uncertainties affecting bottom production are much larger, these contributions are neglected in the rest of the paper.
Before presenting our numerical predictions, we discuss the impact of the two-loop virtual corrections on our NNLO results. The IR finite part of the two-loop contribution is obtained by using the numerical results of Ref. [54], which are provided through a 80 × 21 dimensional grid in the Born level kinematical variables β and cos θ (β and θ are the heavy-quark velocity and scattering angle in the partonic centre-of-mass frame). The grid is then interpolated by using splines, and the ensuing result is supplemented with the analytical results at small and large β [54,66]. The grid of Ref. [54] has, to date, been used only for top-quark production. Since the bottom-quark mass is significantly smaller than the top-quark mass, one may wonder whether the small-angle (collinear) region is sampled sufficiently well in our calculation. We have studied the effect of the IR finite part of the two-loop contribution on our calculations of the total cross section and several differential distributions for bb production. Moreover, to quantify the impact of having a discrete grid, we have repeated our calculations by using only half of the available grid points. We find that the two-loop virtual contribution [54] is below 1% in all cases. The differences obtained by reducing the number of points in the grid by a factor of two are typically at the per mille level for all the distributions we have considered. We conclude that we can safely use the results of Ref. [54] to carry out our fully differential calculation of bb production.

Total cross section
We start the presentation of our results by considering the total cross section. The total cross section for the production of a pair of bottom quarks is controlled by scales of the order of the bottom mass m b . Accordingly, we will consider the two central scales µ 0 = m b and µ 0 = 2m b .
As discussed in Sec. 2, our NNLO results are obtained through an r cut → 0 extrapolation procedure. The r cut dependence at the Tevatron with √ s = 1.96 TeV and at the LHC with √ s = 13 TeV is shown in Fig. 1 in the case µ 0 = m b . The ensuing NNLO cross section with its extrapolation uncertainty is compared with the corresponding result obtained with the numerical program Hathor [46], and the results of the two NNLO calculations are in good agreement. Similar results are obtained for all scale combinations of the applied 7-point variation, and also separated into partonic channels as in Ref. [47]. Fig. 1 shows that the extrapolation uncertainties are larger than the corresponding uncertainties in the case of toppair production [47], but, remarkably, still at the level of about 0.5%. This larger uncertainty in the NNLO result is simply due to the larger relative size of the O(α 4 S ) contribution when considering a lower quark mass. Nevertheless, this level of uncertainty is perfectly acceptable as it is well below other sources of theoretical (and experimental) uncertainties affecting bottomquark hadroproduction.
The NNLO results at the Tevatron and the LHC obtained with µ 0 = m b and µ 0 = 2m b pp → bb @ 13 TeV, µ0 = mb Figure 1: The dependence on r cut of the total cross section for bb production at different hadron colliders.  Table 1: The NNLO total cross section for bb production at the Tevatron and at the LHC: comparison of our results with those obtained by using Hathor. The quoted uncertainties are obtained through scale variations as described in the text. Numerical errors on the last digit are stated in brackets (and include the uncertainties due to the r cut → 0 extrapolation).
together with their scale uncertainties are reported in Table 1, compared to those obtained using Hathor [46]. For both central-scale choices the agreement is excellent, including scale uncertainties. We point out that the two computations are performed by using fully independent methods.
In Table 2 we report the LO, NLO and NNLO results for µ 0 = m b and µ 0 = 2m b . As in Table 1, the cross sections are presented with their perturbative uncertainties estimated through scale variations.
By inspecting the results presented in Table 2, we immediately see that QCD corrections are very large. In order to quantify their impact, we introduce K-factors, K (N)NLO = σ (N)NLO /σ (N)LO , defined as the ratios of the cross section predictions at two subsequent orders. The NLO K-factor K NLO tends to increase as the collider energy decreases. Specifically, the value of K NLO ranges between 1. 29 Table 1.
large QCD corrections, considerably larger for instance than the ones observed in top-quark production, are associated to the relatively low energy scales involved, which lead to large values of the strong coupling. The poor perturbative behaviour is also reflected by the large scale uncertainties that are observed in Table 2. It is important to remark, however, that in all cases the inclusion of the NNLO contribution allows us to strongly reduce the theoretical uncertainties. In addition, K NNLO is always smaller than K NLO , providing a sign of (slow) convergence of the perturbative series We also note that the values of σ NNLO and σ NLO are consistent within their scale uncertainites for each energy and choice of µ 0 (the values of σ NLO and σ LO are consistent as well).
We now turn to the discussion of the scale dependence of our results. Due to the overall proportionality of the cross section to the factor α 2 S (µ R ) at LO, the cross section typically decreases as µ R increases. On the contrary, the cross section generally increases as µ F increases. This is due to the fact that bottom-quark hadroproduction is sensitive to relatively low momentum fractions of the colliding partons, and in this kinematical region PDF scaling violations are typically positive. When performing scale variations, the dominant effect is given by variations of the renormalization scale µ R .
Comparing the predictions corresponding to the two different scale choices, we observe that the results obtained with µ 0 = 2m b present smaller scale uncertainties. Such difference is evident at LO and NLO, and less noticeable at NNLO, especially at the Tevatron. The choice µ 0 = 2m b leads to slightly smaller NNLO K-factors at the LHC (K NNLO = 1.27(1.26) to be compared with K NNLO = 1.31(1.34) at √ s = 7(13) TeV for µ 0 = m b ), and to larger K-factors at the Tevatron (K NNLO = 1.30 to be compared with K NNLO = 1.25 for µ 0 = m b ). Despite the differences between the results obtained with these two central scales, both choices provide predictions that are fully compatible within scale uncertainties, and are equally acceptable on Table 3: Total cross sections and uncertainties for bb production at NNLO for µ 0 = m b . The numerical errors on the last digits are stated in brackets, as in Table 1.
general grounds. The use of µ 0 = m b as central scale choice leads to larger uncertainties due to the lower values of µ R involved, but we do not observe a particularly worrisome perturbative behaviour for this scale choice, in spite of the low value of m b .
Besides the uncertainties related to missing higher-orders in the perturbative expansion, additional uncertainties on the bottom-quark cross section arise from the errors in the determination of the pole mass m b , the parton distributions, and the strong coupling α S (m Z ). In Table 3 we report our result for the NNLO cross section computed with µ 0 = m b including these additional uncertainties. As for the bottom mass, we follow Ref. [67], and we vary m b between m b = 4.79 GeV and m b = 5.05 GeV, corresponding to m b = 4.92 ± 0.13 GeV. The effect of changing the bottom mass is below 10% and slightly decreases with increasing collider energy. The PDF uncertainties are much smaller and increase with the collider energy. This is not unexpected: indeed, bb production at the Tevatron and the LHC is mainly sensitive to PDFs with momentum fractions . In this region the uncertainty in the gluon distribution increases as x decreases. As for the QCD coupling, we compute the corresponding uncertainty by evaluating the NNLO cross section with NNPDF31 NNLO PDFs obtained with α S (m Z ) = 0.119 and α S (m Z ) = 0.117. Naively, one could expect relatively large effects since the process starts at O(α 2 S ) and small variations on α S (m Z ) induce relatively large variations on α S (m b ). However, it is well known that the value of α S (m Z ) and the gluon density are correlated. Correspondingly, the ensuing effect on the NNLO cross section reported in Table 3 is rather small and roughly independent on the collider energy.
In summary, from the results in Table 3 we conclude that, despite the inclusion of the NNLO corrections and the consequent reduction of the scale uncertainties, the missing higher orders in the QCD perturbative expansion still represent the dominant source of theoretical uncertainty in bb production.

Differential distributions: Tevatron
We start the presentation of our differential results by showing selected distributions in pp collisions at the Tevatron ( √ s = 1.96 TeV). Throughout this paper we always consider observables obtained by averaging the corresponding bottom and antibottom quark distributions. In particular we consider the average of the transverse-momentum, rapidity and pseudorapidity distributions of the (anti)bottom quark, which are denoted as p T,bav , y bav and η bav , respectively. These observables are the parton-level equivalent of the corresponding b-hadron distributions.
In Fig. 2 we present our LO, NLO and NNLO predictions for the transverse-momentum distribution. The calculation is carried out without applying any kinematical cut on the finalstate partons. The transverse-momentum distributions of the bottom and antibottom quark are controlled by hard-scattering scales of the order of the corresponding transverse mass, Therefore, in Fig In addition to these two 'natural' scales, we also show predictions with the dynamic scale µ 0 = H T /2 (right panel in Fig. 2), defined as the average of the two transverse masses: Note that we have H T /2 = m T at LO. Differential results at each order in the perturbative expansion are shown in the upper panels of Fig. 2, while in the lower panels the ratio of each perturbative order to our NNLO prediction is presented. From Fig. 2 we see that the distribution is peaked at p T,bav ∼ 3 GeV. The average transverse momentum of the (anti)bottom quark is 4.6 GeV (at both NLO and NNLO), i.e., as expected, quite close to the value of the bottom-quark mass. As a consequence, the total cross section receives a small contribution from the large-p T,bav region. For instance, the region with p T,bav < 20 GeV gives about 99% of the total cross section.
For all the considered central-scale choices LO and NLO predictions are consistent within uncertainties only in the low-p T,bav region, where the bulk of the cross section is located, while they present very different shapes in the tail of the distribution, where the NLO corrections become more sizeable. In this region the LO and NLO scale uncertainty bands do not overlap. The inclusion of the NNLO corrections leads to a nice stabilisation of the perturbative result, analogously to what we have observed for the total cross section. In particular, the uncertainty at NNLO is smaller than at NLO, and the NLO and NNLO bands overlap in the entire region of transverse momenta.
By comparing the left and central panels in Fig. 2 we see that the choice of a smaller scale µ 0 = m T leads to a better overlap between the NLO and NNLO uncertainty bands, similarly to what was observed for the total cross section. In particular, the scale choice µ 0 = 2m T leads to a significantly worse perturbative convergence in the tail of the distribution. The dynamic scale µ 0 = H T /2 presents similar features to those observed for µ 0 = m T , consistently with the fact that both scales are equivalent at LO, with a good overlap of the NLO and NNLO bands for all values of p T,bav .
Perturbative predictions for bottom-quark production can be compared to experimental data for the inclusive production of b-hadrons. A precise comparison in the region of large transverse momenta of the bottom quark (i.e. p T m b ) would require the resummation of the logarithmically enhanced terms at large transverse momenta [24]. Such higher-order contributions are included up to next-to-leading logarithmic accuracy in the resummed predictions of Refs. [27,28]. The non-perturbative effects of the fragmentation of the bottom quark into the triggered b-hadron should eventually be accounted for by folding the perturbative result with an appropriate non-perturbative fragmentation function. In this paper we limit ourselves to considering perturbative predictions up to NNLO, and a thorough comparison with experimental data is beyond the scope of this work. In Appendix B we present a comparison of our NLO and NNLO results with the FONLL prediction from Ref. [68]. Such comparison shows that in the region of transverse momenta considered in this paper the resummation effects have a limited impact on the NLO transverse-momentum and (pseudo)rapidity distributions, and that our NNLO results have smaller perturbative uncertainties than the FONLL results and can thus be considered more accurate. The impact of non-perturbative fragmentation is typically rather small on p T -inclusive observables, definitely smaller than the NNLO perturbative uncertainties.
Having in mind the above issues and limitations, our fixed-order perturbative predictions can be compared with experimental measurements for inclusive b-hadron production.
In Ref. [5], the CDF collaboration performed a measurement of the J/ψ cross section in the rapidity range |y| < 0.6. The total b-hadron production cross section in the same rapidity region was extracted by using a Monte Carlo simulation of the decay kinematics of b-hadrons to the final states containing a J/ψ hadron. The result is In order to study the possible effect of a rapidity cut in our parton-level calculation, we have repeated the computation of the transverse-momentum distribution reported in Fig. 2   observe that the inclusion of the rapidity cut does not significantly modify the behaviour of the perturbative series, and that the shapes of the distributions remains rather similar.
We have computed the parton-level analogue of the b-hadron cross section in Eq. (4). In order to do so, we compute two independent cross sections, with cuts on the rapidity of the bottom or the antibottom quark, respectively. These two cross sections are then averaged. Since the characteristic scale for this observable is the bottom-quark mass m b , we have chosen µ 0 = m b as the central scale. We find and compare these results with the inclusive predictions presented in Table 2. The cross section turns out to be about 4 times smaller in the presence of this rapidity cut. The NLO and NNLO K-factors are close to those for the total cross section, and also the scale uncertainties are rather similar. This is a consequence of the fact that the impact of QCD radiative corrections is rather uniform in rapidity (see e.g. Fig. 5). Comparing our perturbative predictions with the CDF measurement in Eq. (4), we find that NNLO corrections considerably improve the agreement with the data. Scale uncertainties are, although still sizeable, largely reduced at NNLO, and only at this order they approach the size of the experimental uncertainties.

Differential distributions: LHC
We start the presentation of our differential results for the LHC by considering the transverse momentum distribution of the bottom quark. Using different central scales, such as µ 0 = m T , µ 0 = 2m T and µ 0 = H T /2, we obtain relative differences that are similar to those predicted at the Tevatron (Fig. 2). Therefore, we set µ 0 = m T , and in Fig. 4 we show our LHC results at We see that, as already observed at the Tevatron, LO and NLO predictions are consistent within uncertainties only in the low-p T,bav region, while they present very different shapes in the tail of the distribution, where the NLO corrections become very large. In this region the LO and NLO scale uncertainty bands do not overlap. The inclusion of NNLO corrections leads to a nice stabilisation of the perturbative result, analogously to what we have observed for the total cross section. In particular, the uncertainty band at NNLO is smaller than at NLO, and it overlaps with the latter over the entire region of transverse momenta. At small transverse momenta the NNLO scale uncertainty is larger than at the Tevatron, consistently with our observation for the corresponding total cross sections. On the contrary, at large transverse momenta the NNLO band is smaller at the LHC (note that the plots in Fig. 4 extend to p T,bav = 50 GeV, while the Tevatron result in Fig. 2 is shown up to 25 GeV).
The rapidity distribution of the bottom quark, computed with µ 0 = m b , is shown in Fig. 5 for √ s = 7 TeV (left) and 13 TeV (right). The two distributions at different LHC energies show a similar behaviour. The impact of QCD radiative corrections is uniform in rapidity and, therefore, it does not produce sizeable changes in the shape 3 of the distribution (even from LO to NLO). Since the radiative corrections are relatively flat, they are of the same size as those    for the total cross section. The NNLO results almost completely overlap with the NLO results, and they show smaller scale uncertainties, which are about ±30% over the whole spectrum.
In Fig. 6 we report analogous results to those of Fig. 5, but we consider the pseudorapidity (η) distribution of the bottom quark rather than its rapidity distribution. The most striking effect that we observe in going from Fig. 5 to Fig. 6 is the different shape (independently of the perturbative order) of the rapidity and pseudorapidity distributions. The rapidity cross section is maximal at y = 0 and monotonically decreases as y increases. The pseudorapidity cross section has a maximum at a finite (non-vanishing) value of η and a local minumum at η = 0, where it shows a 'central-pseudorapidity dip' (the value of dσ/dη at η = 0 is smaller than the value of dσ/dy at y = 0). These shape differences are a well-known fact: they have a general kinematical origin and are due to the non-vanishing mass of the produced particle. We recall the origin of these shape differences in Appendix A.
Examining the NLO and NNLO radiative corrections, we note that they have a highly similar effect (at both the qualitative and the quantitative levels) on the rapidity and pseudorapidity distributions of the b quark (see the ratio plots in Figs. 5 and 6). This feature is a consequence of the fact that dσ/dy and dσ/dη mainly probe the same underlying dynamics, and that their relative differences are basically of kinematical origin (see the discussion in Appendix A).
Despite the inclusion of NNLO corrections, the results in Figs. 4 and 5 show that perturbative uncertainties in transverse-momentum and rapidity distributions are still relatively large. A possible attempt to reduce scale uncertainties is to consider normalised distributions, as discussed in the following.
In Fig. 7 we consider the normalised p T,bav distribution at LO, NLO and NNLO. The scale uncertainty bands are obtained as follows: for each scale choice needed to study the 7-point scale variation the corresponding distribution is normalised to unity, and the envelope of each normalised distribution is constructed. The scale uncertainties in the peak region are strongly reduced, whereas at high transverse momenta the NLO and NNLO uncertainty bands are slightly larger than those observed in Fig. 4. This is not unexpected, since the low-p T,bav region gives the dominant contribution to the total cross section. Therefore, the scale uncertainties of the total and differential cross sections are strongly correlated at low p T,bav , and increasingly decorrelated as p T,bav increases.
We repeat the same procedure for the rapidity distribution. In Fig. 8 we show the normalised rapidity distribution of the bottom quark, which is constructed in the same manner as for Fig. 7. In this case, since the impact of QCD radiative corrections is rather uniform, the normalised distribution is quite stable and shows reduced scale uncertainties, except for the large-rapidity region.
An alternative strategy to reduce theoretical uncertainties is to consider ratios between distributions computed at different collider energies. In the context of heavy-quark production, the use of such ratios was proposed in Ref. [39], to the purpose of constraining the gluon PDF. Assuming that the scale variations at different energies are fully correlated (i.e. the same values of the renormalisation and factorisation scales can be used in the numerator and the denominator when constructing the ratio), the ratios of differential distributions at different energies exhibit a smaller sensitivity to scale variations. The validity of this assumption can be tested by comparing the ratios at different perturbative orders, and checking the reliability of the (reduced) uncertainty bands obtained assuming such correlation. In Fig. 9 (left) we show the ratio of transverse-momentum distributions. As the centreof-mass energy increases, the p T,bav distribution becomes (slightly) harder, and, therefore, the ratio increases as p T,bav increases. Considering perturbative uncertainties, we notice a strong reduction of the uncertainty bands: while the width of the NNLO band in the original distributions ranges from about ±30% at low p T,bav to about ±10% in the tail, in the ratio the scale uncertainties are reduced to about ±5%. In Fig. 9 (right) we show the corresponding ratio of rapidity distributions. We see that the ratio increases as the rapidity increases, consistently with the fact that at 13 TeV the bottom quarks have a stronger tendency to be produced in the forward direction. Also in this case the LO, NLO and NNLO bands overlap and their size decreases as the order increases, leading to O(±5%) uncertainties at NNLO. Comparing the results at different perturbative orders for both the transverse-momentum and rapidity distributions, we observe an overlap between the LO, NLO and NNLO uncertainty bands over the full kinematical range. This remarkable stability of the perturbative expansion tends to confirm the approach of computing the ratio by using correlated scale variations at different collider energies.
We close this section with an investigation of bottom-quark production in the forward region. In Ref. [18], the LHCb Collaboration has presented results for the measurement of the b-hadron production cross section at the LHC. This measurement used semileptonic decays of b-hadrons into a charmed hadron associated with a muon. The data were collected at √ s = 7 TeV and 13 TeV, in the pseudorapidity interval 2 < η < 5, and inclusively in the b-hadron transverse momentum. Measurements of the ratio between the 7 TeV and 13 TeV rates were provided as well.
Having in mind the limitations of a comparison between theoretical predictions for bottomquark production and b-hadron production data (as already mentioned in Sec. 3.2), in Fig. 10 we present our perturbative predictions for the pseudorapidity distribution of the bottom quark and compare them with the experimental data of Ref. [18].
We first comment on the theoretical results. As previously noticed, the higher order corrections to the pseudorapidity distributions present similar features to those observed for the rapidity spectrum. The corrections are almost independent of the value of η bav and quantitatively very similar to those affecting the total cross sections at both collider energies. At NNLO the scale uncertainty is smaller than at NLO, and the corresponding uncertainy bands largely overlap. The NNLO K-factor is reduced with respect to its NLO equivalent. The effects due to the hadronization of the bottom quarks into b-hadrons are expected to be relatively small (see Appendix B), at the few-percent level in the current range of pseudorapidity and, therefore, well below the perturbative uncertainties at NNLO.
We now comment on the comparison with data. The measured distributions at both collider energies, shown in Fig. 10, are consistent with our NNLO results in the entire pseudorapidity range. While there is an overlap between NLO prediction and data in almost all cases, the inclusion of the NNLO corrections generally improves the agreement with the central predictions and always strongly reduces the scale uncertainties, making the comparison with data more significant. We note, however, that the NNLO scale uncertanties are still considerably larger than the experimental errors. At both collider energies and independently on the perturbative order, we observe that the shapes of the predicted distributions are different from the measured spectra, 4 while all individual points are compatible with the data: at low pseudorapidity values the measurement is below the predicted central value, whereas in the tail the data are closer to the upper bound of the scale uncertainty band. We also observe that the agreement between the NNLO predictions and the data is slightly worse in the tail of the distribution at 13 TeV, although the theory prediction is compatible with the data in each bin.
In Fig. 11 we show the ratio of the 13 TeV and 7 TeV predictions for the pseudorapidity distribution, computed as for the results in Fig. 9. As expected, the ratio features a strong reduction of the scale uncertainties at NNLO, from about ±30% to only about ±5%. Moreover, the LO, NLO and NNLO uncertainty bands overlap, thereby supporting again the use of correlated variations of the renormalisation and factorisation scales. Considering the data-theory comparison, we observe that the reduced NNLO scale uncertainties of the ratio predictions are significantly smaller than the experimental errors.
Comparing the NNLO prediction with the experimental data, we observe that, with the exception of the first bin, the data are systematically above the central NNLO result. In two bins the data points lie outside of the NLO and NNLO uncertainty bands. We note, however, that the η-independent systematic errors are the dominant source of uncertainty in the measurement of the ratio: while the total (statistical plus systematic) uncertainty ranges  Table 4: Cross section for single-inclusive b-quark production at √ s = 7 and 13 TeV in the pseudorapidity region 2 < η bav < 5. The NNLO results are compared with the LHCb measurement of Ref. [18]. Numerical errors on the last digits are stated in brackets, as in Table 1. The LHCb data are presented with their statistical and systematic uncertainties.
between ±8.2% and ±10.5% for the different bins, the η-independent uncertainty is ±7.4% (see Tables 3 and 4 of Ref. [18]). This η-independent uncertainty is of the same order as the observed data-theory discrepancy. In the data-theory comparison of Fig. 10 we had observed similar shape differences at 7 TeV and 13 TeV, which largely cancel out in the ratio.
The pseudorapidity distribution can be integrated over the range 2 < η bav < 5 to obtain the accepted cross section. In Table 4 the corresponding NNLO predictions at √ s = 7 and 13 TeV as well as their ratio are compared with the LHCb data. The NNLO result is stated for two values of central scales, µ 0 = m b and µ 0 = 2m b , and we find the same ratio of the 13 and 7 TeV predictions for these two scale choices. The result for µ 0 = m b suggests larger (and thus more conservative) perturbative uncertainties. Comparing the NNLO predictions for the integrated cross section to the LHCb measurement, we find that both scale choices lead to predictions that are consistent with the data, the choice µ 0 = m b leads to a better agreement though. The NNLO prediction for the ratio is in good agreement with the experimental measurement.
Considering the bb total cross section at NNLO, in Table 3 we have reported its PDF uncertainty (∆ PDFs ) and scale variation uncertainty (∆ scale ), showing that ∆ PDFs is definitely smaller than ∆ scale . In this paper we do not present a study of PDF uncertainties on differential cross sections. We note that, in the very-forward region, PDF uncertainties on (pseudo)rapidity differential cross sections (and on their ratio at different energies) can be larger than the corresponding uncertainties on the total cross section, and their size can be similar to the size of the scale variation uncertainties that we find at NNLO. We remark on this fact by quoting, for instance, some values of ∆ PDFs from the NLO study of Ref. [39]. The PDF uncertainties tend to slightly increase as √ s increases (see accompanying comments to Table 3). In the case of dσ/dη at √ s = 13 TeV, the NLO value of ∆ PDFs is about ±9% at η ∼ 4 and increases to about ±18% at η ∼ 5 [39] (we recall that we find ∆ scale ∼ ±30% at NNLO, see Fig. 10). In the case of the ratio of dσ/dη at 13 TeV and 7 TeV, the NLO value of ∆ PDFs is of O(1%) at η ∼ < 3 and increases to about ±5% (±8%) at η ∼ 4 (η ∼ 5) [39] (at NNLO we find that ∆ scale varies in the range between ±5% and ±3%, see Fig. 11).

Summary
In this paper we have presented the first fully differential NNLO calculation of bottom-quark pair production at hadron colliders. The calculation was carried out by using the q T subtraction formalism to handle and cancel IR divergences. It extends the corresponding calculation for top-quark pair production [47,48]. The computation has been implemented in the Matrix framework, which allows us to evaluate single-and multi-differential distributions with arbitrary IR safe selection cuts.
We have presented results for the total cross section at the Tevatron and the LHC and compared them with predictions obtained by using the numerical code Hathor, finding excellent agreement. We have studied different sources of theoretical uncertainty and observed that, despite the inclusion of NNLO corrections, perturbative uncertainties are still sizeable and dominant over other sources of uncertainties.
We have presented predictions for single-differential distributions at the Tevatron ( √ s = 1.96 TeV) and at LHC ( √ s = 7 TeV and √ s = 13 TeV). As a general feature, we observed that the inclusion of NNLO corrections suggests a (slow) convergence of the perturbative series, with a good overlap between NLO and NNLO bands and a significant reduction of perturbative uncertainties, which are estimated as usually through scale variations. Since perturbative uncertainties at NNLO are still large also for differential distributions, we have investigated possible ways to reduce them. We considered both normalised distributions and ratios of distributions at different energies, and we assumed scale uncertainties to be fully correlated in the respective ratios. In both cases we showed that the ensuing results are perturbatively stable and that LO, NLO and NNLO uncertainty bands overlap, suggesting that this approach indeed provides reliable predictions with reduced perturbative uncertainties.
We have also compared our predictions with experimental measurements for b-hadron production from the CDF Collaboration at the Tevatron and the LHCb Collaboration at the LHC, finding reasonably good agreement. Further studies in the high-p T region require the resummation of large logarithmic contributions of the form ln p T /m b , while more detailed data-theory comparisons could benefit from the inclusion of fragmentation effects. Such studies are left for future work.
have a quantitative dependence on the actual value of p T (and m), but such dependence is 'monotonic': the differences are maximal if p T m, and they tend to vanish if p T m (since η y in this region). Therefore, the main features that we have just discussed (shifted maximun and central dip of the pseudorapidity cross section) remain valid at the single-differential level after integration over p T .
Such features are clearly visible for bottom-quark production by comparing the singledifferential cross sections dσ/dy and dσ/dη in Figs. 5 and 6, respectively. The pseudorapidity distribution has its maximum value at η ∼ 1.5-2.0, and at central pseudorapidity we have with the value r 0 0.7 (at both NLO and NNLO in the bin where 0 ≤ η ≤ 0.5) that is significantly smaller than unity. In the case of bottom-quark production the p T differential cross sections are highly dominated by their values in the region where p T,b/b is of the order of m b (see Figs. [2][3][4]. Therefore, relations between single-differential distributions can be roughly obtained from Eqs. (7) and (9) by performing the replacements dσ/(dp T dx) → dσ/dx (with x = η, y) and J → J , where J is the 'average' (typical) value of J that is obtained by the p T,b/b integration over the region with p T,b/b ∼ m b . In particular, in Eq. (9) this replacement leads to Eq. (10) with r 0 = p T,b/b /m T,b/b 0 , where p T /m T 0 is the average value of p T /m T with respect to the p T integration of dσ/(dp T dy) at y = 0. Such expression is consistent with the numerical value r 0 0.7 (we recall that the average transverse momenta of the b quark are 5.5 GeV and 5.9 GeV at the two LHC energies that we are considering) that we find in Eq. (10) from the results in Figs. 5 and 6.
We remark on the fact that the relation in Eq. (9) between double-differential distributions has an entirely kinematical origin. However, the relation between the single-differential cross sections dσ/dy and dσ/dη (e.g., Eq. (10)) also includes some dynamical effects, since it is controlled by the typical value of p T /m, which does depend on the production dynamics of the observed particle.
To qualitatively highlight these dynamical effects we can compare, e.g., bottom-quark and pion (π) production. In both cases, dσ/dη has a maximum at a finite value of |η| and a centralpseudorapidity dip. However, the shape differences between dσ/dy and dσ/dη are much less pronounced in the case of pion production since p T,π /m π is sizeably larger than p T,b/b /m b . At LHC energies we have p T,π ∼ 400 MeV ∼ 3m π and (consequently) r 0 ∼ 0.95 in Eq. (10). obtained by resumming them up to next-to-leading logarithmic accuracy, and the resummed result is combined and consistently matched to the NLO result to avoid double-counting of perturbative contributions up to O(α 3 S ). The resummation of the large logarithmic terms is carried out by folding the partonic cross section for the inclusive production of a high-p T massless parton i (the heavy quark is also considered to be effectively massless at high p T ) with the perturbative function D i→Q (z, µ F ) [25] that describes the fragmentation of the massless parton at scale µ F into the massive heavy quark Q at a scale of the order of m Q . The perturbative fragmentation function depends on the factorization scale µ F and on the momentum fraction z that is transferred to the produced heavy quark Q.
In the following, we present a comparison of our NLO and NNLO results with the FONLL result, which is obtained from the code in Ref. [68]. For this comparison we limit ourselves to considering b-quark production at the LHC with √ s = 7 TeV. The FONLL results are obtained by using the NNPDF30 NLO PDFs [65] with n f = 5 massless-quark flavours and α S (m Z ) = 0.118. The value of m b is fixed to m b = 4.75 GeV, and the central scale is µ 0 = m T for both the renormalisation and factorisation scales. Our NLO and NNLO results are obtained with the same parameters by using the NNPDF30 NLO and NNLO PDFs with n f = 4 masslessquark flavours.
We find that the p T integral of dσ/dp T increases by about 7% in going from NLO to FONLL accuracy. We note that this increase of the total cross section 5 is not a significant prediction of the FONLL calculation (since the total cross section does not have collinear logarithmic contributions to be resummed), and it has to be regarded as a (higher-order) systematic effect due to the arbitrariness of the procedure that is used to combine the NLO and resummed calculations in the low-p T region. We also note that such systematic effect is 'acceptable', since its size is definitely smaller than the size of the NLO perturbative uncertainty from scale variations (see Table 2).
The comparison of the p T distributions at different perturbative orders is shown in Fig. 12. We start our discussion by comparing the FONLL and NLO results. The high-p T resummation has a non-trivial effect on the shape of the perturbative NLO result. At small transverse momenta (p T ∼ < 10 GeV) the resummation effect is at the few-percent level (this is the p T region that mostly contributes to the value of the total cross section), and it increases to about 25% at p T ∼ 20-30 GeV, decreasing again as p T increases, and it is about 10% at p T ∼ 50 GeV. The region of higher values of p T (which is not shown in the figure) is the most relevant for resummation: here the effect of resummation becomes negative [27] and the FONLL prediction eventually undershoots the NLO result. This behaviour is due to the fact that multiple radiation of collinear partons tends to make the p T spectrum softer at high values of p T . in the low-p T region (p T ∼ < 20 GeV) the FONLL and NLO uncertainty bands are of comparable size, and the effect of resummation is smaller than the scale uncertainties. In this region, the FONLL and NLO predictions are fully consistent, and their differences are due to higher-order effects that are not predicted by resummation. At higher values of p T (30 GeV ∼ < p T ∼ < 50 GeV) the main effect of the resummation is a reduction of the scale uncertainties. We now comment on the comparison between the FONLL and NNLO results. At small transverse momenta the NNLO K-factor increases from about 1.2 to about 1.35 at p T ∼ 15 GeV and then slightly decreases to about 1.3 at p T ∼ 50 GeV. Over the entire region of transverse momenta shown in Fig. 12, the NNLO central result is systematically higher than the FONLL result and has a smaller uncertainty, especially in the low-and intermediatep T regions. The FONLL and NNLO uncertainty bands overlap. Owing to these features, we conclude that the NNLO result is a reliable prediction in the entire region of transverse momenta considered in Fig. 12. For higher values of p T , the impact of resummation will start to be important to reduce the perturbative uncertainties. Eventually, for very high tranverse momenta, a resummed calculation will be needed to obtain reliable predictions.
Comparing the FONLL and NLO results for the rapidity and pseudorapidity distributions (Fig. 13), we observe that the FONLL result is systematically higher than the NLO result. Such enhancement is consistent with the increase of the corresponding total cross sections. For the rapidity distribution the FONLL/NLO ratio is rather flat and about 1.07, while for the pseudorapidity distribution, the FONLL/NLO ratio is about 1.08 in the central region and decreases as |η av | increases. In particular, in the region relevant for LHCb data (see Figs. 10 and 11), the FONLL result is about 6% (3%) higher than the NLO result at η av = 2 (η av = 5). Comparing the scale uncertainty bands at NLO and FONLL accuracies, we see that they have a very similar size, which is definitely larger than the difference between the NLO and FONLL results at central scales. The relatively small and uniform impact of the resummation on the (pseudo)rapidity distribution is not unexpected. The perturbative resummation implemented in the FONLL calculation deals with ln(m T /m b ) terms, and thus it should not significantly affect the shape of such distributions. The FONLL relative effect of O(10%) with respect to NLO is due to the lower-p T region in Fig. 12 where resummation cannot improve the predictivity of fixed-order calculations. In Fig. 13 the NNLO K-factor ranges between 1.2 and 1.3, and the NNLO effects are larger than the FONLL effects by O(20%). The scale uncertainty bands of the FONLL and NNLO results do overlap, and the NNLO result has smaller uncertainties over the entire (pseudo)rapidity region.
Within the FONLL framework, cross sections for the inclusive production of a single B hadron can be computed by supplementing the corresponding b-quark cross section with nonperturbative effects that describe the fragmentation of the b quark into the B hadron. These effects are implemented in the double-differential cross section with respect to the transverse momentum p T and rapidity y of the produced particle (either b quark or B hadron) by performing a convolution of the b-quark cross section with the non-perturbative fragmentation function of the b quark. The convolution acts on the fraction of the transverse momentum that is transferred from the b quark to the B hadron, while the rapidity is kept fixed (i.e., the rapidities of the b quark and B hadron are set to be equal). The multiplicity of B hadrons from the non-perturbative fragmentation function is set to be equal to unity. Using this procedure, the non-perturbative fragmentation produces no effects on the single-differential cross section dσ/dy and, consequently, on the total cross section. Non-perturbative fragmentation effects are instead introduced in other differential distributions.
Using the FONLL code of Ref.
[68], we have computed the effects of non-perturbative fragmentation on the single-differential cross sections with respect to the transverse momentum and to the pseudorapidity of the produced particle. The calculation [68] uses the same nonperturbative fragmentation function as used in Ref. [30], which is extracted [33] from B-hadron production data in high-energy e + e − collisions at LEP (since the extraction is based on the resummed FONLL calculation, the direct use of this non-perturbative fragmentation function in the context of fixed-order calculations, at either NLO or NNLO, is questionable). In the following we comment on the results of the computation.
In the case of the FONLL calculation of dσ/dp T , the non-perturbative fragmentation effects soften the p T spectrum. They decrease the value of dσ/dp T in the high-p T region and, consequently (since the total cross section is unchanged), they increase the value of dσ/dp T in the low-p T region. The crossover point is at p T ∼ 4.5 GeV, where dσ/dp T is almost unchanged. Specifically, we find that the p T differential cross section decreases by about 25% at p T ∼ 50 GeV, about 20% at p T ∼ 20 GeV and about 13% at p T ∼ 10 GeV. In the low-p T region, the p T differential cross section increases by about 10% at p T ∼ 2 GeV.
The non-perturbative fragmentation effects slightly modify the shape of the pseudorapidity cross section (see Fig. 6). Including the non-perturbative fragmentation function, we find that the value of dσ/dη decreases by about 3% at η = 0, is almost unchanged in the region close to the peak at η ∼ 1.8 (the position of the peak is slightly shifted forward by ∆η ∼ 0.1), and increases by about 4% at η = 5. Such effects on the shape of dσ/dη are qualitatively and quantitatively consistent with our expectations. Indeed, as we have discussed in Appendix A (see Eqs. (7)-(10) and accompanying comments), the shape of dσ/dη 'kinematically' follows from the shape of dσ/(dp T dy) and is 'dynamically' controlled by the size of the typical value of p T /m T . The inclusion of the non-perturbative fragmentation function softens the p T spectrum and it slightly reduces (by few percent) the average value of p T , therefore producing the (fewpercent level) effects on dσ/dη that we find.
The results of the FONLL calculation can be used to roughly estimate the size of nonperturbative fragmentation effects on fixed-order calculations. We expect effects of O(10%) on p T -dependent distributions (e.g., dσ/dp T ) in the region from low to intermediate values of p T (say, p T ∼ < 20 GeV), and effects at the few-percent level on p T -inclusive distributions (e.g., dσ/dy and dσ/dη). Comparing these non-perturbative fragmentation effects with perturbative scale uncertainties at NNLO, we see that the fragmentation effects are smaller (although of comparable size) for p T -dependent distributions (see Figs. 4 and 12) and definitely smaller for p T -inclusive distributions (see Figs. 5, 6 and 13).