First NNLO fragmentation functions of $K_S^0$ and $\Lambda/\bar{\Lambda}$ and their uncertainties in the presence of hadron mass corrections

The current paper presents a determination of $K^0_S$ and $\Lambda/\bar{\Lambda}$ fragmentation functions (FFs) from QCD analysis of single-inclusive electron-positron annihilation process (SIA). Our FFs determinations are performed at next-to-leading order (NLO), and for the first time, at next-to-next-to-leading order (NNLO) accuracy in perturbative Quantum Chromodynamics (pQCD) which is designated as {\tt SAK20} FFs. Each of these FFs is accompanied by their uncertainties which are determined using the `Hessian' method. Considering the hadron mass corrections, we clearly investigate the reliability of our results upon the inclusion of higher-order QCD correction. We provide comparisons of {\tt SAK20} FFs set with the available analysis from another group, finding in general a reasonable agreement, and also considerable differences. In order to judge the fit quality, our theoretical predictions are compared with the analyzed SIA datasets. {\tt SAK20} FFs at NLO and NNLO accuracy along with their uncertainties are made available in the standard {\tt LHAPDF} format in order to use for predictions of present and future measurements in high-energy collisions such as LHC and RHIC.


I. INTRODUCTION
In perturbative Quantum Chromodynamics (pQCD), unpolarized fragmentation functions (FFs) are necessary ingredients to calculate the cross-section of inclusive single hadron production in hard scattering processes [1][2][3][4][5][6][7][8][9][10][11]. In perturbative QCD, Collinear FFs D h i (z, µ 2 F ) can be expressed as a probability for a parton i at the factorization scale µ F to fragment into a hadron h which carries the fraction z of the parton momentum.
In addition to study the z dependence of FFs, we can study the FFs dependency on transverse momentum, P hT which are called the transverse momentum dependent fragmentation functions (TMD FFs) [12][13][14][15][16]. From the factorization theorem [17], the leading twist term of single hadron inclusive production measurements can be interpreted as the convolution of universal FFs with partonic cross-sections of real partons, to account for any hadrons in the final state.
The main motivation to improve our understanding of the details of the subsequent hadronization process is provided by the fact that the FFs along with their associated uncertainties play an important role in several applications in hard scattering processes for the present or future hadron colliders such as LHC, LHeC and RHIC [18][19][20][21][22][23].
To begin with, FFs represent one of the dominant theoretical uncertainties at the LHC measurements. FFs along with their uncertainties also affect the productions of light and heavy hadrons at LHC [24,25]. A second example is the precise measurement of SM parameters at hadron colliders such as LHC, and future high-energy LHC (HE-LHC) and proposed post-LHC particle accelerator in which called Future Circular Collider (FCC) [26][27][28].
Several Collaborations provide regular updates of their light and heavy hadrons FFs sets with uncertainties, see arXiv:2009.08139v1 [hep-ph] 17 Sep 2020 for example [1, 3-6, 11, 29] and references therein. For the K 0 S and Λ/Λ FFs which are the main aim of this paper, the results by BKK96 [30], BS [31], AKK05 [32] and AKK08 [33] Collaborations are available in the literature. In Ref. [30], the authors presented new sets of FFs for neutral kaons both at leading order (LO) and NLO accuracy. The inclusive K 0 production in electron-positron annihilation taken by Mark II at SLAC PEP and by ALEPH at CERN LEP have been used. BS Collaboration [31] has calculated unpolarized FFs for the octet baryons by including some data on proton and Λ production in unpolarized DIS [34,35] in addition to octet baryons production in e + e − annihilation. In addition, AKK05 [32] obtained FFs for K 0 S and Λ at NLO accuracy by a QCD analysis using the data from electronpositron collisions. In order to separate the light quark flavor FFs, they have also included for the first time the quark tagging probabilities from OPAL Collaboration [36]. Finally, AKK08 [33] updated their previous study on K 0 S and Λ FFs, AKK05 [32], and also pion, kaon and proton FFs have been determined in this paper, by adding the inclusive hadron production measurements from proton-proton collisions at PHENIX, STAR, BRAHMS and CDF to their data sample of SIA. They also considered the hadron mass effects in their QCD analyses. Actually the last QCD analysis for fragmentation functions of K 0 S and Λ have been done by AKK08. There is a range of differences between the K 0 S and Λ/Λ FFs determined in the mentioned studies and the QCD analyses done in this paper, arising for example at the level of the selection of the input fitted experimental data, methodological choices for the parameterization of FFs, the detailed estimate, and propagation of the FFs uncertainties and finally the presence of high order perturbative QCD corrections.
The FFs presented in this study introduce some methodological and theoretical improvements over previous determinations available in the literature. The main aim of this paper is to extract the FFs of K 0 S and Λ/Λ along with their uncertainties from a QCD analysis of single-inclusive electron-positron annihilation process (SIA). It should be noted here that the FFs uncertainties for K 0 S and Λ/Λ are calculated for the first time in this paper. In addition, this analysis has been done for the first time, at next-to-next-to-leading (NNLO) accuracy in perturbative QCD. All other determinations of K 0 S and Λ/Λ FFs in the literature are restricted to the NLO accuracy in perturbative QCD without determination of their uncertainties.
In order to achieve a reliable estimate of the uncertainties of K 0 S and Λ/Λ FFs, we use the Hessian approach developed in Refs. [37,38]. We discuss the fit quality, the perturbative convergence upon inclusion of higherorder QCD corrections, and the effect arising from the hadron mass corrections. Finally, we compare our FFs determined in this study to other recent sets of FFs available in the literature. Although, in general, we find reasonable agreements, some important differences are also seen. The effect arising from the hadron-mass corrections on the FFs are carefully investigated and discussed in the text.
The following paper is organized as follows: In Sec. II, we discuss in detail the SIA experimental data along with their corresponding observables and the kinematic cuts which are imposed to determine the FFs for K 0 S and Λ/Λ. In Sec. III, we present the theoretical details of the SAK20 determination for K 0 S and Λ/Λ FFs including the evolution of FFs and the hadron mass corrections. SAK20 parameterizations and our assumptions are are discussed in detailed in Sec. IV. Then, Sec. V deals with the χ 2 minimization and the method for calculation of FFs uncertainty. The main results and findings that emerged from this study are presented and discussed in detail in Sec. VI. We first turn to discuss the SAK20 FFs sets for the K 0 S and Λ/Λ. Then, we compare our FFs set NLO and NNLO with other results in the literature. In Sec. VII we also present comparisons between all analyzed SIA data and the corresponding theoretical predictions obtained using the SAK20 FFs. Finally in Sec. VIII, we study the impact of hadron mass corrections at NNLO accuracy for both K 0 S and Λ/Λ FFs. Sec. IX presents our summary and conclusions.

II. EXPERIMENTAL OBSERVABLES
In this section we discuss in details the experimental data used for determination of K 0 S and Λ/Λ FFs. First, we present the datasets from different experiments and their references. Then, the kinematical cuts applied to the datasets at small range of z will be explained. Finally, we report the χ 2 for all experimental collaborations individually both at NLO and NNLO accuracy.
In our analysis, the data included to extract FFs for K 0 S and Λ/Λ are correspond to the inclusive e + e − annihilation and single hadron production which cover the several range of center-of-mass energies. The K 0 S production datasets included in SAK20 analysis is summarized in Table. I. We specify the name of the experiments, the corresponding references, the measured observables, and the number of data points included in the fit. The values of the χ 2 per data point for both the individual and the total datasets are also reported in this  Finally, in order to determine the well-constrained light and heavy quarks FFs, the (u, d, s)-, c-and b-tagged data from SLD collaboration at √ s = M Z are also added to our data sample.  [49] also included in our data sample. Finally, in order to separate the individual quark flavors, we use the (u, d, s)-, c-and b-tagged data from SLD Collaboration at √ s = M Z [48]. The datasets analyzed in our Λ/Λ QCD fit are listed in Table. II. In this table, the experimental Collaborations and the corresponding published reference, the observable and the center-of-mass energies are listed. The table also include the values of the χ 2 per data point in the individual and total datasets extracted at both NLO and NNLO accuracy.
All the experimental data which we used in this analysis in the (z, Q) plane are shown in Fig. 1 for the K 0 S production and in Fig. 2 for the Λ/Λ production in SIA processes. The applied kinematic cut z < 0.05 is illustrated by the vertical dotted lines in the plots. The range of Q for both hadrons varies from the low energy TASSO data with Q = 14 GeV to the high energy Q = 189 GeV from DELPHI Collaboration. As can be seen, a large number of data points are available for the small z region (z < 0.6).
Our baseline determinations of K 0 S and Λ/Λ FFs are based on the data points described above. Since toward the small z region, soft gluon effects lead the DGLAP evolution equation becomes unstable, then the models fall down the experimental data. Hence, all the theoretical models restrict their analyses to the data points with z ≥ z min in which z min indicates to the low-z cut. The QCD analyses for K 0 S and Λ/Λ in Refs. [30][31][32] excluded the low z regions with z < 0.1 and the hadron mass corrections were not considered in their studies. We should stress here that the recent studies have shown that the mass corrections have an important key role in the small z region. Like for the analysis by AKK08 [33], we consider the hadron mass corrections to be able to include more low-z data points by imposing a kinematic cut at the small values of z; z min = 0.05. Hence, we restrict our data sample to the data points with z ≥ z min = 0.05 for both K 0 S and Λ/Λ analyses. The number of data points after the mentioned kinematical cut for our K 0 S and Λ/Λ analyses are 224 and 137, respectively.

III. THE QCD FRAMEWORK AND HADRON MASS CORRECTIONS
QCD formalism allows us to express the hard scattering cross-section in the term of a convolution of the perturbative partonic cross-sections and non-perturbative distribution functions. The scale dependence of nonperturbative FFs can be obtained by the time-like DGLAP evolution equation in z-space. The computation of the cross-section for the SIA processes along with the DGLAP evolution equations is publicly available up to NNLO accuracy via the APFEL package [50].
The cross-section for the single inclusive electron-positron annihilation in production of strange particles K 0 S and Λ/Λ, e + e − → h(K 0 S ; Λ/Λ) + X, can be given in terms of time-like structure functions, F T (z, Q 2 ) and F L (z, Q 2 ) which can be written in terms of convolutions of non-perturbative unpolarized FFs D h i and perturbative partonic cross sections C i . Hence, the differential cross-section is given by, where the SIA differential cross-section was normalized to the total cross-section σ tot as, Here K (i) QCD show the QCD corrections and have been known yet up to O(α 3 s ). The scaling variable is defined as z = 2P h .q/q 2 with hadron four-momentum P h and γ/Z four momentum q. For the structure functions F T (z, Q 2 ) and F L (z, Q 2 ) presented in Eq. (1), the Wilson coefficient C i functions can be written as expansions in term of strong coupling constant. It reads, ji (z, where j = T, L. These coefficient functions are calculated up to the NNLO accuracy in Refs. [51][52][53]. Note that the non-perturbative universal function, D h i (z, µ 2 ) describes density for fragmenting unpolarized parton i into the unpolarized hadron h which carry fraction z of the longitudinal momentum of the incoming parton. In order to calculate the parton FFs at the different scales of energy µ 2 > µ 2 o , the perturbative QCD corrections lead to use the time-like DGLAP evolution equations [54][55][56][57] which is given by, where P ji (x, α s (µ 2 )) are the time-like splitting functions and describe splitting process i → j + X. These functions can be written as perturbative expansions in term of strong coupling constant which have been calculated up to NNLO accuracy in Refs. [58,59]. In the AKK08 [33] analysis the hadron mass effects are studied for π ± , K ± , p/p, K 0 S and Λ/Λ. In their analysis the hadron masses considered as independent parameters in fit procedure. In addition, the hadron mass effects  have been investigated in Refs. [5,60] for charmed meson and proton productions in SIA processes. We follow the strategy presented in Refs. [5,60] to consider such corrections in our analyses.
In the presence of hadron mass effects with the parameter m h as a hadron mass, the scaling variable need to be modified from z = 2E h / √ s to a specific choice of scaling variable η defined as a light-cone scaling. It is given by, Consequently, the differential cross section in the presence of hadron mass effects for a SIA process need to be modified as The values of the hadron masses used in Eqs. (5) and (6) are considered to be m K 0 S = 0.4976 and m Λ = 1.115 GeV [61]. The above equation for the differential crosssection is applied in our K 0 S and Λ/Λ analyses to consider the hadron mass effects. Eq. (6) indicates that including the hadron mass corrections and the effects arising from that, strongly depend on the hadron mass m h , and hence, the kind of hadron.
We should mentioned here that, for the numerical calculations of the time-like DGLAP equations we use the modified Minimal-Subtraction (M S) factorization scheme. We also used the zero-mass variable-flavornumber scheme (ZM-VFNS) which is implemented at open source framework, APFEL [50]. This scheme assumes that quark mass is set to zero. We applied some modifications in APFEL to take into account the hadron mass corrections. We choose the heavy flavor  masses m c =1.51 GeV and m b =4.92 GeV, and we take µ r = µ f = Q for the QCD renormalization and factorization scales. In our analyses, the QCD running coupling constant is fixed to α s (M Z ) = 0.118 [61]. This selection for the strong coupling constant is consistent with very recent determination of the α s (M Z ) reported by the NNPDF3.1 Collaboration [62].
As a short summary, in this section, was briefly review the pQCD framework for the electron-positron annihilation process, the QCD factorization, and the time-like evolution equation up to NNLO accuracy. We refer the reader to the Refs. [1,33] for more details on the QCD framework. Since our aim in this analysis is to investigate the effect of higher order perturbative corrections, we will clearly discuss in Sec. VI the improvements to the fit quality, the value of χ 2 for each dataset, and the total χ 2 /d.o.f at NLO and NNLO accuracy. The comparison of the central values and error bands of K 0 S and Λ/Λ FFs in these two perturbative orders will be presented in Sec. VI as well.

IV. THE TECHNICAL FRAMEWORK AND OUR PARAMETRIZATION
In the following section, we are in a position to describe our methodology, the input parametrization, and the assumptions that we consider to extract the K 0 s and Λ/Λ FFs from perturbative QCD analysis to the available SIA experimental data.
Since the main goal in this analysis is to investigation of the FFs of K 0 s and Λ/Λ, we should parameterized the light and heavy quark FFs at initial scale Q 0 = 5 GeV which should be above the bottom mass threshold in ZM-VFNS. Then QCD evolution will help us to achieve it at arbitrary scale Q. In this analysis we use the most flexible parametrization form for the K 0 S and Λ/Λ FFs with n f = 5 active flavor in which widely used in the analysis of different hadrons in the literature [33,63,64]. Most recently, we also considered such parametrization for the determination of unidentified light charged hadron and pion FFs analyses [4,65]. This parameterization is given by, where the free parameters are N i , α i , β i , γ i , and δ i and A i is the normalization factor. In above parameterization form N i and A i are not independent, N i is the second moment of the parton fragmentation function and A i is the normalization factor which can be computed to be: where B[a, b] is the Euler Beta function. Note that by including only the SIA data in the QCD fit, it is not possible to separate the quark and anti-quark FFs, then we use the quark combination q + = q +q in the parameterization form. In equation (7), i indicates to the d + , u + , s + , c + , b + and g. Now we are in a position to discuss our assumptions for the parametrization of K 0 S FFs. Considering the quark content of the K 0 S (ds), we assume asymmetry between light quarks u, d, s, and parametrize them separately as like heavy quarks and gluon. Since statistically the number of experimental data points from e + e − annihilation to determine the K 0 S is rather limited, all the parameters can not be well constrained by these datasets. The parameters γ and δ are free for d + , b + and g and they need to determine from the QCD fit. For the rest, we fix γ u + ,s + ,c + and δ u + ,s + ,c + to zero. Consequently, the remaining 24 free independent fit parameters of FFs are determined by a standard χ 2 minimization method. Bestfit parameters for the fragmentation of partons into K 0 S obtained through our NLO and NNLO analyses are listed in Table. III.
The Λ/Λ baryon contains the (uds) quarks. Hence, we define separate parametrization for all light quarks and we do not assume SU(2) or SU(3) symmetry for Λ/Λ. Since the number of available data for the Λ/Λ production in SIA process is not enough to constrain all independent fit parameters in Eq. (7), we prefer to consider a simple form of parametrization and fix γ = 0 and δ = 0 for all flavors except for the gluon density. The total number of free parameters for Λ FFs is 20. The Best-fit parameters for the quarks and gluon FFs of Λ/Λ at NLO and NNLO accuracy are presented in Table. IV.

V. χ 2 MINIMIZATION AND METHOD OF ERROR CALCULATION
Our fitting methodology and χ 2 minimization and the uncertainty estimation have been described at length in our previous publications [4][5][6] on the same subjects. In this section, we briefly review methodology specific to the K 0 S and Λ/Λ FFs determinations. We first discuss the minimization strategy to optimize the independent fit parameters, and then we present the uncertainty estimations.
As we mentioned earlier, we perform our QCD analysis using the standard functional form at the initial scale of µ 0 , then we evolve the FFs from the initial scale up to arbitrary scale using the DGLAP evolution equation [54][55][56][57] to calculate the physical observable. By comparing the theoretical prediction with the corresponding experimental data in the full kinematic range, we determine the unknown FF parameters by constructing a global χ 2 function using the experimental measurement and theoretical prediction for i th data point. The χ 2 function is minimized using the CERN MINUIT package [66] and is constructed as follows: where the weight factor w i allows us to apply separate weights to different experimental datasets which in this analysis we take it to be unity. The index i sums over all experimental datasets. For each dataset, the index j sums over all data points. F Theory 1,j is the theoretical prediction for j th bin, ∆O Exp 1,j is included the statistical and systematic errors which we combine in quadrature, and O Exp 1,j is the measured value of the i th data point. In the above, the normalization shifts K i for each experiment, are fitted at the first step of our procedure, and then keep fixed in our analysis. Note that the associated normalized uncertainty ∆K i is obtained by experimental setup.
We now describe the methodology that we applied in this study for the estimation of the K 0 S and Λ/Λ FFs uncertainties. Large amount of QCD analyses use the 'Hessian' method to calculate the FFs uncertainties which is based on tolerance parameter T . They consider ∆χ 2 = T 2 which ensures that each dataset is described within the desired confidence level (CL). The standard error propagation which is given by the statistical error on any given quantity q, is defined as: To calculate the fully 1-σ error bands for the FFs, one could use the Hessian matrix definition, H α,β = 1 2 ∂ 2 χ 2 /∂p α ∂p β which is inverse of the covariance matrix C = H −1 that is obtained at the χ 2 minimum. In this analysis, we adopt the standard parameter fitting criterion by choosing the T=1, which corresponds to the 68% CL, i.e., 1-σ error bands. The details of the 'Hessian method' are fully addressed in Refs. [37,38], and we refer the reader to these published works for more details.

VI. THE RESULTS OF FFS ANALYSIS
The following part of this paper describes in greater details the results of SAK20 K 0 S and Λ/Λ FFs with their uncertainties. First, we present the best-fit parameters for the fragmentation of partons to K 0 S and Λ/Λ. Second, we present SAK20 FFs at NLO and NNLO accuracy and compare them with each other. Next, we quantify the perturbative convergence of the SAK20 FFs upon the inclusion of higher-order QCD corrections. Finally, we compare SAK20 FFs with the corresponding results by AKK08 FFs Collaboration [33]. We show that there are conspicuous differences between SAK20 and AKK08 FFs, specially for light quarks and gluon. As we mentioned before, we will divided our analysis into two separate fit. The first one is perform a fit with all data of Table. 1 to extract the K 0 S partonic FFs using the K 0 S meson production in SIA process, and we identify this as SAK20 K 0 S FF, and the second one is perform a fit with all data of Table. 2 to extract the Λ/Λ partonic FFs using the Λ/Λ baryon production in SIA process, and we identify this as SAK20 Λ/Λ FFs. In the following, we will present the resulting FFs along with their uncertainties. We present the results of K 0 S FFs and their uncertainties in Sec. VI A, and in Sec. VI B, the results of Λ/Λ FFs along with their error bands will be discussed in details.

A.
The results of K 0 S FFs and their uncertainties To initiate the discussions of the K 0 S FFs, we first, present the optimal set of K 0 S FFs parameters which have been derived by minimizing the χ 2 as defined in Eq. 7 by comparing to the measured SIA data presented in Table I. The details of the fit are summarized in Table III which shows the best fit values of the free parameters based on Eq. (7).
We now discuss the overall statistical quality of the fit as measured by the total χ 2 per degree of freedom for the SAK20 K 0 S fit. We find that the total χ 2 /d.o.f. of 1.161 and 1.124 for our NLO and NNLO QCD analysis, respectively, indicating a good quality of fit. Furthermore, as one can see from Table. I, the inclusion of higher-order QCD correction leads to a smaller value for χ 2 /d.o.f. which indicates that our NNLO analysis improve the fit quality.
In the following, we now turn to discuss the K 0 S FFs and their uncertainties obtained from the global fit. As we explained in more detail in Sec. IV, the present analyses adopt the most traditional fitting framework at NLO and NNLO accuracy assuming a very flexible functional form to parameterize the FFs at an initial scale. The SIA data analyzed in this study could not discriminate between quark and antiquark FFs. Hence, we display in Fig. 3 the SAK20 results for zD K 0 S i (z, Q), i = d + , u + , s + , c + , b + and g at the initial scale of Q 0 = 5 GeV. To investigate the effect of higher order correction, we prepare a comparison between NLO and NNLO fit results and their uncertainties in Fig. 3. Although there is no significant change in size for quarks and gluon FFs, a small difference between NLO and NNLO can be observe for zD K 0 S c + (z, Q) and zD K 0 S g (z, Q). As we mentioned before, the uncertainty bands of FFs presented for the choice of tolerance T = ∆χ 2 = 1 for the 68% (one-sigma) confidence level (CL) obtained using Eq. (10). Fig. 3   shows that the errors for the gluon and u + FFs are large in both NLO and NNLO, which means that they are not well determined particularly at small value of Q.
In order to investigate the impact of the inclusion of higher order corrections in more details, in Fig. 4, we show the ratios of NNLO SAK20 FFs (magenta bands) to the corresponding NLO results (green bands) at Q = M Z . As can be seen, the NLO and NNLO uncertainties presented in this figure are similar in size showing that the improvements of FFs uncertainty upon inclusion of higher-order QCD corrections are not significant when going from NLO to NNLO. However, the total χ 2 /d.o.f. that we obtained indicates that the inclusion of NNLO QCD corrections slightly improves the overall fit quality as well as the description of the data.
In Fig. 5, the obtained zD K 0 S i (z, Q), i = d + , u + , s + , c + , b + and g from our analysis as function of z are presented at NNLO accuracy at Q = M Z , and compared them with the results obtained by AKK08 FFs Collaboration [33]. The shaded bands correspond to the uncertainty estimation based on Hessian approach for ∆χ 2 = 1. Concerning the shapes of the K 0 S FFs, several interesting differences between SAK20 and AKK08 can be seen from the comparisons in Fig. 5. Compared to the AKK08 FFs, one can see weak agreements between two results except for the c + distribution. Fig. 5 shows that the SAK20 FFs for u + , b + and gluon densities are smaller than AKK08 for the medium to small value of z. For the case of d + FF, one can see that, the AKK08 is smaller than our result for the whole range of z. The origin of the differences among the SAK20 and AKK08 over the whole z range, is likely to be mostly due to the inclusion of inclusive hadron production measurements from proton-proton collisions data in AKK08 while SAK20 is limited to the SIA data only.

B. The results of Λ/Λ FFs and their uncertainties
In this section, first we mentioned the best fit values of Λ/Λ free parameters presented in Table. IV which have been derived from QCD fit from SIA data. Then we consider the total χ 2 /d.o.f and they are equal to 1.601 and 1.602 at NLO and NNLO, respectively, and indeed we can not see the noticeable improvement in χ 2 /d.o.f value at NNLO in comparison to NLO accuracy.
In the following, we now turn to discuss our results and findings for the determination of Λ/Λ FFs and their uncertainties. In order to study the perturbative convergence of the Λ/Λ FFs upon inclusion of higher order QCD corrections, we first compare our NLO, and NNLO deter-minations among each other at the input scale. The same comparison will be presented as ratios of NNLO SAK20 FFs to the corresponding NLO results. Then, we compare our best-fit NNLO Λ/Λ FFs to their counterparts in the AKK analysis at the scale of M Z .
In the following, we show the FFs results and their uncertainties at NLO and NNLO accuracy, focusing on their perturbative convergence upon inclusions of higherorder QCD corrections. We display the six FFs combinations parameterized in our QCD fits for the Λ/Λ hadrons, and their 1-σ uncertainties in Fig. 6. For each partonic species, the FFs are shown at NLO and NNLO as functions of z at our input scale Q 0 = 5 GeV. A noticeable aspect of the comparisons in Fig. 6 are related to the shape and the size of the FF uncertainties. The u + , d + and b + FFs are similar in shape, while for other FFs, a small differences are observed. For both NLO and NNLO FFs, the s + and gluon FFs turn to zero for small values of z < 0.2. Overall, the differences between NLO and NNLO FFs are slightly small. This is consistent with the perturbative convergence of the global χ 2 that we discussed in Sec. II, (see the χ 2 /d.o.f presented in Table. II).    Table. III, but for Λ/Λ FFs.
In order to judge the effect arising form the these corrections, we present in Fig. 7 the ratios of NNLO SAK20 FFs to the corresponding NLO results at the scale of Q = M Z as functions of z. While the b + and d + uncertainties are seem to be similar in size, the s + , u + and gluon FFs uncertainty bands are in general smaller at NNLO accuracy, particularly those of the s + and u + FFs. The uncertainty band for the c + NNLO FFs is visibly larger at small value of z and smaller at large values for z. Although the total χ 2 /d.o.f presented in Table. II indicate that there is no improvement from NLO to NNLO accuracy, the findings in this figure show the reduction of error bands for three partons at NNLO in comparison to the NLO.
We now compare SAK20 FFs to the recent determination available in the literature, namely the AKK08 [33] FFs set. Their analysis is performed only at NLO accuracy. Such a comparisons are shown in Fig. 8 at Q = M Z for all partonic species. Concerning the shapes of these FFs, several interesting differences between these two sets can be seen from the comparisons in Fig. 8. As can be seen, for the b + FFs, the SAK20 and AKK08 results are in good agreement. For other parton species, the differences in shape among these two FF sets are more marked than in the case of b + FFs and relatively large differences are observed. The origin of differences among these two FF sets, at small to medium values of z for most of the quark FFs and the gluon FF, is likely to be mostly due to the inclusive hadron production measurements from protonproton reactions that AKK08 included in their analysis but it is not considered in the SAK20 FFs sets.

VII. FIT QUALITY AND COMPARISON TO THE SIA DATA
This section deals with our global fit results in term of fit quality and detailed comparison to the SIA experimental measurements. First, we compare our NNLO theory predictions with the K 0 S production data analyzed in this study. Then, we present the comparison of our theory predictions with the Λ/Λ production cross-section measurements. In Tables. I and II, we report the χ 2 per point for each datasets and total χ 2 per degree of freedom included in the SAK20 FFs analysis. The values are shown at NLO and NNLO accuracy for both K 0 S and Λ/Λ FFs determinations. Concerning the fit quality of the total SIA datasets analyzed in SAK20, the most noticeable feature is that the χ 2 /d.o.f value for K 0 S shows almost 3% improvement on total χ 2 when the higher order correction is considered. However, for the case of Λ/Λ the values of χ 2 /d.o.f are almost the same at NLO and NNLO accuracy.
Comparison between the K 0 S production dataset in SIA process analyzed in this study from different experiments and the corresponding theoretical predictions using SAK20 best-fit at NNLO accuracy are shown in Figs. 9, 10 and 11. We show the comparisons as the data/theory ratios. The error bands indicate to the 1-σ FF uncertainties. In Fig. 9, comparisons are displayed for the TASSO data for different center-of-mass energies. In Fig. 10, we show the same comparison for all the inclusive experimental data analyzed in this study, except the SLD. Deviation between the theory and the data can be seen for the large value of z for HRS, DELPHI 183 and DELPHI 189. These findings consistent with the χ 2 values listed in Table. I. The data/theory ratios for the SLD measure-  Fig. 11. One can see a good agreement between SLD inclusive, uds-and b-tagged data and theory predictions at NNLO accuracy obtained from our QCD fit. However, the agreement between SLD c-tagged data and our theory is not as well as the others SLD data.
In Fig. 12, the comparisons have been shown between the theoretical predictions using SAK20 and the TASSO Λ/Λ production for different center of mass energy. According to Fig. 12 and Table. II, the TASSO 34.8 dataset is in good agreement in low and medium z regains with theoretical prediction obtain from our QCD analysis. This plot also reveals that our theoretical prediction could not describe the latest data point in large z region. Comparison between the Λ/Λ production and our results is shown in Fig. 13 for all the inclusive experimental data except SLD. Finally, a comparison with the SLD data in inclusive, uds-, c-, b-tagged for Λ/Λ production are shown in Fig. 14. In general, an overall good agreements between the data from all experiments and SAK20 NNLO theoretical predictions are achieved, which consistent with the individual χ 2 values reported in Tables. I and II. Remarkably, our theoretical predictions and the data are in good agreements from small to and large values of z.

VIII. IMPACT OF HADRON MASS CORRECTIONS
In order to study the effects arising from the hadron mass corrections, we perform the analyses in which we do not consider hadron mass corrections for K 0 S and Λ/Λ and repeat our QCD fits just at NNLO accuracy. We entitle this analysis to massless and our main analyses by considering hadron mass corrections is called massive. Finally, in order to investigate in details the effects arising from the inclusion of such corrections on our K 0 S and Λ/Λ FFs determinations, we compare the results of our massless and massive analyses. It is worth mentioning here that such corrections could affect the small z regions [1].
In Fig. 15 we show the comparison between our massless and massive results as ratio, to investigate the effects of mass correction for all parton species as function of z at Q = M Z . As can be seen, the central values of d + , s + and c + are not affected noticeably by considering mass corrections, the central values for u + , b + and g change remarkably in medium to large z regions. There are differences for the error bands between two analyses, specifically at small values of z. The reduction of uncertainties in the region z < 0.1 can be seen for all partons in massive analysis. However, the error bands increase for the z > 0.1 in massive analysis for all the partons, except c + . Also the rise of error uncertainties for u + and b + are dramatic. According to Fig. 1, statistically the K 0 S experimental data included in our analysis for z < 0.6 is more than for z > 0.6. Consequently, more data points at large values of z could constrain the FFs.
We are now in a position to compare the SAK20 FFs for Λ/Λ with and without the hadron mass effects at NNLO accuracy. Such comparisons are shown in Fig. 16 as a ratios to the SAK20 FFs without such corrections. Concerning the FF uncertainties upon inclusion of hadron mass corrections, we observe that for the quark and gluon distributions, including such corrections significantly affect the uncertainty bands. The smaller uncertainties of the SAK20 FFs in the presence of hadron mass effects as shown in Fig. 16 may be due the fact that such corrections affect the shape and the uncertainty bands at small values of z. As can be seen from Fig. 16, the hadron mass corrections significantly affect the central values of d + , s + , c + and b + more than u + and g FFs. The hadron mass corrections decrease the uncertainties for u + and g in all range of z, and for s + in the range of z < 0.3. The behavior of error bands treat differently for d + , c + and b + for whole range of z.
As a conclusion, applying the hadron mass corrections in our analyses for K 0 S and Λ/Λ can generally decrease the error bands of FFs in small z region. As expected from Eq. (6) including the hadron mass corrections strongly depend on the hadron mass m h , and hence, the kind of hadron. Then this corrections affected the Λ/Λ more than K 0 S , see Figs. 15 and 16. According to Figs. 1 and 2 the number of data from SIA process for K 0 S and Λ/Λ production are being poor at z > 0.6. Hence, one can see the large uncertainties in this region.

IX. SUMMARY AND CONCLUSIONS
As a summary, the main goal of the current study is to present the new sets of FFs for K 0 S and Λ/Λ from QCD analyses of single-inclusive electron-positron annihilation process (SIA) at NLO, and for the first time, at NNLO accuracy in pQCD. These analyses are based on com- prehensive experimental datasets in which include the precise measurements of K 0 S and Λ/Λ production crosssections in the SIA process. Well-established QCD fitting methodology used to determine K 0 S and Λ/Λ FFs. In this analysis, we have introduced several methodological improvements to determined the FFs. As a first improvement, we have performed the first QCD analysis at NNLO approximation to clarify the role of the higherorder QCD corrections on the description of the data. The related results are clearly presented in this paper and compared to our NLO study and the corresponding results from AKK08 [33]. As a next improvement, the calculations of these new FFs for K 0 S and Λ/Λ are along with the determination of uncertainties. The 'Hessian' method is used to provide a faithful representation of the experimental uncertainties. In addition, we apply the K 0 S and Λ/Λ mass corrections in our analyses to include the low-z data points. We also study the hadron mass effects on the different species of patrons and their error bands. As a final point, the analysis reported in this paper represents the first step of a broader project, and hence, a number of improvements are foreseen. SAK20 analyses are based only on the SIA measurements for the K 0 S and Λ/Λ production measurements. Although the SIA data is the cleanest process for FFs determination, it may not possible to consider the flavor separation. It is also a little sensitive to the gluon FF. Hence, the FFs presented in this paper could be improved by adding the data from proton-proton (pp) collisions in our data sample. Gluon FFs could be well-constrained by the hadron collider data. A further improvement would be the inclusion of heavy-quark mass corrections in which could improve the description of the data, especially at the lowest center-of-mass energy. The SAK20 FFs for the K 0 S and Λ/Λ presented in this work are available via the standard LHAPDF interface [67].

ACKNOWLEDGMENTS
Authors thank School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM) for financial support of this project. Hamzeh Khanpour also is thankful the University of Science and Technology of Mazandaran for financial support provided for this research. FIG. 9: (color online). Comparison between the K 0 S production cross-section analyzed in this study from different experiments by TASSO Collaboration [39,40]. The comparisons have been shown as data/theory ratios. The error bands indicate to the one-σ FF uncertainties. The results shown in this plots are correspond to SAK20 NNLO fit in the presence of hadron-mass corrections. . Same as Fig. 9 but for the SIA data from HRS [41], TPC [42], MARK II [43], CELLO [44], TOPAZ [45], ALEPH [46], DELPHI [47,49], OPAL [36] and SLD [48] Collaborations.    Fig. 12 but for HRS [41], MARK II [43], CELLO [44], ALEPH [46], DELPHI [47,49] and OPAL [36] Collaborations datasets.