First determination of $D^{*+}$-meson fragmentation functions and their uncertainties at next-to-next-to-leading order

We present, for the first time, a set of next-to-next-to-leading order (NNLO) fragmentation functions (FFs) describing the production of charmed-meson $D^{*+}$ from partons. Exploiting the universality and scaling violations of FFs, we extract the NLO and NNLO FFs through a global fit to all relevant data sets from single-inclusive $e^+e^-$ annihilation. The uncertainties for the resulting FFs as well as the corresponding observables are estimated using the Hessian approach. We evaluate the quality of the {\tt SKM18} FFs determined in this analysis by comparing with the recent results in literature and show how they describe the available data for single-inclusive $D^{*+}$-meson production in electron-positron annihilation. As a practical application, we apply the extracted FFs to make our theoretical predictions for the scaled-energy distributions of $D^{*+}$-mesons inclusively produced in top quark decays. We explore the implications of {\tt SKM18} for LHC phenomenology and show that our findings of this study can be introduced as a channel to indirect search for top-quark properties.

VI. Discussion of SKM18 fit results 6 A. SKM18 FFs and comparison with other FF sets 7 B. Discussion of fit quality and data/theory comparison 8 C. Importance of theoretical uncertainty at NNLO 9 VII. Energy spectrum of inclusive D * ± -mesons in polarized and unpolarized top decays 9 VIII. Summary and Conclusions 11
Studies over the past two decades have provided valuable important information on the structure of hadrons. The FFs and PDFs extracted from these studies encode the long-distance dynamics of the interactions among quarks and gluons so the partonic cross sections encode the short-distance dynamics of the interactions. The process-independent FFs, D H i (z, µ 2 F ), describe the probability for a parton i at the factorization scale µ F to fragment into a hadron H carrying away a fraction z of its momentum. The specific importance of FFs is for their model-independent predictions of the cross sections (or decay rates) at colliders where the observables involving identified hadrons are detected in the final state [3,10,11]. The knowledge of FFs has received quite some interest in, for example, tests of quantum chromodynamics (QCD) such as theoretical calculations for recent measurements of inclusive production in protonproton collisions at the Relativistic Heavy Ion Collider (RHIC) and the CERN LHC, and in investigating the origin of the proton spin to describe the internal structure of nucleons.
In the naive parton model, the nonperturbative FFs are independent of the factorization scale µ F , but in the QCD-improved parton model their scaling violations are subject to the perturbatively computable Dokshitzer-Gribov-Lipatov-Alteralli-Parisi (DGLAP) evolution equations [12]. In this respect, the PDFs and FFs are on the same footing. Like PDFs, the FFs must be extracted from experimental data through global QCD analyses, possibly from a variety of hardscattering processes [1][2][3][13][14][15]. The FFs are included in hadrproduction processes such as electron-positron Single-Inclusive Annihilation (SIA), lepton-hadron Semi-Inclusive Deep-Inelastic Scattering (SIDIS) and hadronhadron scattering processes. Information obtained from SIDIS multiplicities and from hadron-hadron collisions is particularly useful in order to achieve a complete flavor decomposition into quark and antiquark FFs along with a direct determination of the gluon FF. Among them, SIA remains the theoretically cleanest process to extract the fragmentation densities since its interpretation does not require the simultaneous knowledge of PDFs [11]. Recent progresses in the extraction of FFs have been focused on light charged mesons (pions and kaons), for which data are more numerous as they dominate the identified hadron yields. For example, in Ref. [16] we have determined the π ± and K ± FFs, both at leadingorder (LO) and next-to-leading order (NLO) accuracy in perturbative QCD through a global QCD fit to the SIA data and the SIDIS asymmetry data from HERMES and COMPASS. There, we have broken the symmetry assumption between the quark and antiquark FFs for favored partons by using the asymmetry data. In [17], we have also presented sets of proton FFs with their corresponding uncertainties through a global fit to all relevant data sets of SIA, taking the finite-mass effects of proton into account. Other well-known groups which have determined sets of NLO FFs for these two mesonic species are: AKK [18], LSS [19], DSEHS [20], DSS [21,22], HKNS [23] and JAM [24] collaborations who have used different phenomenological models and variety of experimental data. In these studies, main focus was put on quantifying the effects of the inclusion of new measurements on the FFs, although in the JAM and HKNS fits these were limited to SIA data. Apart form the different data set used in their determinations, these collaborations have also introduced some methodological and theoretical improvements over previous determinations, specifically, in order to achieve a more reliable estimate of the uncertainties of FFs. For example, the iterative Hessian procedure [25] has been used in the DSEHS analysis, while the iterative Monte Carlo approach [26] has been developed in the JAM analysis. Recently, in Refs. [1,27] using data from a comprehensive set of SIA measurements, authors have presented a new determination of the FFs of charged pions, charged kaons, and protons/antiprotons at nextto-next-to-leading order (NNLO) in perturbative QCD. They resulted that the systematic inclusion of higherorder QCD corrections significantly improves the description of the data, especially in the small-z region.
Generally, the theoretical treatment of heavy quarks provides a unique laboratory to test perturbative QCD. In fact, heavy flavor cross sections which have been mea-sured both at very high energies at the LHC and at various low energy experiments poses unique challenges to our deep understanding of QCD. Specifically, charm production cross sections are applied, for instance, to further constrain the gluon PDF at small-x [28], and they play a vital role in neutrino astrophysics and cosmic-ray [29,30]. Another importance of charm production cross sections concerns the modification of heavy flavor yields in heavyion collisions where highly energetic partons can traverse the quark-gluon plasma thereby attaining valuable information about the properties of the QCD medium. For more detail, see Ref. [31].
In [32], using the SIA data from the Belle, CLEO, ALEPH and OPAL Collaborations, authors have determined nonperturbative charmed-meson FFs at NLO in the modified Minimal-Subtraction (M S) factorization scheme. To extract the results, the General-Mass Variable-Flavor-Number scheme (GM-VFNs) was adopted. However, the most recent Belle data which were included in their analysis have been removed due to an unrecoverable error in the measurement. Recently, in [11] authors have calculated the FFs of charged D * -meson at next-to-leading order accuracy using the available data for single-inclusive D * -meson production in e + e − annihilation, hadron-hadron collisions, and in proton-proton scattering. While quark-to-hadron FFs can be relatively well constrained from SIA data alone, the gluon FF can be well constrained via proton-proton scattering data. It should be noted that, to compute the NLO FFs authors have used the Zero-Mass Variable-Flavor-Number scheme (ZM-VFNs) [33] where the heavy quark masses are set to zero in all partonic cross sections so the heavy quarks are treated as the active, massless partons in the proton and the nonzero values of the charm-and bottomquark masses only enter through the initial conditions of the FFs. This scheme is applicable for sufficiently high energies and transverse momenta.
In the present work, using the ZM-VFN scheme we focus on the hadronization of charm-and bottom-quarks into D * ± -mesons, which are of particular relevance in the era of the LHC. We will provide the first global QCD analysis of (c, b) → D * ± FFs at NNLO accuracy through a global QCD fit to SIA data from ALEPH [34] and OPAL [35] collaborations at LEP. Our analysis which is named as "SKM18" from now on, is limited to SIA data only due to the lack of other single-inclusive particle production cross sections at NLO and NNLO accuracy. We will also present an attempt to estimate the uncertainties of the extracted FFs as well as the resulting normalized total cross sections, for which we adopt the Hessian method [36][37][38].
As it is well-known, in the Standard Model of particle physics the top quark has a short lifetime so it decays before hadronization takes place. In fact, at the lowest order the top quark decays as t → W + b followed by b → X + Jets, where X refers to the detected hadrons in the final state. Therefore, at the LHC the study of energy spectrum of produced hadrons through top de-cays is proposed as a new channel to indirect search for the top quark properties. As an example of a possible application, the extracted FFs in our analysis are used to make the theoretical predictions for the energy distributions of D * ± -mesons in polarized and unpolarized top decays. We shall compare our results with the ones obtained through the NLO FFs extracted in [32].
This paper is organized as follows. In Sec. II we explain the QCD analysis of hadronization process in electronpositron annihilation by introducing FFs. We describe our formalism and parametrization form at the initial scale for D * ± FFs in Sec. III. In Sec. IV, we discuss the related experimental data and their effects in our global analysis. In following, the minimization method and the Hessian uncertainty approach to calculate the errors in SKM18 analysis are present in Sec. V. In Sec. VI, the full results for the D * ± -FFs and their uncertainties are listed. We also present a comparison of SKM18 results with experimental data and the other models in this section. The theoretical uncertainties due to the variation of the renormalization and factorization scales at the NLO and NNLO accuracies are studied at the end of this section. Our predictions for the energy spectrum of charmed mesons produced in top quark decays are presented in Sec. VII. Our result and conclusion are summarized in Sec. VIII.

II. QCD ANALYSIS FRAMEWORK UP TO NNLO ACCURACY
The optimal way to determine the D * -FFs is to fit them to experimental date extracted from the singleinclusive e + e − annihilation processes where X stands for the unobserved final jets. Since, this has less contributions from background processes in comparison with the hadron collisions, then one does not need to deal with the uncertainties introduced by PDFs. In the following, we explain how to evaluate the cross section of the above process in the parton model of QCD within the ZM-VFN scheme. Note that, in the ZM-VFNS (or massless scheme), the number of active quark flavors also increases with the flavor thresholds. If we denote the four-momenta of the exchanged virtual gauge boson and the D * meson by q and p D , so that s = q 2 and p 2 D = m 2 D , the scaling variable x D is defined as x D = 2(p D · q)/q 2 , similarly to the Bjorken-x variable in deep inelastic scattering. In the center-of-mass frame, this variable is reduced to x D = 2E D / √ s which refers to the energy of D * -meson in units of the beam energy.
According to the factorization theorem of the improved QCD-parton model [39], the differential cross section of process (1) can be written as a convolutions of perturbatively calculable partonic cross sections i stands for one of the fragmenting partons i = g, u,ū, · · · , b,b. This convolution reads The renormalization (µ R ) and factorization (µ F ) scales are arbitrary quantities, however, a choice often made consists of setting µ 2 = µ 2 F = Q 2 and following Ref. [11] we shall adopt this convention in this work. The variable x i is also defined in analogy to x D as x i = 2(p i · q)/q 2 , where p i is the four-momentum of parton i. To compare our theoretical calculations with experimental data, we normalized the differential cross section dσ/dx D to the total cross section up to NNLO for e + e − annihilation into hadrons (σ tot ) which reads where α and α s are the electromagnetic and the strong coupling constants, respectively, and the coefficients K  [40]. It should be noted that, to calculate the NLO or NNLO FFs one needs to have the partonic cross sections dσ i /dx i (Wilson coefficients) at each desired order. The NNLO Wilson coefficients are analytically presented in Refs. [41][42][43].

III. PHENOMENOLOGICAL PARAMETRIZATION UP TO NNLO
In order to choose the best parametrization for SKM18 global analysis of D * ± -FFs, we tested different models such as Peterson [44,45] and Bowler [32] and other simple forms applied for light hadron FFs (i.e. pion, kaon and proton) [23]. Finally, we adopted the Bowler parametrization form because of its best fit to low number of data for D * ± -meson. Therefore, we parametrize the z distributions of the c(c) and b(b) quark FFs at their starting scales µ 0 as suggested by Bowler, i.e.
with three free parameters; N, α and β. Our fitting procedure is going as follows. At the scale µ 0 , the charm and bottom quark FFs are taken to be of the Bowler form as in Eq. (4), while the FFs of gluon and light quarks (q = u, d, s) are set to zero, i.e.
Then, these light and gluon FFs are evolved to higher scales using the DGLAP evolution equations [12] at NLO or NNLO. To proceed, we set the initial parametrization scale as µ 2 0 = 18.5 GeV 2 which is a little grater than the b-quark threshold, i.e. Q 2 = m 2 b = (4.3) 2 GeV 2 . The advantage of taking this value for the initial scale is due to the fact that the time-like matching conditions are currently known only up to NLO accuracy and with this input scale the heavy-quark thresholds should not be crossed in the QCD evolution.
Since all hadrons in electron-positron annihilation originate from the produced qq pair, the multiplicities for D * + and D * − are the same and the experimental data for charged hadrons are usually presented for their sum, i.e.
According to the parton structure of D * − , the FFs of D * − can be obtained as and for the gluon FF, it reads It seems that one should determine six free parameters for the FFs of charm and bottom quarks into the D * ± -meson using experimental data. Since, N c can not be constrained by data, then one should fix it from the beginning so the actual parameters which should be determined are five parameters. In the next section, we will discuss how further data collections are required to determine all six parameters described above.

IV. EXPERIMENTAL INPUT
Most of experimental data for D * ± in e + e − annihilation is reported by ALEPH [34], OPAL [35], CLEO [46] and Belle [47] Collaborations. An overview of the data included in SKM18 global analysis of FFs is presented in Tables. I and II for the total, c-tagged and b-tagged SIA cross sections from ALEPH [34] and OPAL [35] Collaborations. For each dataset we include the corresponding published reference, the number of data points in the SKM18 NLO and NNLO D * ± -FFs determinations, and also the center-of-mass energy. The total number of data points for the SKM18 FFs determination is 59 at NLO and NNLO. ALEPH and OPAL Collaborations at LEP present their experimental data at Q = M Z , the mass of the Z boson, while Belle and CLEO provide their data in lower energy, i.e. Q = 10.5 GeV. In this range of energy, all D * ± in e + e − annihilation coming from bottom decays are excluded because they are bellow the b-quark mass threshold.
In spite that the Belle data has been published in Ref.
[48], but results are not yet available on the HEP-DATA web page. While KKKS08 group [32] have used the Belle data in their analysis of FFs, on the HEP-DATA web page it is explained that "due to an unrecoverable error in the measurement, the data for this record have been removed by the request of the authors" [48]. Therefore, we are not able to include this data set while KKKS08 group applied them. This makes the KKKS08 analysis unreliable.
Note that, CLEO present the data as distributions dσ/dx p for hadron H where the scaled momentum is defined as where ρ H = 4m 2 H /s in which m H stands for the hadron mass and the allowed values of x p are 0 ≤ x p ≤ 1. Also the conversion formula for differential cross section reads . Ignoring the hadron mass leads to (dσ/dx p )(x p ) = (dσ/dx)(x). This assumption leads to a large difference between theory and experimental results. When we apply the CLEO data in our analysis, it disorganizes our results and also the χ 2 increases. While this data set has been included in the extractions of D * -FFs in [32,49], it has not been included in recent analysis by AKSRV17 group [11] and authors found noticeable tensions between the CLEO and ALEPH data. In conclusion, we do not include the CLEO low-energy SIA data in SMK18 analysis.
The ALEPH and OPAL experimental data sets [34,35] present the total, charm and bottom flavor tagged for the cross section distributions normalized to the total hadronic cross section and include the branching fractions of the decays used to identify the D * ± mesons, namely D * + → D 0 π + followed by D 0 → K − π + . Therefore, we divide these two data sets by the respective decay branching fractions given by [50], i.e. B(D * + → D 0 π + ) = (67.7 ± 0.5)% and B(D 0 → K − π + ) = (3.93 ± 0.04)%. Unfortunately, the c-tagged experimental data from ALEPH Collaboration are presented in graphical form, and we can not reach the numerical values [34]. The normalized c-tagged cross sections from OPAL experiment [35] are included in the SKM18 NLO and NNLO QCD fits.
AKSRV17 group [11] include electron-positron annihilation, hadron-hadron collisions and hadron-in-jet data in their global analysis of D * -FFs at NLO. Since, using different hadronization processes increase the number of data, the c-tagged cross section can be further constrained. Also, the AKSRV17 analysis allows for a nonzero gluon FF at the initial scale in order to get a good global fit of all different hadronization processes data. Note that the Wilson coefficient functions are only known at NNLO just for SIA process [41][42][43]. Since our global analysis is performed at NNLO, which is the main purpose of this study, we limit the potential of global determination of FFs to e + e − annihilation.
In our analysis, we apply the massless scheme where we ignore the finite-mass effects of heavy quarks as well as the charmed-meson mass effects. The cuts applied to SIA data are as follows. We exclude some of the SIA data in our global analysis. The LEP data (ALEPH and OPAL) taken at Q = M Z are used only in the interval 0.1 < z < 0.95.

V. THE CALCULATION METHOD OF ERRORS
In SKM18 global analysis, the total χ 2 is calculated in comparison with the experimental data for D * ± production in e + e − annihilation. For calculating the total χ 2 , the theoretical functions should be obtained at the same experimental z and µ 2 points. As we explained, the µ 2 evolution is calculated by the well-known DGLAP evolution equations [12] and the total cross section for e + e − annihilation into hadrons (σ tot ) is obtained from Eq. (3). The simplest and fastest method to calculate the total χ 2 ({η i }) for independent sets of unknown fit parameters {η i } reads where O data i stand for the experimental values related to the desired observables and T theory i are the corresponding theoretical values of D * ± productions at the same experimental z and µ 2 points. The experimental errors σ data i are calculated from statistical and systematic errors added in quadrature, (σ data ) 2 = (σ sys ) 2 + (σ stat ) 2 . The unknown parameters of FFs, D D * ± i (z, µ 2 0 ), presented in Eq. (4) are determined so as to obtain the minimum χ 2 . The optimization of a given function is done by the CERN program MINUIT [51].
In order to illustrate the effects arising from the use of D * ± production data sets in our analyses, in Tables. I and  II, we show the χ 2 for each data sets. These tables clearly illustrate the quality of our QCD fits at NLO and NNLO accuracy in terms of the individual χ 2 -values obtained for each experiment. The total χ 2 /d.o.f for the resulting fit to the D * ± productions are 1.31 and 1.27 for our NLO and NNLO fits, respectively.
Since most of experiments usually come with an additional information on the fully correlated normalization uncertainty, ∆N n , the simple χ 2 definition presented in Eq. (11) needs to be corrected in order to account for the corresponding normalization uncertainties. In that case and in order to determine the best fit parameters, we need to minimize the χ 2 global ({η i }) function in respect to the free unknown parameters. Considering the discussion mentioned above, the χ 2 global ({η i }) function can be written as, where W n is a weight factor for the n th experiment and in which n exp refers to a given individual experimental data set and N data n corresponds to the number of data points in each data set.
The normalization factors ∆N n in Eq. (13) can be fitted along with the fitted parameters ({η i }) of Eq. (4) and then can keep fixed. The obtained normalization factors ∆N n are also presented in Tables. I and II for our NLO and NNLO analyses.
The determination of nonperturbative FFs through QCD fits to the data is a statistical procedure that necessarily implies a variety of assumptions. The most important one is the input parameterization functions of the charmed-meson FFs and the propagation of the experimental uncertainties into them [1][2][3][52][53][54]. The assessment of uncertainties of PDFs and the corresponding observables has seen significant progress in recent QCD analyses [6,7,27,55]. Among the different approaches in literature, the "Hessian method" [36,37,56], the Lagrange multiplier (LM) technique [38] and newly Neural Network (NN) in which the NNPDF uses [57][58][59][60] are the most reliable ones.
The well-known and practical method is "Hessian method", which has been widely used to extract the uncertainties of the PDFs, polarized PDFs and nuclear PDFs as well as the corresponding observables in our previous analyses. [61][62][63][64][65][66][67][68]. Although, the technical details of the Hessian method are described in the mentioned references, a short summary of the method is explained here.
As we already discussed, our method is the χ 2 n ({η i }) fitting procedure in the global QCD analyses. For the determination of FFs uncertainties, we have used the well-known "Hessian" or error matrix approach. This approach confirms that the SKM18 fitting methodology used in this QCD analysis can faithfully reproduce the input charmed-meson FFs in the region where the D * + productions data sets are sufficiently constraining. The fit parameters are denoted as η i (i = 1, 2, ..., N ), where N refers to the total number of the fitted parameters. One can expand the χ 2 around the minimum χ 2 point, i.e.η, as where, H ij are the elements of the Hessian matrix which can be written as The confidence region is normally can be given in the parameter space by supplying a value of ∆χ 2 . In a standard parameter-fitting criterion, the errors are given by the choice of the tolerance T = ∆χ 2 = 1 in Eq. (14). For the number of one fitted parameter (N = 1), the confidence level (C.L.) is 68% for ∆χ 2 = 1. For the general cases and for more number of fitted parameters (N > 1), the ∆χ 2 should be calculated to determine the size of the resulting uncertainties. The correct determination of the tolerance T should indicates that SKM18 fitting methodology as well as the uncertainties determination does correctly propagate the experimental uncertainty of D * + production data into the uncertainties of fitted charmed-meson FFs.
The determination of the size of uncertainties can be done by applying the "Hessian method" based on the correspondence between the number of fitting parameters with the C.L. P and χ 2 . The C.L. P can be written as follows where Γ is the well-known Gamma function. The value of ∆χ 2 is taken so that the C.L. becomes the one-σerror range, namely P = 0.68. Similarly, for the 90 th percentile, one can use P = 0.90. Using the equation above, the value of ∆χ 2 is numerically calculated.
The Hessian matrix (or error matrix) is accessible by running the subroutine the CERN program library MINUIT [51]. The uncertainty on a given observable O, which is an attributive function of the input parameters of charmed-meson at the input scale µ 2 0 , is obtained by applying the "Hessian method". Having at hand the value of ∆χ 2 , and derivatives of the observables with respect to the charmed-meson fitted parameters, the "Hessian method" gives the uncertainties of a given observable O. It reads, where C j,k is the inverse of the Hessian matrix, C j,k = H −1 jk . For estimation of uncertainties at an arbitrary scale µ 2 , the obtained gradient terms are evolved by the DGLAP evolution kernel [12], and then the charmedmeson FF uncertainties as well as the uncertainties of other observables, such as total cross section for e + e − annihilation into hadrons σ tot , are calculated.
In the next section, we present the main results of this work, namely the "SKM18" set of charmed-meson FFs at NLO and NNLO approximations. First, we discuss the resultant D * ± -FFs and their uncertainties. Then, we show the quality of the fits and compare the SKM18 theoretical predictions to the D * ± -meson production data sets.

VI. DISCUSSION OF SKM18 FIT RESULTS
Now we turn to SKM18 numerical results for the global analysis of D * ± -FFs from SIA data. First, we present the SKM18 NLO and NNLO FFs for the charm and bottom densities in ZM-VFN scheme and investigate the dependence of obtained FFs on various data sets. Second, we quantify the size of uncertainty bands due to the NNLO corrections. Next, SKM18 theoretical predictions will be compared to all SIA data and other well-known phenomenological groups in literature. We finalize this section with comments on the impact of collider data on determination of gluon FFs and the reduction of FFs uncertainties. Finally, in a detailed study we will consider the theoretical uncertainties due to the variation of the renormalization and factorization scales. Our study shows that the NNLO theoretical predictions will indeed be more stable than the NLO ones.
As was mentioned, the SKM18 includes for the fist time the NNLO QCD corrections of SIA cross section in order to achieve a high percentage accuracy. After extracting the NLO and NNLO D * ± -FFs, we will compare them with the NLO ones reported by KKKS08 [32] and the recent ones from AKSRV17 [11]. In addition, we will compare our theoretical predictions for the normalized cross sections with the analyzed experimental data. We use the public code APFEL [69] to compute the NLO and NNLO DGLAP evolution of the FFs and the SIA cross sections in the MS factorization scheme in z-space. In this package the numerical solution for the time-like evolution equations are provided up to NNLO accuracy in QCD. Some of recent analyses on the FFs of light hadrons (pion, kaon and proton) have also been done by using this public code [1,70].
Our results of the optimum fits for the charm and bottom FF parameters are presented in Tables. III and IV for the NLO and NNLO fits, respectively. As was already mentioned, the parameter N c is basically unconstrained by data and comes with a relatively large error. For this reason, we calculate it from the first fit along with six free parameters and then we fix it in the second fit with five free parameters. In Tables. I, II the normalization shifts N i for each data set and χ 2 values are reported. We include 59 data points from single-inclusive hadron production in our global analysis of D * ± -FFs. To see how much the NNLO approximation of D * ± -FFs improves the NLO ones, note that, for the NNLO global fit the value of χ 2 is 1.27 which is better than the NLO one, χ 2 = 1.31. Normally, the quality of fit can be judged from the obtained fitted parameters and the χ 2 values. According to the best fit parameters presented in Tables. III and IV, one can conclude that our parameters are well determined.

A. SKM18 FFs and comparison with other FF sets
The SKM18 D * ± fragmentation densities and their uncertainties are presented in Figs. 1, 2 and 3. The results of this analysis include the one-σ uncertainty bands. They are compared to the central value from the KKKS08 analysis [32] and very recent analysis of AKSRV17 [11]. Now, we take this opportunity to discuss the SKM18 D * ± -FFs in more detail. In Fig. 1, the c and b FFs including their uncertainty bands are shown at NLO (solid lines) and NNLO (dashed lines) at the initial scale µ 2 0 , in which we consider µ 2 0 = 18.5 GeV 2 both for the charm and bottom quark FFs. A sightly small difference is evident between the NLO and NNLO results. A basic question is now that, what impact the higher order QCD corrections can have on the reduction of FFs uncertainties, specifically the NNLO corrections in comparison with the NLO ones? To this end, we show in Fig. 1 the SKM18 NNLO FFs with their error bands. From this figure one can conclude that the NLO and NNLO uncertainties are slightly similar in size.
The resulting SKM18 FFs are depicted in Fig. 2 for higher values µ 2 > µ 2 0 . The NLO (solid lines) and NNLO (dashed lines) results have been shown for the c-and bquarks and also for the gluon at µ 2 = 100 GeV 2 . In this figure, we compared our results with the NLO ones from KKKS08 [32] (dot-dashed lines) [71] and with very recent fit of AKSRV17 (short dashed lines) [11]. In comparison to the AKSRV17 analysis, the SKM18 charm-quark FF is similar to the one from AKSRV17 while the contribution of gluon FF extracted in AKSRV17 is significantly larger than the SKM18 and KKKS08 one. A main reason for this difference is due to the fact that the AKSRV17 allows for a nonzero gluon FF at the initial scale due to inclusion of hadron-hadron collision data in their fit. We will show that a similar behavior is observed for the charm and gluon FFs at µ 2 = M 2 Z . Obviously, the SKM18 bquark FF behaves similar to the one from the AKSRV17 and KKKS08 analysis, although, the result of KKKS08 for the c-quark FF is not compatible with the one from SKM18 and AKSRV17.
In Fig. 3, we repeat the same study as before so the c, b and gluon fragmentation densities are shown at µ 2 = M 2 Z at NLO (solid lines) and NNLO (dashed lines). They are also compared with the KKKS08 (dot-dashed lines) and the AKSRV17 analyses (short dashed lines). As one can conclude from Fig. 3, the most difference between our results and the KKKS08 and the AKSRV17 analysis is in the gluon FF, as was already discussed its reason. While the bottom FF in the analysis of SKM18, AKSRV17 and KKKS08 are in agreement, the charm one behaves differently. One of the differences between the SKM18 and AKSRV17 analysis and the KKKS08 analysis is the charm tagged cross section D * ± from OPAL collaboration which has not been included in the KKKS08 fit. The charm FF is constrained by the charm-tagged data which is included both in the SKM18 and AKSRV17 analysis.
Let us now discuss on the resulting FFs and their uncertainties by focusing on their perturbative convergence as well as the effects arising due to the inclusion of higher-order QCD corrections. Concerning the SKM18 fit quality of the total dataset, the most noticeable feature is the improvement upon inclusion of higher-order corrections. Although, the improvement of the χ 2 /n.d.f is rather marginal when going from NLO to NNLO. This finding demonstrates that the inclusion of the NNLO QCD corrections slightly improves the description of the data. A further noticeable aspect of the comparisons in Figs. 1, 2 and 3 is related to the size of the SKM18 FFs uncertainties which show that the NLO and NNLO uncertainties are slightly similar in the size. One can also conclude from the presented plots that the differences between the NLO and NNLO FFs are slightly small. We should stress here that, from Ref. [1] the same pattern for the NNPDF FFs can be seen at NLO and NNLO. As the NNPDF collaboration are mentioned in their recent paper, while in some cases the LO and NLO distributions can differ by more than one standard deviation, the differences between the NLO and NNLO FFs are relatively    Fig. 4 we present the NNLO/NLO ratio for the FF of c, b and gluon at µ 2 = M 2 Z . From Fig. 4 one can now judge the reduction of error band widths as well as the size of improvements upon inclusion of higher-order corrections.
As was mentioned in Sec. III, at the initial scale µ 0 the FFs of gluon and light quarks (q = u, d, s) are set to zero and they are evolved to higher scales using the DGLAP equations. The light quark contributions to produce the D * ± -mesons are very small in comparison with the heavy quarks so one expects them to have main contributions in production of light hadrons such as pion, kaon and proton. For this reason, in Figs. 2 and 3 we restricted our results to the charm, bottom and gluon FFs. As can be seen from these figures, the differences between the SKM18, AKSRV17 and KKKS08 analysis are considerable for the charm and gluon FFs compared to the bottom quark one. Generally, this is most likely due to the different experimental data which are included in these analysis. Individually, we list some possible explanations for this finding. Firstly, as we mentioned in Sec. IV, the KKKS08 fit [32] have included both the CLEO and withdrawn Belle data at the low energy Q = 10.5 GeV. Since, these data are bellow the bottom threshold, they affect on the charm and anti-charm FFs, directly. On the other hand, because of using CLEO and Belle data, the KKKS08 analysis have included the mass effects of D * ± -meson and heavy quarks into their analysis, according to Eqs. (9) and (10), which might be the source of some differences. Secondly, charm-tagged cross section D * ± from OPAL collaboration which has not been included in the KKKS08 fit are used both in the analysis of SKM18 and AKSRV17. Finally, the AKSRV17 collaboration have used the data from hadron-hadron collisions and jet fragmentation in the proton-proton scattering in addition to the electron-positron data. The collider data is directly sensitive to the gluon FFs. The collider data is taken into account in the AKSRV17 analysis so could carry a large amount of information on the gluon FFs, and may also account for the differences observed between these results. The b-FFs plotted in Figs. 2 and 3 are quite similar to the ones obtained by the KKKS08 and AKSRV17 analysis, because in both fits the bottomtagged data from ALEPH and OPAL [34,35] are included and the bottom FFs are significantly constrained by these data sets.
It should be pointed out that, the authors of AK-SRV17 [11] analysis claimed who applied the online FFs generator [71] to extract the KKKS08 FFs. We are particularly puzzled since there is inconsistency between their KKKS08 curves and the results which can be obtained by the FF generator. We do not judge in this respect.

B.
Discussion of fit quality and data/theory comparison After our detailed discussion on the SKM18 FFs in comparison to the results in literature, we now turn to present our theoretical predictions for the SIA cross sections. In Figs. 5 and 6, the theoretical calculations based on the SKM18 global QCD analysis are compared to the analyzed available data. The ALEPH [34] and OPAL [35] normalized total cross sections at µ 2 = M 2 Z are shown in Fig. 5 along with our calculations for the SIA cross section at NLO (solid line) and NNLO (dashed line) accuracies. According to these results we find that our NLO and NNLO QCD fits are in good agreement with the experimental data and our theoretical predictions describe the data well for the intermediate and large values of z.
Towards small-z region, our model fall down the data because the evolution equation become unstable in smallz region. The splitting functions lead to negative FFs for z << 1 in evolution equations, and additionally mass corrections play important role in this region. So we exclude regions where mass corrections and the singu- lar small-z behavior of the splitting functions are effective. As we mentioned in Sec. IV, we exclude the regions with z < 0.1 and z > 0.95 in our global fit. Also we compare our theoretical results with the KKKS08 model in Figs. 5 and 6. The AKSRV17 theoretical predictions are not shown because they are not available to the authors. According to [32], the KKKS08 analyses are done at both ZM-VFN and GM-VFN schemes. On the other hand, they calculate the free parameters from Belle/CLEO (at 10.5 GeV), ALEPH/OPAL (at M Z ) and global data fits separately. Since our analyses base on the ZM-VFN scheme and we only include the ALEPH and OPAL data, then we compare our results with the corresponding KKKS08 set of FFs in the ZM-VFNS which are calculated with ALEPH/OPAL data fit. As one can see in Figs. 5 and 6, the KKKS08 set is not compatible with data at the small z region. In Fig. 6 the ALEPH [34] and OPAL [35] normalized charm and bottom tagged cross sections at µ 2 = M 2 Z are shown and also our fit results are presented for this cross section at NLO (solid line) and NNLO (dashed line). The most visible conclusion is the very considerable agreements of SKM18 with the analyzed data. From Fig. 5 one can conclude that the KKKS08 overestimate the data for the small values of z. This pattern is also seen in Fig. 6. In the next section, we aim to present our predictions for the energy spectrum of D * ± -mesons produced in polarized and unpolarized top quark decays.

C. Importance of theoretical uncertainty at NNLO
The possible sources of uncertainties on the FFs classify into the experimental errors on the data, and the theoretical or phenomenological assumptions in the global QCD fit. The theoretical uncertainties include, for exam-ple, higher order QCD effects in the calculation of cross sections, the assumption of flavor or charge conjugation symmetries, the parametrization forms of the FFs at an arbitrary initial scale, and so on.
One source of uncertainties on the theoretical results is the values selected for the renormalization (µ R ) and factorization (µ F ) scales, so the former associated with the renormalization of the strong coupling constant. In principle, one can use two different values for these scales, however, a choice often made consists of setting µ F = µ R = µ and we adopted this convention for the results presented.
In this section, we will examine the importance of these scales at the NLO and NNLO approximations and the variation effects of these scales on the theoretical results. The scale dependence of FFs is shown in Fig. 7, where the results of NLO and NNLO accuracy (shaded bands) are shown for Q/2 ≤ µ ≤ 2Q. As is seen, the theoretical calculations depend on the choice of the scale µ and it is, indeed, one of the most important source of the theoretical uncertainties. We expect that this kind of theoretical uncertainties shrinks progressively when higher order corrections are included. Therefore, the NNLO predictions will indeed be more stable than the NLO ones. This is exactly what one can conclude from the SKM18 results presented in Fig.7.

VII. ENERGY SPECTRUM OF INCLUSIVE D * ± -MESONS IN POLARIZED AND UNPOLARIZED TOP DECAYS
In this section, we turn to apply the extracted D * ± -FFs to make our phenomenological predictions for the energy spectrum of charmed mesons produced in top quark Our results are also compared with the KKKS08 (dot-dashed lines) [32] and the AKSRV17 (short dashed lines) [11] results at NLO. decays through the following process where X stands for the unobserved final state. The study of energy distributions of produced mesons through top quark decays might be introduced as an indirect channel to search for the top quark properties. In this channel, both the b-quark and the gluon may hadronize into the charmed-meson so the gluon contributes to the real radiations at NLO and higher orders. Note that, the contribution of gluon FF can not be discriminated so that this contribution is being calculated to see where it contributes to the energy spectrum of produced meson. Therefore, this part of calculation is of more theoretical relevance, however in the scaled-energy of hadrons, as an experimental quantity, all contributions including the b-quark and gluon contribute. Ignoring the b-quark mass and by working in the ZM-VFN scheme, to obtain the energy distribution of D * ±meson in the process (18), we employ the factorization theorem (2) as where, following Ref. [72], we define the scaled-energy fraction of D * ± as x D = 2E D /(m 2 t − m 2 W ) and dΓ/dx i are the parton-level differential decay rates of the process t → i + W + (i = b, g). The analytical expressions for the differential decay widths dΓ/dx i at the QCD NLO approximation are presented in Ref. [72] for the unpolarized top quark decays. The NLO differential decay rates for the polarized top quarks are analytically calculated in [73,74]. Although, in the above relation the factorization (µ F ) and renormalization (µ R ) scales are arbitrary but we here set them to µ R = µ F = m t = 172. FFs of (b, g) → D * + extracted in Ref. [32], we also studied the same quantity (dashed line) in Fig. 8. As is seen, in comparison with the KKKS08 result, our FFs show a reduction for the energy spectrum at the peak position. In Fig. 9, we have presented our NLO predictions for the energy distribution of D * ± -mesons produced through the polarized (dashed line) and unpolarized (dots) top decays.
To have a detailed insight for the size of NNLO corrections, in this plot we have also studied the same distribution at NNLO for the polarized top decays (solid line) and unpolarized ones (dot-dashed). However, these NNLO predictions are not completely correct because for a reliable prediction at this order of perturbative QCD one also needs to have the NNLO parton-level differential decay rates dΓ/dx i to convolute with the NNLO FFs, see Eq. (19). These quantities have not been yet calculated. The study of the x D distribution (i.e. dΓ/dx D ) of the dominant decay mode t → D * ± W + + X at the LHC, as a formidable top factory, will enable us to deepen our understanding of the nonperturbative aspects of D * ±meson formation by hadronization and to pin down the b → D * ± and g → D * ± FFs. By measuring the x D distribution of polarized top decays, the b/g → D * ± FFs can be constrained event further.

VIII. SUMMARY AND CONCLUSIONS
Let us now come to our summary and conclusions. We have determined the nonperturbative FFs of partons into the D * ± -meson at NLO perturbative QCD and, for the first time, at NNLO one from global analyses of singleinclusive electron-positron annihilation. Our analyses are based on the ZM-VFN scheme in which all quarks are treated as massless partons. Our phenomenological analyses (called as SKM18 analyses) to achieve the FFs of D * ± -meson is significant in, at least, two major respects. Firstly, we applied all SIA experimental data as much as possible including most of the data from ALEPH and OPAL Collaborations. Secondly, we considered the NNLO accuracy in our global fit using the public APFEL code. According to the χ 2 values extracted from our optimum fits, increasing the accuracies of theoretical QCD analysis up to NNLO slightly decreases the value of χ 2 .
In addition, we found that the experimental uncertainties for the D * ± -FFs and SIA cross sections are similar in size both for the NLO and NNLO approximations. We have also studied the variation effects of the renormalization and factorization scales considering Q/2 ≤ µ ≤ 2Q where we set µ r = µ f = µ. The obtained results show that our calculations at NNLO come with much smaller theoretical uncertainties relative to the NLO calculations which reflects the stability of our NNLO analysis. These findings are significantly in agreements with previous results reported in the literature.
As is well-known in global QCD analysis of FFs at NLO, some phenomenological collaborations have included single-inclusive D * ± -meson production in electron-positron annihilation, hadron-hadron collisions, and in-jet fragmentation in proton-proton scattering. Unfortunately, the NNLO expressions of Wilson coefficients for the processes of hadron-hadron collision and proton-proton scattering are not yet available. Therefore, we restricted our analyses to the SIA experimental data. Despite to a few number of experimental data available for D * ± -meson, the SKM18 theoretical predictions are in reasonable agreements with the experimental data.
We also applied the SKM18 FFs to make our predictions for the energy spectrum of D * ± -mesons produced in polarized and unpolarized top quark decays. At the LHC, this study can be considered as a channel to indirect search for the top properties and, specifically, it enables us to deepen our understanding of the nonperturbative aspects of D * ± -meson formation and to pin down the (b, g) → D * ± FFs.
As was mentioned, including the mass effects of heavy quarks and meson in the theoretical QCD analysis leads to the significant information about our understanding of FFs, specially at the small-z region. The analyses including the meson and quark masses are reserved for future work. More detail on these effects will be given in our upcoming paper. Although the current study is based on a few sample of datasets, the findings are in good agreements with all experimental data as well as the results in literature. We hope that our research will serve as a base for future studies on D * ± -meson FFs. However, we should emphasize that a number of future  [34] and OPAL [35] are also shown in this scale. The shaded bands refer to our uncertainty results at NLO (green band) and NNLO (yellow band).
studies using the experimental data from colliders are strongly recommended.  [34] and OPAL [35] are also shown in this scale. The shaded bands refer to our uncertainty results at NLO (green band) and NNLO (yellow band).  FIG. 8: xD distribution of dΓ NLO /dxD for the polarized top quark decay at µF = mt. We also show the same distribution using the FFs presented by KKKS08 [32] at NLO. We also show the same results at NNLO.