  The last two decades have seen a growing trend towards the experimental
efforts at the electron-proton collider HERA in which have collected high
precision data on the spectrum of leading neutron (LN) and leading-proton (LP)
carrying a large fraction of the proton's energy. In our recent study [Phys.
Rev. D \textbf{95} (2017), 074011], we have proposed an approach based on the
Fractures Functions (FF) formalism and have extracted the neutron Fracture
Functions (neutron FFs) from a global QCD analysis of LN production data
measured by H1 and ZEUS collaborations at HERA. We have shown that considering
the approach based on the framework of Fracture Functions, one could
phenomenologically parametrize the neutron FFs at the input scale. In order to
access the uncertainties for the obtained neutron FFs as well as the LN
structure functions and cross section, associated with the uncertainties in the
data, we have made an extensive use of the "Hessian method". Our theory
predictions based on the obtained neutron FFs are in satisfactory agreement
with all LN data analyzed, for a wide range of $\beta$ and $x_L$.
