First global QCD analysis of charged hadron fragmentation functions and their uncertainties at next-to-next-to-leading order

In this paper, we present ${\tt SGK18}$ FFs, a first global QCD analysis of parton-to-{\it unidentified} charged hadrons fragmentation functions (FFs) at next-to-next-to-leading order (NNLO) accuracy in perturbative QCD. This analysis is based on single-inclusive charged hadron production in electron-positron ($e^-e^+$) annihilation. The uncertainties in the extraction of ${\tt SGK18}$ FFs as well as the corresponding observables are estimated using the"Hessian"technique. We study the quality of the {\tt SGK18} FFs determined in this analysis by comparing with the recent results in literature. We also show how {\tt SGK18} FFs results describe the available data for single-inclusive unidentified charged hadron production in $e^-e^+$ annihilation. We demonstrate that the theoretical uncertainties due to the variation of the renormalization and factorization scales improve when NNLO QCD corrections are considered. We find that the resulting {\tt SGK18} FFs are in good agreement with all data analyzed and the inclusion of NNLO corrections tends to improve the data description with somewhat smaller uncertainty.

Like PDFs, fragmentation functions (FFs) also plays an important role in our understanding of certain high energy processes with identified hadrons in the final state [42]. According to the asymptotic freedom of QCD, fragmentation functions (FFs) relates to the longdistance dynamics of the interactions among quarks and gluons which cause to their hadronization in a hardscattering process [43,44]. The experimental observables of single-inclusive hadron production involve identified hadrons in the final state. In order to obtain theoretical predictions for such processes, it is written down a factorization formula and FFs are convoluted with partonic cross-sections. FFs plays an important role in understanding of non-perturbative QCD dynamics. Like for the case of PDFs, FFs are determined from global QCD analysis of experimental measurements particularly in hadronization processes. These processes include singleinclusive hadron production of electron-positron (e + e − ) annihilation (SIA), semi-inclusive deep-inelastic leptonnucleon (ℓ ± N) scattering (SIDIS), and single-inclusive hadron production in proton-proton (pp) collisions. According to the QCD-improved parton model, FFs and PDFs scaling violations are subjects to the perturbatively computable Dokshitzer-Gribov-Lipatov-Alteralli-Parisi (DGLAP) evolution equations [45][46][47][48].
In the past few years, several progresses have been done to determine FFs for light and heavy mesons which performed at next-to-leading order (NLO) and next-tonext-to-leading order (NNLO) accuracy in perturbative QCD [42,49,50]. At NNLO, only experimental data from electron-positron annihilation can be used in a QCD analysis while the calculations for the hard processes in SIDIS and pp collisions at NNLO are not accessible yet. Electron-positron annihilation provides the most cleanest and appropriate data sets to access to the FFs, where the final state quarks and gluons fragment into hadrons. Compared to the SIDIS and pp collision, the FFs in SIA processes are the only non-perturbative functions in the calculation of the cross section. Since the SIA experimental measurements are not sensitive to the separation between quark and antiquark FFs, the extraction of favored and disfavored fragmentations is difficult in this processes. In order to allow one to determine quark from antiquark FFs, the data where hadrons of different electrical charge are identified in the final state needs to be taken into account [51]. In addition, the gluon fragmentation density is not exceedingly well constrained by SIA data, since the subleading NLO and NNLO corrections for e − e + annihilation are too weak to determine it. Including the data from SIDIS and pp collision in the extraction of FFs could increase the statistics and also provide a much more complete picture of the fragmentation processes. In SGK18 FFs we restrict our analysis to the SIA data, since the QCD framework for FFs at NNLO only accessible for the e − e + annihilations.
Historically, the knowledge of FFs and their determination through the global analysis of the experimental data has undergone many developments both experimentally and theoretically. For example, in the analysis of Ref. [52], the authors used simultaneously the SIA and SIDIS asymmetry data from the HERMES [53] experiment at HERA and COMPASS [54,55] experiment at CERN to determined the pion and kaon FFs both at leading order (LO) and NLO approximation. In the analysis of Ref. [56], the authors considered the finite-mass effects of the proton to calculate the proton FFs by including SIA data at LO and NLO accuracies. Recently, pion, kaon and proton FFs have been extracted by various groups such as the DEHSS [57,58], HKKS [59], JAM [60], and also by the NNPDF Collaboration [42] using the iterative Monte Carlo method. For the case of charmed-meson D * FFs, we refer the reader to the very recent AKSRV17 [61] and SKM18 [50] global analyses. It should be note that the later one has been done at for the first time at NNLO approximation by including the SIA data.
In this paper, we perform for the first time a comprehensive QCD analysis to obtain a set of unidentified charged hadron FFs and their uncertainties at NNLO. In order to perform our global analysis for determining the FFs of the unidentified charged hadrons at NNLO, we have to limit the potential of global determination of FFs to the SIA measurements. We show that the inclusion of higher order QCD correction could describe the data well, including those data points at rather smaller hadron momentum fraction z, z < 0.02. We extensively discuss the theoretical and phenomenological methodology of the SGK18 analysis, including the exprimental description of the e − e + annihilation exprimental observables in term of the SGK18 FFs, parameterizations and the fitting procedure in next sections of this paper.
Previously available analyses of inclusive charged hadron FFs sum up the pion, kaon and (anti)proton results and ignore the contributions of possible heavier charged hadrons. For example, in Ref. [62], the inclusive charged hadron experimental data have been excluded from the analysis and only the sum of charged pion, kaon and protons FFs obtained from the fit has been compared to the inclusive charged hadron data. While in Ref. [63], the FFs for unidentified charged hadrons has been extracted. The NNPDF Collaboration, after having extracted trustworthy FFs for pion, kaon and proton, as the lightest and most copiously produced charged hadrons, has recently calculated the FFs for unidentified charged hadron up to NLO accuracy [64]. In addition, the analysis by DSS07 [65] included the electron-proton annihilations, SIDIS and proton-(anti)proton collisions experimental data sets. They obtained the contributions from the residual charged hadrons as well as pions, kaons and (anti)protons to the unidentified charged hadron FFs up to NLO accuracy. In this work, we extend the extraction of FFs for charged hadrons for the first time up to NNLO by including the inclusive charged hadron experimental data from e − e + annihilation.
In order to assess the uncertainties of the resulting unidentified charged hadron FFs at NLO and NNLO accuracies as well as the corresponding observable, associated with the uncertanties in the analyzed data, we have applied the "Hessian" method. This paper is organized as follows: In Sec. II, the datasets included in SGK18 FFs analysis, along with the corresponding observables and kinematic cuts are presented. We discuss the QCD analysis of hadronization process in electron-positron annihilation by introducing FFs and their evolution in Sec. III. We describe our formalism, input parametrization at the initial scale for the determination of unidentified charged hadron FFs in Sec. IV. In Sec. V, the minimization strategy and the "Hessian" uncertainty approach to calculate the errors of SGK18 FFs analysis are presented. In Sec. VI, we present the obtained results for the D h -FFs and their uncertainties. We also perform a comparison of SGK18 results with the analyzed experimental data and other available FF sets in this section. The theoretical uncertainties, fit quality and the stability due to the variation of the renormalization and factorization scales are studied at the end of this section. Finally, we conclude and summarize the results in Sec. VII.

II. EXPERIMENTAL OBSERVABLES
We begin this section with discussing the measurements of charged hadron production in e + e − annihilation, collected by a variety of experiments [66][67][68][69][70][71][72][73] at CERN, SLAC and HERA. Our aim is to include all available data sets which help to constrain the resulting charged hadron FFs, and more importantly, provide additional consistency checks of the fitting procedure.
In this analysis, the FFs are determined by including a wide range of the experimental data from electronpositron annihilation into an unidentified charged hadron h and the unobserved jets which are produced along with the detected hadron h. This process is given by: The DIS process is space-like, while the above process is time-like and the related scaling variable is z = 2p h .q/Q 2 , in which the four-momenta of the intermediate gauge boson and hadron h have been denoted by q and p h , respectively, with q 2 = Q. In the center-of-mass energy frame where √ s = Q, the scaling variable can be written as z = 2E h / √ s. In this analysis, the analyzed data sets are based on SIA differential cross-sections for the unidentified charged hadron h = h + +h − . These data sets are differential with respect to the scaling variable z or p h . Actually, the format of the experimental data are different among the various experiments. In Table. I, the SIA cross sections included in SGK18 analysis have been listed for different experiments. The kinematical variables are as follows: scaling variable z = 2E h / √ s, the observed hadron h energy that scaled to the beam energy, and the hadron three-momentum p h . The scaled momentum x p is given by The relation between scaled momentum x p and z is defined as where ρ h = 4m 2 h /s and m h stands for the hadron mass. Note that ignoring the hadron mass leads to z = x p .
In Table I we have listed all analyzed flavor-untagged and tagged measurements used in our analysis which are reported by different experiments. These data sets include the ALEPH [66], OPAL [67,68] and DELPHI [69,70] experiments at CERN; the TPC [71] and SLD [72] experiments at SLAC; and TASSO [73] experiment at DESY.
As one can see from Table I, the measured observables are different for these data sets. Most of experimental collaborations have reported total inclusive and tagged cross sections, while the ALEPH, DELPHI and OPAL have reported longitudinal inclusive and bottom tagged cross section data. Separation of light and heavy quark flavor FFs is provided by the light and heavy flavor tagged experimental data. The longitudinal cross section data are proportional to the longitudinal structure function F L and implemented in SGK18 analysis to put further constraints on the gluon fragmentation function. The gluon coefficient functions were already available from several years ago at LO O(α s ). The NLO O(α 2 s ) coefficient functions have been also used in several analyses, for example in Refs. [74,75]. However, there is no analysis to determine the FFs of the unidentified charged hadrons including the coefficient functions at NNLO. As we mentioned the determination of unidentified charged hadrons at NNLO is the aim of the present paper.
Another point should be mentioned here is on the kinematic cuts applied on the data sets in SGK18 FFs analysis. We study the SIA data in potentially problematic low-z region, and hence, kinematic cuts are chosen consistently. To be on the safe side, we exclude the data points below the scaling variable of z min = 0.02 for the data sets at √ s = M Z , and z min = 0.075 for √ s < M Z . The data points with z max = 0.9 are not included in SGK18 QCD fit. The number of data points which are included in SGK18 fits are shown in the fifth column of Table I for each data sets separately. Moreover, the quality of our fits to SIA data for unidentified charged hadron at NLO and NNLO accuracy in term of the individual χ 2 values for every data set are also reported in the last two columns. The total χ 2 /d.o.f obtained from SGK18 best fits can also be found at the bottom of this table which are equal to 1.64 and 1.62 for NLO and NNLO analyses, respectively. Using the total 474 data points, we determine the 20 free parameters describing SGK18 unidentified charged hadron FFs D h i (z, Q 2 0 ). The details of SGK18 analysis on unidentified charged hadron FFs at NLO and NNLO will be discussed in details in Sec. IV.

III. THE QCD FRAMEWORK OF SGK18 FFS ANALYSIS
In the present SGK18 FFs analysis, we work in the well established pQCD framework for the electron-positron SIA process at the NLO and NNLO accuracy in pQCD. We make an extensive use of the x-space DGLAP evolution implemented in publicly available APFEL code [76] in which developed for a fast computation of the NLO and NNLO cross section of e − e + annihilation. For a clear review, we refer the reader to the Ref. [51,76,77] for further technical details of the QCD framework.
In this section, we review the factorization theorem of the cross-section and fragmentation structure functions in the electron-positron SIA process. We also discuss  the time-like DGLAP evolution of FFs. The differential cross-section for the single-inclusive e + e − annihilation involving a hadron h in the final state, with integrated over the production angle, and at a center-of-mass framework energy of s, is given by: where F h T (z, Q 2 ) and F h L (z, Q 2 ) are the transverse and longitudinal structure functions, respectively.
In the case of multiplicities, the total cross section for the electron positron annihilation into hadrons normalized to the differential cross section up to NNLO is written as: where the coefficients K (i) QCD relate to the QCD perturbative corrections that are currently known up to O(α 3 s ) [78]. Note that we have integrated over the scattering angle θ of the hadron h, and the cross section can be decomposed into transverse (T) and longitudinal (L) parts. Then F h T and F h L are called the time-like structure functions or fragmentation structure functions. The NNLO QCD corrections to the fragmentation structure functions can be expressed in factorized form of fragmentation functions D h i (z, Q 2 ) and calculable coefficient functions C S,NS k,l (z, α s (Q)) as follows The coefficient functions C S,NS k,l with k = T, L and l = q, g have been calculated in Refs. [74,75,79]. The factorization scale µ F and the renormalization scale µ R are set to be equal to the center-of-mass energy of the colli- q is the total cross section for quark production q at LO and σ (0) tot is the corresponding sum over all active flavors n f , σ In this equation, symbol ⊗ also denotes the standard convolution integral defined as The FFs, D h i (z, Q 2 ), which are non-perturbative but universal functions, parametrize the hadronization of massless partons, i = q,q, g, into the observed hadron h which carry fraction z of the hadron momentum. The scale dependence of the FFs which are governed by the renormalization equations are calculable in pQCD using the DGLAP evolution equation. The quark singlet (S) FF D h S (z, Q 2 ), non-singlet (NS) FFs D h NS (z, Q 2 ) as well as the gluon-to-hadron FF D h g (z, Q 2 ) are used in Eq. 6, and the singlet and non-singlet FFs are defined as: and The DGLAP evolution equations [45][46][47][48] evaluate the FFs with the energy scale Q 2 as where i, j = q,q, g and P ji are the time-like splitting functions [80][81][82]. According to the different FFs as nonsinglet, singlet and gluon FFs, one can rewrite Eq. 10 as a decoupled DGLAP equation for the non-singlet FFs and two coupled equations for the singlet and gluon FFs as The coefficient functions in Eq. 6 and the splitting functions in Eqs. 11 and 12 are defined as a perturbative expansion in powers of the α s , where i, j = q, g; k = T, L and a s = α s /(4π). In the MS scheme, the SIA coefficient functions have been computed up to NNLO for the C S,NS T,i . The longitudinal coefficient functions C S,NS L,i vanish at O(a 0 s ) and have been reported up to NLO accuracy in Refs. [74,75,79,83,84]. We should note here that, since C  [85].
In the next section, we briefly highlight the main feature of the SGK18 FFs analysis, specifically discussing SGK18 choice of parameterizations of the unidentified charged hadron FFs at the input scale and the heavy flavor mass scheme. The parameters describing the NLO and NNLO FFs also presented in the next section as well.

IV. OUTLINE OF THE SGK18 FFS ANALYSIS
In this section, we present the methodology of SGK18 FFs analysis, the input functional form and our assumptions we use in this analysis. As we mentioned, determination of individual fragmentation functions D h i for all quark flavors i as well as gluon into unidentified charged hadron at NLO and NNLO is the main aim of the present global analysis. We are also interested in studying the general features of NNLO corrections. As we discussed, the QCD framework for the NNLO corrections are only available for the single electron-positron annihilation among the hard scattering processes and only in the ZM-VFNS. We follow the same flexible functional form to parametrize the non-perturbative input FFs at initial scale Q 0 used in the series of DSS global QCD analyses [49,51,58,65]. In view of this fact, and in order to account the light quark decomposition q +q, we assume the following general initial functional form for SGK18 FFs analysis at a given input scale: where B[a, b] is the Euler Beta function which is used to normalize the parameter N i . We should notice here that the standard electronpositron annihilation data sets only provide information on the certain hadron spices summed over the charge, and hence, they are only sensitive to flavor combinations of q +q, i = u +ū, d +d, s +s, c +c, b +b and g. Since the observables for the unidentified charged hadron are usually presented for the sum dσ h = dσ h + + dσ h − , we only parameterize D h in our analysis. According to the charge conjugation D h + q(q) = D h − q(q) , we can separate quark and antiquark contributions as Since SIA data is sensitive to the D d+s , in SGK18 FFs analysis, we assume the symmetric fragmentation functions for d and s quark as D h d+d = D h s+s . Moreover, since these data sets are not sensitive to all parameters for the charm and bottom FFs, we assume γ c+c = γ b+b = 0 and δ c+c = δ b+b = 0. Hence, we choose the most simple functional form for the heavy charm and bottom FFs as follows, The parameter γ g for the gluon FFs is basically unconstrained by the analyzed data sets, and in order to get the best fit, we decided to keep it fixed at γ g = 70 for both SGK18 NLO and NNLO analyses. We discuss in section VI that the gluon FF obtained in SGK18 analysis is slightly different from the DSS07 analysis which used the SIDIS and hadron collider data. The protonantiproton data from CDF [86,87] experiment at SLAC, the proton-proton data from CMS [88,89] and ALICE [90] experiments at CERN carry a large amount of information on the gluon FF and could constrain it well enough. However, the data from single-inclusive charged hadron production in e − e + annihilation is the major source of exprimental data in our analysis.
We should mention here that in SGK18 FFs analyses, the initial scale for input parametrization is Q 0 = 5 GeV for all parton species. Since the value of bottom mass in our analysis is m b = 4.3 GeV, this initial scale is above bottom threshold. In addition, this value for Q 0 is below the lowest center-of-mass energy of analyzed data sets, s = √ 14 GeV. Since time-like matching conditions are unknown at NNLO, with this value for Q 0 , it is not require heavy quark threshold as well as the matching in the evolution between the initial scale and the data scale. Therefore, in our analysis the number of active flavor keep fixed to the n f = 5.

V. χ 2 MINIMIZATION AND CALCULATION METHOD OF ERRORS
The parameters describing the unidentified charge hadron FFs presented in Eqs. (14) and (16) are determined using a standard χ 2 minimization method. The total χ 2 is calculated in comparison with the single-inclusive charged hadron production data sets in electron-positron annihilation for the unidentified charge hadron FFs. In order to calculate the χ 2 , the theoretical predictions should be obtained at the same experimental z and µ 2 = Q 2 points. As we mentioned, the µ 2 = Q 2 evolution is calculated by the well-known DGLAP evolution equations.
In order to calculate the total χ 2 ({η i }) for independent sets of fit parameters {η i }, one can use the following standard χ 2 definition: where E i is the measured value of a given observable and T i is the corresponding theoretical estimate for a given set of parameters {η i } at the same experimental z and µ 2 = Q 2 points. The experimental errors associated with this measurements are calculated from systematic and statistical errors added in quadrature, (δE i ) 2 = (δE sys i ) 2 + (δE stat i ) 2 . The optimization is done by the CERN program MINUIT [91].
Since most single-inclusive charged hadron production data in e − e + annihilation come with additional information on the fully correlated normalization uncertainty, the above simple χ 2 definition need to be modified in order to account for such normalization uncertainties. Hence, the modified function is given by, where n exp corresponds to the individual experimental data sets for the n th experiment, and N data n refers to the number of data points in each data set. The normalization factors ∆N n in above equation can be fitted along with the fitted parameters ({η i }) of Eqs. (14) and (16) and then keep fixed. In order to illustrate the effects arising from the use of the different single-inclusive charged hadron production data sets, in Table. I, we have shown the obtained χ/n data for each data sets at NLO and NNLO accuracy. This table illustrates the quality of SGK18 NLO and NNLO QCD fits to single-inclusive charged hadron production data in terms of the individual χ 2 -values obtained for each experiment. The total χ 2 /N pts for the SGK18 fits can be found in this table as well. We obtained 1.64 and 1.62 for our NLO and NNLO analyses, respectively.
This section also focuses on the uncertainties of the parameters in Eqs. (14) and (16) to judge the quality of SGK18 QCD fits. In order to determine the uncertainties of unidentified charged hadron FFs as well as the corresponding observable, we apply the "Hessian" method by choosing a particular value of ∆χ 2 = 1. This will provide a clear and comprehensive picture of the uncertainty characteristic of resulting FFs.
The determination of the size of uncertainties using the "Hessian" method is based on the correspondence between the confidence level (C.L.) P and χ 2 with the number of fitting parameters N . The C.L. is given by, where Γ is the Gamma function. The value of ∆χ 2 in Eq. (19) is taken so that the C.L. becomes the one-σerror range, namely P = 0.68. The value for the ∆χ 2 is then numerically calculated by using this equation.
Having at hand the value for ∆χ 2 and the derivatives of given observables with respect to the fitted parameters {η i } (i=1, 2, ..., N ), the Hessian approach provides the uncertainties of desired observables O as, where C m,n is the inverse of the Hessian matrix which can be obtained by running the CERN program library MINUIT [91].
For estimation of uncertainties at an arbitrary Q 2 which is an attributive function of the input parameters, the obtained gradient terms are evolved by the well-known DGLAP evolution kernel. In next section we show that the SGK18 FFs uncertainties determination as well as the fitting methodology can correctly propagate the experimental uncertainty of the single-inclusive charged hadron production data into the uncertainties of the SGK18 FFs.

VI.
SGK18 FIT RESULTS In this section, we present the SGK18 numerical results for the unidentified charged hadron D h FFs obtained from the global analysis of SIA data. Firstly, we present the parameters of the optimum QCD fits describing the unidentified charged hadron and then we present the SGK18 FFs results for different partons at NLO and NNLO accuracy in pQCD. Then, SGK18 results for FFs are compared to the DSS07 FFs for unidentified charged hadron. Secondly, the uncertainty bands at NLO and NNLO accuracy are compared and the improvement of the FFs calculations due to the inclusion of NNLO QCD corrections are discussed. Next, SGK18 theoretical predictions for the total cross sections and all different tagged cross sections are compared with the analyzed SIA experimental data sets. Finally, we present our theoretical uncertainties from the variation of the renormalization and factorization scales.

A. SGM18 FFs and comparison with DSS FF sets
Even though we mainly interested in a precise extraction of unidentified charged hadron FFs D h at NNLO accuracy, we also present the results of our analysis at NLO approximation. As we will discuss in this section, the significantly better NNLO uncertainty highlights the importance of higher order correction in our QCD analysis. In addition our NLO results can be used in calculation of observable which are limited to the NLO corrections.
The 21 best fit parameters describing the optimum NLO and NNLO unidentified charged hadron D h FFs are given in Tables II and III. The SGK18 D h fragmentation functions and their uncertainties have been presented in Fig. 1 at NLO (solid lines) and NNLO (dashed lines) accuracy for Q 2 = 25 GeV 2 . The resulting NLO and NNLO SGK18 D h fragmentation functions and their uncertainties evolved to the scale of 100 GeV 2 and M 2 Z have also been illustrated in Figs.2 and 3. We have included the one-σ uncertainty bands in our analysis.
In these figures, we have also compared SGK18 FFs to the central value of DSS07 FFs analysis [65] (dot-dashed lines) at NLO. The uncertainty bands of the DSS07 results are not shown in our analysis because they are not available. Recently, a preliminary determination of the unidentified charged hadron FFs has been reported in Ref. [64] at NLO but the grid files are not available. The impact of higher order QCD corrections on the reduction of FFs uncertainties at NNLO accuracy in comparison with the NLO analysis can be seen from these figures. As one can see from Figs. 1, 2 and 3, the uncertainty bands for all quark flavors as well as gluon decrease remarkably at NNLO which indicates that the higher order perturbation QCD corrections increase the precision of the calculations.
Let's discuss the results in more details. Focusing on Fig. 1, one can clearly see that the results of SGK18 for light and heavy flavor FFs are similar to the ones obtained from DSS07 analysis, especially at larger values of the momentum fraction z. However, some noticeable differences can be seen in the gluon FFs. As one can see, the contribution of gluon FF in our analysis is significantly larger than the DSS07 one. A main reason for this difference is base on the fact that in SGK18 analysis only SIA data have been included, while the DSS07 included both the electron-positron SIA and proton-proton collision data in their analysis. Since the collider data could directly effect the determination of gluon FF, this clear difference observed between our gluon FF result and DSS07 one is expectable. In addition, we should mentioned this fact that the initial scales used in these two models are different. The SGK18 initial scale has been chosen to be Q 2 0 = 25 GeV 2 , while the DSS07 initial scale is Q 2 0 = 1 GeV 2 for light quarks and gluon, and Q 2 0 = m 2 c for charm quark FF as well as Q 2 0 = m 2 b for the bottom quark FF.
In addition to the points mentioned above, we have excluded the experimental data below z min = 0.02 for data at center-of-mass energy of √ s = M Z and z min = 0.075 for √ s < M Z , while the DSS07 excluded the data points below z = 0.1 for all analyzed data sets. Consequently, the number of SIA data points in DSS07 is 236, but in our analysis there are 474 data points. Because of this difference for the kinematic cuts at small z, the most discrepancy is seen at z < 0.1. In general, the differences become larger towards smaller values of z, z → 0.01, which has been already observed for all parton species.
The SGK18 FFs at higher values of Q 2 = 100 GeV 2 and Q 2 = M 2 Z , have been displayed in Figs. 2 and 3, respectively. It can be conclude from these two figures that as the scale of energy increases, the difference between the SGK18 gluon FF and DSS07 is decreased. Furthermore, the behavior of our light and heavy FFs and the DSS07 ones are slightly in good agreement.
Let us turn to the discussion on the resulting FFs and their uncertainties by focusing on the inclusion of higher order QCD corrections. According to the SGK18 global fit of SIA data, there are some noticeable features that improve global fit at NNLO in comparison to the NLO. First, as one can see from Table. I, the improvement of the χ 2 /d.o.f from NLO to NNLO is slightly better. Actually, in our analysis the value of χ 2 /d.o.f reduces from 1.64 at NLO to 1.62 at NNLO approximation. Second, the size of the SGK18 FFs uncertainties remarkably decrease at NNLO in comparison to our NLO analysis. One can see the reduction of uncertainties for all determined quark flavors as well as gluon at all three scales of energy shown in Figs. 1, 2 and 3. One can conclude that the effects arising due to the inclusion of higher order QCD corrections significantly decrease the obtained error bands.

B.
Discussion of SGK18 fit quality and data/theory comparison In this section, we present the SGK18 NLO and NNLO theoretical predictions for the total and tagged SIA cross sections. We also compare in details our results with all single-inclusive D h charged hadron production in e − e + annihilation data analyzed in this study. The values of the χ 2 for both individual and total data sets included in SGK18 analysis have been reported in Table I  NLO and NNLO accuracy. As one can see, the χ 2 for inclusive and bottom tagged longitudinal data sets remarkably decrease at NNLO approximation.
In order to judge the quality of the fits of SGK18 FFs analysis, we compare the experimental data to their corresponding NLO and NNLO theoretical predictions calculated using the NLO and NNLO FFs obtained from the QCD fits. Fig. 4 shows a comparison between the normalized total cross sections from the ALEPH, DELPHI, OPAL and SLD measurements of unidentified charged hadron and our NLO and NNLO predictions. The uncertainty bands of the predictions, due to the one-σ FF uncertainties, have also been shown in this figure. Moreover, the same comparisons have been performed for light (DELPHI, OPAL and SLD); charm (OPAL and SLD); and bottom (DELPHI, OPAL and SLD) tagged cross sections. Finally, we have shown the same comparison for the inclusive and b-tagged longitudinal cross sections from ALEPH, DELPHI and OPAL data sets in Fig. 5.
Overall, the results obtained demonstrate a good agreement between the SGK18 theoretical predictions and analyzed experimental data. Considering the exclusion of small z data points, the SGK18 results are also in reasonable agreement with data in the small and large z regions for all data sets. In Figs. 4, the NLO and NNLO predictions are in a satisfactory agreement in comparison to the total inclusive, light, charm and bottom tagged data for all range of z. , OPAL [67,68], DELPHI [69,70] and SLD [72] at the scale of Q = MZ. The shaded bands refer to our uncertainty results at NLO (green band) and NNLO (yellow band) and shaded areas indicate the kinematic regions excluded by our cuts.
In order to investigate the fit quality of the total datasets at NLO and NNLO, as a next step, we discuss the size of uncertainty bands at NLO and NNLO. As one can see, the NLO uncertainty bands are slightly larger than NNLO one as presented in Fig. 4. As it is seen from Fig. 5, the SGK18 NNLO theoretical predictions show more consistency with the data in comparison to the NLO ones for the inclusive and b-tagged longitudinal cross sections. The NLO theoretical predictions tend to be overshoot by ALEPH, DELPHI and OPAL longitudinal experimental data for z < 0.1 and DELPHI b-tagged one for z < 0.2. The data sets for the longitudinal inclusive and tagged cross sections have important effect on the determination of gluon FF, because they are nonvanishing already at LO O(α s ) contribution. According to the absence of precise data for wider range of Q 2 , the longitudinal data could help to constrain the gluon FF.
We should notice here that, in spite of the exclusion of small z < 0.02 data points, our NNLO theory predictions are in good agreement with the excluded region for ctagged in Fig. 4 and all other data sets in Fig. 5.

C.
The improvement of NNLO accuracy in theoretical uncertainty The fragmentation function uncertainties have different sources that classify into the experimental data errors and the theoretical uncertainties caused from the phenomenological assumptions in any global QCD fits. The possible sources of theoretical uncertainties may include, for example, higher order correction effects in calculation of cross sections, the phenomenological form for the FFs parametrizations at an arbitrary initial scale and differ- In this section, we present our results for studying the residual dependence on the choice of scale of energy µ. The most important source of theoretical uncertainties is dependence on the choice of the scale of energy µ. We expect to shrink progressively when we include higher and higher order corrections. It is exactly what we find in our study.
In Fig. 6, the best fits of the NLO and NNLO analyses for unidentified charged hadron have been used to demonstrate the residual theoretical uncertainties due to the variations of the renormalization (µ R ) and factorization (µ F ) scales. According to this figure, the SIA cross section depends on the scale of energy and the results have been shown at NLO and NNLO accuracies (shaded bands) for µ R = µ F = µ = Q/2 and µ = 2Q which are normalized to our default choice of µ = Q. It is obvious that the theoretical calculations depend on the scale of µ. Note that our results have been presented for three scales of energy, Q = 10 GeV, Q = 30 GeV and Q = M Z . According to the results presented in Fig. 6, one can clearly conclude that the NNLO predictions are more stable than the NLO ones and come with much smaller theoretical uncertainties.

VII. SUMMARY AND CONCLUSIONS
In this paper, a new determination of unidentified charged hadron FFs at NLO and for the first time at NNLO accuracy in perturbative QCD are presented. The flavor-untagged and the tagged SIA data in e − e + annihilation are included in this analysis that are reported by CERN (ALEPH [66], OPAL [67,68], and DELPHI [69,70]), SLAC (TPC [71] and SLD [72]) and DESY (TASSO [73]). The heavy flavor contributions are considered in the ZM-VFNS in z-space in the framework of publicly available APFEL code.
We illustrate the quality of the SGK18 FFs at NLO and NNLO and show that the results presented in this analysis are in good agreement with the results in literature and all exprimental data analyzed in this study. We have presented the uncertainties for the D h fragmentation functions and the corresponding theory predictions using the "Hessian" approach.
The most striking remarkable improvements to emerge from SGK18 FFs analysis are as follows: As a first improvement, this study is the first step towards enhancing our understanding of parton-to-unidentified charged hadrons FFs by analyzing flavor-untagged and the tagged SIA data considering the NNLO accuracy in perturbative QCD. As a second improvement, we use smoother kinematical cut for the small z regions than other analyses in the literature such as DSS07. Consequently, SGK18 FFs analysis uses a wider range of exprimental data points in the fitting procedure.
As a third improvement, we have presented the perturbative stable QCD fits and observed a reduction of uncertainties for our FFs as well as theory predictions at NNLO with respect to NLO. Finally, we have chosen Q 0 = 5 GeV as an initial scale in SGK18 fits and then the number of active flavor is always fixed to n f = 5. This choice improve the fit because time-like matching conditions are unknown at NNLO. Within this choice of initial scale, the heavy quark threshold as well as the matching condition don't need to be taken into account.
As a final improvement, by using our fit result at NNLO for FFs, the agreement between our predictions for the inclusive and b-tagged longitudinal cross sections have improved in comparison with the NLO analysis which suggest that the inclusion of higher order corrections could improve the fit quality. We hope that our research will serve as a base for future studies on the determination of unidentified charged hadrons FFs from wide range of exprimental observables at CERN, HERA and SLAC. However further works need to be carried out to establish a framework to consider the SIDIS and hadron collider data into the analysis.
A FORTRAN package, which evaluates the SGK18 NLO and NNLO unidentified charged hadron FFs as well as the theory predictions presented in this study can be obtained from the authors upon request via electronic mail.