$B$-hadron fragmentation functions at next-to-next-to-leading order from global analysis of $e^+e^-$ annihilation data

We present nonperturbative fragmentation functions (FFs) for bottom-flavored ($B$) hadrons both at next-to-leading (NLO) and, for the first time, at next-to-next-to-leading order (NNLO) in the $\overline{\mathrm{MS}}$ factorization scheme with five massless quark flavors. They are determined by fitting all available experimental data of inclusive single $B$-hadron production in $e^+e^-$ annihilation, from the ALEPH, DELPHI, and OPAL Collaborations at CERN LEP1 and the SLD Collaboration at SLAC SLC. The uncertainties in these FFs as well as in the corresponding observables are estimated using the Hessian approach. We perform comparisons with available NLO sets of $B$-hadron FFs. We apply our new FFs to generate theoretical predictions for the energy distribution of $B$ hadrons produced through the decay of unpolarized or polarized top quarks, to be measured at the CERN LHC.


I. INTRODUCTION
For a long time, there has been considerable interest in the study of bottom-flavoredhadron (B-hadron) production at hadron and e + e − colliders, both experimentally and theoretically. Historically, the first measurements were performed more than three decades ago by the UA1 Collaboration at the CERN SppS collider [1] operating at center-of-mass (CM) energy √ s = 630 GeV.
In the framework of the parton model of QCD, the description of the inclusive single production of identified hadrons h involves fragmentation functions (FFs), D h a (x, Q 2 ). At leading order, their values correspond to the probability that the colored parton a, which is produced at short distance, of order 1/ Q 2 , fragments into the colorless hadron h carrying the fraction x of the energy of a. Given their x dependence at some scale Q 0 , the evolution of the FFs with Q 2 may be computed perturbatively from the timelike Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [2][3][4]. The timelike splitting functions P T a→b (x, α s (Q 2 )) appearing therein are known through NNLO [5][6][7]. In the case of e + e − annihilation, the hard-scattering cross sections for the inclusive production of parton a, to be convoluted with D h a (x, Q 2 ), are also known through NNLO [5,[8][9][10][11]. This allows one to interpret e + e − data of the inclusive single production of hadron h at NNLO and thus to extract FFs at this order [12][13][14][15]. Owing to the factorization theorem, the FFs are independent of the process by which parton a is produced. This allows one to transfer experimental information from e + e − annihilation to any other production mechanism, such as photoproduction, leptoproduction, hadroproduction, and two-photon scattering.
Of all these processes, e + e − annihilation provides the cleanest laboratory for the extraction of FFs, being devoid of nonperturbative effects beyond fragmentation itself. Presently, there is particular interest in hadroproduction at the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC) due to ongoing experiments.
The parton model of QCD implemented in the modified minimal-subtraction (MS) factorization scheme with n f massless-quark flavors, the so-called zero-mass variableflavor-number scheme (ZM-VFNS), can also be applied to the open production of heavy flavors, such as D and B hadrons, provided the hard energy scale characteristic for the production process is sufficiently larger than the heavy-flavor mass. This is certainly the case for all the applications here, because M B M Z , m t . Recently, D hadron FFs have been provided at NNLO in Ref. [16]. Here, we perform the first NNLO determination of B-hadron FFs.
In Refs. [17,18], B-hadron FFs were extracted at NLO in the ZM-VFNS by fitting to the fractional-energy distributions dσ/dx B of the cross section of e + e − → B + X measured by the ALEPH [19] and OPAL [20] Collaborations at the CERN Large Electron Positron Collider (LEP1) and the SLD Collaboration [21] at the SLAC Linear Collider (SLC).
In the meantime, also the DELPHI Collaboration have reported a similar measurement at LEP1 [22]. In the present work, these data are, for the first time, included in a B-hadron FF fit. On the other hand, we are not aware of any other such data from e + e − annihilation.
In want of NNLO hard-scattering cross sections for inclusive single B-hadron production from other initial states, we do concentrate here on e + e − annihilation. We also go beyond Refs. [17,18] by performing a full-fledged error estimation, both for the FFs and the resulting differential cross sections, using the Hessian approach [23].
The LEP1 experiments [19,20,22]  Section III, we explain the minimization method in our analysis and our approach to error estimation. In Section. IV, our NLO and NNLO results are presented and compared with the experimental data fitted to. In Section. V, we present our NLO predictions for the normalized-energy distributions of B hadrons from decays of (un)polarized top quarks.
Our conclusions are given in Section VI.

II. QCD FRAMEWORK FOR B-HADRON FFS
As mentioned in Sec. I, we fit nonperturbative B-hadron FFs to measured x B distributions of the cross section of where X refers to the unobserved part of the final state. In the following, we explain how to evaluate the cross section of process (1) at NLO and NNLO in the ZM-VFNS. Denoting the four-momenta of the virtual gauge boson and the B hadron by q and p B , respectively, we have s = q 2 , p 2 B = m 2 B , and x B = 2(p B · q)/q 2 . In the CM frame, s is the energy of the B hadron in units of the beam energy. In the ZM-VFNS, we have where i = g, u,ū, . . . , b,b runs over the active partons with four-momenta p i , dσ i (x i , µ R , µ F )/dx i is the partonic cross section of e + e − → i + X differential in x i = 2(p i · q)/q 2 , D B i (z, µ F ) is the i → B FF, and µ R and µ F are the renormalization and factorization scales, respectively.
The latter are a priori arbitrary, but a typical choice is µ F = µ R = √ s. In the CM frame, z = x B /x i is the fraction of energy passed on from parton i to the B hadron. It is customary in experimental analyses to normalize Eq. (2) by the total hadronic cross section, where α and α s are the fine-structure and strong-coupling constants, respectively,ẽ i is the effective electroweak charge of quark i, and the coefficient K QCD ≈ 1.411 [24]. The z distribution of the b → B FF at the starting scale µ 0 is a genuinely nonperturbative quantity to be extracted from experimental data. Its form is unknown, and an educated guess is in order. The selection criterion is to score a minimum χ 2 value as small as possible with a set of fit parameters as minimal as possible. As in Refs. [17,18], we adopt here the simple power ansatz [25], with fit parameters N b , α b , and β b , and choose µ 0 = m b = 4.5 GeV. This ansatz was found to enable excellent fits [17,18]. The i → B FFs for the other quarks and the gluon are assumed to be zero at µ F = µ 0 and are generated through the DGLAP evolution to larger values of µ F . We take α s (M Z ) to be an input parameter and adopt the world average value 0.1181 for n f = 5 [26] both at NLO and NNLO.
As mentioned in Sec. I, we fit to ALEPH [19], DELPHI [22], OPAL [20], and SLD [21] data. These data sets reach down to very small x B values, which fall outside the range of validity of our fixed-order approach. In fact, in the small-x B limit, both the timelike splitting functions and the hard-scattering cross sections develop soft-gluon logarithms that require resummation. At the same time, finite-m b and finite-M B effects become relevant there.
We leave the implementation of these refinements for future work, and instead impose appropriate minimum-x B cuts for the time being. Specifically, we only include ALEPH data points with x B ≥ 0.25, DELPHI data points with x B ≥ 0.36, OPAL data points with x B ≥ 0.325, and SLD data points with x B ≥ 0.28. This enables acceptable fits within the fit range, at the expense of certain deviations in the small-x B range, of course. The fixed-order approach is also challenged in the large-x B limit, by the emergence of threshold logarithms, which also require resummation. In practice, however, these effects do not jeopardize the quality of our fits, so that we refrain from imposing maximum-x B cuts.

III. DETERMINATION OF B-HADRON FFS AND THEIR UNCERTAINTIES
We now explain our fitting procedure. For a given set ters, the goodness of the overall description of the experimental data by the theoretical predictions is measured by the global χ 2 value, where n labels the N exp = 4 experimental data sets, w n are their weight factors [27,28], which we take to be unity, and is the χ 2 value of data set n. On the experimental side, F exp n,i is the central value of (1/σ tot )dσ/dx B measured in bin i out of the N data n bins in data set n, ∆F exp n,i is its individual error obtained by combining statistical and systematic errors in quadrature, N n is the unknown overall normalization factor of data set n to be fitted, and ∆N n is its error as quoted by the experimental collaboration. On the theoretical side, F theo n,i (p) is the respective NLO or NNLO prediction.
We determine the fit parameters p by minimizing Eq. (5) with the help of the Monte Carlo package MINUIT [29] from the CERN program library. We adopt a two-step procedure. In the pre-fitting stage, we determine the four values N n by fitting them simultaneously with the three fit parameters p. In the main fitting stage, we then refine the determination of p with large statistics keeping N n fixed. In the evaluation of χ 2 global (p)/d.o.f., we take the number of degrees of freedom to be the overall number of data points fitted to minus three for the proper fit parameters p. We find the APFEL library [30] to be a very useful FF fitting tool.
We now describe our methodology for the estimation of the uncertainties in the Bhadron FFs. We adopt the Hessian approach to the propagation of uncertainties from the experimental data sets to the FFs, which has proven of value in global analyses of parton distribution functions and has been frequently applied there. For definiteness, we ignore additional sources of uncertainties, which are mostly of theoretical origin and are negligible against the experimental uncertainties taken into account here. In the following, we briefly review our procedure. For more details, we refer to Ref. [31].
In the Hessian approach, the uncertainty bands on the B-hadron FFs, ∆D B b (z), may be obtained through linear error propagation, where p i (i = 1, 2, 3) are the free parameters in Eq. (4),p i are their optimized values, and H −1 (p) is the covariance matrix, which is a default output of the MINUIT program [29]. In Eq. (7), we have suppressed the label µ F for the factorization scale, which we take to be µ 0 .
The error bands ∆D B (z) are subject to DGLAP evolution in µ F along with the central values D B (z). The confidence level (C.L.) is controlled by T 2 = ∆χ 2 global . We adopt the standard parameter fitting criterion by choosing T = 1, which corresponds to the 68% C.L., i.e. the 1σ error band. In Sec. IV, the uncertainty bands thus determined are presented both for the B-hadron FFs and for the physical observables evaluated with them.

IV. RESULTS AND DISCUSSION
We are now in a position to present our results for the B-hadron FFs both at NLO and NNLO and to compare the resulting theoretical predictions with the experimental data fitted to, so as to check directly the consistency and goodness of our fits. We also compare our B-hadron FFs with the NLO ones presented by Kniehl, Kramer, Schienbein and Spiesberger (KKSS) [18].
In Table I, for each of the four experimental data sets n, from ALEPH [19], DELPHI [22], OPAL [20], and SLD [21], the number N data  Table. II.
In Fig. 1, the z distributions of the NLO and NNLO b → B FFs at the initial scale µ F = µ 0 are compared with each other. The NLO and NNLO results agree in shape and position of maximum, but differ in normalization. This difference is induced by the O(α 2 s ) correction terms in the hard-scattering cross sections and in the timelike splitting functions, and it is compensated in the physical cross sections to be compared with the experimental data up to terms beyond O(α 2 s ). The error bands determined as described in Sec. III are also shown in Fig. 1. They are dominated by the experimental errors, which explains why they are not reduced by passing from NLO to NNLO. In Fig. 1, the KKSS b → B FF [18] is included for comparison. It somewhat undershoots our NLO b → B FF, which we attribute to the impact of the DELPHI data [22], which were not available at the time of the analysis in Ref. [18].
In Fig. 2(a), the analysis of Fig. 1 is repeated for µ F = M Z , the CM energy of the experimental data fitted to. Our NLO and NNLO b → B FFs are now closer together, the remaining difference being entirely due to the O(α 2 s ) corrections to the hard-scattering cross sections. On the other hand, the difference between the NLO b → B FF and the KKSS one is hardly affected by the DGLAP evolution from µ 0 to M Z , as it is due to a difference in the collection of experimental data fitted to. Figure 2 Table I. The failure of the theoretical descriptions in the small-x B regime is, of course, a direct consequence of the small-x B cuts applied.
For better visibility, we present the information contained in Fig. 3 as data over theory plots in Fig. 4, one for each experiment. Specifically, the experimental data are in turn normalized to the NLO and NNLO central values. As already explained above, the NLO and NNLO uncertainty bands are very similar. As already visible in Fig. 3  as its degree of polarization in a given production mode, which includes single and pair production. We thus consider both unpolarized and polarized top quarks.
We work in the rest frame of the top quark. The partial width of the decay t → BW + +X, differential in the scaled B-hadron energy x B and the angle θ P enclosed between the topquark polarization three-vector P and the B-hadron three-momentum p B is given by  [19], DELPHI [22], OPAL [20], and SLD [21]. The uncertainty bands (green and red hatched areas) stem from those of the B-hadron FFs.
where P = | P | is the degree of polarization. In the ZM-VFNS, we have where dΓ unpol i /dx i and dΓ pol i /dx i refer to the parton-level decay t → iW + + X, differential in the scaled energy x i of parton i = b, g. In the top-quark rest frame, we have  [32][33][34] and Refs. [35][36][37][38], respectively. In Ref. [38], θ P is taken to be enclosed between P and the W-boson three-momentum p W .
Although a consistent analysis is presently limited to NLO, we also employ our NNLO B-hadron FF set to explore the possible size of the NNLO corrections.
In our numerical analysis, we use m b = 4.5 GeV, m W = 80.379 GeV, and m t = 173.0 GeV [26], and choose µ R = µ F = m t . In Fig. 5(a), we present the NLO predictions of dΓ unpol /dx B and dΓ pol /dx B , evaluated with our NLO B-hadron FF set. For comparison, the evaluations with our NNLO B-hadron FF set are also included. We observe from Fig. 5(a) that switching from the NLO B-hadron FF set to the NNLO one slightly smoothens the theoretical prediction, decreasing it in the peak region and increasing it in the tail region thereunder.
At the same time, the peak position is shifted towards smaller values of x B . The change in normalization is of order 5% at most. These effects should mark an upper limit of the total NNLO corrections because the as-yet-unknown NNLO corrections to dΓ unpol i /dx i and dΓ pol i /dx i are expected to give rise to some compensation if FF universality is realized in nature. In Fig. 5(b), the results for dΓ pol /dx B in Fig. 5(a) are compared to the evaluation with the KKSS B-hadron FF set [18]. As in Figs. 1, 2(a), and 2(b), the NLO result is somewhat reduced by switching to the KKSS B-hadron FF set.

VI. SUMMARY AND CONCLUSIONS
In this paper, we determined nonperturbative FFs for B hadrons, both at NLO and NNLO in the ZM-FVNS, by fitting to all available experimental data of inclusive single B-hadron production in e + e − annihilation, e + e − → B + X, from ALEPH [19], DELPHI [22], OPAL [20], and SLD [21]. We then applied these B-hadron FFs to provide NLO predictions for inclusive B-hadron production by top-quark decay, t → BW + + X, both for unpolarized and polarized top quarks.
Our analysis updates and improves similar ones in the literature [17,18] in the following respects. We included the DELPHI data [22], which had not been available then.
For the first time, we advanced to NNLO in a fit of B-hadron FFs. We performed a careful estimation of the experimental uncertainties in our B-hadron FFs using the Hessian approach.
We adopted the simple power ansatz of Eq. (4) and obtained for the three fit parameters appearing therein the values listed in Table II. The goodness of the NLO and NNLO fits turned out to be excellent, with χ 2 /d.o.f. values of 1.485 and 1.104, respectively (see Table I).
As expected on general grounds, the fit quality is improved by ascending to higher orders of perturbation theory.
We encourage the LHC Collaborations to measure the x B distribution of the partial width of the decay t → BW + + X, for two reasons. On the one hand, this will allow for an independent determination of the B-hadron FFs and thus provide a unique chance to test their universality and DGLAP scaling violations, two important pillars of the QCDimproved parton model of QCD. On the other hand, this will allow for a determination of the top-quark polarization, which should depend on the production mode. t → BW + + X [34], and t(↑) → BW + + X [45], have all been worked out in the GM-VFNS at NLO, but not yet at NNLO. Finite-m B effects may be conveniently incorporated using the approach of Refs. [46][47][48]. The implementation of such theoretical improvements reaches beyond the scope of the present analysis and is left for future work.