Resummation effects in the bottom-quark fragmentation function

We compute the perturbative component of the fragmentation function of the $b$ quark to the best of the present theoretical knowledge. The fixed-order calculation to order $\alpha_s^2$ of the fragmentation function at the initial scale is matched with soft-emission logarithm resummation to next-to-next-to-leading logarithmic accuracy, so that order-$\alpha_s^2$ corrections are accounted for exactly, and logarithmically enhanced contributions from loops of $b$ quarks are included. This requires the calculation of the Mellin transform of the order-$\alpha_s^2$ result in the whole complex plane for the Mellin variable, which we provide for the first time for all the fragmenting partons. Evolution is performed to next-to-next- to-leading log accuracy, and mixing with the gluon fragmentation function is taken into account. The perturbative fragmentation functions are made available via LHAPDF grids.


Introduction
The production of heavy quarks (charm, bottom, top), possibly in association with other particles, is a particularly important class of processes at colliders. Not only they provide key information for advancing our understanding of strong and electroweak interactions in the Standard Model (SM), but their very distinctive experimental signatures typically enter as background in many SM measurements and beyond-the-SM searches. From the theoretical point of view, reliable and accurate predictions for these processes that match the current and foreseen experimental precision, are necessary. This is, however, a tall order. Fixed-order predictions are, in general, not accurate enough and terms appear that hamper the convergence of the perturbative series and need to be included at all orders. A class of such terms involve mass logarithms, i.e. powers of log Q 2 m 2 , where m and Q are the heavy-quark mass and the typical energy scale of the process, respectively. When these terms are large, an all-order resummation, which can be achieved by means of Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution of fragmentation functions (FFs), becomes necessary. At past colliders and at the LHC this is important for the bottom and the charm quark. For top quarks such effects become relevant at transverse momentum above a few TeV, a regime which is marginal at the LHC but possibly of interest at future high-energy colliders. FFs make it possible to organise terms of order α p s log q Q 2 m 2 and to systematically resum them via the DGLAP evolution, thus obtaining physical predictions at a given logarithmic accuracy. Nowadays, the DGLAP evolution is implemented in several computer codes, such as QCDnum [1], ffevol [2], APFEL [3], MELA [4], EKO [5], up to next-to-next-to-leading order (NNLO).
The fragmentation function of a b quark is a special case: not only its dependence on the hard scale Q is under control in perturbation theory via DGLAP equations, but also the initial condition for evolution, typically given at a scale of the order of the b-quark mass, can be computed perturbatively [6][7][8][9][10]. Such perturbative calculation becomes unreliable in the kinematical regime where the produced b quark carries most of the available energy, and therefore the recoiling radiation is soft, giving rise to a further class of large logarithms. The formalism for the resummation of soft emission logarithms is outlined in [8], where resummation is explicitly performed to next-to-leading logarithmic (NLL) accuracy.
We stress that having an excellent perturbative description of FFs represents only part of the solution to the problem of providing an accurate description of heavy-quark related processes. The perturbative result has to be complemented with a non-perturbative part which parametrises the transition from the heavy quark to the corresponding flavoured hadron. This contribution cannot be computed in perturbation theory, and is usually extracted by a direct comparison to data. Once the perturbative part is known, the extraction of the non-perturbative contribution to the FF is a well-defined and independent phenomenological task. Though relevant to provide physical predictions to compare with experiment, it is not addressed in this work. The interested reader can find applications to the case of heavy-quark production in leptonic collisions in [7,8,11,12], of the bottom quark in top decays in [13][14][15], and of the Higgs decay to heavy quarks in [16].
In this work, we improve the perturbative description of the heavy-quark FF, using the framework to compute the coupled evolution of the b-quark fragmentation with the gluon and the other partons at NNLO accuracy described in [17] which builds upon MELA [4]. We asses the impact of resummation of soft logarithms on the initial condition for the evolution of the b fragmentation function D b (x, µ 2 0F ) by implementing the resummation formalism of [8] up to next-to-next-to-leading logarithmic accuracy (NNLL in the following). Resummation at NNLL of the initial condition was already performed in [12]. In that work the resummed result was matched to the exact perturbative calculation at next-to-leading order (NLO). We improve on this result in three ways. First, we compute the process-dependent component of the resummed initial condition in such a way to reproduce the exact order-α 2 s result of [9] to NNLL accuracy, including the contributions originating from loops of heavy quarks. Note that such contributions had not been accounted for in previous works [12,18]. Second, we match the NNLL resummed result to the exact NNLO initial condition provided in [9]. This requires the calculation of the Mellin transform of the full NNLO initial condition, which was not available until now in an analitically-continued version defined in the whole complex plane. The calculations presented here fill this gap. Finally, we include both the heavy-quark and the gluon and light quark FFs in their coupled evolution.
The work is organised as follows. In Sect. 2 we summarise the formalism and the relevant equations for the heavy-quark fragmentation function computation starting from a perturbative initial condition. In particular, we focus on the dependence of the results on the scheme used to deal with heavy flavours, and we present our calculation of the initial bottom fragmentation function with soft logarithms resummed at NNLL and matched with the NNLO expression. In Sect. 3 we present selected numerical results, commenting on the impact of resummation as well as of the inclusion of the NNLO initial conditions. We present our conclusions in Sect. 4. Appendix A contains some details about the calculation of the resummed initial condition, while in Appendix B we illustrate the calculation of the relevant Mellin transforms.

The heavy quark fragmentation function
In the fragmentation-function formalism, the cross section for the production of a heavy quark Q is given by (2.1) where x is the fraction of the available energy carried away by the produced heavy quark, σ j is a hard cross section (possibly including parton distribution functions for the initial state) and D jQ the fragmentation functions of partons j into the heavy quark Q.
The fragmentation functions at the scale µ F are obtained by evolving suitable initial conditions, given at a reference scale µ 0F , through Altarelli-Parisi equations. In the following, we will be interested in the initial condition in the case j = Q = b, which can be computed perturbatively at µ 0F ∼ m. We will denote by D b the b fragmentation function for simplicity. We will also omit the explicit dependence of the fragmentation function on the b-quark mass m.

Fixed-order calculation of the initial condition
The initial condition for the heavy quark fragmentation function can be computed perturbatively, as illustrated in [7,8], at a scale µ 0F of the order of (but not necessarily equal to [19]) the heavy quark mass. The order-α s correction is given in [7], while the order-α 2 s coefficient was computed in [9]. The result is and n = n f − 1 is the number of massless flavors. The functions F are given explicitly in [9]. The calculation in [9] is performed in the MS renormalization scheme for ultraviolet divergences, which means that both massless and massive flavors take part in the evolution of the running coupling.
For the purpose of comparison with the resummed result, it will be useful to rewrite the result of [9] as an expansion in powers of α s = α (n ) s . Furthermore, for greater generality it will be convenient to choose a renormalization scale µ 0 different from µ 0F . To order α 2 s , we have The first term in curly brackets originates from the change of renormalization scheme, while the second one arises from evolution from µ 0 to µ 0F . Replacing Eq. (2.5) in Eq. (2.2) and expanding again in powers of α s (µ 2 0 ) to order 2 we get The calculation of the Mellin transform of this expression, as a function of the complex variable N can be performed as illustrated in Appendix B; this is one of the main results of this paper.

Soft gluon resummation
It was observed in [8] that logarithmic corrections to the initial condition for the heavy quark fragmentation function, arising at all orders because of soft emission, may play a relevant role, at least conceptually, in the large-z region. It is therefore useful to resum such contributions to all orders, to some given logarithmic accuracy, and to assess their impact on the fragmentation function itself and on related observables. Soft gluon resummation of the b fragmentation function at the initial scale was performed in [8] to NLL accuracy, and subsequently improved to NNLL in [12]. Here we recompute the process-dependent component of the resummed initial condition, by including a term that was omitted in previous results, and improve on them by matching the resummed expression to the orderα 2 s result of [9]. Soft gluon resummation is performed in Mellin space, as illustrated in Ref. [8]. For a generic function g(z), with 0 ≤ z ≤ 1, we define (2.9) The resummed fragmentation function at the initial scale for evolution µ 2 0F takes the familiar form of an exponential, The functions A(α s ) and H(α s ) have perturbative expansions in powers of α s starting at order 1: (2.13) NNLL accuracy is achieved including A(α s ) up to order α 3 s and H(α s ) up to order α 2 s . A(α s ) is determined by the Altarelli-Parisi splitting functions, while H(α s ) is process dependent, and must be extracted from the fixed-order calculation. The coefficients H i must be obtained by matching the expansion of the resummed expression Eq. (2.10) to the relevant perturbative order with the fixed-order calculation.
In Eq. (2.10) the strong coupling is computed with n = 4 active flavour; this is the natural choice, because the integration over µ 2 in the exponent ranges between m 2 (1−z) 2 m 2 and µ 2 0F ∼ m 2 . Its evolution from µ 0F up to a hard scale µ F m is however performed with n f = n + 1 massless flavours. 1 It is interesting to discuss the dependence of the resummed initial condition, Eq. (2.10), on the initial factorization scale µ 0F . The fragmentation function evolved up to a generic hard scale µ F is given by is the Altarelli-Parisi evolution kernel at large N . We expect D res b (N, µ 2 F ) to be independent of the initial scale µ 0F , which implies is the relevant MS anomalous dimension in the large-N limit. As already observed, the evolution between an initial scale µ 0F , of the order or the heavy quark mass m, and a hard scale µ F m, is determined by n f massless flavours; the anomalous dimension is therefore naturally expressed as an expansion in powers of α (n f ) s , whose coefficients A k also depend on n f for k ≥ 2. On the other hand, neglecting O(N 0 ) terms, the logarithmic derivative of Eq. (2.10) reads m 2 . After performing the two integrations in Eq. (2.11), the exponent F is expressed as a power expansion in α s at fixed λ. To NNLL accuracy The functions g 0 , g 1 , g 2 , g 3 are given explicitly in Appendix A, together with some details on their derivation, in terms of the coefficients A 1 , A 2 , A 3 and H 1 , H 2 . The functions g 1 (λ) and g 2 (λ) were already given in [8], and we find full agreement with those expressions. The function g 3 (λ) was computed in [12], where NNLL resummation was performed, but not explicitly provided. The extraction of H 2 from the order-α 2 s calculation deserves some comments. It can be obtained by expanding Eq. (2.10) to order α 2 s and comparing the result with the fixed-order calculation in the large-N limit, Eq. (2.2), after the replacement Eq. (2.5). We obtain which depends on log µ 2 0F m 2 as expected. We can check explicitly that, with H 2 given in Eq. (A.15), the resummed initial con- where we have taken into account that H 1 is µ 0F -independent, and we have displayed the dependence of A 2 on the number of flavours n , see Eqs. (A.12,A.15). Since Furthermore, from Eq. (A.12), as announced.
Our result, Eq. (2.20), differs from the value of H 2 given in [18] by the last term, proportional to C F T f and µ 0F -dependent. This extra term is different from zero for all choices of µ 0F . We have seen that a dependence of H 2 on the initial scale µ 0F is expected on general grounds, and can be read off the Altarelli-Parisi anomalous dimension in the large-N limit. We have also shown that our result Eq. (2.20) is consistent with expectations. In [18], H 2 was extracted under the assumption that the b quark appears only as an external line, i.e. no loops involving b quarks was considered; here, we do not make this assumption. The value of H 2 obtained in [18] was later employed in [12].

Matching resummed and fixed-order initial condition
Both the fixed-order initial condition for the b fragmentation function Eq. (2.8), accurate to NNLO, and the resummed initial condition Eq. (2.10), accurate to NNLL, are now expressed as functions of α s (µ 2 0 ) = α (n ) s (µ 2 0 ), and of the complex variable N , Mellin conjugate to z. Their combination requires the subtraction of the resummed result expanded to order α 2 s to avoid double counting. Our final result is therefore

Results
In this Section, we present results for the b-quark fragmentation function D b (z, µ 2 F ) and its Mellin transform D b (N, µ 2 F ). Our aim is to assess the impact of the NNLL soft-gluon resummation for the initial condition, as well as the size of O(α 2 s ) terms in the latter. We will always employ NNLL DGLAP evolution, and we will set the following values for the bottom mass and the factorisation scale: For all results presented in this section we will set µ 0F = µ 0 . Specifically, we will consider two values for the initial scale µ 0 , namely µ 0 = m (displayed in the left panels) and µ 0 = 2m (in the right panels). The Mellin-transfomed initial condition D b (N, µ 2 0 ) displays a logarithmic branch cut on the positive N real axis for λ > 1 2 , or originated by the Landau singularity of the running coupling. The resummed initial condition is meaningless for N too close to N L (µ 0 ). We find  In Fig. 1 we show results for D b (N, µ 2 F ). We include soft resummation at different logarithmic accuracies on top of the NNLO initial conditions: the dashed black curve shows the predictions without resummation, i.e. the same result presented in Ref. [17], which is seen to become negative at N ∼ 20. LL, NLL and NNLL resummation are included in the dashed teal, dot-dashed magenta and solid orange curves, respectively. It can be appreciated that the NNLL resummed prediction is the only one which remains positive over the whole considered N range, while all other predictions become negative for sufficiently large values of N . Also, at large N (N > 10), effects of resummation remain quite important at all considered orders. Finally, when the initial scale µ 0 is increased by a factor 2, differences among the predictions turn larger, a fact related to the smaller impact of the DLGAP evolution, which is common to all predictions, in favour of that of the initial conditions.
The impact of the NNLO initial condition can be better appreciated in Fig. 2, where we show the same predictions as in Fig. 1, but obtained including only NLO initial conditions. The different predictions are displayed as symbols, with the same color code as in Fig. 1 for the corresponding logarithmic accuracy, and compared to the NNLO+NNLL one. We observe that i) resummation effects play an even more important role in this case; otherwise stated, the inclusion of NNLO initial condition has a stabilising effect in this respect; ii) with NLO initial conditions, NLL resummation is sufficient to get positive-definite predictions up to N ≈ 50; iii) if NNLL resummation is included, the effect of finite terms in the NNLO initial conditions is small for N < 20, while it amounts to several 10%s for larger values of Turning to predictions in the physical space of momentum fraction z, Fig. 3 provides a global view of the effects, where the same patterns as in Figs. 1 and 2 are employed for the corresponding predictions. The most visible effect of resummation is to reduce the value of the FF at the large-z peak, slightly moving its position z P to the left, and slightly raising the tail for z < z P . Unlike in N space, predictions generally reflect the hierarchy of the resummation orders, both when NNLO and NLO initial conditions are employed (central and lower inset respectively). The hierarchy is violated only for very large values of z (well right of z P ), where the differences at large N are reflected. Considering the effects of including NNLO initial conditions, by comparing the central and lower inset we see that their inclusion reduces the impact of resummation, as already observed in N space. We also notice that, when NNLL resummation is employed, genuine O(α 2 s ) effects are generally at the percent level (although not uniformly in the z range), about 5% for µ 0 = m and z < z P , and twice as much when µ 0 is increased to 2m.

Conclusions
In this work, we have obtained results for the b-quark fragmentation function including NNLL soft-gluon resummation for the perturbative initial conditions matched with the NNLO result.
Following the approach of [8], we have performed the resummation of soft emission logarithms for the Mellin-transformed initial condition. First, we have identified a term in the process-dependent component of the NNLL resummed initial condition that was omitted in previous works. Note that the inclusion of this term makes the resummed result fully consistent with the exact NNLO perturbative calculation and with the DGLAP evolution. Moreover we have computed the Mellin transform of the full NNLO initial condition, and obtained its analytic continuation to the complex plane of the variable N , Mellin-conjugate to the momentum fraction z, which was not available in the existing literature. Finally we have considered the coupled evolution of the quark and gluon FFs.
We find that NNLL resummation improves the predictions obtained by evolving initial conditions at NNLO. Their behaviour is regular across the relevant range of the Mellin variable N , similarly to what happens in the NLO+NLL case. The inclusion of NNLO initial conditions provides a sizeable stabilisation of resummed predictions in z space, and has a moderate, yet appreciable effect on the evolved fragmentation function. Our work makes it possible to include the NNLO+NNLL initial conditions, convolved with a suitable short-distance cross section, in a wide range of phenomenological applications related to heavy-hadron production at colliders and in particular at the LHC. The perturbative fragmentation functions are available via LHAPDF grids and can be requested to the authors. An easy-to-use computer code implementing the NNLO+NNLL initial conditions supplemented by NNLO DGLAP evolution is in preparation.

A Calculation of the resummed initial condition for D b
In this Appendix we describe in detail the calculation of the resummed b fragmentation function at the initial scale µ 0F , Eq. (2.10): to NNLL accuracy. In this expression, α s = α (n ) s is defined in the n = n f − 1 renormalization scheme: The renormalization scale µ 0 is taken to be different in general from the initial factorization scale µ 0F . Both are taken of the order of the heavy quark mass m. We first consider the exponent F . We have F α s (µ 2 0 ), To NNLL accuracy, and (A.14) where The µ 2 integration of the A term is conveniently performed in the variable α = α s (µ 2 ): To NNLL accuracy, The integrals in eq. (A.18) are easily computed, and the result is an analytic function of λ z : The functions I p can be computed by taking derivatives of a generating function: (A.26) In the last step we have used the Stirling approximation for the Γ function at large values of its argument, which is the limit we are interested in. We obtain We now observe that Replacing eq. (A.29) in (A.27) we obtain Note that the sum has been extended to ∞ because all derivatives of order p + 2 and higher vanish. We have therefore The z integration is conveniently performed in the variable λ z : An even simpler expression is obtained by separating off the first term in the sum, which is the only one which requires an integration: To NNLL accuracy, only the terms k = 0, 1 are relevant in the sum. We therefore obtain our final expression (A.36) which can be cast in the form of an expansion in powers of α s (µ 2 0 ) at fixed λ: for simplicity.) Constant terms, i.e. N -independent terms, are not controlled by Sudakov resummation. Note that, because the expansion of F starts at order α s , we have So, N -independent terms in F first appear at NNLL. We remove them by replacing 3 . (A.41) We find in agreement with the results presented in ref. [8]. We also find This is a new result: NNLL resummation was performed in [12], but an explicit expression of g 3 (λ) is not given there.
We now turn to the pre-exponential factor g 0 in eq. (2.10). The function g 0 is defined in such a way that constant terms are correctly taken into account up to order in α 2 s . It can therefore be read off the result of Ref. [9]:

B Analitically-continued Mellin transforms for the O (α 2 s ) initial conditions
In this Appendix, we report on the computation of the analitically-continued Mellin transforms of the O(α 2 s ) (NNLO) initial conditions, for all fragmenting partons. These are necessary for the inversion of the expressions obtained in Mellin space back to z space. In our case, such an inversion is performed using the Talbot-path method [20]. We remind the reader that the NNLO initial conditions are taken from [9,10], where they are reported in z space. Analitically-continued Mellin transforms for most terms can be computed e.g. by employing the expressions tabulated in [21][22][23][24][25]. However, other terms exist for which the analitically-continued Mellin transform cannot be found in literature. For these terms we will provide details in Sects. B.1-B.6. Finally, in Sect.B.7, we discuss the validation of our results.

B.1 Terms including a factor log(1+z) k 1+z
For those terms which contain the factor log(1 + z) k /(1 + z) we apply the expansion shown in Eqs. (52), (56) of [23], either with an in-house implementation or by using the ancont code [22,23].
For example, we have obtained where the a (p) k coefficients are tabulated in Ref. [23]. We point out that Eq. (62) of [23] has some typographical errors. In the convention of that paper (note the exponent N in the Mellin transform), it should read: where S 1 (N ) is the first-order harmonic sum [26].

B.2 Terms including a factor 1 (1+z) p
Since in the results of Refs. [9,10] one finds terms containing f (x)/(1+x) 3 and f (x)/(1+x) 4 , we generalise Eqs. (52), (56) of Ref. [23] to these cases, and quote the relevant coefficients. Specifically, assuming p > 1: where we have expanded Relevant cases appearing in the O(α 2 s ) initial conditions are e.g.: In Tab. 1 we quote the coefficients b (B.14) By using the first equation reads which agrees with Ref. [22]. For S 1,2 (−z), we have instead performed the Mellin transform of log 2 (1 + z) using the expansion shown in Ref. [23]. Finally, concerning the Mellin transform of S 1,2 (z 2 ), we have exploited the relation

B.4 Terms involving functions of |2z − 1|
The functions appear in the initial conditions of Refs. [9,10], specifically in d (2) g . The Mellin transform of A 1 can be easily computed: . (B.20) Despite its rather simple form, a problem arises when the inverse Mellin transform is considered, (B.21) From Eq. (B.21), in particular the second term on the r.h.s., one immediately sees that the integration contour cannot be closed on the left part of the complex plane if z > 1/2, rather it must be closed on the right, only for this specific term. This poses a practical problem when the inversion is performed numerically, as in our case. However, since the constant c must be larger than the real part of the left-most pole of the integrand (in this case c > 0), when the integration contour is closed on the right it does not encircle any pole, hence its contribution is zero. Therefore, a practical solution is to keep the second term on the r.h.s. of Eq. (B.21) only when z < 1/2.
Moving on to A 2 , its Mellin transform can be obtained semi-analytically. First, while both A 1 and A 2 have a discontinuous derivative at z = 1/2, their difference, is smooth for z ∈ (0, 1), while at the endpoints it has logarithmic singularities. Subtracting the singularities (and imposing that the remainder vanishes for z = 1/2), one gets the function A reg (z) = A diff (z) + 1 2 (log(z) + log(1 − z)) + log (2) , which is now smooth z ∈ [0, 1], is symmetric around z = 1/2, and it is normalised such that A reg (1/2) = 0. A reg (z) can now be fitted with the functional form The values of the coefficient c k , obtained with k max = 7, are reported in Tab. 2. In this way the Mellin transform of each term entering the sum can be expressed in terms of the Euler Beta function: To summarise, we rewrite the function A 2 (z) in the form where the Mellin transform for all terms on the r.h.s. can now be easily computed.  We will consider the Mellin transforms of the following expressions: We will focus on the case for F 1 . This function is divergent at z → 0, therefore we proceed like in sec. B.4, exposing the singular behaviour in this region. However, at variance with what we did there, the remainder is bounded but its derivative is still divergent both at z → 0 and z → 1. These terms have also to be subtracted, in order to have a regular function F 1reg to fit. We call F 1 0 sing the singular part of F 1 at z → 0, and F 1 0 d−sing , F 1 1 d−sing the contributions which make the derivative of F 1 − F 1 0 sing divergent in z → 0 and z → 1 respectively 2 . Thus we have For what concerns the functions F 2 and F 3 , the terms in Eqs. B.31-B.33 are simple enough so that the Mellin transform of their product with either log z or log(1 − z) can be computed with standard methods. The Mellin transform of the corresponding regular part is also trivial, since it amounts just to the Mellin transform of log z or log(1 − z) multiplied by a power of z. The combination appears in the gluon initial condition. Exploiting its symmetry around z = 1/2, we can proceed in a similar way as in Sects. B.4 and B.5. The function G is smooth for z ∈ (0, 1), thus we proceed by extracting the singular parts of the function, and of the derivative of the remainder, when z → 0 and z → 1: The singular parts, written so that the remainder satisfies s G reg (1/2) = 0, read G 0 sing (z) = log(z) 6 π 2 + log 2 (z) + L 2 π 2 6 − L 2 + L where k max = 7 and the coefficients e k are given in Tab. 2.

B.7 Validation
We conclude this Appendix by presenting some validation plots for the analytically-continued Mellin transforms of the O(α 2 s ) initial conditions. For all the cases presented in this Appendix, the fitting parameters have been chosen so that the fitted function differs from the exact one by no more than few parts in 10 5 . Our validation is both in N and in z space, where X = b, b, g, q, i.e. the relative difference between the numeric and analytic Mellin transforms. From this panel we appreciate that the difference between the numeric and analytic Mellin transforms never exceeds 10 −4 , with the exception of D b at very large N . The validation in z space, shown in the right plot of Fig. 4, presents a very similar layout, with the only relevant change being that the numeric (analytic) Mellin transforms of the initial conditions are replaced with their exact expression (inverse Mellin transform) in x space. Again, the relative difference does not exceed 10 −4 , except for those points where the initial conditions, hence the denominator in Eq. (B.42), vanish, or in the region z → 1. Here it has to be stressed that we do not plot the contribution coming from Dirac delta's at z = 1 or the endpoint of plus distributions, which live exactly at z = 1, which we deem responsible for such a discrepancy. Indeed, we have observed similar discrepancies also in the case of very simple distributions, e.g. log(1−z) 1−z + .