  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$.
