NNLO charmed-meson fragmentation functions and their uncertainties in the presence of meson mass corrections

The main aim of this paper is to present new sets of non-perturbative fragmentation functions (FFs) for $D^0$ and $D^+$ mesons at next-to-leading (NLO) and, for the first time, at next-to-next-to-leading order (NNLO) in the $\overline{\mathrm{MS}}$ factorization scheme with five massless quark flavors. This new determination of FFs is based on the QCD fit to the {\tt OPAL} experimental data for hadron production in the electron-positron single-inclusive annihilation (SIA). We discuss in detail the novel aspects of the methodology used in our analysis and the validity of obtained FFs by comparing with previous works in literature which have been carried out up to NLO accuracy. We will also incorporate the effect of charmed meson mass corrections into our QCD analysis and discuss the improvements upon inclusion of these effects. The uncertainties in the extracted FFs as well as in the corresponding observables are estimated using the"Hessian"approach. For a typical application, we use our new FFs to make theoretical predictions for the energy distributions of charmed mesons inclusively produced through the decay of unpolarized top quarks, to be measured at the CERN LHC. As a result of this analysis, suggestions are discussed for possible future studies on the current topic to consider any theory improvements and other available experimental observables.

Previous measurements of charmed mesons production at high energy scattering experiments [1] have shown that inclusive heavy quarks production provides an interesting test for the strong interaction dynamics. Historically, the first results for the charm quark photoproduction cross section σ(γp → cc + X) were presented by ZEUS [2] and H1 [3] Collaborations at ep collider HERA. With the advent of high statistics data from the CERN Large Hadron Collider (LHC), one can expect that it allows to test the predictions of QCD with much higher precision, however, this is not only important for testing QCD itself; but a better understanding of heavy quarks productions is also vital for searches of new physics (NP) phenomena and beyond standard model (BSM) researches [4,5].
Over the last few years, several experimental collaborations at pp and ep colliders have presented data on the differential cross section d 2 σ/(dydp T ) for the inclusive production of D 0 , D + and D + s mesons [1,6,7], where y and p T refer to the rapidity and transverse momentum of mesons, respectively. On the theoretical side, fragmentation functions (FFs) describing the transition probability of the initial partons into observed hadrons, are needed as nonperturbative inputs for the calculation of all the relevant cross sections measurements [8][9][10][11][12][13]. In the framework of QCD, FFs encode the long-distance dynamics of the interactions among quarks and gluons which lead to their hadronization in a hard-scattering process. Generally speaking, to obtain theoretical predictions for all observables involving identified hadrons in the final state, FFs have to be convoluted with partonic cross-sections encoding the short-distance dynamics of hadron production process.
In Ref. [14], non-perturbative D 0 , D + and D + s FFs are determined both at leading order (LO) and next-toleading order (NLO) in the modified minimal subtraction (MS) factorization scheme by fitting the fractional energy spectra of these mesons measured by the OPAL Collaboration [15,16] in electron-positron annihilation on the Z-boson resonance at the CERN LEP1 collider. Authors in Ref. [14] have applied the massless scheme or zero-mass variable-flavor-number scheme (ZM-VFNS) where heavy quarks are treated as any other massless parton. This scheme can be applied to the open production of heavy flavors, such as D and B hadrons, provided the hard energy scale characteristic for the production process is sufficiently larger than the heavy-flavors mass. This is certainly the case for our applications, because m D ≪ m Z , m t . In Ref. [17], the authors updated their previous results [14] using newer data on charmed meson production presented by the ALEPH [18], BELLE [19] and the CLEO [20] Collaborations. Compared to the data from ALEPH and OPAL, the data from BELLE and CLEO are located much closer to the thresholds √ s = 2m c and √ s = 2m b of the transitions c → H c and b → H b , where H c and H b stand for generic c or b hadrons, respectively. However, at the present, the BELLE data included in their analyses have been removed due to an unrecoverable error in the measurements.
In the present analysis, using the ZM-VFN scheme we focus on the hadronization of charm-and bottomquarks into D 0 and D + mesons and provide the first QCD analysis of (c, b) → (D 0 , D + ) FFs at NNLO accuracy in perturbative QCD. We perform a QCD fit to the electron-positron single inclusive annihilation (SIA) data measured by OPAL Collaboration. It should be mentioned here that, our analysis is restricted to the SIA data only due to the lack of other theoretical partonic cross sections for single inclusive production of partons at NNLO accuracy. We shall also go beyond Refs. [14,17] by performing a full-fledged error estimation, both for the FFs and the resulting normalized total cross sections, employing the Hessian approach. Furthermore, we include the effects of meson mass on the FFs, which modify the relations between partonic and hadronic variables and reduce the available phase space. Hadron masses are responsible for the low-energy thresholds.
As an application of our analysis, we apply the extracted FFs to make our theoretical predictions for the energy distribution of charmed mesons produced through the top quark decays. In the SM of particle physics, the top quark is the heaviest elementary particle so that its large mass does not allow it to make a constrained state. Therefore, the top quark decays before hadronization takes place, i.e. t → bW + → HW + + Jets, where H refers to the observed colorless hadron in the final state. Therefore, at the CERN LHC, the study of energy distribution of observed hadrons can be considered as a channel to indirect search for detailed study of top quark properties. The study of this energy spectrum can be even considered as a new window towards searches on new physics. In fact, any considerable deviation of energy distribution of mesons from the theoretical predictions in the SM theory can be related to the new physics. In this work, using the obtained results for the charmed meson FFs, we make our predictions for the energy distributions of charmed mesons produced through the decay of unpolarized top quarks at NLO and NNLO accuracies.
The outline of this paper is as follows. In Section II, we shortly describe our theoretical framework to compute charmed meson FFs. In Section III, we describe how one should incorporate the effects of meson mass into the formalism. Section IV includes a detailed description of experimental data sets used in this analysis. The results of this study are presented in Sec. V. The extracted D 0 and D + FFs are discussed in this section and comparisons with other results in literature will be also presented. This section does also include comparisons of analyzed SIA data sets with our NLO and NNLO theoretical predictions for the normalized cross sections along with the one-sigma error bands. In Sec. VI, we will present an example for the application of our analysis. Sec. VII includes the summary and conclusions.

II. THEORETICAL FRAMEWORK
In the present work using the ZM-VFN scheme, where the heavy quarks are treated as other massless partons, we mainly focus on the hadronization of charm-and bottom-quarks into charmed mesons which are of particular relevance in the era of the LHC [21][22][23][24][25]. In particular, we perform a new QCD analysis of OPAL experimental data [15] for D + and D 0 mesons up to NNLO accuracy. As usual the optimal way to determine the FFs is to fit them to experimental date extracted from the single-inclusive electron-positron annihilation processes via a virtual photon (γ) or Z boson, i.e.
where X stands for the residual final state including unobservable jet and H c denotes the charmed meson; D 0 and D + in our analysis. In comparison with the hadron collisions, for example ep or pp, the SIA process has less contributions from background processes so one does not need to deal with the uncertainties due to the parton distribution functions (PDFs). Hence, SIA would provide the cleanest process for the determinations of FFs rather than the ep of pp processes. The differential cross section of hadron production can be calculated theoretically in the frame work of QCD through the factorization theorem [26]. According to this theorem, the differential cross section of e + e − annihilation can be written as a convolution of hard scattering subprocess, i.e. e + e − → iī, and the non-perturbative meson FFs describing i/ī → D, i.e.
Denoting the four-momenta of the virtual gauge boson and the D meson by q and p D , respectively, so that s = q 2 and p 2 D = m 2 D , in (2) we introduced the scaling variable x D = 2(p D · q)/q 2 . In the center-of-mass (c.m.) frame where q = ( √ s, 0) this variable is simplified to as x D = 2E D / √ s where E D refers to the energy of D-meson. Thus, x D represents the energy of Dmeson in units of the beam energy. In Eq. (2), µ F and µ R are the factorization and renormalization scales, respectively. The scale µ F is associated to the divergences due to collinear gluon radiation off a massless primary quark or antiquark. To avoid these divergences which are proportional to ln(s/µ 2 F ) we adopt the usual convention µ F = µ R = √ s [27]. In Eq. (2), y = 2(p i · q)/q 2 where p i denotes the four-momentum of parton i so by working in the c.m. frame it is simplified as y = 2E i / √ s. In fact, y refers to the fraction of the beam energy passed on to the parton i. In Eq. (2), the Wilson coefficients dσ i /dy are the differential cross sections at the parton level which can be calculated perturbatively. Nowadays, these coefficients are known up to NNLO approximation in the perturbative QCD theory, see for example Refs. [28][29][30]. Since, most of experimental data are presented in the form of 1/σ tot ×dσ/dx D , therefore, to be able to compare our theoretical results with experimental data we need to normalize Eq. (2) to the total cross section which reads whereẽ i is the effective electroweak charge of quark i, α and α s are the electromagnetic and the strong coupling constants, respectively. The coefficients K  (1) which are determined phenomenologically. In our analysis, to perform the best parametrization we adopt the optimal functional form suggested by Bowler [32] for the parametrization of c and b quark FFs. It is given by where i = c, c, b, b. This parameterization form includes three free parameters α i , β i and N i which should be determined phenomenologically from the QCD fit to the SIA data sets. In the proposed form, µ 0 is the initial fragmentation scale which is taken to be µ 2 0 = 18.5 GeV 2 for heavy flavors b and c. The advantage of choosing this amount for the initial scale µ 0 , which is a little larger than m 2 b , is due to the fact that the time-like matching conditions are currently known up to NLO accuracy and with this input scale the heavy quark thresholds should not be crossed in the QCD evolution. The FFs of gluon and light flavors (u, d, s) are taken to be zero at the initial scale and are generated at higher scales through the Dokshitzer-Gribov-Lipatov-Alteralli-Parisi (DGLAP) evolution equations [33]. For the value of strong coupling constant, we adopt α 118 at the Z boson mass scale [34]. In this analysis, we employ the publicly available APFEL package [35] in order to the evolution of zD D 0 (z, µ 2 ) and zD D + (z, µ 2 ) FFs as well as for the calculation of the SIA cross sections up to NNLO accuracy. To determine the free parameters, we applied the CERN program MINUIT [36] to minimize the χ 2 global which is defined in our previous work in detail [37,38]. This quantity includes the overall normalization errors of the D-meson production data sets. The uncertainties are also estimated using the standard Hessian error propagation, as outlined in Refs. [39,40], with ∆χ 2 = 1, which corresponds to a 68% confidence level (CL) in the ideal Gaussian statistics. For more detail, see Ref. [41].

III. MASS CORRECTIONS
The factorization theorem described in Eq. (1) is routinely used in the literature for the inclusive production of single hadrons. In proving this formula it is assumed that quarks and hadrons are massless so the non-zero values of the c-and b-quark masses only enter through the initial conditions of the FFs, and the mass of the heavy baryon sets the lower bounds on the scaling variable x H , i.e. x min H = 2m H / √ s. In Ref. [42], we have proved how to enter the hadron mass effects into the inclusive hadron production in e + e − reaction. To incorporate the hadron mass effects we used a specific choice of scaling variables by working in the light-cone (L.C) coordinates. Ignoring the detail of calculations, as a generalization of the massless case, the differential cross section in the presence of hadron mass m H reads The above formula is a fundamental relation for the factorization theorem extended in the presence of hadron mass and would be more effective and applicable when the data are presented in lower energy scales. Among all well-known collaborations, the only collaboration who studied the effects of hadron mass into their calculations is the AKK collaboration [43,44]. Although, due to the data applied in their work, they have identified the scaling variable x p = 2| p|/ √ s, where | p| stands for the threemomentum of produced hadron, while in our analysis we use the energy scaling variable as x H = 2E H / √ s. Therefore, the differential cross section or, equivalently, the extended factorization formula including the hadron mass effects computed in our work is different from the one (dσ/dx p ) presented in Refs. [43,44].
It should be noted that the effect of hadron mass is to increases the cross section dσ/dx H at small x H so this treatment acts inverse for large x H . We incorporated the mass corrections into the publicly available APFEL package [35]. Our results show that the inclusion of hadron mass leads to a change in the value of χ 2 /(d.o.f) up to about 5%. Although, due to the scale of energy given for the data used in our work, i.e. √ s = 91.2 GeV, this improvement is small but having the data at lower energies it is expected to have more effect on the χ 2 /(d.o.f).

IV. EXPERIMENTAL DATA SELECTION
After introducing the theoretical framework of our analysis in the previous sections, in the following we discuss the experimental data which are used in this analysis. The first experimental information corresponding to the D meson production through e − e + annihilation comes from the measurements performed by the OPAL Collaboration [15]. The OPAL results have been reported for D 0 and D + production at the Z-boson resonance. In this process, two mechanisms contribute with similar rate; Z → cc decay followed by c/c → H c fragmentation and Z → bb decay followed by b/b → H b fragmentation and weak decay of the bottom-flavored hadron H b into the charmed-meson via H b → H c + X. Note that, the energy spectrum of H c hadrons originating from decays of H b hadrons is much softer than that due to primary charm production. For separating charmed hadron production through Z → cc decay from Z → bb decay, OPAL used the apparent decay length distributions and energy spectra of the charmed hadrons. The OPAL Collaboration has presented x D -distributions for D 0 and D + samples and for their z → bb subsamples (b-tagged events). They are displayed in the form 1/N had × dN/dx D where N is the number of charmed meson candidates reconstructed through appropriate decay chains. To convert these data into the desired cross section 1/σ tot ×dσ/dx D , one should divide them by the convenient branching fractions of decays used in Refs. [15,16] for the reconstruction of the various charmed mesons, i.e.
These data are displayed in Figs. 3 (for the D 0 meson) and 4 (for the D + meson) and are listed in Tables. I and II along with the corresponding values of χ 2 . The number of data points and the data properties as well as the center-of-mass energy ( √ s) are presented in these tables. We have also presented the total χ 2 divided by the number of degrees of freedom χ 2 /d.o.f for each analysis. These values show the quality of our well-satisfying fit. Another important measurement on the charmed mesons production has been done by CLEO Collaboration [20] in the e − e + annihilation. Excluding the decay products of B mesons, their results are presented for the scaled momentum spectra dσ/dx p at √ s = 10.5 GeV where x p = | − → p |/| − → p | max . As was mentioned, the OPAL Collaboration [15,16] has presented the cross section as a function of the scaled energy x D . The relation between these two variables reads [17] x where ρ H = 4m 2 H /s. Since the variable x D ranges as √ ρ H < x D < 1, consequently x p takes values from 0 to 1. For differential cross section the conversion formula is Although the CLEO data can provide useful information on the D 0 and D + FFs, but there is a problem in using them in the analysis. In fact, as was mentioned, we have chosen the ZM-VFN scheme in our analysis which is reliable only in the region of large transverse momenta or equivalently in high energy scales. On the other hand, the data from CLEO are located much closer to the thresholds √ s = 2m c and √ s = 2m b of the transitions c → H c and b → H b , than those from OPAL. Hence, including the CLEO data sets in the analysis might be a reason for tension. To check this point, we converted the CLEO data to the desired form dσ/dx D using Eqs. (8) and (9) and included them into our analysis. We observed that inclusion of these data increased the value of χ 2 , significantly. Accordingly, we restricted ourselves to the OPAL data sets for extracting the charmed mesons FFs and neglected the CLEO data sets. We remind that our choice of the massless scheme was due to the lack of massive partonic cross sections at NNLO, currently. In next section, we shall present the main findings of our study and discuss on the results. We shall compare our zD D 0 (z, µ 2 ) and zD D + (z, µ 2 ) FFs with other results in the literature. A detailed comparison of our NLO and NNLO theoretical predictions for the SIA cross sections with the analyzed data sets will be also presented.

V. RESULTS AND DISCUSSION
After introducing all ingredients needed for analysis, we are now in a situation to perform a fit with experimental data sets for each hadron species (D 0 and D + ). In this section, at first, the numerical results obtained for the (c, b) → D 0 /D + FFs through the perturbative QCD analysis are reported at the initial scale µ 0 . Next, as an example, we present our results for the FFs along with their uncertainty bands for the gluon, charm and bottom quarks up to NLO and NNLO corrections at µ 2 = 100 GeV 2 . Our results for the FFs are compared to the KKKS08 FFs [17] for D mesons. In addition, our theoretical predictions as well as the error bands for total and b-tagged cross sections are compared with the analyzed SIA experimental data sets.  The best values of the fit parameters for the FFs of charm and bottom quarks into the D 0 and D + mesons are listed in Tables. III and IV for the NLO and NNLO accuracies at the initial scale µ 0 . Since the time-like matching conditions are not known for NNLO correction we chose the input scale as µ 0 = 18.5 GeV 2 [13,45] which is above the b-quark mass threshold. According to the Table. I, the value of χ 2 /d.o.f for the D 0 FFs at NLO and NNLO accuracies are 1.149 and 1.05, respectively. Hence, one can conclude that the inclusion of higher order QCD corrections can improve the fit quality. The same conclusion does also hold for the D + FFs analysis, see Table. II. In fact, a comparison between the values of χ 2 NNLO/d.o.f = 0.69 and χ 2 NLO/d.o.f = 0.75 is a reason for the goodness of fit at NNLO analysis. In Fig. 1, we plotted the FFs of D 0 and D + at NLO and NNLO accuracies at higher energy scale µ 2 = 100 GeV 2 . To do this, we studied the D 0 fragmentation from the c +c and b +b quarks as well as the gluon contribution for both NLO (solid lines) and NNLO (dashed lines) accuracies. For a quantitative comparison, our results are compared with the ones obtained through using the fit parameters presented by the KKKS08 Collaboration (dotdashed lines) [17]. As is seen, our result for the transition of b +b into the D 0 meson is in good agreement with the corresponding result obtained from the KKKS08 analysis. For the FF of c +c, our results are different from the KKKS08 analysis, specifically, in the large z. For the case of gluon FF zD D 0 g (z, µ 2 ), a significant difference is seen for a wide range of z. The error bands in Fig. 1 correspond to the choice of one-unit tolerance ∆χ 2 = 1 obtained using the Hessian approach. The green and yellow bands show the NLO and NNLO error uncertainties, respectively. For the differences between our results and the ones from the KKKS08 analysis, it should be noted that the data sets included in our fit and the KKKS08 analysis are different. In fact, the KKKS08 Collaboration has included the data from BELLE in their analyses as well, which is now removed due to an unrecoverable error in the measurements. Moreover, CLEO data has been also included in the KKKS08 analysis. Since, both the CLEO and Belle dataset have been measured at low energy then they cause a slightly harder FF, as is expected. According to the results presented in tables 1 and 2 of Ref. [17], these two datasets impact the c-quark FF parameters more than the b-quark one. On the other hand, the initial scales of energy are not the same in both analyses so that the KKKS08 has set the initial scale as µ 2 0 = m 2 c while ours is µ 2 0 = 18.5 GeV 2 which is a little larger than m 2 b . In Fig. 2, the same study is done for the D + at the scale µ 2 = 100 GeV 2 for both NLO (solid lines) and NNLO (dashed lines) analyses. They are also compared with the results obtained through the KKKS08 analysis (dotdashed lines). Our result for the b-quark FF is in good agreement with the one obtained by KKKS08 analysis. In the case of charm quark and gluon FFs our results behave differently in comparison to the KKKS08 FFs. Nevertheless, as is seen from Tables. I and II, the inclusion of higher order QCD corrections leads to improve the value of χ 2 /d.o.f slightly when going from NLO to the NNLO accuracy. The central values and the size of uncertainty for the fragmentation of heavy quarks q = c +c, b +b and gluon into the D + meson are found similar in behavior at NLO and NNLO QCD corrections. Concerning the gluon and light quarks FFs, remember that there are no free fit parameters associated to their FFs so their forms were taken to be zero at the initial scale. Available data do not provide constraints on the gluon and light quark FFs so they are generated at higher scales through the DGLAP evolution equations [33]. The small uncertainty bands for the gluon FF in Figs. 1 and 2 are expected because they are only sensitive to the quark FFs uncertainties, directly, through DGLAP evolution. Since, the gluon FF really matters, for example, in the process p + p → D + X at the Tevatron or the LHC with a contribution to the differential cross section easily exceeding 50%, then its accurate determination as much as possible is needed. Hence, in order to well constrain the gluon FF, one needs to includes hadron collider observables which also have a significant impact on the gluon FF through a direct determination. However, this is beyond the scope of our present analysis and left for future work.  Fig. 1, but for the D + meson.

B. Data/theory comparisons
In what follows we mainly focus on our NLO and NNLO theoretical predictions for the D 0 and D + mesons production in e + e − annihilation process. Then, we compare the predictions with the inclusive total and b-tagged data sets to judge the goodness of our fit. Fig. 3 includes our theoretical calculation for the total and b-tagged cross sections by applying the evolved FFs at the scale µ = M Z at NLO (solid line) and NNLO (dashed line) accuracies as well as their 1-σ uncertainty bands. Our results are compared to the OPAL data for cross section measurements at the scale µ = M Z . As is seen, Fig. 3 shows a clear agreement between our NLO/NNLO QCD fits and the data sets which indicates a good fit quality of our QCD analyses. Our NLO and NNLO theoretical predictions for the normalized inclusive total and b-tagged cross sections of D + meson production are also shown in Fig. 4. These theoretical predictions are also in good agreements with the analyzed data sets. With a few accuracy, it is seen that the inclusion of higher order QCD corrections up to NNLO leads to reduction in the size of cross section at the region x D ≤ 0.2 and for larger scaling variable x D the NNLO correction does not change the NLO theoretical predictions, significantly. The error bands of our theoretical predictions for NLO and NNLO corrections are rather the same. Besides the (c, b) → H c FFs themselves, their first two moments are also of phenomenological interest. The first one corresponds to the branching fractions which is defined as The second moment corresponds to the average fraction of energy that the produced meson receives from the parent quark Q, namely where the cut z cut = 0.1 excludes the problematic z range where our formalism is not valid, however, in Figs. 3 and 4 it is shown that there are no experimental data at z < z cut . These two quantities stringently characterize the lineshape in x of the Q → H c FFs at a given value of µ and simplify the comparisons with previous FF sets introduced by other authors. Tables. V and VI contain   A comparison between these experimental values and the theoretical ones presented in Tables. V and VI shows the convenient effect of NNLO corrections which lead to much more consistencies. We note also that at each order of perturbative QCD corrections the values of z (µ) are shifted towards smaller values through the DGLAP evolution in µ, as is expected.

VI. ENERGY SPECTRUM OF CHARMED MESONS IN TOP QUARK DECAYS
As a topical application of our D-meson FFs, we study inclusive single D-meson production at the LHC. These mesons may be produced directly or through the decay of heavier particles, including the Z boson, the Higgs boson, and the top quark. For definiteness, we turn to apply the extracted D 0 /D + -FFs to make our phenomenological predictions for the energy distribution of charmed mesons produced in top quark decays through the following process where X collectively denotes any other final-state particles. The study of energy spectrum of produced mesons through top decays might be considered as an indirect channel to search for the top quark properties. At the parton level in the process (13), both the b-quark and gluon may hadronize into the charmed mesons. The gluon fragmentation contributes to the real radiations at NLO and higher orders. Although, the contribution of gluon fragmentation can not be discriminated but it is considered to have accurate predictions for the energy spectrum of outgoing mesons. Ignoring the b-quark mass and by working in the ZM-VFN scheme, to obtain the energy distribution of charmed mesons through the process (13) we employ the factorization theorem as where, following Ref. [46], we defined the scaled-energy fraction of charmed meson as x D = 2E D /(m 2 t − m 2 W ). In the above relation, dΓ/dx i are the differential decay rates of the process t → i + W + (i = b, g) at the parton level. In Refs. [46,47], the NLO analytical expressions of the differential decay widths dΓ/dx i are presented for the polarized and unpolarized top quark decays. In (14), the factorization (µ F ) and renormalization (µ R ) scales are arbitrary but to remove the large logarithms which appear in the differential decay rate, here, we set them to µ R = µ F = m t = 172.9 GeV.
Adopting m W = 80.39 GeV, in Fig. 5 we plotted the energy distribution of D 0 -meson produced in the unpolarized top quark decay at NLO (dot-dashed line) and NNLO (solid line) accuracies. We have also presented an uncertainty band for the NNLO result. Using the NLO FFs of (b, g) → D 0 from Ref. [17], we have also compared our results with the one from the KKKS08 Collaboration (dashed). In Fig. 6, we have presented our NLO and NNLO predictions for the energy distribution of D +mesons along with the uncertainty band for the NNLO result. As is seen from these figures, in comparison with the KKKS08 result our FFs lead to an enhancement(reduction) in the energy spectrum of D 0 (D + ) meson at the peak position. Study of the x D -distribution (i.e. dΓ/dx D ) of the dominant decay mode t → W + + D 0 /D + + X at the CERN LHC will enable us to deepen our understanding of the nonperturbative aspects of D-meson formation by hadronization and to pin down the (b, g) → D + /D 0 FFs. At the LHC, the study of this energy spectrum can be also considered as a new window towards searches on new physics. In Ref. [48], Fig. 8, the energy spectrum of B-mesons have been studied through the dominant decays of top quark in the general two Higgs doublet model (2HDM), i.e. t → b(→ B + X) + H + . There is shown that there is a considerable difference between the energy spectrum of mesons produced in the Standard Model, i.e. t → b(→ B + X) + W + , and 2HDM. In other words, any considerable deviation of energy distribution of produced mesons from the theoretical predictions in the SM theory can be related to the new physics, see also [49,50].

VII. SUMMARY AND CONCLUSIONS
The first set of nonperturbative FFs for bottomflavored hadrons (B-hadrons), both at NLO and NNLO in the ZM-VFNS, was presented in our previous work [51] by fitting to all available experimental data of inclusive single B-hadron production in e + e − annihilation, e + e − → B + X, from ALEPH [52], OPAL [53], SLD [54], and DELPHI [55]. In Ref. [45], we have also provided the first QCD analysis of (c, b) → D * ± FFs up to NNLO accuracy through a QCD fit to SIA data from ALEPH [18] and OPAL [16] collaborations at LEP. In both works, the effects of hadron mass were ignored. In the present study, a new determination of nonperturbative FFs for D 0 and D + mesons are presented at NLO and, for the first time, at NNLO in perturbative QCD considering the hadron mass effects. Of all the available data sets for D 0 and D + productions in e + e − annihilation, i.e. OPAL, CLEO and BELLE, the D 0 /D + -FFs are determined by fitting the data taken by the OPAL Collaboration at LEP, because the data from BELLE is now removed due to an unrecoverable error in the measurements and including the CLEO data sets in the analysis might be a reason for tension. The uncertainties of extracted FFs and SIA cross sections are also estimated using the standard 'Hessian" approach. The determination of uncertainty includes the error estimation for the tolerance criterion of χ 2 = 1.
We judged the quality of our fit by comparing with the prior result presented by the KKKS08 Collaboration [17] ® ® dG/dx D [GeV] x D and also showed how they describe the available data for charmed meson productions. Considering the χ 2 values, we have shown that the inclusions of higher order QCD corrections and hadron mass effects could slightly improve the fit quality. As is seen from Figs. 1 and 2, our results for the fragmentation of b +b into the D 0 and D + mesons are in good agreement with the ones presented by KKKS08. In comparison to the KKKS08's results, there is a considerable difference between the KKKS08 FFs for the transition c +c → D 0 /D + and ours at large values of z. This difference holds for the gluon FFs over the whole range of z. The reason for this inconsistency might be due to the different initial energy scales chosen, the different SIA data sets used and different schemes applied. To check the validity of our results, we have also calculated the branching fractions of mesons and the average fractions of energy carried by the outgoing mesons and compared them with the available experimental data. A good consistency between them, specifically at NNLO, guarantees the validity of our analysis.
As a typical application of extracted FFs, we employed the NLO and NNLO FFs to make our theoretical predictions for the scaled-energy distributions of D 0 /D +mesons inclusively produced in top quark decays.
To describe the novelty and innovation of our work, it should be noted that this study introduces the following improvements. It is one of the first attempts to contain higher order QCD corrections into the D 0 /D + mesons FFs considering all available data. As was explained, the results of Refs. [14,17], as the only references focused on the D-mesons FFs, have included the BELLE data in their analyses as well which are now removed from the Durham database due to some errors reported by the authors. Our work also includes the effects of hadron mass corrections into the FF analysis. Hence, the findings of this paper could make several contributions to the current literature. For future study, it would be interesting to include other sources of experimental observables which carry more information on quark and gluon FFs. Apart from the data, this study is possible if other theoretical partonic cross sections for single inclusive production of partons at NNLO accuracy are available. Another possible area of future research would be to investigate the electromagnetic initial-state radiation (ISR) as well as the effects of heavy quarks mass corrections, however, the latter is of minor importance. The bulk of ISR corrections is due to the effect of electromagnetic radiation emitted from the incoming electrons and positrons. As was explained in Ref. [17], the impact of ISR on the determination of FFs is found to be non-negligible for the analysis of the BELLE and CLEO data which are not included in our work.
The NLO and NNLO D 0 and D + FFs sets presented in this work are available in the standard LHAPDF format [56] from the author upon request.