Heavy Quark Fragmentation in e + e − Collisions to NNLO+NNLL Accuracy in Perturbative QCD

: Fragmentation of heavy quarks into heavy-flavoured hadrons receives both per-turbative and non-perturbative contributions. We consider perturbative QCD corrections to heavy quark production in e + e − collisions to next-to-next-to-leading order accuracy in QCD with next-to-next-to-leading-logarithmic resummation of quasi-collinear and soft emissions. We study multiple matching schemes, and multiple regularisations of the soft resummation, and observe a significant dependence of the perturbative results on these ingredients, suggesting that NNLO+NNLL perturbative accuracy may not lead to real gains unless the interface with non-perturbative physics is properly analysed. We confirm previous evidence that D ∗ + experimental data from CLEO/BELLE and from LEP are not reconcilable with perturbative predictions employing standard DGLAP evolution. We extract non-perturbative contributions from e + e − experimental data for both D and B meson fragmentation. Such contributions can be used to predict heavy-quark fragmentation in other processes, e.g. DIS and proton-proton

Fragmentation of heavy quarks, namely charm and bottom, into heavy-flavoured hadrons has been extensively studied in the past.Just a few years after the charm quark discovery the fragmentation function of heavy quarks was understood to be much harder than that of light quarks and peaked at large hadron energy fractions [1,2], while phenomenological models trying to describe its functional form soon appeared (see e.g.[3][4][5]).
Perturbative calculations of heavy quark fragmentation were also soon considered, and resummation of collinear logarithms up to next-to-leading order accuracy, as well as of soft logarithms to leading logarithmic accuracy, was performed in [6,7].Refs.[8,9] extended soft resummation to next-to-leading logarithmic accuracy, while [10] performed a comprehensive phenomenological analysis at NLO+NLL level.The interplay between a massive quark and the resummation of soft gluons has been investigated e.g. in [11][12][13][14].
A number of papers considered the interplay of perturbative and non-perturbative physics in this process using an array of techniques: heavy quark effective theory [15,16], renormalons [17,18], effective strong coupling [8,19,20], soft and collinear effective theory and boosted heavy quark effective theory [21].
Despite fixed order calculations to next-to-next-to-leading-order accuracy having been available since quite some time [22][23][24][25][26][27], it is only quite recently that they were implemented into codes and tools useful for studying them numerically and for doing heavy quark fragmentation phenomenology.Notably, bottom quark fragmentation at next-tonext-to-leading order (NNLO) plus next-to-next-to-leading logarithmic (NNLL) collinear and soft accuracy has been considered in [21,[28][29][30], but the implications of the higher order perturbative accuracy have probably not been fully explored, and charm quark fragmentation has not been studied.
The purpose of this paper, which considers both charm and bottom fragmentation, and also addresses top fragmentation at very large energies, is to provide a phenomenological and user-friendly implementation of heavy quark fragmentation in e + e − collisions at NNLO+NNLL level.The paper also aims to study more comprehensively the impact of NNLO+NNLL corrections, and their interface with the non-perturbative phase describing the hadronisation of the heavy quark into the heavy flavoured hadron.Finally, by comparing to experimental data the paper extracts a non-perturbative component which can be reused to calculate heavy hadron production in hadronic collisions.
Section 2 sketches the setup and establishes nomenclature and notations.Section 3 presents numerical results, while section 4 introduces comparisons to experimental data and extracts non perturbative components by comparing to them.All details and formulas, with the exception of the NNLO fixed order results of [22,[24][25][26][27], too lengthy to be reproduced here, are given in a series of Appendices.

Perturbative Setup
We denote the cross section for the full, observable-level process for the production of a heavy-flavoured hadron H in e + e − collisions with an energy fraction x = E H /E beam , i.e. a one-particle inclusive fragmentation function (FF), as where Q = 2E beam is the total-centre-of-mass energy of the collision, here implicitly understood to be that of electron-positron beams 1 .This cross section will be modeled by splitting it into a convolution2 of a perturbative and a non-perturbative component, 2) The non-perturbative component, D np Q→H (x, {params}), is here understood to integrate to f Q→H , the fraction of heavy quarks Q hadronising into the given hadron H, whereas the perturbative component dσ Q /dx(x) is the cross section for producing the heavy quark Q with energy fraction x.
To simplify the notation, we can rewrite eq.(2.2) as and this form allows one to read σ H (•), σ Q (•), and D np Q→H (•) either as differential distributions in x (i.e.considering • = x, and σ actually being dσ/dx), as above, or as Mellin moments 3 of the original distributions, with • = N and the convolution operation ⊗ becoming simply a product in moments space.
In this expression we have explicitly shown that the perturbative cross section for production of a heavy-quark Q depends on its mass, and that the non-perturbative component D np Q→H (•, {params}) can depend on a set of parameters that will be made explicit later.Of course, neither σ Q nor D np Q→H are physical observables, and the parameters of the latter will in particular depend on the details of the definition of the former.

Fixed order
In this paper we are concerned with evaluating σ Q up to next-to-next-to-leading order (NNLO) accuracy in QCD, and complementing the fixed-order calculation with up-to nextto-next-to-leading-logarithmic (NNLL) accuracy resummation of quasi-collinear and of soft emissions.
To do so, we first introduce the approximation, pioneered by Mele and Nason [6,7], of factorising the perturbative contribution into a process-specific short-distance cross section σi and a universal decay function 4 with p ≥ 1.In this expression i denotes a partonic channel (quarks and gluon, or singlet, gluon and non-singlet, depending on the choice of the basis) and µ F a factorisation scale 5 , taken to be of the order of Q ≫ m.Both the σi and the D i→Q are calculable in perturbation theory (the latter as long as the heavy quark mass is sufficiently larger than Λ QCD ), and they are known up to NNLO accuracy [22,[24][25][26][27].Note that the choice of a factorisation scheme is implicit in these formulas: neither the σi nor the D i→Q are physical objects, and depend on the scheme.In this paper we will use the MS scheme, and leave to future work the study of alternative factorisation schemes.
We shall denote the fixed order accuracy as LO, NLO or NNLO for a perturbative expansion up to order α 0 s , α 1 s or α 2 s respectively.Higher order calculations also introduce a running strong coupling α s (µ R ) and corresponding renormalisation scale µ R , and a beyondaccuracy dependency on the number of active flavours n f used in the renormalisation procedure.

Collinear resummation
In order to avoid the appearance of large unresummed logarithms α n s log k (µ F /m) of quasicollinear origin, Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations can be used to evolve the perturbatively calculated decay functions D i→Q from an initial scale µ 0F ≃ m up to a final scale µ F ≃ Q.We will use the DGLAP evolution implemented in MELA [33] to evaluate where the E ij (µ F , µ 0F ) represent the evolution factors stemming from the solution of the DGLAP evolution equations, leading to the resummation of collinear logs to leading [34] (α n s log n (µ F /µ 0F )), next-to-leading [35,36] (α n s log n−1 (µ F /µ 0F )), or next-to-next-to-leading [37,38] (α n s log n−2 (µ F /µ 0F )) logarithmic accuracy. 4The nomenclature 'decay function' was originally introduced in [31,32] to denote the fragmentation of a parton i into a hadron H, time-like counterpart to the 'parton distribution functions' (PDFs).While the use of the 'PDFs' nomenclature has remained current, the one of 'decay functions' has mostly run out of favour, 'fragmentation functions' being often preferred.In this paper we will refrain from using the term 'fragmentation function' for an unphysical object, preferring to reserve it for the observable-level cross section defined in eq.(2.2), and we will rather use the term 'initial conditions' or 'evolved initial conditions' or 'perturbative fragmentation functions' -all found in the literature -for the decay functions. 5Here and in the following we list in the arguments of the functions only the 'real' scales on which they explicitly depend, while we avoid writing also the 'residual' scales, on which a quantity can still depend beyond a given perturbative order.
The use of the decay functions D i→Q (µ 0F , m), perturbatively calculated at O(α s ) (nextto-leading order) [6,7] and O(α 2 s ) (next-to-next-to-leading order) [26,27] at µ 0F ≃ m, as inputs for the DGLAP evolution justifies the name of 'initial conditions' that is often employed to describe them and that we will use from here on.

Soft resummation
Both the coefficient functions and the initial condition D Q→Q exhibit large-logarithmic behaviour of soft origin in the x → 1 region (or, equivalently, large-N in moment space, α n s log k N with k ≤ 2n).Such logarithms can be resummed to all orders, and results are known up to NNLL accuracy for the initial conditions [9,19,29,30] and up to N3LL for the coefficient function [9,39,40].Explicit expressions are given in Appendix A.3 and A.4.
The resummed expressions exhibit a Landau pole related to the onset of non-perturbative phenomena.In order to push the prediction above this point (i.e. to large N or, equivalently, until x = 1), such a pole must be regularised.Various options are available, they are detailed in Appendix B, and we will study the quantitative impact that they have on the final result.
We denote a soft-resummed and regularised initial condition as where 'res' denotes soft-resummation up to a given logarithmic accuracy, and 'reg' the regularisation chosen for the Landau pole.D res,reg i→Q is defined in more detail in eq.(A.36), with the constants d(1) and d (2) given in eqs.(A.37) and (A.38).An equivalent expression for the coefficient function is defined in eqs.(A.23), (A.24) and (A.25).
Finally, we match the resummed expression to the fixed order (f o) one, so as to avoid double counting, and denote it as D f o+res,match,reg i→Q . (2.8) We can perform either a simple additive matching, match = add, or a so-called log-R one [41], match = logR (2.10) In these expressions, the notation [...] α p s denotes expansion and truncation up to order α p s , where p = 1 for NLO accuracy and p = 2 for NNLO accuracy.The '(,reg)' label in the expanded terms [...] α p s is a reminder that whatever regularisation procedure may have been employed in the resummed expressions, it can also be applied to the expanded version, even if it is not necessary because the latter does not exhibit a Landau pole.The choice of using or not the regularisation in both terms amounts to a freedom contributing to an uncertainty stemming from this procedure.
An equivalent notation is used for the coefficient functions.

Perturbative accuracy and scheme choices
The numerical value of the overall prediction for a heavy quark fragmentation function will depend on a number of choices, some discrete and some continuous ones: 1. the accuracy of the fixed order calculation used in the coefficient functions and in the initial conditions.We denote it as 'fo' and it can take the values 'LO', 'NLO' and 'NNLO' 2. (a) the logarithmic accuracy of the collinear resummation.We denote it as 'res coll ' and it can take the values 'LL', 'NLL' and 'NNLL' (b) the logarithmic accuracy of the soft resummation.We denote it as 'res sof t ' and it can take the values 'LL', 'NLL' and 'NNLL' Unless otherwise noted, we usually synchronise the soft and collinear resummation accuracies, res coll = res sof t , and collectively denote them through the value of 'res'.
3. the scheme used for the matching of the fixed order and of the soft-resummed calculation.We denote it as 'match' and can take the values 'add' (for additive) and 'logR' (for log-R matching) 4. the choice of regularisation of the Landau pole.We denote it as 'reg' and it can take the values 'noreg' (for no regularisation), 'CNO' (with an additional value for a parameter f ) [10], 'CGMP' [30].We will also briefly consider a modified CNO, to be denoted 'CNOmod'.Appendix B details the definitions of these regularisations.
5. the values of the renormalisation and factorisation scale around m and around Q, respectively µ 0R , µ 0F and µ R , µ F , can also affect the numerical value of the result.
We do not include in this list, for instance, the choice of the active flavour number scheme, and various options in the solution of the DGLAP equations as offered by MELA, for which we make a single choice, to be detailed in Appendix C.
A final result for σ Q can hence be labeled as7 Our default ("best") choice, to be justified below, will be (2.12)

Numerical Results
We produce numerical results setting the pole masses of the heavy quarks to m c = 1.5 GeV, m b = 4.75 GeV and m t = 172.5 GeV.The value of the strong coupling is set at a scale Q = 91.2GeV as α s (Q = 91.2GeV) = 0.118, and run from this initial condition to the scale where it is needed using an NNLO-accurate evolution equation in the MS scheme with a variable number of active flavours.
We evolve numerically all partons and the strong coupling using the package MELA [33].Details are given in Appendix C.
We set the thresholds for changing the number of active flavours in a variable flavour number evolution scheme (VFNS) at the same values as the quark masses, i.e. we use 3 flavours below m c , 4 flavours between m c and m b , 5 flavours between m b and m t , and 6 flavours above the top mass. 8  No other physical numerical parameter enters the perturbative predictions, other than the value of the centre-of-mass energy Q.

Initial conditions
Before considering the full, physical fragmentation functions, we study the initial conditions which, because of the low scales ≃ m at which they are evaluated (this is the case for the charm and the bottom quark), are the ingredient which is the most sensitive to the details of the perturbative setup.
We start by considering the bottom quark case, for easier comparison with the existing recent literature [29,30].In particular, we study in some detail the D b→b (N, µ 0F , m) component.
Figure 1 displays, as a function of the Mellin moment index N , the moments of the function D b→b (N, µ 0 = m, m) for the full range of perturbative orders, logarithmic accuracies of the soft resummation, matching schemes and regularisations of the Landau pole considered in this paper. 9 8 Since the initial conditions are always used in the n l = n f − 1 scheme, when they are evaluated with an α (n f ) s value -which can happen when one varies the renormalisation scale µ0R slightly above m -a compensating factor for the value of αs and the number of active flavours is introduced.The numerical impact of this procedure is anyway minimal. 9The discontinuities or irregularities in some of the curves denote the presence of the Landau pole for a specific value of N = N L 0 in the resummed expressions, with the initial condition becoming infinite at that point.As expected from the expression N L 0 = exp(1/(2b0αs(m)) (see Appendix B), this value is N L 0 ≃ 30 for bottom quarks.Beyond N L 0 the function becomes complex-valued, and what is plotted is its real part.
-7 - The initial condition for the bottom-to-bottom fragmentation function, without any evolution.'CNO' is always used here with the parameter f = 1.25.The plot with a grey background contains the curve that will be used as a reference, i.e. 'NNLO+NNLL logR CNO'.bottom initial condition, m=4.75 GeV The initial condition for the bottom-to-bottom fragmentation function, without any evolution, normalised to the curve 'NNLO+NNLL logR CNO' (the red curve in the plot with a grey background in figure 1).'CNO' is always used here with the parameter f = 1.25.
Figure 2 displays, for better readability, the ratio of all the curves in figure 1 to the curve 'NNLO+NNLL logR CNO', i.e. the red curve in the plot with a grey background in figure 1.
The main take-away from this survey is that no clear and consistent hierarchy of perturbative orders (i.e.LO > NLO > NNLO) can be observed, nor a systematic convergence can be seen as the perturbative accuracy increases, either in the fixed order or in the resummed cases, the one exception seemingly being the CGMP regularisation of the Landau pole, with a log-R-type matching of the fixed-order and resummed results.This also holds true for the charm and the top quark cases (see figures 12 and 13 in Appendix D.1).
While it is tempting to conclude, on the basis of this observation, that CGMP is an ideal choice for the regularisation 10 , we will see that it makes it impossible to describe charmed meson data with the same simple non-perturbative functional forms that instead work well for the bottom meson data.This lack of universality will lead us to prefer using the CNO regularisation in this work.It remains of course worth exploring if it is possible to combine the CGMP regularisation, and the better perturbative convergence it seems to lead to, with an effective non-perturbative parametrisation allowing for good fits to all the data.
In order to better appreciate the impact of different non-perturbative regularisations of the Landau pole on the otherwise perturbative initial conditions, we show in figure 3 the NNLO+NNLL curve with log-R matching with different regularisations and for all three heavy quarks, both in N and in x space.We see that the resulting curves can be very different, especially in the large-x region.The non-perturbative contributions that will have to be added in order to fit the experimental data, where available, will consequently also be very different.In some cases, the shape of the perturbative fragmentation function can make it impossible to fit the experimental data with simple non-perturbative forms: we will see that this is notably the case for charmed hadron data and CGMP regularisation, which leads to a flat and null initial condition in a large region x ≳ 0.8 (this characteristic remains after evolution and convolution with the coefficient function).
In figure 3 we have also introduced a new regularisation for the Landau pole, dubbed 'CNOmod', which is meant to be CNO but implemented as it is presumably done in ref. [30], i.e. rescaling only the resummation term but not the subtraction one (which is devoid of a Landau pole, since it is expanded at a fixed order in α s ) in the matching (see appendix B).We can see that the feature of the initial condition with CNO not vanishing for x → 1, pointed out in ref. [30] as being problematic (cfr.figure 2 there and the text below) may in fact stem from the specific implementation of CNO that that work adopts. 11When CNO is implemented as in [10], rescaling the resummed term but also its expansion, the initial condition appears to vanish quickly when x → 1 (see also Appendix B).
It is also worth noting how, in figures 1 and 2, the plots employing additive matching and CNO regularisation of the Landau pole with NLO and NNLO fixed order appear to be outliers to some extent, with the moments becoming very rapidly small and eventually negative at large N .This behaviour can eventually be traced back to the feature of the x = 1 improves, but at the expense of a less effective screening of the Landau pole.
-11 - CNO regularisation of rescaling both the resummation term and the subtraction one: this leaves the (unrescaled) logarithms present in the fixed order partially unsubtracted, leading to the observed patterns.Employing the CNOmod regularisation instead of the CNO one, i.e.only rescaling the resummation term, prevents this behaviour, as shown in figure 4.
In this work we have nevertheless chosen CNO as our default regularisation of the Landau pole, in part because it better preserves the exact fixed order calculation at small N , and in part because of easier comparison with ref. [10].When employing a log-R rather than an additive matching, which will also be our default, the difference is minimal anyway.

Full e + e − fragmentation function
The full fragmentation function in e + e − collisions is obtained convoluting the coefficient functions with the evolved initial conditions, summing over all parton channels, as shown in eq.(2.5).We present here numerical results for fragmentation to bottom quarks, σ b , and we study theoretical uncertainties stemming from missing higher perturbative orders by varying the renormalisation and the factorisation scales, both around the initial scale µ 0 , set equal to the heavy quark mass, and around the final scale Q, the centre-of-mass energy of the collision.More specifically, we will consider for the coefficient functions, and for the initial conditions, resulting in the customary 7-point rule in the first case, and in a 5-point rule in the second.Figure 5 shows the uncertainty bands, in N -and in x-space, due to µ 0R and µ 0F variation in the region detailed above, for four different perturbative approximations, NLO, NNLO, NLO+NLL resummation (collinear and soft) and NNLO+NNLL resummation (collinear and soft), the latter two with log-R matching and CNO or CGMP Landau pole regularisation.In the N -space (left) plots, one can see some degree of improvement (i.e.narrower bands) as the accuracy of the perturbative series increases.However, on one hand the NNLO+NNLL band (red) is not significantly narrower than the NLO+NLL (blue) one, and on the other hand the two bands are hardly or not at all overlapping (especially in the CNO regularisation case), indicating poor convergence of the series.This poor convergence, and the significant differences between the CNO and the CGMP cases, are consistent with the observations made in section 3.1.One last consideration that we make is that Scheme with n flavours from a starting scale below the mass of the n th flavour.While this could in principle be fixed, the remaining variations already give a sufficiently realistic picture of the uncertainty, and in most practical cases the initial scales will be fixed at m anyway.Lifting this limitation was therefore considered not to be a priority.the shape and the position of the bands in x-space do not tell us much about convergence.Moreover, they can drastically depend on the chosen regularisation for the Landau pole, again as observed in section 3.1.
In figure 6 we repeat the same exercise, but varying this time the scales µ R and µ F around Q. The bands are much narrower, due to the value of α s being significantly smaller around the large scale Q = 91.2GeV than around the mass scale m = 4.75 GeV.Because of this smaller value of α s , the expected perturbative hierarchy is better respected, with the bands becoming generally narrower as the perturbative precision increases.As was the case in figure 5, bands in x-space are instead not easily interpretable in terms of behaviour of the perturbative series.
Analogous figures for the charm and top quark cases can be found in Appendix D.2.

Comparisons to Experimental Data
In this section we study how the fragmentation functions that we have calculated, possibly complemented by a non-perturbative component, fare in describing D and B meson fragmentation data measured in e + e − collisions.We consider the following datasets: Bmesons production data at the Z peak (Q = 91.2GeV), as measured by the ALEPH [42], OPAL [43], SLD [44] and DELPHI [45] collaborations 13 ; D-mesons production data at the Z peak, as measured by the ALEPH [46] collaboration; D-mesons production data near the Υ(4S) peak (Q = 10.6 GeV), as measured by the CLEO [47] and BELLE [48] collaborations. 14Electromagnetic initial-state radiation (ISR) effects are accounted for in the D-mesons production data near the Υ(4S) peak, according to the procedure detailed in section 3 of [10].
We first show how perturbative curves for different perturbative orders (NLO+NLL and NNLO+NNLL) and different prescriptions for the regularisation of the Landau pole (CNO, CGMP) compare to representative data points for B-meson production in figure 7 and D-meson production in figure 8.We present curves both in Mellin and in direct space, respectively on the left and the right panels.
As already stressed in section 3, for the bottom quark the behaviour of CNO and CGMP is similar at NNLO, and both results can be easily used as perturbative contri- 13 While the ALEPH set refers specifically to B-mesons, the SLD, DELPHI and OPAL data are for all weakly decaying bottom-flavoured hadrons and include about 10% of baryons. 14The BELLE data from ref. [48] have been removed from HEPData in 2014 'at the request of the authors due to an unrecoverable errror in the measurement' (see https://www.hepdata.net/record/ins686014).However, since the corresponding paper does not seem to have been formally retracted, and the CLEO data appear compatible where comparable, we have decided to keep using these BELLE data in this paper, in part to facilitate comparisons with ref. [10].butions to fit a non-perturbative component to the data.For the charm quark instead, CNO and CGMP behave qualitatively and quantitatively differently, both at NLO and at NNLO.In direct space (figure 8, right panel), we see that both CNO with f = 1.25 and CGMP feature a sharp peak around x ∼ 0.7-0.8 and a steep drop around x ∼ 0.8-0.9.
In moment space (figure 8, left panel), such behaviour translates to the NNLO curves dipping below the experimental data points at N ≳ 6.This fact immediately prevents one from fitting the data using the NNLO result and simple non-perturbative functional forms: their Mellin moments being normalised to values smaller than one (because these non-perturbative forms are meant to degrade the quark momentum when it fragments into a hadron), they can only lower the value of the theoretical prediction, but not raise it.A simple practical solution to this problem is to increase the value of the free parameter of the CNO regularisation to f = 2. Figure 8 shows that, in this case, for the charm quark the NNLO+NNLL curves remains above the data in Mellin space up to larger values of N (left panel), suggesting that a simple non-perturbative form may be able to describe the data.The right panel in the same figure shows that, indeed, the x distribution is closer to the shape of the data.

Single-moment and single-parameter fits
Here we try to describe single moments of experimental data using a very simple nonperturbative functional form, dependent on a single parameter, on top of our perturbative description.We use the form introduced by Kartvelishvili et al. [49] whose Mellin transform is  In the following we present our results for the values of the parameter α for B and Dmesons data.We determine the size of the non-perturbative component in Mellin space by dividing the moments of the experimental data by the perturbative result according to eq. (2.3).
The result for B mesons is shown in the upper panel of figure 9, for values of N between 2 and 8.This range has been chosen because, if one were to use the results of these fits in hadronic collisions, where steeply falling transverse momentum distribution ∼ 1/p 5  T would have to be convoluted, the moments at N ≃ 5 would provide the dominant contribution to the result.
In the lower panel of figure 9, for each moment N , we show the value of α determined by inverting eq.(4.2).We notice that both at NLO+NLL and NNLO+NNLL the extracted values of α are fairly consistent at different values of N , showing that this simple nonperturbative form manages to describe fairly well the shape of the whole distribution.We also observe that the larger value of α extracted at NNLO than at NLO is consistent with the comparison shown in figure 7, since the NNLO curve is closer to the data and needs therefore to be supplemented with a smaller non-perturbative component (i.e. a harder We consider similar fits to D meson data.The results are shown in figure 10 where we consider experimental data for D * + and D * 0 production15 .At NLO+NLL, the data can be described quite well with a single value of α within uncertainties, independently of the value of f used.Instead, and at variance also with the single-moment B meson results of figure 9, this is not the case at NNLO+NNLL, where a large dependence of α on N can be observed in the lower panel of figure 10.We also observe that at NNLO level it is mandatory to choose f = 2.0, as shown by the perturbative result/data comparison in figure 8.This points to a stronger dependence of the non-perturbative component of the D meson fragmentation on the details of the Landau pole regularisation applied to the NNLO+NNLL perturbative result. Numerical results for D np and for the α parameter of eq.(4.1) for B and D mesons are collected in tables 1 and 2. At NLO+NLL level they are not identical to those presented in table 4 of [10], but we have checked that this is exclusively an effect of having used different perturbative parameters (essentially the value of α s ).This is of course legitimate, the non-perturbative component being unphysical by itself and dependent on the setup used in the perturbative part.

Charm ratio
We now consider an observable which is derived from experimental data, but which is directly comparable to perturbative predictions, namely the ratio of Mellin moments of charmed hadrons data measured at LEP, at centre-of-mass energy Q f = 91.2GeV, and at BELLE and CLEO, at As shown by the arguments in the equation above, the moments at the two energies Q i and Q f depend individually on the mass of the heavy quark and on the parameters that describe the non-perturbative component that must be added to the perturbative part.However, these 'low scale' dependencies drop out in the ratio 16 , leaving R D (N, Q i , Q f ) to depend only on the two 'high' scales Q i and Q f .Since moments at these two scales are bridged by the DGLAP evolution, this ratio essentially measures the degree of correctness (or lack thereof) of the collinear factorisation picture that underlies the collinear resummation.Deviations from this picture, for instance under the form of power corrections larger than expected, can show up as a discrepancy between the experimental ratio and the perturbative prediction.
This comparison had already been performed in [10] using the NLO+NLL perturbative prediction for the fragmentation function and the D * + fragmentation data from ALEPH [46] and from BELLE [48] (see footnote 14).A visible discrepancy with the data had been observed.We update that comparison to NNLO+NNLL in figure 11.We see that the increased perturbative accuracy shrinks the uncertainty bands but, as expected and already mentioned in [10], it does not help in reconciling the prediction with the data, the discrepancy remaining largely the same.

Conclusions
We have implemented in a C++ code, that will eventually be publicly released, the NNLO+NNLL calculation for heavy quark fragmentation in e + e − collisions.The code can calculate charm and bottom production at B factories and LEP energies, as well as top production at a hypothetical higher energy collider.A non-perturbative component, describing the fragmentation of the heavy quark into the heavy hadron, can be modeled with a simple function and convoluted with the perturbative result.-21 - The first important observation that we make is that increasing the perturbative accuracy from NLO+NLL to NNLO+NNLL does not necessarily improve the quality or accuracy of the prediction.The need to regularise the soft-gluon resummation to avoid the non-perturbative Landau pole means that the result depends on this procedure, whose choice is not unique.We observe that, depending on the choice, the numerical value and the quality of convergence of the perturbative series can vary widely.
Another observation is that, while the description of all heavy quarks (top, bottom and charm) is in principle perturbatively identical, important qualitative and quantitative differences can appear.
Top production is usually very well behaved.Due to all scales involved being large, and the value of the strong coupling consequently small, one observes good convergence and limited dependence on the choice of Landau pole regularisation.Unfortunately, this is also the least phenomenologically relevant case.
Bottom production is quite dependent on the choice made for the Landau pole regularisation, but the fragmentation function calculated in perturbative QCD has perturbative uncertainties that are still under control, and it provides a sufficiently good baseline: a simple model for the non-perturbative component is usually sufficient to provide a numerically small correction and achieve a reasonably good description of the data.
For charm production, due to the very low scales and large value of the coupling, all the sensitivities mentioned above are magnified, as one can expect.There are significant perturbative uncertainties, and significant quantitative differences between regularisation schemes.Moreover, in many cases the perturbative description provides a very poor de-scription of the data, to the extent that a simple non-perturbative model cannot compensate sufficiently to describe them.
We conclude by observing that the time-honoured 'perturbative QCD convoluted with a simple non-perturbative function', that has served us well up to NLO+NLL accuracy, appears to be challenged by the upgrade to NNLO+NNLL accuracy.The upgrade to a higher perturbative level leads to larger sensitivity to the details of the regularisation at the boundary with non-perturbative phenomena, largely washing out -at the phenomenological level and especially for charm and bottom quarks -any gain that one may have hoped to glean from the higher perturbative accuracy.We find it likely that the achievement of a really improved theoretical predictivity will necessitate a more thorough study and modeling of the interface between perturbative and non-perturbative regions.
We can rewrite eq.(A.1) as where the sum runs over n f quark flavours.This form (obtained after summing over transverse and longitudinal components) makes it explicit the contribution of the leading order electroweak cross sections for production of a quark in e + e − collisions, the σ i , and that of the non-singlet, gluon and singlet QCD coefficient functions, C N S q , C g and C S q , described in Appendix A.3.
The electroweak cross sections σ (0) i are given by 17 In the above and Γ Z mass and width of the Zboson 18 respectively and e i quark's i electric charge.The electroweak couplings of the charged lepton ℓ with charge e ℓ = −1 and quarks are with θ W Weinberg angle and sin 2 θ W = 0.23.
We also introduce the total cross section for production of a heavy quark Q = c, b, t, with n f number of active flavours (including the heavy quark) and QCD colour factors 17 As in [10], we complement the electroweak cross sections with a threshold factor for the heavy quarks and antiquarks, The numerical impact is however usually negligible at the combinations of heavy quark type and centre-ofmass energy Q considered in this paper. 18We adopt the values of MZ = 91.2GeV and ΓZ = 2.5 GeV in our numerical implementation.
The D i→Q functions in eq.(A.1) are the perturbative fragmentation functions at the hard scale of the process, namely initial conditions D j→Q (µ 0F , m) evolved through DGLAP equations from the initial factorisation scale µ 0F ∼ m to the final one µ The initial conditions and their soft-gluon resummation are described in detail in Appendix A.4.Their DGLAP evolution is performed using the MELA package as described in Appendix C.
The relations between the partonic initial conditions D i→Q and the singlet and nonsinglet combinations that enter eq.(A.2) are 19 Finally, we can rewrite eq.(A.2) in a partonic basis (i.e.recover the form of eq.(A.1)) as with At NLO the singlet and the non-singlet coefficient function are equal and the first term cancels.At NNLO the difference C S q − C N S q is the pure-singlet (PS) contribution, which vanishes in the large-N limit.In the partonic-channel basis, the quarks coefficient functions contain the quark's charge, leading to a distinction between (anti)up-type and (anti)downtype.
The coefficient functions and the soft-gluon resummation of the C i in eq.(A.11) are described in detail in Appendix A.3.

A.2 Perturbative parameters
We list here the parameters that enter the perturbative results and the soft-gluon resummation up to NNLL accuracy.
The strong coupling α s (µ) evolves as dα s d ln µ 2 = β(α s ) , (A.12) 19 Since we are inclusive over the angle of emission of the identified hadron with respect to the electron beam direction in the center-of-mass frame, we have no contribution from the asymmetric component defined in [22].
and the coefficients of the QCD β-function in a theory with n massless flavours read b In practice, we use the evolution of the strong coupling α s (µ) as implemented in MELA, see Appendix C.
We define The coefficients of the cusp anomalous dimensions A[α s ] up to O(α 3 s ) and of the functions [9,19,29,40] are For the resummation of the coefficient functions up to NNLL we use {b and for the resummation of the initial conditions up to NNLL {b with n l = n f − 1.

A.3 QCD corrections to the coefficient functions
The coefficient functions are process-dependent objects, computed in perturbative QCD with n f massless partons.In this section, it is understood that α s = α . The (sum of transverse and longitudinal) singlet and non-singlet combinations up to O(α 2 s ) read Since the two basis described in Appendix A.1 are degenerate at NLO, the â(1) i coefficients are insensitive to the basis chosen.The coefficients â are given e.g. in [7,50] at NLO and in [22] (x-space), [24,25] (N -space) at NNLO.The factors 1/4 are meant to make contact with the notation of [25], where the expansion is given in terms of a s = α s /(4π), and whose â(2) expressions are used in the Fortran implementation that we integrate into our code.The terms proportional to α 2 s ln(µ and C S g account for the µ R ̸ = µ F scale-dependence from the NLO running of the strong coupling, with the special functions S i (N ), written in terms of the polygamma functions ψ m (N ).
To account for soft-gluon effects, the quark coefficient function C i of eq.(A.11)20 must be resummed to all orders in α s .The Sudakov-resummed coefficient function, up to NNLL accuracy is where ā(1) q (Q, µ F ) and ā(2) q (Q, µ F ) are the constant (i.e.N -independent) terms of the large-N limit of â( 1) We had extracted the constant ā q from the non-singlet coefficient function presented in [24], and we have found it to coincide with the one successively published in [30].
The Sudakov factor C sud q takes the exponential form [9,30,40] ln where we have defined λ = α s (µ R )b 0 ln N and the functions A[α s ] and B[α s ] are given in Appendix A.2.The function D[α s ] can be neglected at NNLL accuracy [9,40].
The functions g k up to at least k = 3 are available in the literature [9,30,40].We have recalculated them and found perfect agreement.They read For completeness, we give here the expansion in α s of the Sudakov factor up to O α 2 s , to be used in the matching of the resummed result to the fixed order one: (A.32)

A.4 Initial conditions
As first introduced in [7], the initial conditions D i→Q for heavy hadron production are perturbative objects, and therefore admit a perturbative expansion in α s .The O(α 2 s ) corrections were computed in [26,27] in x-space and more recently transformed to Nspace [29].The calculations of [26,27] were performed in an MS renormalisation scheme for ultraviolet divergences where both massless and massive flavours contribute to the evolution of the running coupling (i.e. an n f -scheme).To match the fixed order results with the resummed one, where only massless flavours contribute the evolution of the running coupling (an n l scheme), the initial conditions must be written in the n l -scheme.In Nspace they read, expanded in powers of α (n l ) s , i.e. renormalised with n l = n f − 1 massless flavours, where Q = c, b or t, and q are the flavours lighter than the quark identified as heavy, and therefore considered massless.The terms proportional to α 2 s ln(µ 2 0F /µ 2 0R ) account for the renormalisation-scale dependence at NNLO from the running of the coupling at NLO.The terms proportional to α 2 s ln(m 2 /µ 2 0R ) in eq.(A.33) account for the change of scheme, as described in [29].For sake of clarity, in eq.(A.33) b shows explicitly the number of flavours used in its definition. 21he fixed order NLO coefficients, see e.g.[7,50], are given by while the NNLO coefficients d i→Q can be found in refs.[26,27].To account for soft-gluons effects, the heavy quark initial condition D Q→Q must be resummed to all orders in α s .The Sudakov-resummed heavy quark initial condition reads (A.36) The constants d( 1) Q→Q (µ 0F , m) and d( 2) Q→Q (µ 0F , m) from the large-N limits of d Q→Q (N, µ 0F , m) and d (2) Q→Q (N, µ 0F , m) are given in [9,26].We copy them here as written in [29] for completeness: The Sudakov factor is where we have defined λ 0 = b s (µ 0R ) ln N , and used we can write explicitly the g ini k functions (that we have independently verified) as [9,29,51] For completeness, the expansion of the Sudakov factor up to O α 2 s is with coefficients having used again eq.(A.40).

B Landau pole prescription
As a manifestation of non-perturbative phenomena and eventual failure of perturbation theory, the Sudakov factors of initial condition eq.(A.39) and coefficient function eq.(A.26) present singularities at large-N values.The latter manifest for specific values of the parameters λ and λ 0 at the so called Landau poles, sensitive to the values of α s .
For the coefficient function the Landau pole N L is at and for the values of α s used appears at very large values of N (> 100), having no impact on the perturbative prediction.For the initial conditions instead the Landau pole N L 0 is at and with typical values of α s and n f we find that for the bottom N L 0 ∼ 30 while for the charm N L 0 ∼ 7 as can be seen in figures 1 and 12.The behaviour of the cross section σ Q (N ) for production of a heavy quark in the large-N region is dominated by the Landau pole of the initial condition and, for N > N L 0 , the moments can become negative and unphysical.In order to produce sensible phenomenological results, a prescription must be applied to regulate the pole, and many choices can be made.We notice that a reasonable Landau prescription should not introduce power corrections that are larger than generically expected for the process in question [52].
A prescription that respects this condition is the so called CNO prescription (Cacciari-Nason-Oleari), that regulates the pole via a tunable parameter f .It amounts to the following remapping of N in the initial conditions and analogously in the coefficient functions replacing N L 0 with N L .The Landau singularity is regulated when f ≥ 1 and we have chosen f = 1.25 for the bottom and f = 1.25 or f = 2.0 for the charm.
In general it is ambiguous whether one should apply such a prescription also to the expansion of the Sudakov factors when matching to the fixed order, since it is free from the Landau singularity.Both choices are acceptable and regulate successfully the Landau pole, but give different results.As in section 3.1, we denote as CNO the prescription in which the re-mapping is applied to the Sudakov factor and its expansion, and with CNOmod the prescription in which the remapping is applied only to the Sudakov factor.Since the resummed result and its expansion are both rescaled and cancel well, the CNO choice ensures that the fixed-order result is well preserved at small N .On the other hand, at large N the (unrescaled) fixed-order result and the (rescaled) expansion of the resummation will not cancel perfectly.While this could be seen as spoiling the resummation, it mainly happens around and beyond the Landau pole, where phenomenology is anyway expected to be mainly driven by empirical modelling of non-perturbative effects.CNOmod ensures the large-N cancellation, but it does not preserve the fixed-order calculation at small N .Finally, we note that a third choice may in principle be legitimate, namely rescaling all three components (fixed-order, resummed, and expansion) of the matched result.We refrain from considering this option, which has the same properties of CNOmod.
As mentioned in section 3.1, we have decided to adopt the CNO prescription for our comparisons with the data since the CGMP prescription yields poorly behaving results for the charm.This can be seen both in N -space, where the CGMP-regularised cross section becomes smaller than the data at relatively small N -values, and in x-space, where the perturbative cross section has a steep fall to zero around x ∼ 0.8 (see figure 8), making it impossible to describe the data with any of the simple non perturbative prescriptions used in [10] 22 .

C MELA Setup
In the following we detail the choices made in the usage of the evolution library MELA [33].
The obvious choice for our process is the evolution of type 'TIME', setting to 'OFF' the polarization mode.We fix the method used in the DGLAP solution to 'TRN' and the quark masses to be of 'POLE' type.We observe negligible discrepancies under the choice of the DGLAP solution method.Our inputs for the values of quarks masses and flavour thresholds are given in section 3.
For the value of α s at a given scale in the initial conditions and coefficient functions we have always used the internal routine of MELA for α s in a VFNS starting from a reference value of α s (91.2 GeV) = 0.118.Where needed, we apply a correction at the threshold for the scheme used in the calculation of the fixed order ingredients.
For the evolution, the 'VFNS' mode is always our default.Given that the matching conditions at thresholds [52] are not available in the time-like case at NNLO accuracy, the VFNS of MELA can normally only be used with a NLO evolution.With a small modification in the source code we have allowed for a NNLO evolution with only NLO matching thresholds.Such modification is now available in the branch https://github.com/vbertone/MELA/tree/timelike_thresholds.The numerical impact of the unknown and missing NNLO matching thresholds is expected to be negligible. 23

D.1 Initial conditions
For completeness we show here, in figures 12 and 13, the analogous of figure 2 for charm and top quarks.The pattern of convergence, or lack thereof, of the perturbative series and of the resummation is consistent with the bottom quark case, though respectively more and less extreme for charm and top quarks respectively.This is due to the relevant scaleequal to the quark mass -being smaller in the charm case and larger in the top case, and the corresponding strong coupling consequently larger and smaller.Again for completeness, we show here the scale dependence of the full e + e − fragmentation function into charm (at Q = 91.2GeV) and into top (at Q = 1000 GeV) quarks.Note that, at variance with the bottom and top quark cases, for the charm quark we only vary the renormalisation scale in the range m < µ 0R < 2m, to avoid reaching too deep into a barely perturbative region.These figures are to be compared with the ones for bottom quarks, in figures 5 and 6.
One can see that in the charm quark case (figures 14 and 15) the uncertainty bands are wider than in the bottom quark case.In the case of scale variations around the initial scale m, and despite the reduced range of variation of the renormalisation scale, any perturbative hierarchy is essentially absent, and the bands are so large that they make perturbative predictivity quite weak.Sensitivity on the choice of regularisation of the Landau pole, clearly seen in the x-space plots, also appears to be even larger than in the bottom case.For the top quark case in figure 16 instead, thanks to the small value of α s at all scales involved in the process (m = 172.5 GeV and Q = 1000 GeV), the perturbative behaviour is much better.The bands in N -space display the expected hierarchy, and in x-space all curves are very similar, except very close to x = 1, which is the only region shown in the plot.While the top case is presently not of phenomenological interest, it serves as a cross-check that the 'confusing' behaviour observed in the bottom quark and -to an even larger extent -in the charm quark case can be ascribed to the large value of α s in play there.We do not show results for variations of the final scales around Q for the top case as they result in bands even narrower than in the bottom case and are of no pedagogical value.

2 2
Full e + e − fragmentation function 4 Comparisons to Experimental Data 4.1 Single-moment and single-parameter fits 4Full e + e − fragmentation function References 1 Introduction

Figure 1 .
Figure1.The initial condition for the bottom-to-bottom fragmentation function, without any evolution.'CNO' is always used here with the parameter f = 1.25.The plot with a grey background contains the curve that will be used as a reference, i.e. 'NNLO+NNLL logR CNO'.

Figure 3 .
Figure 3. N -space (left column) and x-space (right column) representation of the perturbative initial conditions for the three heavy quarks, evaluated at NNLO+NNLL accuracy with log-R matching, with different kinds of Landau pole regularisations.In all cases µ 0F = µ 0R = m, and α s (91.2 GeV) = 0.118.CNO and CNOmod are always implemented with f = 1.25.

Figure 5 .
Figure 5. Bottom fragmentation function.Uncertainty bands due to variation of the µ 0R and µ 0F scales around m. Left plots in N -space, right plots in x-space.Upper plots with CNO regularisation (f = 1.25), lower plots with CGMP regularisation.

Figure 6 .
Figure 6.Bottom fragmentation function.Uncertainty bands due to variation of the µ R and µ F scales around Q. Left plots in N -space, right plots in x-space.Upper plots with CNO regularisation (f = 1.25), lower plots with CGMP regularisation.

Figure 9 .
Figure 9.Single-point fits to B-meson data from experiments at the Z peak, with the oneparameter non-perturbative expression of Kartvelishvili et al. in eq.(4.1).Upper pane: value of normalised non-perturbative moment, fitted to a single, specific moment.Lower pane: corresponding value of α.NLO stands for 'NLO+NLL logR', while NNLO stands for 'NNLO+NNLL logR'.

Figure 11 .
Figure 11.ALEPH/BELLE and ALEPH/CLEO ratios for D * + data, compared to theoretical predictions.The BELLE and CLEO data are corrected for initial state electromagnetic radiation effects[10].

Figure 12 .
Figure 12.The initial condition for the charm-to-charm fragmentation function, without any evolution, normalised to the curve 'NNLO+NNLL logR CNO'.'CNO' is always used here with the parameter f = 1.25.

Figure 13 .
Figure 13.The initial condition for the top-to-top fragmentation function, without any evolution, normalised to the curve 'NNLO+NNLL logR CNO'.'CNO' is always used here with the parameter f = 1.25.

Figure 14 .
Figure 14.Charm fragmentation function.Uncertainty bands due to variation of the µ 0R and µ 0F scales around m. Left plots in N -space, right plots in x-space.Upper plots with CNO regularisation (f = 1.25), lower plots with CGMP regularisation.Note that, at variance with the bottom and top quark cases, here we only vary the renormalisation scale in the range m < µ 0R < 2m.

2 NLO 2 NLO 2 NLO 2 NLOFigure 15 .
Figure 15.Charm fragmentation function.Uncertainty bands due to variation of the µ R and µ F scales around Q. Left plots in N -space, right plots in x-space.Upper plots with CNO regularisation (f = 1.25), lower plots with CGMP regularisation The little spikes that one can observe on the curves in x space, here and possibly elsewhere, are artefacts of the numerical Mellin inversion.

Table 1 .
Values of the first eight moments of the non-perturbative fragmentation function D np (N ), as extracted from the experimental data.

Table 2 .
Values of the parameter α of the non-perturbative fragmentation function D np K (N ), corresponding to the values of the first eight moments extracted from the experimental data and displayed in table 1.