Phenomenology of leading nucleon production in $ep$ collisions at HERA in the framework of fracture functions

In recent years, several experiments at the $e^-p$ collider HERA have collected high precision deep inelastic scattering (DIS) data on the spectrum of leading nucleon carrying a large fraction of the proton's energy. In this paper, we have analyzed recent experimental data on the production of forward proton and neutron in DIS at HERA in the framework of a perturbative QCD. We propose a technique based on the fractures functions framework, and extract the nucleon fracture functions (nucleon FFs) ${\cal M}_2^{(n/p)} (x, Q^2; x_L)$ from global QCD analysis of DIS data measured by ZEUS collaboration at HERA. We have shown that an approach based on the fracture functions formalism allows us phenomenologically parametrize the nucleon FFs. Considering both leading neutron as well as leading proton production data at HERA, we present the results for the separate parton distributions for all parton species, including valence quark densities, the anti-quark densities, the strange sea distribution, and the gluon distribution functions. We proposed several parameterizations for the nucleon FFs and open the possibility of these asymmetries. The obtained optimum set of nucleon FFs is accompanied by Hessian uncertainty sets which allow one to propagate uncertainties to other observables interest. The extracted results for the $t$-integrated leading neutron $F_2^{\rm LN(3)} (x, Q^2; x_L)$ and leading proton $F_2^{\rm LP(3)} (x, Q^2; x_L)$ structure functions are in good agreement with all data analyzed, for a wide range of fractional momentum variable $x$ as well as the longitudinal momentum fraction $x_L$.

These days, a complete understanding of the nucleon and nuclear structures in terms of the underlying partonic and nucleonic constituents is one of the outstanding challenges in hadron physics.
High energy lepton-proton and proton-proton scattering provide the most powerful tools to investigate the hadron structures. In such processes, contributions to the measured nucleon F 2 and nuclear F A structure functions can be expressed in terms of the parton distribution functions (PDFs), nuclear PDFs and spin-dependent PDFs of the nucleon. Precise understanding of PDFs will be a key ingredient in searches for new physics at the LHC through, for example, top-quark and Higgs-boson coupling measurements [1,2]. In consequence, reliable extraction of information on the unpolarized PDFs [3][4][5][6][7][8][9], helicity-dependent PDFs [10][11][12][13][14][15][16][17] and global nuclear PDFs fitting efforts [18][19][20][21] from global QCD analyses of DIS data as well as all related studies , provides deep understanding of the structure of hadrons in term of their quarks and gluon constituents.
HERA as a e ± p collider and unique particle physics data sets collected by the H1 and ZEUS experiments, have provided opportunities to study high-energy electron-proton collisions beyond the electroweak scale [2]. The main process in HERA which is the DIS, probes the internal quark structure of the proton via exchanged virtual photons. The point-like nature of the virtual photon with Q 2 ≫ Λ QCD ensures that the photon can successfully probes the inner structure of the nucleon.
The high center-of-mass energy at HERA has provided searches for the rare processes and physics beyond the Standard Model (SM).
The advent of the HERA collider made it possible to explore a much wider region in momentum fraction x and photon virtuality Q 2 than that previously accessible at fixed target experiments. By applying QCD factorization [45,46] and employing the well-known DGLAP [47][48][49][50] parton evolution scheme, the HERA PDFs can extracted from the high-precision H1 and ZEUS measurements across a large range in x. The extracted PDFs can be used as input to calculate predictions for the Large Hadron Collider (LHC), Large Hadron Electron Collider (LHeC) [51] as well as Future Circular Hadron Collider (FCC-he and FCC-hh) [52] at much higher values of Q 2 .
In addition to the points mentioned above, the productions of energetic neutrons and protons in electron-proton collisions have been extensively studied with the H1 and ZEUS detectors at HERA [53][54][55][56] as well as the related phenomenological studies [57][58][59][60][61]. The H1 and ZEUS experiments at HERA have studied the production of forward protons, neutrons and photons which carry a large fraction of the longitudinal momentum of the incoming proton [53][54][55][56][62][63][64][65]. These analyses have demonstrated that models of DIS are able to reproduce the forward nucleons measurements if contributions from different production mechanisms are considered, such as one pion exchange (OPE) [56], diffractive dissociation, elastic scattering of the proton [56,64] as well as the fracture function formalism [66]. The measurements of leading nucleons also confirm the hypothesis of limiting fragmentation [67,68], according to which, in the high-energy limit, the cross section for the inclusive production of particles in the target fragmentation regions is independent of the incident projectile energy [53,56].
In this paper, we have extracted the nucleon FFs M (n/p) 2 (x, Q 2 ; x L ) from global QCD analysis of DIS data at next-to-leading order (NLO). We have also shown that the fracture functions approach, works well in describing the deep inelastic leading nucleon data measured by the ZEUS collaborations at HERA. We argued that the fracture functions could open some new possibilities for studying hadron structure and open a new window to predict a variety of hard processes at hadron colliders.
The remainder of the paper is organized as follows: We begin in Sec. II by reviewing the formalism for the leading nucleon production including a summary on the several models in which can be used to describe such processes. We will focus on the fracture functions formalism in Sec. II A in which our analysis is based on. The semi-inclusive cross section and the corresponding leading nucleon structure functions have been discussed in Se. II B. The singlet and gluon evolution equations for the nucleon FFs are presented in Sec. II C. The methodology underpinning our global QCD analysis is presented in Sec. III, where we describe in details the parametrizations employed highlighting several improvements in the methodology compared to that introduced originally in the global analysis of the SKTJ17 neutron FFs [57]. The ZEUS-02 [53] leading neutron, and ZEUS-06 [55] and ZEUS-09 [64] leading proton productions data sets analyzed in this study are summarized in Sec. III B. The χ 2 minimization as well as the treatment of uncertainties are presented in Sec. III C.
Our analysis results have been discussed in details in Sec. IV. We compare our theory predictions with the fitted leading neutron and leading proton observables finding overall a good agreement with all data analyzed. Finally, in Sec. V, we summarize our findings and preview future extensions of the present analysis.

II. THEORY SETUP
In the framework of perturbative QCD, the study of leading nucleon production in lepton-proton scattering represents an important field of investigation. Indeed leading nucleon carry a significant fraction of the initial momentum x L > 0.2 and have low transverse momentum p T < 0.7 GeV, therefore covering kinematic regions of the phase space are not accessible from other processes.
However, due to the difficulty of detecting the leading particles in high energy physics experiment, the data available are scarce.
From the theoretical point of view, much successful phenomenological models have been developed to explain the leading nucleon production mechanism. An alternative model to describe the leading nucleon production is based on the the Fracture Functions (FFs) formalism, where the leading particles production is described in terms of structure functions of the fragmented nucleon [61,66,69]. Another picture is given by the Regge formalism, in which the leading nucleon are produced via the exchange of a particle mediating the interaction. Regge theory [70] gives a good description of soft hadronic interactions and can explain the leading nucleon production mechanism [61]. The leading nucleon production can be explained by the exchange of the Reggeon trajectory as well as a Pomeron exchange which produces only leading proton and dominates at x L ∼ 1 [71,72].
Successful descriptions of the available data on the charge-exchange processes, p → n, in hadronhadron and lepton-hadron interactions has been obtained using the exchange of virtual particles with the quantum numbers of the π and ρ mesons. In these kind of processes, the pion, due to its small mass, dominates the p → n transition amplitude. One-Pion-Exchange (OPE) model, therefore, can describes by the process in which a leading neutron is produced. For the leading neutrons production in DIS as well as in dijet photoproduction, the dominating mechanism for x L > 0.6 is the One-Pion-Exchange (OPE) [56]. Based on the assumption that at high x L the leading-neutron production is dominated by the pion exchange mechanism, the measurements of can provide an important information on the pion structure function F π + 2 [56].
It is worth noting here that, the present global QCD analysis is based on the fracture functions approach which can provide a QCD-based description of semi-inclusive DIS in the target fragmentation region [61,66,69]. This formalism, where the leading particles production is described in terms of structure functions of the fragmented proton, has been successfully used to describe leading-nucleon production from the H1 and ZEUS collaborations [57][58][59][60].

A. The Fracture Functions
The fracture functions approach [61,66,69] has been introduced to extend the usual QCD- The first term known as σ current in which arising from target structure function, and only knowledge of the perturbatlvely uncalculable fragmentation function D h i is needed. Such contributions have been widely discussed in the literature (see, e.g, [17,22,73,74] for a clear review).
The second term in Eq. (1) which known as σ target requires a new non-perturbative (but measurable) quantity, a fragmentation-structure or nucleon fracture "functions" (FFs). The σ target (x L ) is given by [66,69,75,76] Here n refers to the detected leading nucleon in the final state and M n/h i (x, Q 2 ; x L ) are the ordinary nucleon fracture functions which cannot be computed in perturbative QCD and like PDFs one need to determine it from fit to leading nucleon data. They are expected to give the dominant contributions to the cross sections for the production of leading hadrons and they can be thought as an ingredient of the perturbative QCD treatment. They refers to the non-perturbative parton distributions of a incoming hadron fragmented into a forward nucleon. Hence, they contains information about partons inside the detected leading nucleon in the final state. σ hard i (x, Q 2 ) in Eq. (2) is the particular hard lepton-parton cross-section we are interested in. Indices i in Eq. (2) stands for the active parton, which interact with the incoming virtual photon, therefore the other partons in the incoming hadron would be the spectators. As we mentioned, the kinemtical variable x is the momentum fraction of struck parton, and x L is defined as the longitudinal momentum fraction in which carried by the final state leading nucleon. It it worth noting here that, 1 − x L is the maximum available fractional momentum of the parton participating in the hard scattering, and hence, in Eq. (2), one need to consider the maximum available fractional momentum in the integration.
We should stress here that the Eq. (2) for the semi-inclusive DIS cross section at low transverse momentum, valid up to power corrections in QCD and to all orders in perturbation theory [75][76][77][78].
On the theory side, there has been great progress in the last few years for the factorization hypothesis in lepton-nucleon collisions. There were also great progress for applying the perturbative QCD to the description of the semi-inclusive DIS, in particular, for the QCD factorization for SIDIS process. Unlike the inclusive DIS, SIDIS is more involved because of additional hadron measurement in the final state. For the SIDIS process, one need to integrate out the transverse momentum (p T ) of the final state hadrons, and similar to the case of inclusive DIS, a collinear factorization is applicable.
The cross section for the SIDIS process can be written as a convolution of the integrated PDFs and the hard partonic cross sections which can be calculated from perturbative QCD. For the low transverse momentum hadron production, in which a collinear factorization approach may not be applicable because the transverse momentum (p T ) of the final state hadron is small compared to the hard scale quantity, Q, one need to introduce a new factorization theorem, involving the transverse momentum (p T ) dependent PDFs. It is instructive to demonstrate this factorization, but it is beyond the scope of our analysis. We refer the reader to the detailed discussions in which can be found in Ref. [75,77]. shown that the fracture functions can well produce the HERA diffractive structure functions [60] as well as leading neutron data [57,58], thus convalidating a common perturbative QCD approach to these particular classes of semi-inclusive processes. Hence, the fracture function approach provides an alternative tool in the framework of QCD to describe the leading nucleon production mechanism.

B.
Leading-nucleons structure functions and observables In comparison to the DIS and diffractive DIS (DDIS) formalism, the cross section of leading nucleons production d 4 σ LB (4) dxdx L dQ 2 dp 2 T can be expressed in terms of leading-Baryons structure functions F LB(4) 2 [71,72] d 4 σ LB(4) dxdx L dQ 2 dp 2 where x L ∼ E B Ep is the energy fraction carried by the produced Baryons or equivalently the longitudinal momentum fraction of the detected Baryons (B) in the final state, and ∆ LB takes into account the effect of the longitudinal structure functions F The transverse momentum of the nucleons in given by p T ≃ x L E p θ n = 0.656x L GeV for ZEUS-02 experiment [53], T < 0.5 GeV 2 for ZEUS-06 [55] and ZEUS-09 [64] experiments. In the OPE approach, the so-scaled fractional momentum variable β is defined by and x is the Bjorken scaling variable. The quantity β may be interpreted as the fraction of the exchange object's momentum carried by the gluon or quarks interacting with the virtual photon.
In term of these variables, the squared four-momentum transfer from the target proton is given by In the measurements presented by H1 and ZEUS collaborations, p 2 T is not measured. The integration over measured range of dp 2 T up to maximum experimentally accessible range of θ max n , corresponding to a p max T , gives It is sometimes more convenient to discuss the measurements in term of the reduced e + p cross section σ LB(3) r which can be written as where F

LB(3) 2
is the leading baryons transverse and F is the longitudinal structure functions [56,58]. We have to mention here that, based on hard scattering factorization and like for the case of diffractive DIS [79], the leading-baryons structure functions can be written in terms of the fragmentation-structure or nucleon "fracture" function and hard-scattering coefficient functions [58,77] as, The index i runs on the flavour of the interacting parton and the Wilson coefficient functions, C q and C g , are the same as in fully inclusive DIS [80]. The p T -unintegrated leading nucleon FFs appearing in Eq. (8) obey the standard DGLAP evolution equations [58]. As one can see from Refs. [53,54,56,64,81] as well as we discussed in section II, in order to describe the leading Baryons production, two more variable in addition to the DIS kinemtical variables are also needed which are x L and t = (p p − p B ) 2 . One can conclude that the M i in Eq. (8) is the t or p Tdependent fracture functions of Eq. (2). In the limit of t ≪ Q 2 , the dominant mechanism is target fragmentation and this expansion holds up to corrections suppressed by powers of 1 Q 2 [77]. The relation between ordinary fracture functions M h/A i (x, Q 2 ; x L ) and the extended fracture functions The above ordinary fracture function has been obtained by integral over p T up to a cut-off of order Q 2 , e.g. ǫQ 2 . In the leading nucleon production, the p T of the leading nucleon is integrated up to some p T,max . We will return to this issue in section III B.

C. Evolution of the nucleon FFs
As we mentioned earlier, fracture functions can provide a QCD-based description of semiinclusive DIS in the target fragmentation regions. This approach has been used to describe the leading-nucleon production data [59,60] and to extract the neutron FFs [57,58]. Like for the DIS structure function, QCD can not predicted the shape of the nucleon FFs. As for the global analysis of the parton distribution functions (PDFs), the non-perturbative nucleon FFs can be parameterized at a given initial scale Q 2 0 . Having at hand the input functional form for the nucleon FFs and the QCD-evolution as well as the corresponding observable, one can extract the nucleon FFs from a global analysis of the DIS leading-nucleon production data. The extracted nucleon FFs can be used to describe the leading nucleon production. These kind of studies can open a new window to predict a variety of hard processes at hadron colliders in which suitable hadronic triggers are used. The presence of leading nucleons in the final state indicates to the long distances processes while the short distances are probed by high photon virtuality Q 2 . In consequence, the study of leading nucleon productions provides information on the relationship between the soft and hard aspects of the strong interaction. These statements and the corresponding discussions can be found in Refs. [72,82].
In Refs. [53,56,81] have been shown that the leading-neutron structure function F LN (3) 2 and inclusive DIS structure functions F p 2 have a similar (x, Q 2 ) behavior as expected from the hypothesis of limiting fragmentation [67,68] in which states that the production of of leading-neutron in the proton fragmentation region is independent of x and Q 2 . The results demonstrated that the similarity between the Q 2 evolution of F LN (3) 2 and Q 2 evolution of F p 2 . Therefore, as in the inclusive case, the Q 2 evolution of the semi-inclusive multi-particle distributions can be predicted in perturbative QCD. Since the scale dependence of the cross section in forward particle production in DIS can be calculated within perturbative QCD [66], therefore the nucleon FFs also obey the standard DGLAP evolution equations [57,58,69,83]. The evolution equations of nucleon FFs are easily obtained by the DGLAP evolution equations [57,66,69] as where M n Σ/P (x, x L , Q 2 ) and M n g/P (x, x L , Q 2 ) are the singlet and gluon distributions, respectively. P Σ and P g are the common NLO contributions to the splitting functions which are perturbatively calculable as a power expansion in the strong coupling constant α s . They are the same as in the case of fully inclusive DIS [84,85]. We have parametrized these non-perturbative distributions, nucleon FFs, at an input scale Q 2 0 < m 2 c . Their evolution to higher scale, Q 2 > Q 2 0 , have been described by using the DGLAP evolution equation within the zero-mass variable flavour number scheme (ZM-VFNS) at the NLO accuracy of perturbative QCD.

III. OUTLINE OF THE ANALYSIS
In this section, we present the method of global QCD analysis of STKJ17 nucleon FFs. The need for well-constrained nucleon FFs has recently been highlighted in global analyses of leading neutron production [57,58]. Differences between these nucleon FFs determinations come from a variety of sources, including different data sets used in the analyses, the choice of parametrization for the nucleon FFs, assumptions about nucleon FFs that are not well constrained by data, or even the method of minimizations. Most of the analyses to date have been performed at leading-order (LO) [60] and NLO [57,58] accuracy in the strong coupling constant.
Note that our aim in this analysis is a definitive determination of nucleon FFs and to explore the application of the fracture function approach to determine the maximal information that can be extracted from the HERA leading nucleon production processes. We will try to propose several parameterizations to consider the asymmetric parton densities in which require the inclusion of most of the possible processes that have sensitivity to nucleon FFs.

A. Input parametrizations
In this section we describe STKJ17 nucleon FFs at the input scale Q 2 0 = 2 GeV 2 . In choosing a functional form for the nucleon FFs at the input scale, it is important to note that the current leading nucleon production observables are sensitive only to the singlet and gluon distributions.
In our previous analysis we therefore seek only to extract the singlet xM n Σ/P and gluon xM n g/P distributions, and do not attempt to separate quark and antiquark FFs [57]. This would require another kind of parametrizations as well as additional data to provide a filter on the quark and antiquark flavors.
It is worth mentioning to notice that the quark distributions for the nucleon FFs at large values of β = x 1−x L could show valence-like structures for some quark-flavour combinations. Although, the accessible values of β or equivalently x in the analyzed experimental data on leading nucleon production are quite low. We should emphasized again that, using hypothesis of limiting fragmentation [67,68], the structure functions for the leading Baryons production is given by [81] As one can see, the first term relates to Baryon variables, {x L , p T }, and the second term related to Lepton variables, {x, Q 2 }. In view of this fact, and in order to account the light quark decomposition, one can assume the following general initial functional form for the singlet xM N Σ/P (x, where n refers to the nucleon and the singlet W q (x L ) and gluon W g (x L ) weight functions define as The label of Σ/P and g/P correspond to the singlet and gluon distributions inside proton, respectively. The x L dependence of the nucleon FFs is encoded in the weight function W i (x L ). We parametrized only the W q (x L ) and W g (x L ) which are related to the Baryon variables. For the xf q and xg(x, Q 2 ) we have used the GJR08 PDFs [86], and hence, this method allows the consideration of the asymmetry of leading nucleon FFs. The initial scale, Q 2 0 = 2 GeV 2 , is chosen at the lowest possible value where a perturbative QCD description can be applied. xf GJR08 q (x, Q 2 ) and xg GJR08 (x, Q 2 ) are the GJR08 parton densities for singlet and gluon distributions [86]. We have selected GJR08 PDFs set because we have used the inclusive structure functions F p 2 (x, Q 2 ) which are obtained from the GJR08 parameterization to account the ZEUS data for the cross section ratio . This selection help us to avoid inconsistencies between the theoretical formalism and the data.
We have also consider a "MRST-type" [87] or old "GRV-type" [88] instead of GJR08, and found consistence results. We have also tested our analysis by considering an independent parametrization for the x dependence of xf q (x, Q 2 ) and xg(x, Q 2 ). The functional forms in Eqs. (13) and (14) state that the dependence of the input parametrizations on the lepton variables, x and Q 2 , should be independent of the leading neutron variables x L .
Considering Eqs. (13) and (14), one can also apply quark and antiquark asymmetries and present the results for the separate parton distributions for all parton species, including valence quark densities, the anti-quark densities, the strange sea distribution functions, and the gluon distribution.
We proposed the following parameterizations for the nucleon FFs and open the possibility of these asymmetries. STKJ17 nucleon FFs at input scale reads where n refers to the nucleon and the weight functions W q (x L ) and W g (x L ) are now defined as where in our analysis. Therefore, our strange densities can be obtained from this assumption and directly sensitive to theū-type fromd-type distributions. We should mentioned here that, one could also consider a standard "MSTW-type" parameterizations for the total light sea quark distribution as xS = 2ū + 2d + s +s as [89] However, in our analysis we prefer to consider the GJR08 NLO parton sets.
As can be seen from Eq. (16), we consider a standard "GJR-type" parameterizations for the up and down valence quark distributions. Since we have used both the leading proton (LP) and leading neutron (LN) data in our analysis with a possible separation of the u and d densities, we believe that these data sets can provide the possibility to determination of valence quark distributions.
The distributions of nucleon FFs in the quark sector at large x may show valence-like structures for some quark-flavour combinations [58]. Since the accessible values of x and x L in the experimental data are quite wide, these data sets can constrain the up and down valence quark distributions well enough. The kinematical coverage of our analyzed data sets can be found in Figs. 1, 2 and 3. Finally, note that the nucleon FFs parametrizations adopted in this analysis is intrinsically more flexible than those used in our previous global analysis. We will show that the above parametrizations will provide us the enough flexibility to extract the non-singlet, valence-like and gluon densities.
The effects of heavy quark masses (m c , m b ) in hard processes and the appropriate definition of parton densities for these quarks have been very actively studied in recent years. When heavy quarks participate in hard processes, the simplest and easiest choice is to treat these particles as massless throughout the calculation, rather than appeal to the more conventional massive parton approximation. However, in most modern global QCD analysis of data, the so called general mass variable flavor number scheme (GM-VFNS) [89] for parton densities are used which is the actual relevance of heavy quark mass corrections. For the case of leading nucleon production at collider DIS as well as hadron collider, a little attention has been paid to the heavy quark contributions. In our previous analysis [57], we have used the leading neutron productions data from ZEUS-02 [53] and H1 [56] collaborations. The H1 semi-inclusive DIS structure function F (β, Q 2 ; x L ) to the inclusive F p 2 (x, Q 2 ) presented in H1 measurement [56] show that these two structure functions have a similar (x, Q 2 ) behavior. This result suggests the validity of hard scattering factorization presented in Eq. (8) and scale dependence of the nucleon FFs in Eq. (10).
As we mentioned, our analysis is based on leading proton ep → epX and leading neutron ep → enX productions in semi-inclusive DIS at HERA. In Table I, we list all data sets included in STKJ17 global analysis along with the corresponding references, the kinematical coverage of x L , x and Q 2 , the number of data points, and finally the normalization shifts N n obtained in the QCD fit. Table I: List of all the leading neutron and leading proton productions data points above Q 2 = 2.0 GeV 2 used in STKJ17 global analysis. For each dataset we provide the references, the kinematical coverage of x L , x and Q 2 , the number of data points and the fitted normalization shifts N n obtained in the fit.

Experiment
Observable In next section, we briefly summarize the selection of data sets and the cuts imposed on them, and the treatment of experimental normalization uncertainties.

ZEUS-02 data on leading neutron production
As we summarized in Table I, we have used three data sets in our QCD fit. For the leadingneutron production in e + p collisions, we have used the data from ZEUS collaboration [53]. At ZEUS experiment, a Forward Neutron Calorimeter (FNC) detected leading neutrons with the energy resolution of σ/E = 65%/ E(GeV), covering 0 ± 0.8 mrad. In order to reduce the systematic uncertainties, the neutron-tagged cross section, ep → enX, is measured relative to the inclusive DIS cross section, ep → eX. Leading neutron production has been studied for x L > 0.2, production angle of the neutron θ n < 0.8, and the photon virtuality up to Q 2 ∼ 10 4 GeV 2 . The transversemomentum acceptance of the leading neutron in ZEUS experiment has been set to p T ≈ x L E p θ n = 0.656 x L (GeV).
In order to minimize the systematics uncertainties, ZEUS collaboration at HERA studied the relationship between the leading neutron production and the inclusive e + p scattering in terms of the ratio of the cross sections. The corresponding ratio r LN (3) in bins of width ∆x L is then defined where the correction term ∆ is small and can be neglected, and F p 2 (x, Q 2 ) is the inclusive DIS structure function and neutron-tagged structure function F LN (3) 2 (x, Q 2 ; x L ) is obtained by the integration over p 2 T from Eq. (6). In our analysis, we have used the inclusive DIS structure function of GJR08 at NLO approximation [86]. Multiplying the r LN (3) ratios by a fit to the inclusive In our global analysis, we have considered both leading neutron and leading proton production data. For ZEUS-02 [53] experiment, we have used the semi-inclusive cross section σ Our theory setup also includes the negligible contributions from the longitudinal structure functions . Since in ZEUS-02 experiment, the cross section for the production of leading neutrons has been determined as a ratio r LN (3) (x, Q 2 ; x L ) relative to the inclusive neutral current cross section, we have used Eq. (18) to account this data. The values of inclusive structure functions F p 2 (x, Q 2 ) are obtained from the GJR08 parameterization of the parton density [86].
In Fig. 1, we plot the nominal coverage of ZEUS-02 leading neutron data sets used in STKJ17 QCD fits for one representative bin of x L = 0.24. The plot nicely summarizes the universal x, x L , and Q 2 dependence of the leading neutron production at HERA.

ZEUS-06 data on leading proton production
In addition to the ZEUS-02 data on leading neutron production, we also include the ZEUS-06 data on leading proton production [55]. Events of the type e + p → e + pX with a final-state proton with x L > 0.6 have been studied in e + p collisions at HERA using the ZEUS detector [54,55], emphasizing the non-diffractive region. The high-energy leading protons with low transverse momentum carrying at least 60% of the incoming-proton momentum and measured in the ZEUS leading proton spectrometer (LPS). The ZEUS-06 data were taken during 1994-1995 using the  The ratio r LP (3) (x, Q 2 ; x L ) is given by The contributions from the longitudinal structure function F L is assumed to be the same for the inclusive DIS and the proton-tagged reactions. For the inclusive DIS structure function, F p 2 (x, Q 2 ), we have used the results of GJR08 [86]. Fig. 2 shows the distribution of the leading proton production data in the Q 2 and x plane for two representative bins of x L = 0.575 and 0.89. We can see that the ZEUS-06 data set provides a rather wide coverage in the (x; Q 2 ) plane.
As we already mentioned, the ZEUS-06 leading proton production has been measured in several x and Q 2 bins for p 2 T < 0.50 GeV 2 which is different from ZEUS-02 leading neutron production which is measured in neutron transfer momentum of p 2 T < 0. b(x L ) = (16.3x L − 4.25) GeV −2 . For the leading proton production at ZEUS experiment [54,55], the values of the slope-parameters b is independent of Q 2 and x L and mean value of b is b = 6.6 ± 0.6 (stat.) ± 0.8 (syst.) GeV −2 for 0.6 < x L < 0.97.

ZEUS-09 data on leading proton production
In addition to the mentioned data sets, we also have used the recent ZEUS data on leading proton production in semi-inclusive reaction e + p → e + pX [64]. The leading proton production cross section and its ratio to the inclusive DIS cross section were studied with the ZEUS detector at HERA with an integrated luminosity of 12.8 pb −1 . In this experiment, leading protons carried a large fraction of the incoming proton energy, x L > 0.32, and its transverse momentum squared satisfied p 2 T < 0.5 GeV 2 . An approximately 24% of DIS events in this measurement have a leading proton. In Fig. 3 we show the kinematic coverage in the (x; Q 2 ) plane of the ZEUS-09 datasets included in our global analysis. The plot nicely summarizes the universal x, x L , and Q 2 dependence of the leading proton productions using ZEUS detector at HERA collider.

C. χ 2 analysis and uncertainties of nucleon FFs
The total χ 2 is calculated in comparison with the leading nucleon data for the nucleon FFs in Eq. (15). The theoretical functions should be obtained at the same experimental x, x L and Q 2 points for calculating χ 2 . As we explained in Sec. II C, the Q 2 evolution is calculated by the DGLAP evolution equations of Eq. (10). The simplest method to calculate the total χ 2 ({ζ i }) for independent sets of unknown fit parameters {ζ i } is as follow where D data Since most experiments come with additional information on the fully correlated normalization uncertainty ∆N n , the simple χ 2 definition in Eq. (20) need to be modified in order to account for such normalization uncertainties. In this case and in order to determine the best fit parameters of Eq. (15), we need to minimize the χ 2 global ({ζ i }) function with the free unknown parameters. This function is given by, where w n is a weight factor for the n th experiment and where n exp correspond to the individual experimental data sets and N data  Table. I. In order to illustrate the effects arising from the use of the different data sets, in Tables. II, III and IV we  The determination of PDFs, polarized PDFs, nuclear PDFs as well as nucleon FFs through a QCD fit to the experimental data is a procedure that necessarily implies a variety of assumptions, mostly concerning their input parameterization and the propagation of the experimental uncertainties into them [1]. In recent years, the assessment of uncertainties has seen significant progress in the QCD analyses. Among the different approaches the Hessian method, the Lagrange multiplier technique (LM) and Neural Network (NN) are the most reliable ones. The "Hessian method" in which use have used to extract the uncertainties of the neutron FFs as well as the corresponding observables were estimated in Refs. [7,8,18,57,89,[96][97][98][99]. Although technical details of the Hessian method are described in these references, outline of this method is explained here because it is used in our analysis.
As we discussed, our method is the χ 2 n ({ζ i }) fitting procedure used in the global QCD analysis and for the determination of the uncertainties, we have used the well-known "Hessian" or error matrix approach. This method confirms that the STKJ17 fitting methodology used in this QCD analysis can faithfully reproduce the input nucleon FFs in the region where the leading nucleon All data sets 388.713 300 Table II: The values of χ 2 /n data for the ZEUS-02 dataset [53] included in the STKJ17 global QCD analysis. More detailed discussion of the description of the individual data sets, and the definitions of   Table IV: The values of χ 2 /n data for the ZEUS-09 dataset [64] included in the STKJ17 global QCD analysis, see details in Table. II.
where H ij are the elements of the Hessian matrix which is given by The confidence region normally is given in the parameter space by supplying a value of ∆χ 2 .
In the standard parameter-fitting criterion, the errors are given by the choice of the tolerance T = ∆χ 2 = 1. It is also known that the confidence level (C.L.) is 68% for ∆χ 2 = 1 if the number of the fitted parameters is one N = 1. It is important to know that for the general case in which the number of fitted parameters is N > 1, the ∆χ 2 value needs to be calculated to determine the size of the uncertainties. This indicates that our fitting methodology as well as the uncertainties determination correctly propagate the experimental uncertainty of the data into the uncertainties of the fitted nucleon FFs.
The determination of the size of uncertainties have been done by applying the Hessian method based on the correspondence between the confidence level P and χ 2 with the number of fitting parameters N . The confidence level P is given by, where Γ is the Gamma function. The value of ∆χ 2 in above equation is taken so that the confidence level (C.L.) becomes the one-σ-error range, namely P = 0.68. Similarly, for the 90th percentile, we have P = 0.90. The value of ∆χ 2 is then numerically calculated by using equation above.
The Hessian matrix or error matrix can be obtained by running the subroutine the CERN program library MINUIT [95]. The uncertainty on any observables O, which is an attributive function of the input parameters obtained in the QCD analysis at the input scale Q 2 0 , is obtained by applying the Hessian method. Having at hand the value of ∆χ 2 , and derivatives of the observables with respect to the fitted parameters, the Hessian method gives the uncertainties on any observables O as, where C j,k is the inverse of the Hessian matrix, H −1 jk . For estimation of uncertainties at an arbitrary Q 2 , the obtained gradient terms are evolved by the DGLAP evolution kernel, and then the nucleon FFs uncertainties as well as the uncertainties of other observables, such as leading nucleon structure function or cross sections, are calculated. We should notice here that a set of uncertainties due to the theoretical method have been included in our analysis. Specially we considered 5% uncertainty due to GJR08 PDFs in which we have used in our definition for the leading nucleon FFs at the input scale, Q 2 0 = 2 GeV 2 . In the next section, we present the main results of this work, namely the "STKJ17" set of nucleon FFs at NLO approximation. First we discuss the resulting nucleon FFs and their uncertainties.
Then, we show the quality of the fits and compare the STKJ17 predictions to the fitted leading nucleon production data sets.

IV. THE RESULTS OF GLOBAL QCD ANALYSIS
In this section we will present and discuss in depth the main results of STKJ17 QCD global analysis of nucleon FFs. First, we present the optimum fit parameters and the constraints applied to control the nucleon FFs parameters. Next, the newly obtained nucleon FFs and their uncertainty estimates are shown. The quality of the fit to ZEUS-02, ZEUS-06 and ZEUS-09 datasets and potential open issues and tensions among the different sets of data are illustrated and discussed in this section.
STKJ17 fitted parameters in the NLO approximation at the input scale Q 2 0 = 2 GeV 2 obtained from the best fit to the combined ZEUS-02 leading neutron data, and ZEUS-06 and ZEUS-09 leading proton datasets are listed in Table. V. All these datasets could provide sensitivity to the flavor separation of the nucleon FFs that was not available in the our previous analysis [57]. In this analysis, it was difficult to determine all unknown parameters of the valence quark densities, the anti-quark, the strange sea and gluon functions of Eq. (15), so that we decided to fix some parameters. It indicates that the data are not sensitive to all parton species at this stage even in the NLO analysis. For valence quark and W ∆ (x L ) densities, we prefer to set C i and D i to zero.
The values without errors in Table. V have been fixed after the first minimization since the data do not constrain these unknown parameters well enough. With the available leading nucleon data sets, for example, the fit can not constrain distinct strange-quark fracture functions. We therefore assume a symmetric strange-quark distributions. As we discussed in section. III A, we considered Fig. 4 we present our results for the W uv (x L ) and W dv (x L ) as a function of x L at the input scale Q 2 0 = 2 GeV 2 . The shaded bands correspond to the uncertainty estimates at 68% confidence level (C.L.) for ∆χ 2 = 1. As can be seen from the plot, both W uv (x L ) and W dv (x L ) have similar pattern and pick at x L ≈ 0.65. For the leading nucleon production at HERA, the shape of the x L spectra is independent of the kinematic variables x and Q 2 confirming the hypothesis of limiting fragmentation [67,68]. The nucleon FFs xM i (x, Q 2 ; x L ) for all parton species resulting from our QCD analysis are shown in Fig. 5 and 6 at the input scale, which is taken to be Q 2 0 = 2 GeV 2 . The shaded bands correspond to the uncertainty estimates at 68% confidence level (C.L.). In terms of uncertainties, the gluon FFs is less well constrained by the leading nucleon production data. As can be inferred from the figures, the error bands for light quark FFs at x L = 0.8 are bigger than those of x L = 0.5.
In Fig. 7 we present a detailed comparison of our theory prediction for the r LP (3) (x, and its uncertainties at 68% C.L. with the ZEUS-06 data already included. The ratio is estimated by dividing the tagged-proton structure function F is approximately constant over a large kinematic range in the leptonic variables is in good agreement with ZEUS-06 data [55]. Our results show that, for nucleon in DIS, hypothesis of limiting fragmentation works well at medium and high values of x L . with the ZEUS-02 leading neutron data [53]. First and foremost, these results demonstrate a good global fit of data taken at different energies Q 2 and kinematic ranges of x with our universal set of neutron FFs.
In order to investigate the validity of our QCD analysis and our obtained theory predictions, in Fig. 9, we present a detailed comparison of our results with the ZEUS-02 leading neutron data [53] for a higher values of x L = 0.61. Overall, one can conclude that our obtained results are in good agreement with all data analyzed, for a wide range of fractional momentum variable x as well as the longitudinal momentum fraction x L .  and its uncertainties at 68% C.L. in comparison with the ZEUS-06 data [55].    [53].

V. SUMMARY AND CONCLUSIONS
Events with a high energetic leading nucleon carrying a large fraction of the proton beam energy have been observed in ep scattering at HERA [53,55,56,64]. In this paper, we have shown that a complete description of semi-inclusive hard processes in perturbabative QCD needs the introduction of new factorizable quantities, known as fracture functions. We have also shown that the fracture functions formalism, in which offers a general theoretical framework for a QCD-based study of leading baryon physics [57-60, 66, 100], works well in describing the deep inelastic leading-nucleon data measured by the ZEUS collaborations at HERA. We argued that the fracture functions could open some new possibilities for studying hadron structure and open a new window to predict a variety of hard processes at hadron colliders. Our analysis is based on the fracture function approach, in which in this framework, semi-inclusive cross sections of leading nucleon production may be written in terms of perturbatively calculable hard-scattering coefficient functions convoluted with appropriate sets of non-perturbative but universal input nucleon FFs constrained by data.
Such a picture was presented here, together with the results of a NLO QCD global analysis of leading nucleon data where it was implemented. We discuss novel aspect of the methodology used in the present analysis, namely an optimized parametrization of nucleon FFs. We compare STKJ17 nucleon FFs set to available leading neutron and leading proton data finding in general a reasonable agreement. We have shown that the partial separation of the nucleon FFs for the various quark flavors has been possible because of the existence of tagged-neutron structure function data from ZEUS-02, and cross section ratio from ZEUS-06 and ZEUS-09 experiments as well as a proposed parametrizations. To further decompose the quark and antiquark neutron FFs, and better constrain the gluon density, additional information will be needed from leading neutron and proton productions in proton-proton collisions. We obtained a relatively good overall description of the leading neutron and proton productions data at both low and high values of x, x L and Q 2 .