Determination of neutron fracture functions from a global QCD analysis of leading neutron production at HERA

TL;DR

This paper performs a global QCD analysis to determine neutron fracture functions from leading neutron production data at HERA, providing a detailed NLO fit and uncertainty estimation that aligns well with experimental results.

Contribution

It introduces a novel NLO global QCD fit of neutron fracture functions using HERA data, including uncertainty quantification with the Hessian method.

Findings

Neutron fracture functions are successfully extracted from data.

Predictions agree well with experimental measurements across kinematic ranges.

Uncertainty estimates are provided for the neutron FFs and related observables.

Abstract

In this article, we present our global QCD analysis of leading neutron production in deep inelastic scattering at H1 and ZEUS collaborations. The analysis is performed in the framework of a perturbative QCD description for semi-inclusive processes which is based on the fracture functions approach. Modeling the non-perturbative part of the fragmentation process at the input scale Q$_0^2$, we analyze the Q$^2$-dependence of the leading neutron structure functions and obtain the neutron fracture functions (neutron FFs) from next-to-leading order (NLO) global QCD fit to data. We have also performed a careful estimation of the uncertainties using the "Hessian method" for the neutron FFs and corresponding observables originating from experimental errors. The predictions based on the obtained neutron FFs are in good agreement with all data analyzed, at small and large longitudinal momentum fraction $x_L$ as well as the scaled fractional momentum variable $\beta$.

Tables

Table 1. Table 1: List of all the leading neutron production data points above Q 2 superscript 𝑄 2 Q^{2} = 1.0 GeV 2 used in SKTJ17 global analysis. For each dataset we provide the corresponding references, the kinematical coverage of x L subscript 𝑥 𝐿 x_{L} , x B subscript 𝑥 𝐵 x_{B} and Q 2 superscript 𝑄 2 Q^{2} , the number of data points and the fitted normalization shifts 𝒩 n subscript 𝒩 𝑛 {\cal{N}}_{n} .

Experiments	[$x_L^{\min },x_L^{\max }$]	[$x_B^{\min },x_B^{\max }$]	Q2 GeV2	Number of data points	${𝒩}_n$
H1 Aaron:2010ab	[$0.365$–$0.905$]	[$1.5\times {10}^{-4}$ – $3.0\times {10}^{-2}$]	7.3–82	203	0.9922
ZEUS Chekanov:2002pf	[$0.240$–$0.920$]	[$1.1\times {10}^{-4}$ – $3.2\times {10}^{-2}$]	7–1000	300	1.0033
Total data				503

Table 2. Table 2: The values of χ 2 / N pts . superscript 𝜒 2 subscript 𝑁 pts \chi^{2}/N_{\rm pts.} for the data sets included in the SKTJ17 global QCD analysis. More detailed discussion of the description of the individual data sets and the definitions of χ 2 superscript 𝜒 2 \chi^{2} are contained in the text.

Experiment	Data set	${\chi }^2$	Npts
	$x_L$ = 0.365	24.13	29
	$x_L$ = 0.455	25.62	29
	$x_L$ = 0.545	19.36	29
H1	$x_L$ = 0.635	19.28	29
	$x_L$ = 0.725	17.33	29
	$x_L$ = 0.815	13.23	29
	$x_L$ = 0.905	10.15	29
All data sets		130.05	203
	$x_L$ = 0.240	24.84	25
	$x_L$ = 0.310	9.68	25
	$x_L$ = 0.370	11.68	25
	$x_L$ = 0.430	48.45	25
ZEUS	$x_L$ = 0.490	21.11	25
	$x_L$ = 0.550	24.84	25
	$x_L$ = 0.610	17.74	25
	$x_L$ = 0.670	28.18	25
	$x_L$ = 0.730	6.61	25
	$x_L$ = 0.790	3.02	25
	$x_L$ = 0.850	7.88	25
	$x_L$ = 0.920	15.02	25
All data sets		219.10	300

Table 3. Table 3: Parameter values { p i } subscript 𝑝 𝑖 \{p_{i}\} for SKTJ17 QCD analysis at the input scale Q 2 0 superscript subscript absent 0 2 {}_{0}^{2} = 1 GeV 2 from the combined H1 and ZEUS data sets. The values without errors have been fixed after the first minimization since the data do not constrain these unknown parameters well enough. The details of the χ 2 superscript 𝜒 2 \chi^{2} analysis and the constraints applied to control the neutron FFs parameters are contained in the text.

Parameters	$\beta M_{\Sigma /P}^N(\beta ,Q_0^2,x_L)$	$p_i\pm \delta p_i$	$\beta M_{g/P}^N(\beta ,Q_0^2,x_L)$	$p_i\pm \delta p_i$
$a$	$a_q$	$0.116\pm 0.031$	$a_g$	${0.0}^{\ast }$
$b$	$b_q$	${0.260}^{\ast }$	$b_g$	${4.884}^{\ast }$
$c$	$c_q$	${0.523}^{\ast }$	$c_g$	${9.969}^{\ast }$
$𝒩$	${𝒩}_q$	$0.245\pm 0.023$	${𝒩}_g$	$0.130\pm 0.027$
$A$	$A_q$	${0.0}^{\ast }$	$A_g$	${0.201}^{\ast }$
$B$	$B_q$	$1.430\pm 0.092$	$B_g$	$1.740\pm 0.117$
$C$	$C_q$	$12.071\pm 2.270$	$C_g$	${29.865}^{\ast }$
$D$	$D_q$	$5.307\pm 0.390$	$D_g$	$6.733\pm 0.646$

Equations

$$Q^2=-q^2,x=\frac{Q^2}{2p.q},y=\frac{p.q}{p.k},$$

$$\frac{d^4\sigma (ep\to e^{\prime }BX)}{d\beta d{Q}^{2}d{x}_{L}dt}$$

$$x_L\simeq \frac{E_B}{E_p},\beta =\frac{x}{1-x_L},$$

$$\frac{d^3\sigma (ep\to e^{\prime }BX)}{d\beta d{Q}^{2}d{x}_L}$$

$$t_{min}$$

$$t_0$$

$${\sigma }_r^{LN(3)}$$

$$Q^2\frac{\partial M_{\Sigma /P}^B(\beta ,Q^2,x_L)}{\partial Q^2}$$

$$Q^2\frac{\partial M_{g/P}^B(\beta ,Q^2,x_L)}{\partial Q^2}$$

$$\beta M_{\Sigma /P}^N(\beta ,Q_0^2,x_L)$$

$$\beta M_{g/P}^N(\beta ,Q_0^2,x_L)$$

$${\mathcal{A}}_q(x_L)={\mathcal{N}}_q{x}_L^{A_q}(1-x_L{)}^{B_q}(1+C_q{x}_L^{D_q}),$$

$${\mathcal{A}}_g(x_L)={\mathcal{N}}_g{x}_L^{A_g}(1-x_L{)}^{B_g}(1+C_g{x}_L^{D_g}).$$

$$\frac{d{\sigma }^{{\gamma }^{\ast }p\to Xn}}{d{p}_T^2}\propto e^{-b(x_L)p_T^2},$$

$${\chi }^2(\{p_i\})=\sum _{n=1}^{N^{\exp }}\sum _{j=1}^{N_n^{\mathrm{data}}}{w}_j\frac{\text{leavevmode}({\mathrm{Data}}_j-{\mathrm{Theory}}_j(\{p_i\})\text{leavevmode}{)}^2}{\delta {\mathrm{Data}}_j},$$

$${\chi }_{\mathrm{global}}^2(\{p_i\})=\sum _{n=1}^{N^{\exp }}{w}_n{\chi }_n^2,$$

$${\chi }_n^2(\{p_i\})$$

$$\Delta \beta M(\beta ,Q^2,x_L)$$

$$\Delta {\chi }_{\mathrm{global}}^2$$

$$H_{ij}=\frac{1}{2}\frac{\partial^{2}\chi^{2}_{\rm global}}{\partial p_{i}\,\partial p_{j}}\Bigg{|}_{0}\,.$$

$$\sum _{j=1}^n{C}_{ij}{\upsilon }_{jk}={\lambda }_k{\upsilon }_{ik},$$

$$p_i-p_i^0=\sum _{k=1}^n{e}_{ik}{z}_k.$$

$$\Delta {\chi }_{\mathrm{global}}^2={\chi }_{\mathrm{global}}^2(\{p\})-{\chi }_0^2(\{p^0\})=\sum _{k=1}^n{z}_k^2.$$

$$\sum _{k=1}^n{z}_k^2\le T^2.$$

$$p_i(s_k^{\pm })=p_i^0\pm t\sqrt{{\lambda }_k}{v}_{ik},$$

$$\Delta \mathcal{O}=\frac{1}{2}{\left[\sum _{k=1}^n(\mathcal{O}(s_k^{+})-\mathcal{O}(s_k^{-}){)}^2\right]}^{\frac{1}{2}}.$$

$$\Delta \mathcal{O}={\left[\Delta {\chi }_{\mathrm{global}}^2\sum _{i,j=1}^k\frac{\partial \mathcal{O}}{\partial p_i}{C}_{ij}\frac{\partial \mathcal{O}}{\partial p_j}\right]}^{\frac{1}{2}}.$$

$$P_N({\chi }^2)=\frac{({\chi }^2{)}^{\frac{N}{2}-1}}{2^{\frac{N}{2}}\Gamma (\frac{N}{2})}{e}^{\frac{-{\chi }^2}{2}}\text{leavevmode},$$

$${\int }_0^{\Delta {\chi }^2}{P}_N({\chi }^2)d{\chi }^2=P(\approx 0.68)\text{leavevmode},$$

Full text

1. I Introduction

2. II Theoretical framework

3. III NLO QCD analysis of neutron FFs and parameterization

4. IV Leading neutron production data

5. IV.1 H1 data

6. IV.2 ZEUS data

7. V The method of $\chi^{2}$ minimization and neutron FFs uncertainties

8. V.1 $\chi^{2}$ minimization

9. V.2 Neutron FFs uncertainties

10. VI Results and discussions

11. VI.1 SKTJ17 neutron FFs and their uncertainties

12. VI.2 Comparison to leading neutron data

13. VII Leading-baryons production at the LHC

14. VIII Summary and Conclusion

I Introduction

Over the past decade, our knowledge of the quark and gluon substructure of the nucleon has been extensively improved due to the high-energy scattering data from fixed target experiments, the precise data from electron-proton collider HERA Abramowicz:2016ztw ; Abt:2016zth ; Abramowicz:2015mha ; Abramowicz:2014jak ; Andreev:2013vha ; Aaron:2009aa ; Aaron:2009kv ; Aaron:2009bp , the data from high energy proton-proton scattering at the Tevatron TevatronElectroweakWorkingGroup:2016lid ; Abazov:2016kfr ; Abazov:2016tba ; Aaltonen:2016pos ; D0:2016ull ; Aaltonen:2008eq ; Abulencia:2005yg ; Abazov:2008ae ; Abulencia:2007ez ; Abbott:2000ew and up-to-date data from LHC CMS:2017tvp ; Sirunyan:2017azo ; Aad:2013lpa ; Aad:2011fc ; Chatrchyan:2013tia ; Chatrchyan:2011wt ; Aad:2013iua ; Aaij:2012mda ; Aaij:2015vua ; Chatrchyan:2011jz ; Aad:2011dm ; Chatrchyan:2012xt . Deep inelastic scattering (DIS) data as well as data from hadron colliders has been successfully used in many Global QCD analyses to extract the unpolarized parton distribution functions (PDFs) Bourrely:2015kla ; Ball:2012cx ; Harland-Lang:2014zoa ; Ball:2014uwa ; Martin:2009iq ; Gao:2013xoa ; Alekhin:2012ig , polarized PDFs Shahri:2016uzl ; Ball:2013lla ; Jimenez-Delgado:2014xza ; Sato:2016tuz ; Nocera:2014uea ; Leader:2014uua ; Nocera:2014gqa , nuclear PDFs Khanpour:2016pph ; Eskola:2016oht ; Kovarik:2015cma ; Klasen:2017kwb ; Wang:2016mzo ; deFlorian:2011fp , and related studies Goharipour:2017rjl ; Carlson:2017gpk ; Accardi:2011fa ; Owens:2012bv ; Monahan:2016bvm ; Dahiya:2016wjf ; Jimenez-Delgado:2013boa ; Nocera:2016zyg ; Ball:2016spl ; Goharipour:2017uic ; Ru:2016wfx ; Haider:2016zrk ; Accardi:2016qay ; Armesto:2015lrg ; Frankfurt:2015cwa ; Guzey:2013xba ; Frankfurt:2016qca ; Guzey:2016qwo ; Accardi:2011mz ; Frankfurt:2011cs ; Alekhin:2017kpj . Beside the mentioned data sets, the production of leading neutron and proton in deep inelastic scattering (DIS) opens a new window for the theory of strong interaction in the soft region and provides a probe of the relationship between QCD of quarks and gluons and the strong interaction of hadrons Aaron:2010ab ; Chekanov:2002pf . Consequently, in the framework of perturbative QCD (pQCD), the study of leading-baryon production represents an important field of investigation. In the leading-baryon productions, $ep\to e^{\prime}BX$, the energetic neutron or proton which are produced in the fragmentation of the proton remnant, carry a large fraction $x_{L}$ of the longitudinal momentum of the incoming proton  Aaron:2010ab ; Chekanov:2002pf ; Adloff:1998yg ; Breitweg:2000nk ; Chekanov:2004wn ; Aktas:2004gi ; Chekanov:2004dk ; Chekanov:2007tv ; Chekanov:2008tn ; Chekanov:2002yh ; Andreev:2014zka . These events are measured at small polar angle with respect to the collision axis.

However, due to difficulty of detecting the leading-baryon in high energy physics experiment, the data available are scarce. More recently the H1 and ZEUS collaborations at HERA have measured events, in which neutron is produced in the forward region, obtain sizable contributions of leading neutrons to the DIS cross sections, $\sim 8-10$% Aaron:2010ab ; Chekanov:2002pf . These kinds of processes open a new window to study hard processes in a new kinematical region to obtain information on soft quantum chromodynamics (QCD) dynamics. In that case, the measurements of leading-baryon structure functions can be used as a test of new aspects of QCD. Along with these experimental developments, the fracture functions approach has been developed in the framework of perturbative QCD in order to deal with such kind of forward processes Trentadue:1993ka ; Ceccopieri:2007th ; Szczurek:1997cw .

Fracture functions provides a QCD-based description of semi-inclusive DIS in the target fragmentation region. The formalism of fracture functions, where the leading particles production is described in terms of structure functions of the fragmented nucleon, has been successfully used to describe forward neutron data from the H1 and ZEUS collaborations Ceccopieri:2014rpa ; deFlorian:1998rj ; deFlorian:1997wi . As for DIS structure function, QCD can not predicted the shape of fracture functions. As for parton distribution functions (PDFs), the non-perturbative neutron fracture function (neutron FFs) can be parameterized at a given initial scale Q${}_{0}^{2}$. Fracture functions, the probabilities of finding a parton and a hadron in the target, can be related to the parton distributions of the object exchanged between the initial and final states. In the production of the leading neutron in the target fragmentation region this object is $\pi^{+}$ in the $ep\to enX$ process. In the target fragmentation region, the corresponding cross sections is expressed as a convolution of the fracture functions, $M_{h/p}$, with the point like partonic cross sections. In this paper, we present the results of our QCD global analysis of recent and up-to-date experimental data for the production of neutrons in the forward direction in DIS. As we mentioned, the results obtained in this analysis, is in the framework of fracture functions by modeling the neutron FFs at the input scale, Q${}_{0}^{2}$. We propose a standard parametric form for the neutron FFs at a given initial scale Q${}_{0}^{2}=1$ GeV2 and obtain their parameters by next-to-leading order (NLO) global QCD fit to forward neutron production data measured by H1 and ZEUS collaboration at HERA. We find that our theory predictions are in satisfactory agreements with all data analyzed.

The outline of the present paper is the following. First, in Section II we present the theoretical settings of the analysis. The details of the fitting methodology applied in this work and the functional forms used to extract neutron FFs are presented in Section III. The details of the forward neutron production data from H1 and ZEUS collaboration are discussed in Section IV. Section V provides the method of the $\chi^{2}$ minimization, uncertainties estimation and error calculations. The results of present NLO neutron FFs fits and detailed comparison with available observables are discussed in Section VI. In Section VII, we briefly discuss the present and upcoming experimental data on the production of leading-baryons at LHC and at Jefferson Lab. Finally Section VIII contains the summary and conclusions.

II Theoretical framework

We can now specify the theory settings used for the neutron FFs fits presented in this work. We will restrict ourselves to a brief summary of the theoretical framework relevant for our global QCD analysis of leading neutron structure functions in which we closely follow Ref. Ceccopieri:2014rpa .

We use the NLO theory with $\alpha_{s}(M_{Z}^{2})=0.1184$ in variable flavour number scheme (VFNS) with charm and bottom masses of $m_{c}=1.41$ and $m_{b}=4.50$ GeV. In order to describe the hard scattering DIS process, we use the usual kinematic variables $x$, Q2, $y$ in which are defined as

[TABLE]

where in the DIS process, $p$ is the four-momenta of the incident proton, $k$ is the four-momenta of the incident positron and $q$ is the four-momenta of the virtual photon. The four-fold differential cross section to describe the baryon production processes $ep\to e^{\prime}BX$ can be obtained by semi-inclusive leading-baryon transverse and longitudinal structure functions, $F_{2}^{\rm LB(4)}$ and $F_{L}^{\rm LB(4)}$, which is defined as Aaron:2010ab ; Chekanov:2002pf

[TABLE]

The longitudinal momentum fraction $x_{L}$ and the scaled fractional momentum variable $\beta$ are defined by

[TABLE]

where $x$ is the Bjorken variable, $E_{p}$ is the proton beam energy and $E_{B}$ is the energy of final-state baryon. In Eq. (2), $t$ is the squared four-momentum transfer between the incident proton and the final state neutron. The $t$ integrated differential cross section can be obtained by

[TABLE]

where the integration limits are

[TABLE]

$m_{N}$ is the mass of final-state baryon, $m_{p}$ is the proton mass, and $p_{T}^{max}$ is the upper limit of the neutron transverse momentum used for the $F_{2}^{\rm LB(3)}$ measurement. For the semi-inclusive processes which have final-state proton and neutron, the structure function $F_{2,L}^{\rm LB(3)}$, is denoted by $F_{2,L}^{\rm LP(3)}$ and $F_{2,L}^{\rm LN(3)}$ respectively. In this paper which is correspond to a QCD analysis of forward neutron production, we define the reduced $e^{+}p$ cross section $\sigma_{r}^{\rm LN(3)}$ in term of leading neutron transverse $F_{2}^{\rm LN(3)}$ and the longitudinal structure functions $F_{L}^{\rm LN(3)}$ as Aaron:2010ab ; Chekanov:2002pf

[TABLE]

It is noteworthy to mention here that the leading neutron structure functions in above equations can be written in terms of neutron FFs and hard-scattering coefficients Ceccopieri:2014rpa . The Wilson coefficient functions are the same as in fully inclusive DIS Vermaseren:2005qc .

The well-known DGLAP evolution equations Dokshitzer:1977sg ; Gribov:1972ri ; Lipatov:1974qm ; Altarelli:1977zs which are a set of an integro-differential equations can be used to evolve the polarized and unpolarized parton distributions functions to an arbitrary energy scale, Q2. The solutions of these evolution equations will provide us the valance, gluon, and sea quark distributions inside the nucleon. These equations widely can be used as fundamental tools to extract the deep inelastic scattering (DIS) structure functions of proton, neutron and deuteron which enrich our current information about the structure of the hadrons. Since the scale dependence of the cross section in forward particle production in DIS can be calculated within perturbative quantum chromodynamics (pQCD) Trentadue:1993ka , consequently the neutron fracture functions also obey the standard DGLAP evolution equations Ceccopieri:2014rpa ; Camici:1998bg .

In Refs. deFlorian:1998rj ; Daleo:2003xg ; Daleo:2003jf ; Trentadue:1996jc ; Camici:1998bg ; Trentadue:1994iw have been shown that, in the phenomenological level, the fracture functions well reproduce the leading proton data, thus one can use the common perturbative QCD approach to these particular classes of semi-inclusive processes. So, like for the case of parton distributions functions (PDFs), one can use phenomenological model to describe forward neutron production and extract the neutron FFs from QCD fit to the data Ceccopieri:2014rpa ; Ceccopieri:2007th . The evolution equations of neutron FFs are easily obtained by DGLAP evolution equations Trentadue:1993ka as

[TABLE]

where $M^{B}_{\Sigma/P}(\beta,Q^{2},x_{L})$ and $M^{B}_{g/P}(\beta,Q^{2},x_{L})$ correspond to the singlet and gluon distributions, respectively Ceccopieri:2007th . These non-perturbative distributions in which hereafter indicated by “neutron FFs” need to be parametrize at an input scale, Q${}_{0}^{2}$. Their evolution to higher scale, $Q^{2}>Q_{0}^{2}$, can be described by using the evolution equation given above. $P_{\Sigma}$ and $P_{g}$ in Eq. (II) are the common NLO contributions to the splitting functions governing the evolution of unpolarized singlet and non-singlet combinations of quark densities in perturbative QCD. Splitting functions are perturbatively calculable as a power expansion in the strong coupling constant $\alpha_{s}$. The splitting functions $P_{\Sigma}$ and $P_{g}$ in Eq. (II) are the same as in fully inclusive DIS Vogt:2004mw ; Moch:2004pa ; vanNeerven:2000uj ; vanNeerven:1999ca ; Moch:2014sna

In the next sections, we give a detailed account of the first global analysis of neutron FFs performed in this study which in the following will be referred to as “SKTJ17”. We first discuss in details the parameterization of neutron FFs and then we will present data selection and the determination of the best fit, which we compare to the fitted data. We then focus on the studies of uncertainties using the standard Hessian error matrix approach.

III NLO QCD analysis of neutron FFs and parameterization

In order to obtain a parametrization for the neutron FFs, $\,\beta M^{N}_{i/P}(\beta,Q_{0}^{2},x_{L})$ with $i=\Sigma$ and $g$, at a given initial scale $Q_{0}^{2}$, we select a relatively simple functional dependence in the variables $\beta$ and $x_{L}$ with enough flexibility as to reproduce the data accurately. We assume the following general initial functional form at Q${}_{0}^{2}$ = 1 GeV2

[TABLE]

where ${\cal A}_{q}(x_{L})$ and ${\cal A}_{g}(x_{L})$ define as

[TABLE]

The label of $\,\Sigma/P$ and $\,g/P$ correspond to the singlet and gluon distributions, respectively. The $x_{L}$ dependence of the neutron FFs is encoded in ${\cal A}_{i}(x_{L})$. Since the present leading neutron data are not yet sufficient to distinguish $q(=u,d,s)$ from $\bar{q}(=\bar{u},\bar{d},\bar{s})$, we assume symmetric sea distributions throughout. We will show that these kinds of parametrizations give relatively good initial approximations to the description of the H1 and ZEUS leading neutron data sets Aaron:2010ab ; Chekanov:2002pf , however, their survival seems unlikely in a more precise analysis. The available forward neutron production data are not accurate enough to determine all the shape parameters with sufficient accuracy. Eq. (III) includes 18 free parameters in total in which we further reduce the number of free parameters in the final minimization.

The parameters $\{p_{i}\}$ representing our best global QCD fit of neutron FFs in Eq. (III), henceforth denoted as SKTJ17 are given in Table 3. A few additional remarks will be presented in Sec VI. As we mentioned, the currently available leading neutron data do not fully constrain the entire $\beta$ and $x_{L}$ dependence of $\beta M^{N}_{\Sigma/P}$ and $\beta M^{N}_{g/P}$ imposed in Eq. (III). Consequently we are forced to make some restrictions on the parameter space $\{p_{i}\}$. We will return to this issue in a separate section.

Rather than determining also the strong coupling $\alpha_{s}(Q_{0}^{2})$ in the global QCD fit along with the neutron FFs parameters, we fixed $\,\alpha_{s}(M_{Z}^{2})$ value close to the updated Particle Data Group (PDG) average. The scale dependence of $\alpha_{s}$ is normally computed by numerically solving its renormalization group equation at next-to-leading order accuracy. For the evolution we take $\,\alpha_{s}(M_{Z}^{2})=0.1184$ Olive:2016xmw ; Agashe:2014kda , and we choose to work in the variable flavor number scheme (VFNS) where charm and bottom quark distributions are radiatively generated from their corresponding thresholds Martin:2009iq ; Ball:2014uwa . In the present analysis, all quarks are treated as massless and we fixed the heavy quark masses at $m_{c}=1.41$ GeV, $m_{b}=4.50$ GeV and $m_{b}=175.0$ GeV. Our choice for the VFNS scheme is due to that for all presently available leading neutron observables, heavy quarks play a negligible role. The scale evolution equations for the neutron FFs are solved in $x$-space at next-to-leading order. Likewise, all leading neutron observables used in our QCD fit are computed consistently at next-to-leading order accuracy in the $\overline{MS}$ factorization scheme.

IV Leading neutron production data

Our first physics objective is to establish the set of neutron FFs that gives the optimum theoretical description of the available hard scattering leading neutron production data. In this section, we will present the data sets used in the present analysis. The data sets that we will use is the following: The H1 data on the leading neutron production in DIS scattering Aaron:2010ab as well as the data from leading neutron production in $e^{+}p$ collision from ZEUS collaboration Chekanov:2002pf . The detail of the data sets will be presented in the next section.

IV.1 H1 data

The semi-inclusive cross sections data for the production of leading neutron are taken during the years of 2006 and 2007 by the H1 collaboration at HERA in DIS positron-proton scattering Aaron:2010ab which is correspond to an integrate luminosity of $\cal{L}$ = 122 pb*-1*, much larger than the previous H1 measurement Adloff:1998yg . Better experimental capabilities in this measurement lead to the extension of the kinematical coverage of $x$ and $Q^{2}$ to higher values. This leading neutron structure function $F_{2}^{\rm LN(3)}$ in which has been measured by H1 experiment at HERA covers a large range of kinematics of $Q^{2}$, $\,6\leq Q^{2}\leq 100$ GeV2, and $x$, $\,1.5\times 10^{-4}\leq x\leq 3\times 10^{-2}$, for average $y$ values between $0.05$ and $0.68$, and the upper limit of the neutron transverse momentum of $p_{T}^{\rm max}<200$ MeV. The value of longitudinal momentum fraction $x_{L}$ covers the range from $0.365$ to $0.905$. In order to enhance the relative contribution of pion exchange Holtmann:1994rs ; Kopeliovich:1996iw for these selected of DIS events, the value of $p_{T}^{\rm max}$ in which used for the measurement of $F_{2}^{\rm LN(3)}$ is set to 200 MeV. Considering that the pion exchange mechanism dominates leading neutron production, these data sets can provide constraints on the shape of the pion structure function McKenney:2015xis . In Fig 1, we plot the nominal coverage of H1 data sets used in our QCD fits. The plot nicely summarizes the universal $\beta$, $x_{L}$ and $Q^{2}$ dependence of the forward neutron production at HERA.

IV.2 ZEUS data

The semi inclusive cross section for production of leading neutron measured by ZEUS collaboration are also used in our global QCD analysis. The ZEUS collaboration presented the leading neutron production cross sections for $x_{L}>0.2$ in neutral current electron-proton collisions at HERA Chekanov:2002pf . Positron and proton energies are $E_{e}=27.5$ GeV and $E_{p}=820$ GeV, respectively, correspond to a center of mass energy of $\sqrt{s}=300$ GeV. Similarly to H1 experiment, extensive range of kinematics was covered by the ZEUS data, for $1.1\times 10^{-4}\leq x\leq 3.2\times 10^{-2}$ from photoproduction up to $Q^{2}\sim 10^{4}$ GeV2, with $0<y<0.8$ and neutron scattering angle $\theta_{n}<0.8$ mrad. The HERA magnet apertures limit the FNC (Forward Neutron Calorimeter) acceptance to neutron with the production angle less than $\theta_{n}^{\rm max}=0.8$ mrad, which is corresponding to the transverse momenta of $p_{T}^{\rm max}=E_{n}\theta_{n}^{\rm max}=0.656x_{L}$ GeV. As we already mentioned in the previous section, the distribution of the neutron for H1 data is integrated only up to $p_{T}^{\rm max}=200$ MeV, so the H1 and ZEUS data can only be used in the analysis for the longitudinal momentum fraction of $x_{L}=0.3$, which is correspond to $p_{T}^{\rm max}=197$ MeV. For higher values of $x_{L}$, the ZEUS data should be scaled to account for the smaller $p_{T}$ range measured by H1 collaboration. Of course, in general we would like to maximize the $\beta$, $x_{L}$ and $Q^{2}$ coverages included in the analysis in order to increase the statistics of our fit. Therefore, we have scaled down the ZEUS data to the H1 $p_{T}$-range by using the form of $p_{T}^{2}$ distribution for the fixed values of $x_{L}$ as Chekanov:2002pf :

[TABLE]

where $\sigma^{\gamma^{*}p\to Xn}$ is the the virtual photon-proton cross section for the process $\gamma^{*}p\to Xn$. The slope $b(x_{L})$ can be parameterised as $b(x_{L})=(16.3\ x_{L}-4.25)$ GeV*-2* which is in reasonable accord with the data Chekanov:2002pf ; Chekanov:2002yh . In order to reduce the systematic uncertainties, ZEUS collaboration is measured the neutron-tagged cross section $ep\to e^{\prime}Xn$ relative to the inclusive DIS cross section $ep\to e^{\prime}X$. Considering this ratio as well as proton structure function, one can obtain the $F_{2}^{\rm LN(3)}$ values for various bin of $x$, Q2 and $y$. The kinematic range of ZEUS forward neutron data are shown in Fig. 2. We should notice here that the H1 leading neutron data were collected during the 2006-2007 run by an integrated luminosity about 3 times that of the ZEUS data in the DIS region. In consequence, the statistical uncertainties of the H1 data are much smaller than those for the ZEUS leading neutron spectra.

The measured leading neutron production data points above $Q^{2}$ = 1.0 GeV2 used in the SKTJ17 global analysis are listed in Table 1. For each data set we provide the corresponding references, the kinematical coverage of $x_{L}$, $x_{B}$ and $Q^{2}$, the number of data points and the fitted normalization shifts ${\cal{N}}_{n}$.

V The method of $\chi^{2}$ minimization and neutron FFs uncertainties

In this Section, we outline the details of SKTJ17 analysis. More specifically, we discuss the selection of data sets, treatment of experimental normalization uncertainties, as well as the determination of the parameters by global $\chi^{2}$ minimization. We also briefly present the details of the Hessian matrix method for estimating uncertainties. As we noted before, we have performed a careful estimation of the uncertainties using the “Hessian method”. An advantage of the Hessian technique is that it allows us to produce sets of eigenvector PDFs, which can be straightforwardly used in computations of other observables such as reduced $e^{+}p$ cross section $\sigma_{r}^{\rm LN(3)}$ as well as leading neutron structure function $F_{2}^{\rm LN(3)}$.

V.1 $\chi^{2}$ minimization

Global QCD extractions of PDFs, nuclear PDFs as well as polarized PDFs are implemented around an effective $\chi^{2}$ function that quantifies the goodness of the fit to data for a given set of theoretical parameters in which determines the PDFs at some input scale Q${}_{0}^{2}$. In order to search for optimum PDFs by minimization, the simplest $\chi^{2}$ function is usually given by

[TABLE]

The simple form of $\chi^{2}(\{p_{i}\})$ presented above is appropriate only in the ideal case of data sets with uncorrelated errors. Since most experiments come with additional information on the fully correlated normalization uncertainty $\Delta{\cal N}_{n}$, Eq. (11) need to be modified in order to account for such normalization uncertainties. In order to determine the best fit parameters of Eq. (III), we need to minimize the $\chi^{2}_{\rm global}$ function with the free unknown parameters. $\chi_{\rm global}^{2}(\{p_{i}\})$ quantifies the goodness of fit to the data for a set of independent parameters $(\{p_{i}\})$ that specifies the neutron FFs at the input scale Q${}_{0}^{2}$ = 1 GeV2. This function is expressed as,

[TABLE]

where $w_{n}$ is a weight factor for the $n^{\rm th}$ experiment and

[TABLE]

where $N^{\rm exp}$ correspond to the individual experimental data sets and $N^{\rm data}_{n}$ correspond to the number of data points in each data set. The normalization factors $\Delta{\cal N}_{n}$ in Eq. (V.1) can be fitted along with the fitted parameters $(\{p_{i}\})$.

The $\chi_{\mathrm{\rm global}}^{2}(\rm p)$ function is minimized by the CERN program library MINUIT James:1994vla . From the $\chi_{\rm global}^{2}(\{p_{i}\})$ analysis, an error matrix that is the inverse of a Hessian matrix is obtained. In order to determine the sensitivity of the fit to different values of $x_{L}$ collected by H1 and ZEUS collaborations at HERA, we compute the $\chi^{2}$ values for each data sets. The data sets included in SKTJ17 analysis are listed in Table 2, together with the $\chi^{2}$ values, defined in Eq.(11), corresponding to each individual data set for each of $x_{L}$. This suggests that reasonable fits to the leading neutron cross sections can be obtained within most of the $x_{L}$ values. More detailed discussion of the description of the individual data sets has been given in Section IV.

V.2 Neutron FFs uncertainties

As in the case of standard PDFs, the evolved leading neutron fracture functions are linear functions of the input densities. Let $M(\beta,Q^{2},x_{L};p_{i}|_{i=1}^{k})$ be the evolved neutron FFs at $Q^{2}$ depending on the parameters $p_{i}|_{i=1}^{k}$. Then its correlated error as given by Gaussian error propagation as Bluemlein:2002be

[TABLE]

where $C(p_{i},p_{j})$ are the elements of the covariance matrix obtained in the QCD fit procedure at the input scale $Q_{0}^{2}$. The covariance matrix can be used at any scale of $Q^{2}>Q_{0}^{2}$. The gradients $\partial\beta M/\partial p_{i}$ at this scale can be calculated analytically. Their value at $Q^{2}$ is then calculated by evolution in $x$ space and are used according to Eq. (V.2).

In addition to the method presented above, one can also determine the uncertainties of obtained neutron FFs via well-known Hessian method and diagonalize the covariance matrix and work in terms of the eigenvectors and eigenvalues. Here, we briefly review the important points for studying the neighborhood of $\chi^{2}_{0}$. The basic procedure is provided in Refs. deFlorian:2011fp ; Hou:2016sho ; MoosaviNejad:2016ebo ; Khanpour:2016uxh ; Khanpour:2016pph ; Martin:2002aw ; Martin:2009iq ; Pumplin:2001ct .

As we have mentioned earlier, one can find the appropriate parameter set in which minimize the $\chi^{2}_{\rm global}$ function. We call this neutron FFs set $S_{0}$. The parameters value of $S_{0}$, i.e. {$p_{1}^{0}\dots p_{n}^{0}$}, in which extracted from QCD fit to H1 and ZEUS leading neutron data, will be presented in Sec. VI. As we will mention latter, we simply fix some of the parameters of our input functional from presented in Eq. (III) at their best-fit values, so that the Hessian matrix only depends on a subset of parameters.

By moving away the parameters from their obtained values, $\chi^{2}$ increases by the amount of $\Delta\chi^{2}$

[TABLE]

where $H_{ij}$ is the Hessian matrix which defined as

[TABLE]

Now it is convenient to work in term of the eigenvalues and their corresponding orthogonal eigenvectors of covariance matrix. It is given by

[TABLE]

and we should notice here that $C_{ij}\equiv H_{ij}^{-1}$ is the error (or covariance) matrix. The displacement of the parameter {$p_{i}$} from its obtained minimum $p_{i}^{0}$ can be expressed in terms of the rescaled eigenvectors $e_{ik}=\sqrt{\lambda_{k}}\,v_{ik}$, that is

[TABLE]

Considering the orthogonality of eigenvectors $\upsilon_{ik}$ and putting Eq. (18) in (15), one can write

[TABLE]

The relevant neighborhood of $\chi^{2}$ is the interior of hypersphere with radius $T$. This means that

[TABLE]

Finally the neighborhood parameters can be written as

[TABLE]

with $s_{k}$ is the $k^{th}$ set of neutron FFs, $t$ adapted to make the desired $T^{2}=\Delta\chi^{2}_{\rm global}$ which is the tolerance for the required confidence interval (C.L.) and $t=T$ in the quadratic approximation.

Using the method we mentioned above, we accompany the construction of the QCD fit by reliable estimation of uncertainty. Finally uncertainties of any observables ${\cal O}$, which can be the neutron FFs or reduced cross sections in our case, in the Hessian method can calculate as Martin:2002aw ; Martin:2009iq ; Pumplin:2001ct

[TABLE]

In above equation, ${\cal O}(s_{k}^{+})$ and ${\cal O}(s_{k}^{-})$ are the value of ${\cal O}$ extracted from the input set of parameters $p_{i}(s_{k}^{\pm})$ obtained from Eq. (21). In this paper, we follow the standard Hessian method to calculate the neutron FFs error band as well as the corresponding observables such as the reduced cross sections. The evolved neutron FFs are attributive functions of the input parameters obtained in the QCD fit procedure at the scale $Q_{0}^{2}$, then their uncertainty can be written applying the standard Hessian method

[TABLE]

The $\Delta\chi^{2}$ values determine the confidence region, and it is calculated so that the confidence level (C.L.) $P$ becomes the one-$\sigma$-error range ($P=0.68$) for a given number of parameters ($p_{i=N}$) by assuming the normal distribution in the multi-parameter space. Since the neutron FFs are provided with many parameters, so that the $\Delta\chi^{2}_{\rm global}$ value should be calculated.

Assuming correspondence between the confidence level (C.L.) of a normal distribution in multi-parameter space and the one of a $\chi^{2}$ distribution with $N$ degree of freedom, one can define the probability density function as

[TABLE]

then the confidence level $P$ can be obtain as

[TABLE]

and similarly for the 90th percentile we have $P=0.90$. The parameter number in our analysis is eight ($N=8$), and it leads to $\Delta\chi^{2}$ = 9.27. The uncertainty of a neutron FFs with respect to the optimized parameters $p_{i=N}$ is then calculated using Eq. (23) by using Hessian matrices and assuming mentioned linear error propagation. For the neutron FFs uncertainty estimation, one can analytically calculate the gradient terms in Eq. (23) at the initial scale $Q_{0}^{2}$ = 1 GeV2. For the estimation at arbitrary Q2, each gradient term is evolved by the DGLAP evolution kernel, and then the neutron FFs uncertainties as well as the uncertainties for any other observables such as cross sections are calculated. Here we calculate the neutron FFs uncertainty with $\Delta\chi^{2}=1$ and $9.27$ which is the most appropriate choice. The Hessian method discussed in the present analysis has been used for estimating TKAA16 NNLO polarized PDFs Shahri:2016uzl as well as KA15 nuclear PDFs analysis Khanpour:2016pph . The details of the uncertainty analysis are discussed in details in Refs. Martin:2002aw ; Martin:2009iq ; Pumplin:2001ct .

VI Results and discussions

We are in a position to describe the details and all techniques we used for the parametrizations of neutron FFs in SKTJ17 global analysis. The minimization is carried out with respect to the set of parameters in Eq. (III), $\{p_{i}\}=\{a_{i},b_{i},c_{i},{\cal N}_{i},A_{i},B_{i},C_{i},D_{i}\}$. The neutron FFs are evolved to the scales Q2 > Q${}_{0}^{2}$ relevant in experiment. Like for the case of PDFs parameterization, particular functional form and the value for Q${}_{0}^{2}$ are not too crucial. The parameterization at the input scale should be flexible enough to accommodate all DIS data within their ranges of uncertainties. As we mentioned, our input distributions in Eq. (III) follow the standard form used in fits to DIS data. In addition to our much more flexible input parametrization presented in Eq. (III), we have repeated our QCD fit with alternative parametrizations, some of them even more flexible than the one we choose. For example, we have also included $\sqrt{x}$ terms, both for the singlet and gluon distributions, even allowing the fit to vary them. We have found no significant improvement in the quality of the fit to data or changes of the uncertainty bands. This indicates that the present H1 and ZEUS leading neutron production data is not really able to discriminate between various forms of the input distributions.

As will be seen from our results presented in this section, we found that our input distributions in Eq. (III) could be considered as good parametrizations to the leading neutron production experimental data.

The parameter values $\{p_{i}\}$ of the next-to-leading order input neutron FFs at Q${}_{0}^{2}$ = 1 GeV2 obtained from the best fit to the combined H1 and ZEUS leading neutron data sets are presented in Table 3.

Parameters marked with (∗) are fixed. This is due to that these parameters are only very weakly determined by the fit, consequently we fixed them to their preferred values. For the sea quark density we set $A_{q}$ to [math] and for the gluon density we set $a_{g}$ to [math] in Eq. (III). These only marginally limit the freedom in the functional form. We found that singlet small-$x_{L}$ coefficient $A_{q}$ as well as gluon small-$\beta$ coefficient $a_{g}$ is determined with rather large error and also compatible with zero, so that we fixed them to these values. These are because there are no enough data sensitive to smaller values of $\beta$ and $x_{L}$. Moreover we found that the factor $(1+c_{i}\,\beta)$ in SKTJ17 parametrization provides flexibility to obtain a good description of the data. Thus, we will make use of the $c_{i}$ coefficients for the sea quark and gluon densities. The parameters $B_{q}$ and $D_{q}$ always came out close to $B_{g}$ and $D_{g}$, so one can set them equal. In order to let enough flexibility to the sea quark and gluon densities, we prefer them to vary differently in the QCD fit. In total this leaves us with 8 free parameters in the SKTJ17 QCD fit, (5 for sea quarks and 3 for the gluon density), which we include later on also in our uncertainty estimates. We also tried to relax the imposed constraints discussed above, but found that present leading neutron data are not really sensitive to them. We find ${\chi}^{2}/{\rm{d.o.f.}}=349.16/495=0.705$ which yields an acceptable fit to the experimental data.

VI.1 SKTJ17 neutron FFs and their uncertainties

Our newly obtained singlet and gluon momentum distributions at the input scale Q${}_{0}^{2}$ = 1 GeV2 are shown in Figs. 3 and 4 along with estimates of their uncertainties using the Hessian methods for a tolerance of $\Delta\chi^{2}=1$ and $9.27$. The results presented for three representative bins of $x_{L}$ = 0.24, 0.55 and 0.92. The inner error band is obtained with the standard “parameter-fitting” criterion, by the choice of tolerance T = $\Delta\chi^{2}=1$ for the 68% (one-sigma) confidence level (C.L.) limit while the outer one is obtained with the choice of tolerance T = $\Delta\chi^{2}=9.27$ using Eq. 25. The main conclusion that can be drawn about the gluon and singlet distributions from SKTJ17 analysis is that the distributions are important at large $\beta$. As was stated earlier in Sec.VI, their behavior cannot be precisely determined yet from the available leading neutron production data. In particular, the behavior of the exponent of the $(1-\beta)$ factors in the parametrization, $b_{q}$ and $b_{g}$.

VI.2 Comparison to leading neutron data

In order to check the reliability of the distributions obtained in our analysis, in the following we compare results obtained using our best parametrization in Eq. (III) with the leading neutron production data sets presented by the H1 and ZEUS collaboration in which have been included in the SKTJ17 fit. In Figs. 5, SKTJ17 theory predictions for the reduced cross section $\sigma_{r}^{\rm LN(3)}(\beta,Q^{2},x_{L})$ are plotted as a function of $\beta$ for some selected values of Q2. For better description of the fit quality for different region of $x_{L}$, three representative bins of $x_{L}$ = 0.365, 0.55 and 0.725 are shown. The reduced cross section $\sigma_{r}^{\rm LN(3)}(\beta,Q^{2},x_{L})$ is scaled by a factor of $3^{i}$ for better visibility in the plots. In order to see the fit quality, the leading neutron production data from H1 and ZEUS collaborations Aaron:2010ab ; Chekanov:2002pf also added to these plots. From the figures it is clear that SKTJ17 QCD fit based on hard-scattering formula in Eq. (6) together with the neutron FFs initial conditions in Eq. (III) are in acceptable agreement with the H1 and ZEUS data. The plots also show that our results describe the data well, down to the lowest accessible value of Q2 as well as for different region of $x_{L}$.

In order to study the scale dependence of H1 and ZEUS leading neutron data, we have plotted the reduced cross sections $\sigma_{r}^{\rm LN(3)}(\beta,Q^{2},x_{L})$ as a function of Q2 in Fig. 6 for some selected values of $\beta$ and for three representative bins of $x_{L}$ = 0.365, 0.550 and 0.725. The reduced cross section $\sigma_{r}^{\rm LN(3)}(\beta,Q^{2},x_{L})$ is scaled by a factor of $0.08i$ for better visibility in the plots. One can conclude our results show that the scale dependence induced by the evolution equations of Eq. (II) is perfectly consistent with the leading neutron production data. The results clearly show that one can use the fracture functions approach to describe semi-inclusive hard processes in perturbative QCD at the kinematic region covered by electron-proton collider HERA and hadron colliders.

In Fig. 7, our theory predictions for the reduced cross sections $\sigma_{r}^{\rm LN(3)}(\beta,Q^{2},x_{L})$ shown as a function of $\beta$. The H1 (ZEUS) data correspond to Q2 = 7.3 (7.0) GeV2, and $x_{L}$= 0.365 (0.370) in the left panel and $x_{L}$= 0.725 (0.730) in the right panel. As shown in the plots, we obtain remarkable agreement with the data in the common $x_{L}$ and $\beta$ range. The plots also clearly show that our approach based on the fracture functions formalism allow a unified description of leading neutron deep inelastic cross sections.

For completeness, we finally show SKTJ17 theory predictions as a function of Q2 for the reduced cross sections $\sigma_{r}^{\rm LN(3)}(\beta,Q^{2},x_{L})$ with a representative selection of H1 and ZEUS data in Fig. 8. In the right panel, the H1 (ZEUS) data corresponds to $\beta=7.29\,(7.77)\times 10^{-3}$ and $x_{L}=0.365\,(0.370)$. The H1 (ZEUS) data in the left panel corresponds to $\beta=1.02\,(1.08)\times 10^{-2}$ and $x_{L}=0.545\,(0.550)$. The results demonstrate that SKTJ17 theory predictions can provide good description of the HERA leading neutron spectra at all kinematics.

In this section, we turned to present our perturbative predictions for the reduced cross section and detail comparison with the available leading neutron production data. Summarizing, our analysis provided a good description of H1 and ZEUS data for leading neutron production in DIS, as a function of $\beta$, $Q^{2}$ and $x_{L}$. The analysis results presented in this section enabled us to established the models and parameters which are best able to well describe the existing leading neutron production data from H1 and ZEUS collaborations. In spite of the fact that excellent descriptions of the H1 and ZEUS leading neutron spectra are obtained over the entire range of $\beta$, $x_{L}$ and Q2 covered by the data, new data could enable further constraints on the extracted neutron FFs. The success of the SKTJ17 global analysis performed here, stands for an explicit check of the pQCD framework in the fracture functions approach for the description of the leading neutron production processes.

VII Leading-baryons production at the LHC

Let us here conclude by listing some further possible developments of the present framework as well as experimental efforts. One of the important goals in high energy particle physics is to understand the production of leading-baryons which have large fractional longitudinal momentum $x_{L}\geq 0.3$. Recent measurements of leading proton and neutron spectra in electron-proton collisions by H1 and ZEUS collaborations at HERA have open a new window on this subject. Very recently, H1 collaboration at HERA has been measured for the first time the photoproduction cross section for exclusive $\rho^{0}$ production associated with a leading neutron H1:2015bxa . Since there is no hard scale presented in exclusive $\rho^{0}$ production, one can use a phenomenological approach such as Regge theory or color dipole formalism, to describe these kind of reactions Ducati:2016jdg ; Lebiedowicz:2016ryp ; Goncalves:2015mbf .

Nowadays, our understanding of the hadron structure as well as the QCD dynamics have advanced with the successful operation and precise data at HERA collider. In addition to the HERA collider, the next generation of high energy and high luminosity electron-proton colliders, such as Large Hadron Electron Collider (LHeC) AbelleiraFernandez:2012ni ; AbelleiraFernandez:2012ty ; Acar:2015cva as well as Future Circular Hadron-Electron Collider (FCC-he) Acar:2016rde which are proposed to build on the same site with LHC, could help to study the leading-baryon processes.

Another possibility is the use of the hadronic colliders. One of the important issue which have strong implications in the forward physics at hadron colliders, is the understanding of the leading neutron processes. A very rich program at the Large Hadron Collider (LHC) is being pursued in forward physics with sufficient experimental information N.Cartiglia:2015gve ; Albrow:2006xt . Finally, the upcoming experiment at Jefferson Lab (JLAB) plans to take data on the production of leading protons in the $en\to epX$ process Goncalves:2016uhj ; TDIS-JLAB ; R.A.Montgomery:2017hab . With the help of more and precise upcoming experimental data on such processes, a new era for theoretical understanding of strong interactions in the soft, non-perturbative regime will be open Goncalves:2016uhj .

VIII Summary and Conclusion

In the recent years, several dedicated experiments at the electron-proton collider HERA have collected high-precision data on the spectrum of leading-baryons carrying a large fraction of the proton’s energy. However the experimental information on leading-baryons production in lepton DIS, $ep\to e^{\prime}BX$, is still rather scarce. In addition to these experimental efforts, much successful phenomenology has been developed in understanding the mechanism of leading-baryon productions. The presence of a leading-baryon in the final state of lepton DIS provides valuable information on the relationship between the soft and hard aspects of the strong interaction.

In this work we have presented SKTJ17 NLO QCD analysis of neutron FFs using available and up-to-date data from forward neutron production at HERA Aaron:2010ab ; Chekanov:2002pf . It is shown that an approach based on the fracture functions formalism allows us phenomenologically parametrize the neutron FFs. We also have shown that a standard simple parametric form for this function gives a very accurate description of the available leading neutron production data. Finally, one can conclude that our obtained results based on the fracture function approach agree well with the scale dependence of the leading neutron production data. Completing such a picture is crucial as hadron colliders enter an era of new-generation of experimental data capable of testing this formalism. In order to asses the uncertainties in the resulting neutron FFs and the corresponding observables, associated with the uncertainties in the data, we have made an extensive use of the Hessian method.

A FORTRAN package containing SKTJ17 neutron FFs parameterization as well as the corresponding error set can be obtained via Email from the authors upon request. This FORTRAN package also includes an example program to illustrate the use of the routines.

Acknowledgments

The authors are especially grateful Garry Levman, Katja Kruger, Stefan Schmitt, Mohsen Khakzad from H1 and ZEUS collaborations for many useful discussions and comments. We are also thankful to Federico Ceccopieri, Fatemeh Arbabifar, Muhammad Goharipour and S. Atashbar Tehrani for numerous informative discussions. Hamzeh Khanpour is indebted the University of Science and Technology of Mazandaran and the School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM), to support financially this project. Fatemeh Taghavi-Shahri and Kurosh Javidan also acknowledge Ferdowsi University of Mashhad.