$\Lambda_c^+$ fragmentation functions from pQCD approach and the Suzuki model

Through data analysis, we present new sets of nonperturbative fragmentation functions (FFs) for $\Lambda_c^+$ baryon both at leading and next-to-leading order (NLO) and, for the first time, at next-to-next-to-leading order (NNLO) in the minimal subtraction factorization scheme with five massless quarks. The FFs are determined by fitting all available data of inclusive single $\Lambda_c^+$ baryon production in $e^+e^-$ annihilation taken by the OPAL Collaboration at CERN LEP1 and Belle Collaboration at KEKB. We also estimate the uncertainties in the $\Lambda_c^+$ FFs as well as in the corresponding observables. In a completely different approach based on the Suzuki model, we will theoretically calculate the $\Lambda_c^+$ FF from charm quark and present our result at leading order perturbative QCD framework. A comparison confirms a good consistency between both approaches. We will also apply the $\Lambda_c^+$ FFs to make theoretical predictions for the energy distribution of $\Lambda_c^+$ produced through the top quark decay, to be measured at the CERN LHC.


I. INTRODUCTION
Study of heavy hadrons properties provides a possibility for better understanding the quark-gluon interaction dynamics in the QCD framework. In this regards and due to ongoing experiments there are particular interests in hadron productions at the CERN LHC and the BNL Relativistic Heavy Ion Collider (RHIC). In this work, we study the production mechanism of heavy baryons through the fragmentation process. In a general expression, the fragmentation mechanism describes the hadronization process where a parton carrying large transverse momentum decays to form a jet containing the expected hadron [1]. The hadronization processes are described by the fragmentation functions (FFs) which refer to the probability densities of hadron productions from initial partons and make the nonperturbative aspects of hadroproduction processes. These functions along with the parton distribution functions (PDFs) [2] construct nonperturbative inputs for the calculation of hadroproduction cross sections. In this work, we focus on the FFs of Λ + c baryon through two different approaches. In the first approach which is usually called as the phenomenological approach, a specific form including some free parameters is proposed for the desired FF so that all these free parameters are extracted through experimental data analysis, see, for example, Refs. [3][4][5][6][7][8][9][10][11][12][13][14][15]. In an alternative scheme based on the theoretical models, the heavy hadron FFs might be computed by the virtue of perturbative QCD with limited phenomenological parameters, see, for example, Refs. [16][17][18][19][20][21][22][23] where the heavy hadron FFs are computed by use of the Suzuki model [24]. This elaborate model is related to the perturbative * Electronic address: Delpasand.mahdi@yahoo.com † Electronic address: mmoosavi@yazd.ac.ir ‡ Electronic address: maryam_soleymaninia@ipm.ir QCD framework where all convenient Feynman diagrams at each order of perturbative QCD are considered for the parton level of hadronization process. In this approach, the nonperturbative aspect of hadronization is emerged in the bound state wave function. Detailed description is presented in Section III. Independent of approach used to extract the FFs, when these functions are computed at the initial scale of fragmentation they can be evolved to higher scales using the timelike Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evaluation equations [25]. In this paper, we first use the phenomenological approach to obtain a set of gluon, charm-and bottom-quark FFs into the Λ + c baryon through a global QCD fit to all available e + e − single inclusive annihilation (SIA) data measured by OPAL Collaboration [26,27] at the CERN LEP1 and very recent data from Belle Collaboration [28]. In Ref. [29], the Λ + c FFs were determined both at leading and next-to-leading order in the minimal subtraction factorization scheme (M S) by fitting the fractional energy spectra of Λ + c baryon measured by the OPAL in the e + e − annihilation on the Z-boson resonance. In their work, authors have applied the zero-mass variable-flavor-number scheme (ZM-VFNS) or, in a short expression, the massless scheme, where heavy quarks are treated as massless partons as well. Recently, KKSS20 Collaboration has updated their previous analysis [29] combining the data for e + e − annihilation from OPAL and the recent one from Belle, see Ref. [30]. Their strategy for constructing the Λ + c -FFs is the same as their previous work for D-meson [31]. Their work is restricted to the NLO accuracy and they have not also evaluated the uncertainties for FFs and corresponding observables. In the present work, we focus on the hadronization of gluon, charm-and bottom-quarks into the Λ + c using the massless scheme and provide the first QCD analysis of (g, c, b) → Λ + c FFs at next-to-next-to-leading order (NNLO). Meanwhile, we go beyond Refs. [29,30] and perform a full-fledged error estimation for the parton FFs as well as the resulting differential cross sections. In order to evaluate the error estimations we apply the wellknown Hessian approach [32]. Note that, our analysis is restricted to only three data sets from single inclusive annihilation process due to two reasons: firstly, we are not aware of any other such data from electron-positron annihilation, and secondly, due to the lack of other theoretical partonic cross sections for the production of partons at NNLO accuracy. Although, among all processes producing baryons, the e + e − annihilation process provides the cleanest environment to calculate the FFs, being devoid of nonperturbative effects beyond fragmentation itself. In the following and in a theoretical approach independent of data analysis, we compute the initial scale fragmentation function of charm quark to split into the Swave Λ + c baryon at lowest-order of perturbative QCD. For this approach, we employ the elaborate Suzuki model which contains most of kinematical and dynamical aspects of hadroproduction process. Finally, we shall compare the initial scale FF of c → Λ + c determined in both approaches. Our comparison shows a good consistency between both results.
In the Standard Model (SM) of particle physics, the top quark has very short life time so does not have enough time to form a bound state, then before it decays hadronization takes place. At the lowest order of perturbative QCD and at the parton level, the decay mode t → bW + followed by b → X+Jets is governed. Here, X refers to the detected hadrons in the final state. Thus, at the CERN LHC a proposed channel to indirect search for the top quark properties is to study the energy spectrum of produced hadrons through top decays. In this work, as an example of possible applications of extracted FFs, we make the theoretical predictions for the energy distributions of Λ + c baryons in top quark decays at LO, NLO and NNLO. This prediction will be compared with the one obtained through using the fit parameters reported in Ref. [30]. This comparison does also shows a good consistency between our analysis and the one performed in Ref. [30].
The outline of this paper is as follows: In Section II, we explain the theoretical framework of hadron production in e + e − annihilation in the massless scheme and introduce our parametrization of the c/b → Λ + c FF at the initial scale. We will also describe the minimization method in our analysis and the approach used for determination of error estimation. All applied experimental data will be describe in this section and our LO, NLO and NNLO results will be presented and compared with the data fitted to. In Section III, through the perturbative QCD approach we provide a general discussion of the fragmentation process for the S-wave heavy baryon and determine the fragmentation distribution of c-quark to fragment into Λ + c baryon at lowest-order of perturbative QCD. In Section. IV, predictions for the normalizedenergy distributions of Λ + c baryon produced from top decay are presented. Our conclusions are listed in Section V.

II. PHENOMENOLOGICAL APPROACH FOR
DETERMINATION OF Λ + c FFS AND THEIR UNCERTAINTIES As was mentioned, one of the most common approach to calculate the unpolarized nonperturbative FFs is the phenomenological approach based on the data analysis. In order to present the theoretical predictions for the observables involving cross section of identified hadrons in the final state, considering different hierarchical features is vital which is mentioned in this section. In this regard, we first review the QCD framework including the standard factorization theorem for differential cross section in a hadronization process of single inclusive electronpositron annihilation. We shall also introduce the OPAL and Belle experimental data as the only data sets for Λ + c production in SIA process. Finally, we indicate our theoretical formalisms for determination of Λ + c FFs and describe our strategy to determine the uncertainties of FFs as well as corresponding theoretical cross sections.

A. QCD framework for Λ + c baryon FFs
Our analysis depends on the normalized differential cross section 1/σ tot × dσ/dx Λ of the annihilation process where, X refers to the unobserved hadrons in the final state and Λ + c is the identified hadron. As usual, the scaling variable x Λ is defined as x Λ = 2p Λ · q/q 2 where, p Λ and q refer to the four-momentum of detected baryon and intermediate gauge boson, respectively, so that s = q 2 is the collision energy. In the center of mass (CM) frame, the scaling variable is simplified as x Λ = 2E Λ / √ s where E Λ shows the energy of detected baryon.
The key implement to divide the perturbative and nonperturbative parts of the e + e − annihilation process (1) is the factorization theorem in the QCD-improved parton model [33]. According to this theorem, the differential cross section of process (1) is written as the convolution of differential partonic cross sections dσ i (e + e − → i + X)/dx i , with the Λ + c -FFs which is denoted by D Λ + c i . Here, i = g, u,ū, . . . , b,b runs over the active partons so that the number of active flavors is dependent on the energy scale. In the ZM-VFNS (or zero-mass scheme) where all light and heavy quarks are considered as massless partons, the differential cross section normalized to the total one is written as [33] 1 In the CM frame, the scaling variable x i is also defined as s which refers to the energy of produced parton i in units of the beam energy. In the relation above, the scales µ F and µ R are the factorization and renormalization scales, respectively. Normally, they are arbitrary quantities which appear in each order of perturbation but in order to omit the ambiguous logarithmic terms ln(s/µ 2 F ) in the partonic cross sections, they are chosen to be µ F = µ R = √ s. The experimental data included in our analysis are normalized to the total hadronic cross section for the e + e − annihilation. This cross section reads Here, α s and α are the strong-coupling and fine-structure constants, respectively, andẽ i is the effective electroweak charge of quark i. The QCD perturbative coefficients K (i) QCD are currently known up to NNLO accuracy [34] QCD ≈ 1.411 [35]. In Eq. (2), the nonperturbative part of the process (1) related to the transition i → Λ + c is described by the D Λ + c i (z, µ F )-FF, where the fragmentation parameter z = x Λ /x i indicates the energy fraction passed on from parton i to the Λ + c baryon, i.e., z = E Λ /E i . Since the FFs depend on the factorization scale µ F , they are evaluated to the various scales of energies by the DGLAP evolution equations, i.e.
where, P ij are the splitting functions which have been computed up to NNLO accuracy [36][37][38].

B. Theoretical formalism
According to the phenomenological approach, in order to extract the nonperturbative FFs, the z-distributions of i → Λ + c FFs at the starting scale µ 0 are parametrized from the beginning and the free parameters are constrained from the SIA experimental data. Note that, the selection criterion for the best parametrization form is to score a minimum χ 2 global value as small as possible with a set of fit parameters as minimal as possible. The χ 2 global function is defined in our previous works in more details [8,9]. Following Ref. [30], we parametrize the z-dependence of c → Λ + c and b → Λ + c FFs at the starting scale as suggested by Bowler [39] while the FFs of light flavors are assumed to be zero at the starting scale. Moreover, since the data set included in our analysis is limited to the SIA process we can not constrain the gluon FF at the initial scale, thus its corresponding FF is also set equal to zero at the initial scale, i.e., The FFs of gluon and light flavors will be generated to higher energy scales via the DGLAP evaluation equations. As was mentioned, for the c → Λ + c and b → Λ + c fragmentation the following Bowler parametrization [39] is considered: which includes six free parameters: N c , a c , b c , N b , a a and b b . It is found that the Bowler parametrization enables one to do excellent fits at each order of perturbation, i.e. LO, NLO and NNLO. Here, the initial scale is set as µ 0 = 4.301 GeV which is a little grater than the bottom mass threshold m b = 4.3 GeV. Consequently, we have six free parameters which should be extracted from the best QCD fit on the experimental data. The optimal values of fit parameters for the charm and bottom FFs are reported in Table I at each order of perturbation. Technically, it should be mentioned that for the evolution of D Λ + c i (z, µ 2 )-FFs as well as for the calculation of SIA cross sections up to NNLO accuracy we employed the publicly available APFEL package [40] and the free parameters of FFs are determined by minimizing the χ 2 global function using the CERN program MINUIT [41]. Furthermore, to estimate the uncertainties of D Λ + c i (z, µ 2 )-FFs the data uncertainties are propagated to the extracted QCD fit parameters using the asymmetric Hessian method (or Hessian methodology), as is outlined in [42,43]. More details can be found in our previous analysis [9].

C. Experimental data and fit results
In our analysis, we applied two data sets measured by OPAL Collaboration at the CERN LEP1 Collider [26]. In fact, for the SIA process (1) two different mechanisms contribute with similar rates; direct production through Z → cc decay followed by the fragmentation c/c → Λ + c and the decay Z → bb (b-tagged events) followed by the fragmentation of b (orb) into the bottom-flavored hadron H b , i.e., b/b → H b , where finally the weak decay of H b into the Λ + c -baryon occurs; H b → Λ + c + X. Consequently, the energy spectrum of Λ + c -baryon originating from the decay of H b -hadron is much softer than that due to primary charm production, as is expected. In order to separate charmed hadron production through the decay process Z → cc from the decay Z → bb, the OPAL Collaboration investigated the apparent decay length distributions as well as the energy spectra of charmed hadrons. As is expected, the decay lengths of H b hadrons into Λ + c baryon are always longer than those from prompt production. Note that, the OPAL Collaboration has presented x Λ -distributions for their Λ + c sample and for the b-tagged (Z → bb) subsamples. In addition to the OPAL data which includes only 4 points with rather large uncertainties we have also included a very recent data set measured by Belle Collaboration [28] at √ s = 10.52 GeV. This new data set is much more precise and contains more points.
Belle data dose not have contributions from B-meson decays so that the contribution of b → Λ c is not needed to be taken in the calculation of cross sections. However, the FFs from charm and bottom quarks are coupled through the DGLAP evolution equation. Following Ref. [30], we fix the b → Λ c FF using the values of N b , a b and b b extracted from the OPAL fit. Therefore, the values of N c , a c and b c are yielded from the fit to the Belle data. Here, we describe a technical point over the reconstruction of old OPAL experimental data. OPAL data sets have been displayed in the form (1/N had )dN/dx Λ where N refers to the number of charmed-flavor hadron candidates which are reconstructed through appropriate decay chains. Therefore, in order to convert these data into the convenient cross section (1/σ tot )dσ/dx Λ , it is needed to divide them by the corresponding branching fractions of decays for the reconstruction of charmedflavored baryons. In Refs. [26,27], for the required branching fraction the following value is applied: Since, in our analysis we are including both the old OPAL data and more recent Belle data which are based on the observation of decay Λ + c → pK − π + , then we have to use the same branching ratio to reconstruct both data sets. Since, the Belle analysis has applied the newest branching ratio as Br(Λ + c → pK − π + ) = (6.635)% [44], we therefore rescale the old OPAL data by the factor 0.044/0.0635 = 0.6929.
Another point about the Belle data is that we have to exclude data at small values of the scaling variable z ≤ 0.5, since the theory is not reliable in this range without taking resummation of soft-gluon logarithms into account. Note that, no kinematical cut is taken over the OPAL data. In Table. II, the characteristics of available experimental data including the number of data points N data n are presented. We have also listed the individual values of χ 2 for inclusive and b-tagged data sets at LO, NLO  In order to show the consistency and goodness between the theoretical prediction and the experimental data used in the fits, in Fig. 1 we have plotted the inclusive differential cross section, as reported by Belle, and in Fig. 2 we have plotted the b-tagged and total differential cross sections normalized to the total one evaluated with our respective Λ + c FFs. Both are compared with the Belle and OPAL data sets fitted to. In these figures, the uncertainty bands are also plotted using the Hessian approach. Through this approach, we just considered the uncertainties due to the experimental data sets so that we ignored additional sources of uncertainties. As is seen the quality of fit is improved when passing to higher order corrections. From Fig. 2, it is seen that our theoretical descriptions at LO, NLO and NNLO for both b-tagged and total differential cross sections are in good mutual agreement. The consistency seems to be better for the normalized total differential cross section (shown in lower panel) in comparison with the b-tagged one because our theoretical predictions do not go across one of the b-tagged data points located at x Λ = 0.8. This is why higher values of individual χ 2 occur for b-tagged, see Table. II.
In Fig. 2, it is also seen that the behavior of theoretical predictions in small values of x Λ for the lowest order accuracy is completely different with the ones at NLO and NNLO. Obviously, in one hand, the theoretical cross section at LO goes to infinity when x Λ → 0 and, on the other hand, the LO uncertainty band is anomalously wider than the NLO and NNLO ones in all ranges of x Λ . Therefore, the LO results are not reliable so that higher order radiative corrections need to be considered. In order to show the fragmentation contribution of gluon, charm-and bottom-quark to the production of Λ + c , in Fig. 3 we have plotted these contributions at the scale √ s = M Z . The total differential cross section at NLO is also shown which is obtained by the sum of all contributions. As is expected, the contribution of gluon fragmen-   In Fig. 4, the z-distributions of Λ + c -FFs are plotted at µ F = M Z ; the energy scale of OPAL data sets fitted to. For this purpose, we plotted the (c, b, g) → D Λ + c FFs at LO (solid lines), NLO (dashed lines) and NNLO (dotdashed lines). From this plot, it is seen that the fragmentation of charm-quark is peaked at large-z whereas the bottom fragmentation has its maximum at small-z. This behavior is due to the fact that the fragmentation process b → Λ + c contains two-step mechanism. In Fig. 4, the uncertainty bands of Λ + c -FFs are also presented which are needed to visually quantify the remaining error of analysis. Since, the Belle date does not include the contributions from the b → Λ c fragmentation in the calculation of cross sections, then the uncertainties of b-quark FF are considerably much wider than the charm and gluon ones in each order of perturbation. Moreover, the error bands of all flavors decrease by increasing the order of perturbation. In Fig. 4, the NLO KKSS20 results [30] (dashed-dashed-dot curves) are also plotted. As is seen, our result for gluon fragmentation is in good agreement with the one presented by KKSS20. In comparison to the KKSS20's results, there is a considerable difference between the KKSS20 charm FF and ours in the range z < 0.5. However, the behavior of our bottom-FF is the same as the KKSS20 one in the whole range of z. Unlike our procedure in which we set the same scale for the cand b-quark FFs, the KKSS20 collaboration has selected different initial scales so that in their work the starting scales for the charm-and bottom-FFs were taken to be µ 0 = m c = 1.5 GeV and µ 0 = m b = 5 GeV, respectively. III. THEORETICAL APPROACH FOR Λ + c

BARYON FF
As was mentioned in the Introduction, apart from the phenomenological approaches to determine the nonperturbative FFs there are also some theoretical models to compute them. In fact, it was fortunately understood that for heavy hadron productions these functions can be analytically calculated by virtue of perturbative QCD (pQCD) including limited phenomenological parameters [45,46]. The first theoretical effort to illustrate the procedure of heavy hadron production was established by Bjorken [47], so that in the following, Suzuki [24], Ji and Amiri [48] have applied the pQCD approach considering elaborate models to describe the hadronization process. Since the Suzuki model includes most of the kinematical and dynamical aspects of hadroproduction process it gives us a detailed insight on the hadronization process. Especially, this model is much suitable to consider the spin effects of produced hadron or fragmenting parton which is absent in the phenomenological approach, see for example [49]. It does also enable us to describe the gluon hadronization process which is not well determined in the phenomenological approach, see Refs. [50][51][52]. In this section, using the Suzuki model we focus on the fragmentation of Λ + c -baryon from charm quark for which the respective Feynman diagrams at leading order in α s are shown in Fig. 5. According to this model, the FF for the production of an S-wave heavy bound state M in hadronization of an initial heavy quark Q might be put in the following general relation [24] where, p Q and p f are the momenta of the fragmenting quark and the final particles, respectively. The fragmentation parameter z is as the one introduced in the phenomenological approach (section II.A), i.e., z = E M /E Q which takes the values as 0 ≤ z ≤ 1. Furthermore, µ 0 stands for the initial fragmentation scale which is in order of fragmenting heavy quark mass and s Q refers to the fragmenting quark spin. In the above relation, the quantity T M is the probability amplitude for the hadron production. In the Suzuki model, this amplitude at large momentum transfer is expressed in terms of the hard scattering amplitude T H and the process-independent distribution amplitude Φ M describing the nonperturbative dy-namics of bound state. In fact, the long-distance amplitude Φ M is, in essence, the probability amplitude for constituent quarks to be evolved into the final bound state. Therefore, the amplitude T M is expressed as [48,53] where, x j = E j /E M is the energy fraction carried by the constituent quark j of heavy bound state M . Considering the general definition (8) and the Feynman diagrams shown in Fig. 5, for the production of S-wave Λ + c -baryon from the initial charm quark, the FF of c → Λ + c is written as where, four-momenta are as labeled in Fig. 5, and  The advantage of above scheme is to absorb the soft behavior of produced bound state into the hard scattering amplitude T H [53]. Ignoring the details, the distribution amplitude Φ B is related to the valence wave function Ψ [53]. Following Ref. [24], we adopt the infinite momentum frame where the distribution amplitude Φ B , with neglecting the Fermi motion of constituents, reads [54,55] Φ where, the baryon decay constant f B is related to the nonrelativistic S-wave function Ψ(0) at the origin as f B = 12/M|Ψ(0)|. In Eq. (11), the short-distance amplitude T H can be calculated perturbatively considering the Feynman diagrams shown in Fig. 5, where a charm-quark creates a heavy baryon Λ + c along with two light antiquarksū and d. In the old-fashioned noncovariant perturbation theory, the hard scattering amplitude T H may be expressed as follows to keep the initial heavy quark on mass shell all the time [56] where, D 0 =P 0 + s ′ 0 + t ′ 0 − p 0 is the energy denominator, C F is the usual color factor and Γ i represents an appropriate combination of the propagators and the spinorial parts of the amplitude. In the above equation, the amplitudes Γ i stand for each Feynman diagrams in Fig. 5. Substituting all expressions in Eq. (10) and carrying out the necessary integrations, we find where, In the above relation, using the well-known completeness relations spin u(p)ū(p) = ( p+ m) and spin v(q)v(q) = ( q − m) in the unpolarized Dirac string, one has 2 i,j=1 (15) where, In the above expressions, after using the Dirac algebra and the traditional trace technique the dot products of four-momenta will appear. To proceed we need to specify our kinematics to determine the relevant dot products. Considering the Feynman diagrams shown in Fig. 5, where by ignoring the Fermi motion of quark constituents the Λ + c baryon is replaced by its collinear constituents, the relevant four-momenta are set as where,P L = r L + t L + s L and we also assumed that the produced baryon moves along theẑ-axes (fragmentation axes). According to the definition of fragmentation parameter, i.e., z =P 0 /p 0 , the baryon takes a fraction z of the energy of initial heavy quark (each constituent a fraction of x 1 , x 2 and x 3 ) and two antiquarks take the remaining 1 − z (each one with a fraction of x 4 and x 5 ). Thus, the parton energies can be expressed in terms of the initial heavy quark energy p 0 , as where, the condition x 1 + x 2 + x 3 = 1 holds as well as x 4 + x 5 = 1. Moreover, according to our assumption that baryon moves along theẑ-axes, the transverse momentum of initial quark is carrying by two antiquarks so that in the infinite momentum frame we have s ′ T = x 4 p T and t ′ T = x 5 p T . With the approximation (12), we are postulating that the contribution of each constituent from the baryon energy is proportional to its mass, namely, We also assume that Regarding the kinematics introduced, the dot products of relevant four-momenta read where, G(z) = M 2 −m 2 c −m 2 u −m 2 d −2t ′ ·s ′ +2p·t ′ +2p·s ′ = m ′ (2m c + β), so that for the remaining integrals, according to the Suzuki model and for simplicity, we replace the transverse momentum integrations by their average values, e.g.,ˆd where, according our assumption one has t ′ T = x 5 p T = (m u /m ′ )p T . Substituting all in Eq. (14), the hadronization process c → Λ + c is described by the following function [48]. In the above relation, 2 i,j=1 Γ iΓj is given in (15) which is simplified in terms of dot products of four-momenta after using the Dirac algebra. Due to the lengthy and cumbersome expression for this term we ignore to present analytical result and just show our numerical analysis. Note that, in the Suzuki model the fragmentation function depends on both the fragmentation parameter z and the phenomenological parameter p 2 T . Although, the zdependence of FFs is not yet calculable at each desired scale, but once they are computed at the initial fragmentation scale µ 0 , their µ evolution is determined through the DGLAP equations [25]. In the Suzuki model, the initial scale is the minimum value of the invariant mass of the fragmenting parton. Therefore, the FF presented in Eq. (24) should be regarded as a model for the c → Λ + c transition at the initial scale µ 0 = m Λ + m u + m d .
For our numerical analysis, we adopt the input parameters as m c = 1.43 GeV, m d = 4.67 MeV, m u = 2.16 MeV, f B = 0.25 GeV, and α s (m c ) = 0.38 [57]. The color factor C F is calculated using the simple color line counting rule so we applied C F = 7/6 for our purpose.
In Fig. 6, taking p 2 T = 1 GeV 2 our theoretical prediction for the D Λ + c c -FF at the starting scale µ 0 is shown (dotted line). In Refs. [16,22], it is shown that the choice of p 2 T = 1 GeV 2 is an optimum value for this quantity so that any higher value of this parameter will produce the peak position even at lower-z regions. To check the validity of the Suzuki model, using the parameters presented in Table. I we have also plotted the D Λ + c c (z, µ 0 )-FF at LO (dashed line), NLO (solid line) and NNLO (dot-dashed line). As is seen, there is a considerable consistency between both approaches. This allows one to rely on the Suzuki model to determine the heavy quark FFs. As was mentioned previously, the Suzuki model is much suitable to consider the spin effects of produced hadrons or fragmenting partons which is absent in the phenomenological approach. It does also enable us to describe the gluon hadronization process which is not well determined in the phenomenological approach. It also gives one a detailed insight on the hadronization process because includes most of the kinematical and dynamical aspects of hadroproduction process.

IV.
Λ + c BARYON PRODUCTION BY TOP QUARK DECAY In this section, as a topical application of our baryon FFs, we study the inclusive single production of Λ + c at the CERN LHC. Generally, the Λ + c baryon may be produced directly or via the decay process of heavier particles, including the Higgs boson, the Z boson, and the top quark. At the LHC, the study of energy distributions of observed hadrons through top quark decays might be considered as an indirect channel to search for the properties of top quarks. Since the top quark discovery by the D0 and CDF experiments at Fermilab Tevatron [58], the full determination of its properties has not yet been performed so it has been one of the main aims in top physics theories.
In the Standard Model of particle physics, the top has a very short life time (τ t ≈ 0.5 × 10 −24 s [59]) which is much shorter than the typical time to form the QCD bound states, i.e., τ QCD ≈ 1/Λ QCD ≈ 3 × 10 −24 s, then the top quark decays rapidly before hadronization takes place. Related to the Cabibbo-Kobayashi-Maskawa (CKM) mixing matrix element for which V tb ≈ 1 [60], top quarks almost exclusively decay to bottom quarks via t → bW + so, in the following, produced bottom quarks hadronize by producing final jets. Therefore, a suggestion for a new channel to look for top properties at the LHC is to study the inclusive single Λ + c -baryon produc-tion through the following process where, X collectively represents any other final-state particles. At the parton level, both the bottom quark and the gluon may hadronize into the Λ + c -baryon so that the gluon fragmentation contributes to the real radiations at NLO. According to the factorization theorem of QCD-improved parton model [61], the partial width of the decay process (25) differential in the scaled Λ + c -baryon energy, x Λ , is expressed as where the factorization and the renormalization scales, i.e., µ F and µ R , are arbitrary but to avoid large logarithms appearing in the parton differential decay rates dΓ i /dx i , we set µ R = µ F = m t , as usual. For simplicity, we shall work in the top-rest frame in which the scaling variables are defined as , where E Λ and E i stand for the energies of Λ + c baryon and parton i, respectively. Here, E max b = m t (1 − ω)/2 is the maximum energy of the bottom quark in the process (25), where ω = (m W /m t ) 2 . At present, analytic expressions for the Wilson coefficient functions dΓ i /dx i are only available at NLO accuracy which are computed in Refs. [62,63]. Using our extracted (b, g) → Λ + c FFs at LO and NLO, we make our predictions for the energy spectrum of Λ + c -baryon produced through the unpolarized top quark decay. However, a consistent analysis is presently restricted to NLO approximation, but we also employ the NNLO set of D Λ + c c -FF to probe the possible size of NNLO corrections.
Adopting the input parameters as m t = 173 GeV and m W = 80.379 GeV, in Fig. 7 we studied the energy distribution of Λ + c -baryon in unpolarized top decays at LO (dotted line), NLO (dashed line) and NNLO (solid line). As is seen, switching from the LO Λ + c -baryon FF set to the NLO one slightly smoothens the theoretical prediction, decreasing it in the peak region and the tail region thereunder. In Fig. 7, the results for dΓ/dx Λ are also compared to the evaluation with the KKSS20 Λ + c -baryon FF set [30]. As is seen, there is a good consistency between both results. In comparison with the KKSS20 spectrum, the peak position of our results is shifted towards larger values of x Λ .
At the LHC, the study of energy distribution of Λ + cbaryon may be also considered as a new window towards searches on new physics. Practically, for the energy distribution of produced hadrons any considerable deviation from the SM predictions can be assigned to the new physics. For example, it would be a signal for the existence of charged Higgs bosons produced from t → Λ + c H + in the theories beyond the SM [64][65][66][67][68]. Meanwhile, the study of x Λ -distribution in the decay mode (25)  able us to deepen our understanding of the nonperturbative aspects of baryon formation by hadronization. Moreover, through studying these distributions the b/g → Λ + c FFs can be also constrained event further.

V. SUMMARY AND CONCLUSIONS
Through this work, we determined the nonperturbative unpolarized FFs for the charmed baryon Λ + c in two various approaches; phenomenological analysis and theoretical approach based on the Suzuki model. Employing the theoretical model we computed the Λ + c -FF at lowest-order of perturbative QCD whereas using the phenomenological approach we determined the Λ + c -FFs both at LO, NLO and, for the first time, at NNLO accuracy in the ZM-VFN scheme by fitting to all available data of inclusive single Λ + c -baryon production in e + e − annihilation from OPAL and Belle Collaborations [27,28]. A comparison between both approaches showed a good consistency between results. Note that, the theoretical framework provided by the ZM-VFN scheme is quite appropriate for our data analysis because the characteristic energy scales of annihilation process, i.e., M Z , greatly exceeds the cand b-quark masses, which could thus be neglected. For our data analysis, we adopted the same functional form for the parameterization of charm and bottom FFs with three free parameters, see Eqs. (6). From Fig. 2, it is seen that in the lowest-order approximation the behavior of theoretical cross section is not acceptable at low-x Λ region but it is reasonably improved when higher order radiative corrections are considered. Through the following aspects, our analysis on the Λ + c -baryon FFs improves a similar analysis in previous works [29,30]. Firstly, we increased the precision of calculation to NNLO, however due to few numbers of experimental data for Λ + c -baryon production our results showed that the effect of this correction is not considerable. Secondly, we did an accurate estimation of the experimental uncertainties in the Λ + c -FFs using the Hessian approach. The uncertainties bands of FFs as well as corresponding observables show that the NNLO radiative corrections affect the error band and decrease them considerably. Meanwhile, we have compared, for the first time, the analytical result obtained for the D Λ + c c (z, µ 0 )-FF through the Suzuki model with the one extracted via the phenomenological analysis. A good consistency between both approaches ensure the Suzuki model, see Fig. 6. This model is suitable to consider the spin effect of produced hadron or/and initial parton on the corresponding FFs, a subject absent in data analysis approach. On the other hand, as is well-known, the gluon FF's play a significant role in hadroproduction but they are only feebly constrained by e + e − data. But, through the Suzuki model it would be possible to determine them analytically, see Refs. [18,23,49]. As a topical application of our obtained FFs, we used the LO, NLO and NNLO FFs to make our theoretical predictions for the scaled-energy distributions of Λ + c -baryon inclusively produced in unpolarized top decays. This channel is proposed for independent determination of Λ + cbaryon FFs which provides a unique chance to test their universality and DGLAP scaling violations; two important pillars of the QCD-improved parton model. Furthermore, this study provides a new window towards searches on new physics. For theoretical approach, one can think of other possible improvements including the Fermi motion of constituents. This is done by considering the real aspects of the valence wave function of baryon [53,69], etc. Related to the phenomenological approaches, improvements due to the inclusion of finite quark masses and the resummation of soft-gluon logarithms would be effective. These effects extend the validity of analysis towards small and large values of x Λ , respectively. In this regards, the general-mass variable-flavor-number scheme (ZM-VFNS) [70,71] where the charm-and bottom-quark masses are preserved from the beginning provides a consistent and natural finite mass generalization of the ZM-VFNS on the basis of the MS factorization scheme [33]. The implementation of these improvements reaches beyond the scope of our present analysis and is left for future researches.