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 nonperturbative part of the fragmentation process at the input scale Q20, 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 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 xL as well as the scaled fractional momentum variable β.