Phenomenology of diffractive DIS in the framework of fracture functions and determination of diffractive parton distribution functions
Hamzeh Khanpour







Hamzeh.Khanpour@mail.ipm.ir



Department of Physics, University of Science and Technology of Mazandaran, P.O.Box 48518-78195, Behshahr, Iran    
 
School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM), P.O.Box 19395-5531, Tehran, Iran 










The goal of this study is to determine a set of diffractive parton distribution function (diffractive PDFs) from a QCD analysis of all available and up-to-date diffractive deep inelastic scattering (diffractive DIS) data sets from HERA ep collider, including the most recent H1 and ZUES combined inclusive diffractive cross section
measurements. This extraction of diffractive PDFs, referred to as HK19-DPDF, is performed at next-to-leading (NLO) and next-to-next-to-leading (NNLO) in perturbative QCD. This new determination of diffractive PDFs is based on the fracture functions methodology, a QCD framework designed to provide a statistically sound representation of diffractive DIS processes. Heavy quark contributions to the diffractive DIS are considered within the framework of the FONLL general mass variable flavor number scheme (GM-VFNS) and the “Hessian approach” is used to determine the uncertainties of diffractive PDFs.
We discuss the novel aspects of the approach used in the present analysis, namely an optimized and flexible parametrization of the diffractive PDFs as well as a strategy based on the fully factorization theorem for diffractive hard processes.
We then present the diffractive PDFs, and discuss the fit quality and the stability upon variations of the kinematic cuts and the fitted data sets. We find that the systematic
inclusion of higher-order QCD corrections could improves the description of the data. We compare the extracted sets of diffractive PDFs based on the fracture functions approach to other recent sets of diffractive PDFs, finding in general very good agreements.








 Phenomenology of diffractive DIS in the framework of fracture functions and determination of diffractive parton distribution functions  
    Hamzeh Khanpour
    December 30, 2023
=========================================================================================================================================








§ INTRODUCTION




In the past decades, there has been an increasing interest in attempting to understand the structure of hadron <cit.>.
The high energy processes of interest which contain information on the hadron structure in Quantum Chromodynamics (QCD) mostly include the hadron production in in lepton-nucleon (ℓ p) deep-inelastic scattering (DIS) and in proton-proton (pp) collisions. Information from hadron collisions, especially at the large hadron collider (LHC) at CERN, is particularly useful in order to achieve a precise understanding of non-perturbative QCD dynamics.
Since its start of data taking, the H1 and ZEUS experiments at HERA-I and HERA-II have provided an impressive wealth of information on the
quark and gluon structure of the proton, and hence, a considerable amount of literature has been published over the past decade. Indeed, modern global analyses of parton distribution functions (PDFs) <cit.> as well as diffractive PDFs include a wide range of HERA measurements of electron-proton process. Recent developments in the field of PDFs of the proton have led to a renewed interest in diffractive DIS process to extract the non-perturbative diffractive PDFs from a global analysis. Significant progresses in understanding diffraction processes have been made at the HERA collider, where typically a 27.5 GeV electron (or positron) collides with a 820 or 920 GeV proton. Recent analyses of diffractive DIS data collected by H1 and ZEUS collaborations at HERA have confirmed substantial contributions of perturbative QCD-based effects in diffractive DIS cross sections. In general, in ep collisions at HERA, the hard diffraction contributes a fraction of order 8-10% to the total DIS cross sections.

Future high precision and high energy DIS experiments are expected to reach a wider kinematic range of momentum fraction x and photon virtuality Q^2 which could not be explored previously by ep HERA. Among them are the Large Hadron-electron Collider (LHeC) <cit.> and Future Circular Collider in electron-hadron mode (FCC-eh) <cit.>. LHeC would utilize the 7 TeV proton beam from the LHC and collides it with a 60 GeV electron beam, and would extend the available kinematic range in x and photon virtuality Q^2 by a factor of order 20 and 100, respectively. Beyond the LHeC, the next generation ep collider FCC-eh, utilizing the 50 TeV proton beam from the FCC which would probe DIS processes at center-of-mass energy of 3.5 TeV, much higher than the HERA collider, leading to a better understanding of the proton structure with extremely high precision. Recently, an investigated of the potential of theses high energy and high luminosity machines for the measurement of diffractive DIS cross sections and to constrain the diffractive PDFs has been performed, see the analysis of Ref. <cit.> for details.

On the phenomenological side, a similar strategy to determine the PDFs can also be adapted to the case of diffractive DIS, considering the collinear factorization <cit.> and the validity of proton vertex factorization <cit.>. Diffractive PDFs can be extracted from QCD analyses of diffractive DIS data sets. Similarly to the PDFs, the diffractive PDFs are expected to obey the by the (Dokshitzer?Gribov?Lipatov?Altarelli?Parisi) DGLAP equations. In this picture of proton vertex factorization, the diffractive DIS processes are described by the exchange of colourless object such a pomeron. Recent progress in the determination of diffractive PDFs widely used the diffractive DIS data sets from H1 and ZEUS experiments at HERA (see, for example, Refs. <cit.> for recent reviews). In the last few years, at least three groups have reported sets of diffractive PDFs with uncertainties using the mentioned data sets: H1-2006-DPDF <cit.>, ZEUS-2010-DPDF <cit.>, and the most recent analysis by GKG18-DPDF <cit.>. All these of diffractive PDFs determinations were performed at next-to-leading order (NLO)
accuracy in perturbative QCD. The primary focus of these QCD analyses were put on quantifying the effects of the inclusion of new measurements of diffractive DIS at HERA as well as the diffractive dijet production. The analysis by GKG18-DPDF introduced some improvements over previous determinations. Specifically, in order to achieve a more reliable estimate of the diffractive PDFs uncertainties, the xFitter package have been employed. In addition, GKG18-DPDF used for the first time the most recent H1/ZEUS combined diffractive DIS cross section measurements. Up to now, predictions for diffractive DIS, and in particular for diffractive dijet production, were performed only at next-to-leading order (NLO) accuracy. Recently, predictions for the diffractive dijet production, also is provided at next-to-next-to-leading order (NNLO) in Ref. <cit.>. Although much theoretical work remains to be done for this process, it is already clear that the accumulation of precise data for diffractive DIS processes will greatly deepen our understanding of perturbative QCD (pQCD).

In the present study we construct for the first time a set of NLO and NNLO diffractive PDFs using all available and up-to-date diffractive DIS cross section measurements <cit.> including the most recent H1 and ZEUS combined measurement for the inclusive diffractive DIS cross sections <cit.>. We do so by using a methodology which has been suggested in Refs. <cit.>, and recently used to study leading neutron production <cit.> as well as inclusive diffractive DIS <cit.>: the so called “fracture function approach”.

Similar to the case of ordinary structure functions in totally inclusive DIS, QCD does not predict the shape of the fracture functions unless it is known at
a certain initial scale, Q_0^2. This universal and non-perturbative information has to be determined from a QCD analysis of experimental data sets, and hence, can be parametrized finding inspiration in non-perturbative models. More recently, the fracture functions approach has been successfully applied to describe leading neutron and leading proton productions at H1 and ZEUS experiment, considering a model as non-perturbative input <cit.>.

The main aim of this analysis is to provide a conceptual theoretical framework based on fracture function approach to describe the diffractive DIS processes. 
As a result of the following analysis, we present a parametrization for the fracture function that characterizes the underlying diffractive DIS process at
an initial scale, as extracted from H1 and ZEUS data sets. The obtained parametrization is used to compute other observables measured by H1, not included in our QCD fit, finding also an outstanding agreement with the data. In this analysis we emphasis that an approach in the framework of fracture functions phenomenologically allows a accurate and reliable description of diffractive DIS cross sections. Our results also verify that the scale dependence of the diffractive DIS data agrees well with the one predicted by the use of fracture function formalism. Therefore, this paper provides an important opportunity to advance the understanding of diffractive DIS events measured by H1 and ZEUS collaboration at HERA.

This paper has been organized in the following way: Sec. <ref> begins by laying out the theoretical framework and assumptions of this research. In
particular, we discuss the methodology used for this study including the diffractive structure functions, diffractive PDFs, hard-scattering factorization and heavy quark contributions. Then in Sec. <ref> we present our input for the diffractive PDFs at a given initial scale Q^2 = Q_0^2. In Sec. <ref> we concentrate on the description of the H1 data sets which have been added in our QCD analysis, together with a discussion on the inclusion of H1 and ZEUS combined data sets. 
Sec. <ref> is concerned with the methodology used in this study to determine the diffractive PDFs uncertainty. The results of the global analysis can be found in Sec. <ref>. This section starts with a presentation of our NLO and NNLO diffractive PDFs and their uncertainties, together with the values of the input parameters. Detailed discussions of the main results and comparisons with the analyzed data sets are also presented as well. Finally, the conclusion in Sec. <ref> gives a brief summary and critique of the findings.









§ THEORETICAL FRAMEWORK AND ASSUMPTIONS






Predictions for the diffractive processes in DIS can be obtained in the framework of perturbative QCD (pQCD). As we already mentioned, hard processes in diffractive DIS can be described by a factorization into parton level subprocesses and diffractive PDFs. In this framework, cross sections for diffractive deep inelastic lepton proton (ℓ p) scattering can be computed at NLO and NNLO accuracy in pQCD.












 §.§  Diffractive structure functions 





The kinematical variables to describe the diffractive DIS events can be inferred from the momenta of the incoming lepton (ℓ ), proton (p) and the outgoing lepton. The leading order Feynman diagram for diffractive DIS at HERA is displayed in Figs. <ref> and  <ref> in the picture of proton vertex and collinear and QCD hard scattering
collinear factorizations, respectively. In the diffractive DIS in which belongs to the semi-inclusive lepton proton DIS (SIDIS), 











    ℓ (k)  +  p(P) →ℓ (k^')
    +  p(P^')
    +  X(p_X)   ,



in which in addition to the outgoing leptons, an additional proton (p) can be detected in the final state of the diffractive DIS processes.
In above equation, X indicates to the unobserved part of the hadronic final state. The kinematics of such events is specified by the following kinematical variables: 



    Q^2 = -q^2,
        
    x_B = Q^2/2P.q,
        
     y = P.q/P.k .



The momentum transfered to the target proton  (p)  is given by the momentum q^2 = (k - k^')^2 of the virtual photon γ^*.
Hence, the Q^2 in Eq. (<ref>) is the well-known photon virtuality. x_B is the usual x-Bjorken variable in DIS process which idicate the longitudinal momentum fraction of the proton carried by the struck quarks, and y is refereed to the inelasticity of the DIS scattering.

In the lepton-proton (ℓ p) center-of-mass (CM) system, diffractive DIS processes are interpreted by an outgoing proton (p) with a large momentum fractions of the incident protons as well as quite small values of the  invariant transverse momentum t measured with respect to the collision axis in the “target fragmentation” regions of the incident protons.

In order to the describe such diffractive events, one need to introduce additional kinematics invariants. As illustrated in Fig. <ref>, the longitudinal
momentum fractions, x_ of the colourless exchange with respect to the incoming proton (p) momentum, and the invariant momentum transfer t at the proton vertex t, and β of the struck quarks with respect to the colourless exchange, are also needed to describe the diffractive DIS processes. These two new variables are given by,



    x_ = q.(P - P^')/P.q,     β = Q^2/Q^2 + M_X^2 .



M_X is the invariant mass of the hadronic final state X. The scaled fractional momentum variable β which is also defined by β = x_B/x_ interpreted as the fractional momentum of interacting partons in the protons with respect to the pomeron  fractional momentum x_. 

Now we turn to inclusive diffractive DIS cross sections. The neutral current (NC) measurements in diffractive DIS are often presented in terms of the reduced lepton-proton (ℓ p) cross section, “neutral current diffractive reduced cross section” σ_r^D(4). The t-unintegrated reduced cross section depends on the diffractive transverse and longitudinal structure functions F_2^D(4) and F_L^D(4), respectively. In the one-photon exchange approximation, it is given by:



    σ_r^D(4) (β, Q^2; x_, t) =        F_2^D(4)  (β, Q^2; x_, t)  - 
          y^2/1+(1-y)^2 F_L^D(4)  (β, Q^2; x_, t)  .



In the measurements of inclusive diffractive DIS at HERA, the data were presented in terms of σ_r^D(4) (β, Q^2; x_).
In this analysis, we analyze all available and up-to-date diffractive DIS data sets  <cit.> including the most recent data for combined H1 and ZEUS diffractive DIS cross sections measurements <cit.>.










 §.§  Diffractive parton distributions and hard-scattering factorization 





According to the factorization theorem for diffractive DIS in pQCD <cit.>, if the process is sufficiently hard, the calculation
can be subdivided into two components: the hard partonic cross sections which are calculable within the pQCD in powers of strong coupling constant α_s, which need to be convoluted with soft diffractive PDFs, F_i/p^D  (β, μ^2_F; x_, t) that specify the contributing partons inside the incoming hadrons. Likewise the PDFs, diffractive PDFs are also universal for all diffractive DIS events  with the hardness of the process being ensured by the virtuality Q^2 of the exchanged photons. The structure functions (SF) appearing in Eq. (<ref>) are given by



    F_k^D(4) (β, Q^2; x_, t)
    =       ∑_i∫_β^1d ζ/ζ F_i/p^D
    (β, μ^2_F; x_, t)  
           C_ki (β/ζ,
    Q^2/μ^2_F, α_s(μ^2_R) )
    +  O (1/Q^2)  .



The index i in Eq. (<ref>) runs on the flavors of the interacting partons in the proton <cit.>. The quark and gluon hard scattering coefficient functions, C_q and C_g, which appeared in above equation are perturbatively
calculable as a power expansion in the strong coupling and are the same as in fully inclusive DIS process <cit.>.  As one can find in literature, for example <cit.> as well as our discussion in previous section, in order to describe the diffractive DIS events, two more variable in addition to the usual DIS kinemtical variables are also needed which are the longitudinal momentum fraction x_ and the invariant momentum transfer t.
In the general factorization theorem <cit.> for the diffractive DIS events in the form of above equation holds at fixed values of x_ and invariant momentum transfer t <cit.>.
Hence, the parton content of the colour singlet exchange described by diffractive PDFs, F_i/p^D, is uniquely controlled by the kinematics of the outgoing protons in which carry the fractional momentum of 1 - x_.

The diffractive PDFs, F_i/p^D(β, μ^2_F; x_, t), appearing in Eq. (<ref>) are known as “proton-to-proton fracture functions” <cit.> in the very forward kinematical region.
It can be interpreted as the number density of interacting partons at a given scale of μ^2 with fractional momentum of β conditional to the detection of a final state proton (p) with fractional momentum 1-x_ and invariant momentum transfer t <cit.>. 
As in the inclusive case, the Q^2 evolution of the diffractive DIS events can be predicted in pQCD. Since the scale dependence of the cross section in leading particle production in DIS processes can be calculated within pQCD <cit.>, therefore the diffractive PDFs also obey the standard DGLAP evolution equations <cit.>. The t-unintegrated diffractive PDFs presented in Eq. (<ref>) also obey the standard DGLAP evolution equations <cit.>. 

By integrating Eq. (<ref>) with respect to the momentum transfer t over up to values of order Q^2, the diffractive PDFs obey an inhomogeneous DGLAP-type evolution equations <cit.>.
We should note here that most of diffractive DIS events occur for the small values of t. In the case for the H1/ZEUS combined data, the squared four-momentum transfer at the proton vertex, t, is integrated in the restricted range of 0.09 GeV^2 < |t| < 0.55   GeV^2 . In the case where the diffractive PDFs presented in Eq. (<ref>) are integrated over t in a limited range, they can obey the standard DGLAP evolution equation <cit.>.  Hence, we considered such approximation in our analysis and therefore one can use the standard DGLAP evolution when F_i/p^D  (β, μ^2_F; x_, t) are integrated over t in a limited range <cit.>



    F_i/p^D  (β, μ^2_F; x_) =  ∫_t_ min^t_ max dt   F_i/p^D  (β, μ^2_F; x_, t) ,      t_ max≪ Q^2 .



Hence, like for the case of PDFs, the evolution equations of diffractive PDFs in diffractive DIS are easily obtained by the DGLAP evolution equations <cit.> as



    Q^2 ∂ F^D_Σ/p (β, Q^2; x_)/∂ Q^2        =  α_s(Q^2)/2 π
          ∫_β^1du/u P_Σ^j(u)   F^D_Σ/p (β/u, Q^2; x_) , 
    
    Q^2 ∂ F^B_g/p (β, Q^2; x_)/∂ Q^2        =  α_s(Q^2)/2 π
          ∫_β^1du/u P_g^j(u)   F^D_g/p (β/u, Q^2; x_) ,



where F^D_Σ/p (β, Q^2; x_) and F^B_g/p (β, Q^2; x_) correspond to the singlet and gluon diffractive distributions, respectively.
These non-perturbative distributions need to be parametrize at an input scale Q_0^2 and their evolution to higher scale, Q^2>Q_0^2, can be described by using the evolution equation given above.
P_Σ and P_g in Eq. (<ref>) are the common NLO and NNLO contributions to the splitting functions governing the evolution of unpolarized singlet (S) and non-singlet (NS) combinations of quark densities in perturbative QCD.
Splitting functions are perturbatively calculable as a power expansion in the strong coupling constant α_s. The splitting functions P_Σ and P_g in Eq. (<ref>) are the same as in fully inclusive DIS <cit.>.

We have now, in principle, all the essential ingredients to write down the diffractive DIS cross sections of Eq. (<ref>) in terms of diffractive structure functions.
In the next sections, we give a detailed account of the global analysis of diffractive PDFs performed in this study. We first discuss the heavy flavour contributions, the details parameterization of diffractive PDFs 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.












 §.§  Heavy flavour contributions to the diffractive DIS structure function 




Is is shown in literature that the correct treatments of heavy flavors have an important impact on the PDFs extracted from the global analysis due to the data available for
heavy flavor structure functions, F^h_L (x, Q^2, m_h^2) with h = c, b and k= 2, L, as well as the heavy flavor contributions to the total inclusive structure functions at small values of momentum fraction x <cit.>.
There are several scheme available in literature where heavy quark production can be readily described. Among them are zero-mass variable flavor scheme (ZM-VFNS), fixed-flavor number scheme (FFNS) and general-mass variable flavour number scheme (GM-VFNS). 

The resulting calculations for the heavy quark coefficient functions are more accurate in two different and complementary regimes. The FFNS is more accurate for values of Q^2 comparable to the mass of the heavy quark involved in the calculation, while the ZM-VFNS is instead more accurate for values of Q^2 much larger than the heavy-quark mass, Q^2 ≫ m_h^2. It has been shown in literature that the GM-VFNS is the most appropriate scheme to extract PDFs from a global QCD analysis <cit.>. 
This scheme could extrapolates smoothly from the FFNS at low value of Q^2 to the ZM-VFNS at high Q^2, and therefore, produces a good description of the effect of heavy quark contributions to the structure functions as well as inclusive cross section over the whole range of Q^2. Well known examples of GM-VFNS are the original Aivazis-Collins-Olness-Tung (ACOT GM-VFNS) scheme <cit.>, the Thorne-Roberts (TR GM-VFNS) scheme <cit.> and the so-called FONLL GM-VNFS scheme <cit.>. In our previous diffractive PDFs analysis, GKG18-DPDF <cit.> we adopted the more recent “TR” prescription <cit.>.

In the present analysis, we prefer to use the FONLL GM-VFNS <cit.> which provides the charm structure functions by using the exact x-space O (α_s^2) heavy quark coefficient functions as well as computes neutral- and charged-current DIS observables up to this order. The evolution of diffractive PDFs F^D (β, Q^2; x_) is performed within the FONLL GM-VFNS by using the publicly available APFEL program <cit.>. The QCD parameters are the ones quoted in the analysis by NNPDF Collaboration <cit.>. In particular, we use the heavy flavor masses for charm and bottom as m_c = 1.51 GeV and m_b = 4.92 GeV, respectively. Also the Z-boson mass is chosen to be M_Z = 91.1876 GeV and the the strong coupling is evaluated at two loop setting α_s,  NLO^n_F=5 (M_Z) =  0.1185 <cit.>. This selection of α_s (M_Z) is consistent with the very recent determination of the strong coupling constant reported by the NNPDF3.1 <cit.>, in which for the first time uses the jet production and t t̅ differential distributions at NNLO accuracy.











§  INPUT DISTRIBUTIONS 





In this section, we present the input distributions for the diffractive PDFs at the reference scale of Q^2 = Q_0^2. We also glance ahead to mention some of the main feature of this selection.
As is clear from the discussion in Sec. <ref>, one improvement in the following analysis is to use parameterizations for the input diffractive PDFs based on the fracture function approach. Following the detailed studies in Refs. <cit.>, we take for most PDFs a parameterization of the form



    β  F_q (β, Q_0^2;x_)  =   W (x_)  β  f_q (β, Q_0^2)  ,  
    β  F_g (β, Q_0^2;x_)  =   W (x_)  β  f_g (β, Q_0^2)  ,



where Q_0^2 = 2 GeV^2 is the the input scale, and  W (x_) is the flux factor which we assume that it depends only to the x_ variable. For the diffractive flux factor, we use the standard functional form which can be written as,



    W (x_) = x_^w_1 (1 - x_)^w_2 (1 + w_3 x_^w_4)  .



The global fit determines the values of the set of parameters w_1, w_2, w_3 and w_4 for diffractive quark and gluon PDFs. 
We consider that this flux factor to be the same for the light sea and gluon, as there is not enough data which can constrain these distributions, while leaving all four parameters w_i in the polynomial free leads to instabilities in the fit. The dependence of the diffractive PDFs on the x_ depends very much on the parameterization that one chooses for the flux factor. In Eq. (<ref>), β  f_q (β, Q_0^2) and β  f_g (β, Q_0^2) are the quark and gluon densities at the reference scale Q_0^2. The quark and gluon densities depend on the scaled fractional momentum variable β and photon virtuality Q^2. 



    β  f_q (β, Q_0^2) =  N_q  β^α_q (1 - β)^β_q ( 1 + γ_q √(β) + η_q β^2 ) , 
    β  f_g (β, Q_0^2) =  N_g  β^α_g (1 - β)^β_g ( 1 + γ_g √(β) + η_g β^2 ) .



For the quark density β  f_q (β, Q_0^2), we consider a standard functional parameterization. Not all the parameters in our inputs for the quark density are free. We will return to this issue that due to the lack of enough diffractive DIS data sets, one needs to fixed some of these variables at their best fit values, and therefore, there are potentially less free parameters in the diffractive PDFs fit rather than what we presented in Eq. (<ref>).
Like to the case of quark density, the poorly determined gluon density β  f_g (β, Q_0^2) is also taken to has a simpler input functional form as presented in Eq. (<ref>).
It is worth to mention here that our parametrizations for the input distributions and all free parameters listed there allow the QCD fits a large degree of flexibility. 

Considering the flux factor introduced in Eq. (<ref>) as well as the quark and gluon densities presented in Eq. (<ref>), one can express the diffractive PDFs in terms of the following set of basis functions for the quark and gluon diffractive PDFs,



    β  F_q (β, Q_0^2; x_)        =   W (x_) 
           N_q  β^α_q (1 - β)^β_q ( 1 + γ_q √(β) + η_q β^2 ) ,  
    β  F_g (β, Q_0^2; x_)        =  W (x_) 
           N_g  β^α_g (1 - β)^β_g ( 1 + γ_g √(β) + η_g β^2 )  .



As already emphasized, the introduction of such flexible parameterization for the diffractive PDFs at the reference scale gives a much better description of the diffractive DIS data. 
Eq. (<ref>) states that the variables x_, which related to the loosely scattered proton in diffractive DIS, are factorized from the variables characterizing the diffractive system (β, Q^2). As we mentioned in Sec. <ref>, for a fixed value of x_, the Q^2 evolution of diffractive PDFs β  F (β, Q^2; x_) is given by the DGLAP equations. 












§  DIFFRACTIVE DIS DATA SETS 





We strive to include as much of the available diffractive DIS data as possible in our data sets.
However, one need to apply some certain kinematical cuts in order to ensure that only the data sets for which the available
perturbative QCD treatment is adequate are included in the QCD fit. All the data sets used in the GKG18 diffractive PDFs analysis are also included in this analysis. 
Statistically, most significant data set that we use in our QCD analysis are the HERA measurements of the diffractive DIS “reduced” cross sections. 
An overview of all available and most up-to-date diffractive DIS data set is presented in Table. <ref>, where we indicate, for each data set: the name of experiments and the corresponding references, observables, the kinematic range of β, x_ and Q^2, and finally the number of data points before kinematic cuts.








We now discuss the inclusion of the HERA diffractive DIS data sets into our diffractive PDFs fit. The data sets that are included in our analysis are as follow:
The H1 and ZEUS combined diffractive DIS cross section measurement <cit.>, the H1-LRG-12 <cit.> data and finally the H1-LRG-11 <cit.> data for three different center-of-mass energies of √(s) = 225, 252 and 319 GeV.

Now we are in a position to present the details of these data sets in turn. As we mentioned, in our QCD fit, we use the most recent H1 and ZEUS
combined measurement for inclusive diffractive scattering cross sections <cit.>. Until recently, most of the other groups in literature that have performed global diffractive PDFs analyses do not include this combined data sets. An exception is the analysis of Ref. <cit.> by GKG18-DPDF. In the present work,
and in GKG18-DPDF, we have added the H1/ZEUS combined data in order to present a well-determined diffractive PDFs with a reliable uncertainty. This combined data set contains complete information on diffractive DIS cross sections measurement published by the H1 and ZEUS Collaborations. The kinematic range of these combined data sets is 2.5   GeV^2 < Q^2 < 200   GeV^2 in photon virtuality, 3.5 × 10^-4 < x_ < 9.0 × 10^-2 in proton fractional momentum loss, 1.8 × 10^-3 < β < 0.816 in scaled fractional momentum variable and finally 0.09 < |t| < 0.55 GeV^2 in squared four momentum transfer at the proton vertex. This measurement used samples
of diffractive DIS data at a center-of-mass energy of √(s) = 318  GeV.  The combination by H1 and ZEUS Collaborations is based on the cross sections measured
with the ZEUS LPS 1 <cit.>, ZEUS LPS 2 <cit.>, H1 FPS HERA I <cit.> and finally the H1 FPS HERA II <cit.> data sets. 

Another data sets we have used in our analysis is the recent inclusive measurement of diffractive DIS at HERA by H1 Collaboration, entitled as H1-LRG-12 <cit.>. The measurement is restricted to the phase space region 3 < Q2 < 1600 GeV^2 of the photon virtuality , the square of the four momentum transfer at the proton vertex with |t| < 1.0 GeV^2 and finally the longitudinal momentum fraction of the incident proton x_ < 0.05 carried by the colourless exchange.
These high statistics measurements at HERA cover the data taking periods of 1999-2000 and 2004-2007. This measurement is combined with previously published results by H1 Collaboration <cit.> in order to provide a single data set of diffractive DIS cross sections using the large rapidity gap (LRG) selection method.  Like for the case of H1/ZEUS combined data sets, this measurement is also presented as the reduced diffractive cross section, σ_r^D(4) (β, Q^2; x_, t).

Finally, we have used H1 measurements of the diffractive DIS reduced cross section <cit.> at center-of-mass energies of √(s) = 225, 252 GeV as well as the precise measurement at √(s) = 319 GeV. This reduced cross section measurements in this experiment, entitled H1-LRG-11, is measured in the photon virtuality range of 4.0 < Q^2 < 44.0 GeV^2 and longitudinal momentum fraction of the diffractive exchange x_ of   5 × 10^-4 <  x_ < 3 × 10^-3.

To ensure the validity of the well-known DGLAP evolution equations one needs to impose certain cuts on the data sets we discussed in this section.
We follow the analysis by H1-2006 DPFs <cit.> and continue to use the same cuts on the diffractive DIS data, i.e. Q^2> Q^2_min=8.5 GeV^2 in order to make sure that higher-twist corrections (HTs) are not needed in the analysis. As an aside, we should comment on the very large x_ domain. One can not impose any cut at large x_, although, at present, there are essentially no diffractive DIS data available probing the x_ > 10^-1 domain. In addition to the cut on Q^2, we include the data with M_X > 2 GeV in the fit.

Having discussed the details of kinematic cuts applied on the data sets in this analysis, we are now ready to discuss the corrections need to be taken into account for the H1 and ZEUS combined data. The first correction that one needs to apply on the data sets is due to the proton dissociation background. 
We should notice here that distinct methods have been used by the H1 and ZEUS experiments at HERA, and hence, the measured cross sections are not always presented with
the corrections for proton dissociation backgrounds. In this analysis, the combined H1/ZEUS diffractive DIS cross sections are corrected by a global factor of
1.21 to account for the such contributions.
Another correction one needs to consider comes from the fact that the all H1-LRG data sets used in our analysis are given for the range |t| < 1 GeV^2 while the combined H1/ZEUS diffractive DIS cross sections, which is based upon proton-tagged samples, are restricted to the t range of 0.09 < |t| < 0.55 GeV^2. Hence, one need to consider a global normalization factor to account this correction. The extrapolations for these two distinct range of t have been performed by considering an exponential t dependence of the inclusive
diffractive DIS cross sections, using the H1 default value of exponential slope parameter b = 6 1/GeV^2 <cit.>. Therefore, our analysis has been carried out by considering the above corrections as well as the pre-selection kinematical cuts applied on the analyzed data sets.

It is important to emphasize here that the list of diffractive DIS data sets considered in this work contains enough information to extract diffractive PDFs from a global QCD analysis. Clearly, there are other source of important processes that will provide additional information on the diffractive PDFs in the LHC era. Among these, one could consider diffractive dijet productions <cit.>, providing useful information on the diffractive gluon density and, possibly, on the quark flavor separation respectively. As we discussed in Introduction, 
the future high energy and high luminosity machines such as LHeC and FFC-eh will measure the diffractive DIS cross sections with the highest possible precision leading to the much better constrains for the diffractive PDFs <cit.>. Having discussed the diffractive DIS data set as well as the kinematic cuts that we apply, we are now turn to discuss the method of χ^2 minimizations and our approach to determine the diffractive PDFs uncertainties.

	






§  THE METHOD OF Χ^2 MINIMIZATIONS AND DIFFRACTIVE PDFS UNCERTAINTIES 




In the present analysis, a χ^2 (p) minimizations method to extract the fit parameters and an error calculation based on the “Hessian approach” to determine the uncertainties of diffractive PDFs are applied. As discussed in literature review (e.g. see <cit.>), a precise understanding of uncertainties due to PDFs in a QCD analysis is crucial to precision studies of the standard model (SM) of particle physics <cit.>, as well as to searches for new physics (NP) beyond the standard model (BSM) <cit.>, especially for the future high-luminosity LHC and high-energy LHC.. In turn, new measurements of SM processes from present or future high energy collider experiments such as HERA, LHC, Tevatron, RHIC and LHeC can be used to constrain the uncertainties on the PDFs in a global QCD analysis. The most complete method for obtaining well constraints from the new data sets on the PDFs would be to add the up-to-date data sets into the global QCD analysis.
Currently, the two most commonly used methods to estimate the PDFs uncertainties from a QCD analysis are the Monte Carlo approach <cit.> and the standard Hessian method <cit.>. The details of Monte Carlo sampling techniques and use of Neural Networks (NN) as unbiased parametrization can be found in the analysis by NNPDF Collaboration <cit.>.

The details of “Hessian approach” is fully addressed in Refs. <cit.>.
In the standard “Hessian approach??, smaller number of` “error sets?? are considered to obtain an estimate of the PDFs errors <cit.> in QCD fits. These error sets are correspond to the plus and minus eigenvectors directions in the space of PDFs fit parameters, which finally can be used to estimate the χ^2 function near their global minimum. The Hessian approach relies on a quadratic approximation for the fit parameters dependence of the χ^2({η_i}) minimization function and a linear approximation  for the parameter dependence of the observables in question.  In practice, the Hessian method works quite well for the most high energy observables in which used to determine the PDFs <cit.>. A version of this method, “Hessian profiling”, has been included in the xFitter package <cit.>. Our previous analysis on diffractive PDFs, GKG18-DPDF <cit.>, was based on this package.

Let's discuss the χ^2  ({η_i}) minimization procedure first. The goodness of fit is traditionally determined by the effective global χ_ global^2({η_i}) minimization algorithm that measures the quality of QCD fit between theory predictions and diffractive DIS experiments. This minimization strategy allow an extraction of independent parameters {η_i} which specify the diffractive PDFs at the input scale Q_0^2 =  2 GeV^2. This function is given by,



    χ_ global^2({η_i}) =
    ∑_n=1 ^ N^ exp w_n   χ_n^2({η_i})  ,



where n labels a particular data set and w_n is a weight factor for the n^ th experiment with default value of 1 <cit.>. χ_n^2({η_i}) can be written as



    χ_n^2({η_i}) =       ( 1 - N_n /Δ N_n)^2 + 
          ∑_j=1^N_n^data( N_n    O_j^ exp -  T_j^ theory ({η_i}) / N_n  δ O_j^ exp)^2 .



The minimization of χ_ global^2({η_i}) function presented above is done applying the CERN program library MINUIT <cit.>. 
In Eq. (<ref>), O represents the experimental measurement, and δ O denotes the experimental uncertainty
(statistical and systematic combined in quadrature), and T({η_i}) is the theoretical value for the i^ th experimental data point which depends on the input diffractive PDFs parameters {η_i}. In this equation, N_n is the overall normalization factors for a given data point of experiment n and the Δ N_n is the experimental normalization uncertainty. The (fitted) normalization factors N_n need to be extracted along with the diffractive PDFs parameters. We minimize the χ_ global^2 ({η_i}) value with the final 7 fit parameters. The obtained normalization factors N are presented in Table. <ref> for each experiment.  














		






			






			







In order to show the effects arising from the use of the different diffractive DIS data sets, in Tables. <ref>, <ref>, <ref>, <ref> and <ref>, we present the χ^2 for each bin of x_ for our NLO QCD analysis. These tables illustrate the quality of our NLO and NNLO QCD fits to diffractive DIS cross sections data sets in terms of the individual χ^2 values obtained for each experiment at a certain value of x_. As shown in these tables all the data sets are well-fitted.  We find χ^2/ d.o.f= 312.07/346 = 0.91  ( NNLO) and  χ^2/ d.o.f = 310.83/346 = 0.89  ( NLO) which yield an acceptable fit to the experimental DIS data. The obtained χ^2/ d.o.f shows a slightly improvement in the quality of the fit at NNLO accuracy.

In the following, we briefly discuss how to apply the “Hessian method” to determine the uncertainties of extracted diffractive PDFs. The basic procedure of this method is provided in Refs. <cit.>. The diffractive PDFs, β F (β, Q^2; x_), defined at the initial scale Q_0^2, are parametrized by η parameters. As we mentioned, the determination of the diffractive PDFs is obtained using a χ_n^2 ({η_i}) function given in Eq. (<ref>), which quantifies the discrepancy between the QCD theory predictions and the experimental observables of a global set of experiments, including the experimental errors.

In the well-known “Hessian” method one can diagonalize the covariance matrix and work in terms of the eigenvectors and eigenvalues.
As we have mentioned earlier, the appropriate parameters set can be obtained by minimizing the χ_n^2 ({η_i}) function. We entitled this as diffractive PDFs set
s_0. The parameters values of s_0, i.e. {η_1^0 …η_n^0}, in which extracted from QCD fit to diffractive DIS data, will be presented and discussed in details in Sec. <ref>.

By moving away the parameters from their best fitted values, χ^2 increases by the amount of Δχ^2


    Δχ^2_ global   =   χ^2_ global({η}) -
    χ^2_0({η^0}) 
           = ∑_i, j=1^n (η_i -
    η_i^0)  H_ij  (η_j - η_j^0) ,


where H_ij is the Hessian or error matrix which is given by


    H_ij = 1/2 ∂^2 χ^2_ global/∂η_i  ∂η_j|_ min .


These covariance or Hessian matrix can be obtained by running the CERN program library MINUIT <cit.>.
Having ay hand the derivative of the observable O with respect to each parameter {η}, one can use the following equation


    Δ O^2 = 
    Δχ^2_ global ∑_i, j = 1^k∂ O/∂η_i  C_i j ∂ O/∂η_j ,


and can calculate the diffractive PDFs symmetric error bands as well as the corresponding observables such as the diffractive DIS cross sections for a desired values of confidence interval, T = Δχ^2_ global. We have to mentioned here that, in Eq. (<ref>), C_ij≡ H_ij^-1  is the elements of covariance or error matrix.

We should stressed that it is convenient to work in term of the eigenvalues and orthogonal eigenvectors of covariance matrix which is given by


    ∑_j = 1^n C_ijυ_j k =
    λ_kυ_i k .


The displacement of the parameter {η_i} from its obtained minimum values η_i^0 can be expressed in terms of the rescaled eigenvectors e_i  k = √(λ_k)   v_i  k. It reads


    η_i - η_i^0 =
    ∑_k = 1^n
      e_i k  z_k .


Considering the orthogonality of eigenvectors υ_i k and putting Eq. (<ref>) in (<ref>), one can obtain


    Δχ^2_ global =
    χ^2_ global({η}) - 
    χ^2_0({η^0})
    = ∑_k=1^n  z_k^2 .


The relevant neighborhood of χ ^ 2 is the interior of hypersphere with radius T.
The uncertainty in the diffractive PDFs is then set by the requirement Δχ^2_ global  <  T^2 at some prescribed confidence level (CL). This means that


    ∑_k = 1^n
    z_k^2 ≤ T^2 .


Finally the neighborhood parameters
can be written as


    η_i(s_k^±) =
    η_i^0± t  √(λ_k) 
    v_i k ,


with s_k is the k^ th set of diffractive PDFs, t adapted to make the desired T^2 = Δχ^2_ global which is the tolerance for the
required confidence interval (CL) and t=T in the quadratic approximation.

Using the method discussed above, we accompany the construction of our QCD fit by reliable estimation of diffractive PDFs uncertainty.
Finally uncertainty of a given observables O in the Hessian method, which can be the diffractive PDFs or or related sort of diffractive observables,  can calculate as <cit.>


    Δ O =
    1/2[∑_k = 1
    ^n( O(s_k^+) 
    -  O(s_k^-))^2]
    ^1/2 .


In above equation, O(s_k^+) and O(s_k^-) are the value of observables O extracted from the input set of parameters η_i(s_k^±) obtained from Eq. (<ref>).
The Δχ^2_ global values determine the confidence region, and it is calculated so that the confidence level (CL) becomes the one-σ-error range for a given number of parameters. In this analysis, we calculate the diffractive PDFs uncertainty with Δχ^2_ global = 1 and present their symmetric errors. We should mentioned here that a one-unit tolerance of T=1 would be the most appropriate choice if all the uncertainties
in the data points as well as the theory would be gaussian and perfectly accounted for. Hence, in order to obtain a more conservative error estimate one need to choose a reliable choice of tolerance. 

For the case of η degrees of freedom, the Δχ^2_ global value needs to be calculated to determine the appropriate size of diffractive PDFs uncertainty. Assuming that the Δχ^2_ global follows the χ^2 distribution with η degrees of freedom, we have the confidence level P as <cit.>




    P =
    ∫_0^Δχ^2_ global1/2 Γ(η/2)
    (ξ/2)^(η/2)-1
    e^-ξ/2  dξ ,



where Γ is the Gamma function. For the case of the one-free-parameter fit, we obviously have Δχ^2_ global=1. 
Sine the diffractive PDFs in common QCD fits are considered with several free fit parameters, η>1, the value of Δχ^2_ global should be calculated form Eq. (<ref>). Considering the parameter number in our diffractive PDFs analysis, Eq. (<ref>) leads to the tolerance criterion of Δχ^2_ global = 5.0
































In the next section, we will present full details of our fit quality and the extracted diffractive PDFs.






§  RESULTS AND DISCUSSIONS 




In this section, we discuss the main results of this work, namely the NLO and NNLO diffractive PDFs sets extracted form QCD analyses of diffractive DIS data sets. 
We first discuss the numerical results of present work including the values of fitted parameters. Secondly, we discuss the quality of the QCD fits and compare  our NLO and NNLO predictions to the fitted data sets. Then we show the resulting diffractive PDFs and their uncertainties. We mainly focus on the perturbative convergence upon applying the fracture function framework as well as inclusion of higher-order QCD corrections. Then we compare the extracted diffractive PDFs with the most recent results in literature, especially the results from GKG18-DPDF <cit.> analysis. 

In Table. <ref>, we present the fit parameters obtained in this work for our NLO and NNLO QCD analyses at the initial scale of Q_0^2 = 2  GeV^2. Values marked with (*) are fixed in the fit since the analyzed data sets do not constrain these parameters well enough. These fitted parameters include {N, α, β, γ, η} for the diffractive quark and gluon densities. The parameters for the weight factor W (x_) also presented as well.
As one can see from this table, the values of the physical parameters used in the computation of diffractive DIS cross section and in the evolution of diffractive PDFs are the same as those used in the NNPDF global analysis of unpolarized PDFs and FFs <cit.>. Specifically, we use α_s(M_Z^2) = 0.1185 as reference value for the QCD couplings, and m_c = 1.51 and m_b = 4.92 GeV for the charm- and bottom-quark masses.









In the present fit of diffractive PDFs, the analyzed data sets is much more limited than in a typical global proton PDFs fit. As we discussed earlier, diffractive DIS data allow for
the determination of only two independent combinations of densities, namely the singlet and the gluon. In addition, these distributions at some region of β sill remain unconstrained. This effect is more enhanced for the case of diffractive gluon density, due to the reduced sensitivity of the diffractive DIS data included in our QCD fit to the gluon distribution. Therefore one would expect that the diffractive gluon PDFs is determined with larger uncertainties than the diffractive quark PDFs. This issue deserves a separate comments which we will discuss later. In addition, diffractive DIS is blind to the separation between quark and antiquark distributions. These points confirm that the input functional form and it's sensitivity to the parameter space of {η_i} and {w_i} need to be clearly investigated. As we mentioned, the currently
available diffractive DIS data sets do not fully constrain the entire β and x_ dependence of quark and gluon diffractive PDFs presented in Eq. (<ref>). Consequently, we are forced to make some restrictions on the parameter space of {η_i} and {w_i}.
For the diffractive gluon density, we set η_g to 0 and β_g to 0.5. These only marginally limit the freedom in the input functional form for the gluon density. We fixed the γ_q, η_q and γ_g to their best fit values. The lack of diffractive DIS data at high-x_, mean that the flux factor is not really well determined in this region. Hence, for the flux factor, we set the w_2 to 0 and fixed other variable {w_3} to their best fit values. In total, these leave us with 7 free parameters in our QCD fit (three for quarks, two for the gluon density and two for the flux factor), which we include later on our diffractive PDFs uncertainty estimations. 

In the following, we discuss the overall fit quality. Figs. <ref> and <ref> show the NLO theory predictions based on the extracted diffractive PDFs for the diffractive reduced cross sections x_σ_r^D(3) (β, Q^2; x_) for two representative bins of x_ = 0.01 and 0.003. The uncertainty bands correspond to the choice of tolerance T=Δχ^2_ global=1 also have been shown as well. These error bands represent the
fit uncertainties derived only from the experimental input. In order to judge the fit quality, our NLO theory predictions are compared with the H1-LRG-2012 <cit.> measurements. It would be also interesting to examine the fit quality in comparison to the data sets that we did not include in our fit such as old diffractive DIS data measurements at HERA. Hence, in Figs. <ref> and <ref>, we also compare our NLO theory predictions with the H1-LRG-1997 <cit.> measurements. As one can see, prediction based on our diffractive PDFs sets are in reasonably good agreements with the HERA data. The NLO theory prediction based on the GKG18-DPDF <cit.> also have been shown for comparison. From the theory predictions presented in these figures one can conclude that the scale dependence induced by the evolution equations of Eq. (<ref>) is perfectly consistent with the diffractive DIS data. The results clearly indicate that one can use the fracture functions approach to describe diffractive DIS in perturbative QCD at the kinematic region covered by the ep collider HERA as well as other hadron colliders.



The previous figures show the quality of the description of the H1-LRG-2012 data with our NLO theory predictions. It would be interesting to observe that the fit quality to the inclusive H1 and ZEUS combined measurement for the inclusive diffractive DIS cross sections. In the following we discuss the overall fit quality of this data sets. For completeness, in Fig. <ref>, the NLO theory predictions for the diffractive reduced cross sections x_σ_r^D(3) (β, Q^2; x_) have been presented as a function of Q^2 for some selected values of β and for two representative bins
of x_ = 0.05 and 0.075. The H1/ZEUS combined <cit.> diffractive DIS measurements also has been shown for comparison. The uncertainty bands for HK19-DPDF analysis as well as for the GKG18-DPDF are correspond to the choice of tolerance T=Δχ^2_ global=1 and T=Δχ^2_ global=5.  These results show that the quality of the description between our NLO theory predictions and all the H1/ZEUS combined data points analyzed in this study is quite acceptable.



In the rest of this section, we present the resulting diffractive PDFs at NLO and NNLO accuracy. We also compare the HK19-DPDF results with the most recent analysis of GKG18-DPDF <cit.>. We mainly discuss the main difference of the results as well as the effect arising from including the higher order QCD corrections. 
We now compare the diffractive PDFs β  F (β, Q^2; x_) obtained from this analysis with other results in literature, and discuss how their central values and uncertainties vary. As we presented in Sec. <ref>, this analysis is based on a standard “parameter-fitting” criterion, and we plan to present the uncertainty bands of diffractive PDFs considering the choice of tolerance T=Δχ^2_ global=1 for the 68% (one-sigma) confidence level (CL) uncertainty.

In the following, we now turn to discuss the obtained diffractive PDFs and their uncertainties. In Fig. <ref> detailed comparisons of the NLO and NNLO diffractive quark  β  F_q (β, Q^2; x_)  and gluon β  F_g (β, Q^2; x_) PDFs have been presented at the input scale Q_0^2 = 2 GeV^2 for two different x_ bin of 0.01 and 0.003. The uncertainty bands of diffractive PDFs presented for the choice of tolerance T=Δχ^2_ global=1 for the 68% (one-sigma) confidence level (CL). In order to study the perturbative convergence of the diffractive PDFs upon inclusion of higher order QCD corrections, in Fig. <ref> we also compare our NLO and NNLO determinations among each other. Comparing our results, it can be seen that the NLO and NNLO quark diffractive PDFs β  F_q are similar in size for all range of β. However a small difference can be seen in range of β < 0.1 for the x_ = 0.003. The single most striking observation to emerge from the NLO and NNLO comparisons is for the case of the diffractive gluon PDFs β  F_g. One can see that the diffractive gluon PDFs is affected by including the higher-order QCD corrections. Concerning the shapes of the diffractive gluon PDFs, a number of interesting differences between the NLO and NNLO results can be seen from the comparisons in Fig. <ref>. Significant differences in shape are observed for the small values of β, especially for the range of β < 0.2. In this range, the diffractive gluon PDFs at NLO accuracy have a larger magnitude than the NNLO results.

As one can see the NLO and NNLO uncertainties are similar in size showing that the inclusion of higher-order QCD corrections do not improve the uncertainty.  
As we presented in Tables. <ref> to <ref>, concerning the fit quality of the total diffractive DIS data set, the most noticeable feature is the improvement upon inclusion of higher-order corrections. However, the improvement of the total χ^2/ d.o.f is not significant when going from NLO to NNLO. This demonstrates that the inclusion of the NNLO QCD corrections could slightly improves the fit quality as well as the description of the data.

In Fig. <ref>, we present the quark β  F_q (β, Q^2; x_) and gluon β  F_g (β, Q^2; x_) diffractive PDFs
at the scale of Q^2 = 100 GeV^2  and for the x_ = 0.01. Comparing these two results, it can be seen that the differences between the NLO and NNLO sets are very small, both for central values and uncertainties. 



We are in a position to compare our best-fit NLO diffractive PDFs to their counterparts in the GKG18-DPDF analysis <cit.> which applied the general method used to extract diffractive PDFs from available data by considering a number of assumptions motivated by the Regge phenomenology. In Fig. <ref>, we present the results for the diffractive gluon and quark PDFs at Q^2 = 10 GeV^2 for two selected bin of x_ = 0.01 and 0.003. We see that the overall effects of the new methodology are comparable to those reported by <cit.>, but that they act in different kinematical regions of β and x_ and for different diffractive PDFs. From Fig. <ref>, it is clear that central values move very little while diffractive PDFs uncertainties are slightly increased. For instance, for the light quarks and the gluon diffractive PDFs the impact of the new methodology mainly effect the small regions of β; β < 0.1, where it produces an enhancement for the gluon density and a reduction for the quark density. We should emphasize here that, the QCD analysis of diffractive PDFs motivated by the Regge phenomenology uses the GRV parametrization derived from a fit to pion structure function data <cit.> for the Reggeon parton density, and hence, an additional source of uncertainty need to be taken into account in presenting the diffractive PDFs in the kinematic space of z, Q^2 and x_. Overall, there is satisfactory agreement between the two methodology of diffractive PDFs determinations. Hence, this study strengthens the idea of using the fracture function approach to determine diffractive PDFs from a QCD analysis of diffractive DIS data sets.



Let us conclude this section with a presentation of heavy quark diffractive PDFs obtained in this study. As we discussed in section <ref>, for the calculation of heavy-quark structure functions is performed in the FONLL GM-VFNS <cit.>. In Fig. <ref>, the charm β  c (β, Q^2; x_) and bottom β b (β, Q^2; x_) quark diffractive PDFs obtained from our NLO QCD fits have been shown at Q^2 = 100 GeV^2 and for x_ = 0.01. The error bands shown in these figures ar correspond to the fit uncertainties derived only from the experimental input. The results from GKG18-DPDF analysis <cit.> also presented for comparison. As one can see from these results, the overall agreements are well and only small differences between our NLO results and GKG18-DPDF can be found for all heavy quark diffractive PDFs at lower values of β; β < 0.02.



Through our results presented in section, we have shown that our QCD analyses show good agreements with the results obtained by GKG18-DPDF parameterization. In addition, we found that our theory predictions based on the extracted diffractive PDFs are in satisfactory agreements with the H1 and ZEUS diffractive DIS data sets over a wide range of DIS kinematics. The scale dependence of these data sets is found to be in satisfactory agreement with the one predicted for the fracture function, driven by the DGLAP evolution equations. As a short summary, the results of this research support the idea of extracting the diffractive PDFs using the fracture functions approach. 
The presented QCD-based predictions of diffractive DIS processes also provide an important step towards an improved understanding of such processes and also represent a precise test of the employed theoretical concepts to extract PDFs from a QCD analysis of diffractive DIS observables. In particular, it is observed in this analysis, that for the given kinematical range of the HERA hard diffraction events, higher-order QCD corrections are of crucial importance.










§ SUMMARY AND CONCLUSIONS
 



In summary, in this work, we have presented a set of diffractive PDFs at next-to-leading order (NLO) and next-to-next-leading order (NNLO) accuracy obtained from the most up-to-date diffractive DIS data including the most recent combined data sets from H1 and ZEUS Collaborations. The new combined diffractive DIS data from run II at HERA allow for a very accurate determination of the quark and gluon distributions in a wide range of scaled fractional
momentum β and longitudinal momentum fractions x_. We supplement our best-fit diffractive PDFs parameterizations with the reliable uncertainties obtained according to the “Hessian approach” which allows the experimental uncertainties to propagate to an arbitrary observable such as diffractive DIS cross sections. The theory predictions for the hard-scattering diffractive DIS processes in this analysis maintain the NLO accuracy, and for the first time the NNLO accuracy in QCD and employ the FONLL GM-VNFS, which has been shown to provide an excellent description of the existing diffractive DIS data.

In addition to the points mentioned, the theory framework applied in this analysis features a number of new improvements. We work in the framework of fracture functions, as a new method to extract diffractive PDFs inspired by the fully factorization theorem for diffractive DIS processes, which allows a good description of diffractive DIS cross sections. We have demonstrated that diffractive DIS is consistent within this picture and one can introduced diffractive PDFs accordingly.
We have shown that a simple parameterization form for the diffractive PDFs along with the fully factorized approach for the cross section provide very accurate descriptions of the diffractive DIS data sets measured by H1 and ZEUS collaborations at HERA. The diffractive PDFs extracted for a QCD analysis in the fracture functions are also in satisfactory agreements by other analysis in literature. Finally, our results verify that the scale dependence of the data agrees well with the one predicted by the fracture function formalism. Hence, we can conclude that in the diffractive DIS kinematics analyzed in this study, the fracture functions framework can provide a good understanding of the physical picture of this sort of high energy processes. 

Our analysis can be extended in various different directions. First, our analysis can be extended to the diffractive dijet events at HERA <cit.>. For the future, our main goal is to assess the impact of diffractive dijet productions data on the diffractive PDFs and their uncertainties. More detailed discussions on this new extraction of diffractive PDFs at NLO accuracy will be presented in our next study. Second, in terms of future work, it would be interesting to repeat the analysis described here and present a combined QCD analysis of recent data sets measured by the H1 and ZEUS collaborations at HERA on the diffractive DIS and leading-nucleon productions <cit.>. The LO, NLO and NNLO diffractive PDFs sets β  F (β, Q^2; x_) presented in this work are available in the LHAPDF format <cit.> from the author upon request. 











The author is especially grateful to Muhammad Goharipour, Maryam Soleymaninia and S. Atashbar Tehrani for carefully reading the manuscript and fruitful discussions. Author thanks University of Science and Technology of Mazandaran, and School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM) for financial support of this project.




















 















Gao:2017yyd 
J. Gao, L. Harland-Lang and J. Rojo,
“The Structure of the Proton in the LHC Precision Era,”
Phys. Rept.  742, 1 (2018),

[arXiv:1709.04922 [hep-ph]].








Newman:2013ada 
P. Newman and M. Wing,
“The Hadronic Final State at HERA,”
Rev. Mod. Phys.  86, no. 3, 1037 (2014),

[arXiv:1308.3368 [hep-ex]].









H1:2018flt 
H. Abramowicz et al. [H1 and ZEUS Collaborations],
“Combination and QCD analysis of charm and beauty production cross-section measurements in deep inelastic ep scattering at HERA,”
Eur. Phys. J. C 78, no. 6, 473 (2018),

[arXiv:1804.01019 [hep-ex]].







Ball:2014uwa 
R. D. Ball et al. [NNPDF Collaboration],
“Parton distributions for the LHC Run II,”
JHEP 1504, 040 (2015),

[arXiv:1410.8849 [hep-ph]].










Harland-Lang:2014zoa 
L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne,
“Parton distributions in the LHC era: MMHT 2014 PDFs,”
Eur. Phys. J. C 75, no. 5, 204 (2015),

[arXiv:1412.3989 [hep-ph]].









Accardi:2016qay 
A. Accardi, L. T. Brady, W. Melnitchouk, J. F. Owens and N. Sato,
“Constraints on large-x parton distributions from new weak boson production and deep-inelastic scattering data,”
Phys. Rev. D 93, no. 11, 114017 (2016),

[arXiv:1602.03154 [hep-ph]].










Klein:2018rhq 
M. Klein,
“Future Deep Inelastic Scattering with the LHeC,”

arXiv:1802.04317 [hep-ph].








Mangano:2018mur 
M. Mangano et al.,
“Future Circular Collider : Vol. 1 Physics opportunities,”
CERN-ACC-2018-0056.







Benedikt:2018csr 
M. Benedikt et al.,
“Future Circular Collider : Vol. 3 The Hadron Collider (FCC-hh),”
CERN-ACC-2018-0058.












Armesto:2019gxy 
N. Armesto, P. R. Newman, W. Slominski and A. M. Stasto,
“Inclusive diffraction at the LHeC and FCC-eh,”
arXiv:1901.09076 [hep-ph].









Collins:1997sr 
J. C. Collins,
“Proof of factorization for diffractive hard scattering,”
Phys. Rev. D 57, 3051 (1998)
Erratum: [Phys. Rev. D 61, 019902 (2000)],

[hep-ph/9709499].








Ingelman:1984ns 
G. Ingelman and P. E. Schlein,

Phys. Lett.  152B, 256 (1985).
doi:10.1016/0370-2693(85)91181-5






Helenius:2019gbd 
I. Helenius and C. O. Rasmussen,
“Hard diffraction in photoproduction,”
arXiv:1901.05261 [hep-ph].







Rasmussen:2018dgo 
C. O. Rasmussen and T. Sjstrand,
“Models for total, elastic and diffractive cross sections,”
Eur. Phys. J. C 78, no. 6, 461 (2018),

[arXiv:1804.10373 [hep-ph]].











Aktas:2006hy 
A. Aktas et al. [H1 Collaboration],
“Measurement and QCD analysis of the diffractive deep-inelastic scattering cross-section at HERA,”
Eur. Phys. J. C 48, 715 (2006),

[hep-ex/0606004].








Chekanov:2009aa 
S. Chekanov et al. [ZEUS Collaboration],
“A QCD analysis of ZEUS diffractive data,”
Nucl. Phys. B 831, 1 (2010),

[arXiv:0911.4119 [hep-ex]].








Goharipour:2018yov 
M. Goharipour, H. Khanpour and V. Guzey,
“First global next-to-leading order determination of diffractive parton distribution functions and their uncertainties within the xFitter framework,”
Eur. Phys. J. C 78, no. 4, 309 (2018),

[arXiv:1802.01363 [hep-ph]].








Britzger:2018zvv 
D. Britzger, J. Currie, T. Gehrmann, A. Huss, J. Niehues and R. leb?k,
“Dijet production in diffractive deep-inelastic scattering in next-to-next-to-leading order QCD,”
Eur. Phys. J. C 78, no. 7, 538 (2018),

[arXiv:1804.05663 [hep-ph]].








Aaron:2012zz 
F. D. Aaron et al. [H1 Collaboration],
“Measurement of the Diffractive Longitudinal Structure Function F_L^D at HERA,”
Eur. Phys. J. C 72, 1836 (2012)
[Eur. Phys. J. C 71, 1836 (2011)],

[arXiv:1107.3420 [hep-ex]].








Aaron:2012ad 
F. D. Aaron et al. [H1 Collaboration],
“Inclusive Measurement of Diffractive Deep-Inelastic Scattering at HERA,”
Eur. Phys. J. C 72, 2074 (2012),

[arXiv:1203.4495 [hep-ex]].








Aaron:2012hua 
F. D. Aaron et al. [H1 and ZEUS Collaborations],
“Combined inclusive diffractive cross sections measured with forward proton spectrometers in deep inelastic ep scattering at HERA,”
Eur. Phys. J. C 72, 2175 (2012),

[arXiv:1207.4864 [hep-ex]].









Trentadue:1993ka 
L. Trentadue and G. Veneziano,
“Fracture functions: An Improved description of inclusive hard processes in QCD,”
Phys. Lett. B 323, 201 (1994).










deFlorian:1998rj 
D. de Florian and R. Sassot,
“QCD analysis of diffractive and leading proton DIS structure functions in the framework of fracture functions,”
Phys. Rev. D 58, 054003 (1998),

[hep-ph/9804240].









Ceccopieri:2014rpa 
F. A. Ceccopieri,
“QCD analysis of forward neutron production in DIS,”
Eur. Phys. J. C 74, no. 8, 3029 (2014),

[arXiv:1406.0754 [hep-ph]].









Shoeibi:2017lrl 
S. Shoeibi, H. Khanpour, F. Taghavi-Shahri and K. Javidan,
“Determination of neutron fracture functions from a global QCD analysis of the leading neutron production at HERA,”
Phys. Rev. D 95, no. 7, 074011 (2017),

[arXiv:1703.04369 [hep-ph]].








Shoeibi:2017zha 
S. Shoeibi, F. Taghavi-Shahri, H. Khanpour and K. Javidan,
“Phenomenology of leading nucleon production in ep collisions at HERA in the framework of fracture functions,”
Phys. Rev. D 97, no. 7, 074013 (2018),

[arXiv:1710.06329 [hep-ph]].









Ceccopieri:2016rga 
F. A. Ceccopieri,
“Single-diffractive Drell-Yan pair production at the LHC,”
Eur. Phys. J. C 77, no. 1, 56 (2017),

[arXiv:1606.06134 [hep-ph]].









Collins:2001ga 
J. C. Collins,
“Factorization in hard diffraction,”
J. Phys. G 28, 1069 (2002),

[hep-ph/0107252].









Grazzini:1997ih 
M. Grazzini, L. Trentadue and G. Veneziano,
“Fracture functions from cut vertices,”
Nucl. Phys. B 519, 394 (1998),

[hep-ph/9709452].








Vermaseren:2005qc 
J. A. M. Vermaseren, A. Vogt and S. Moch,
“The Third-order QCD corrections to deep-inelastic scattering by photon exchange,”
Nucl. Phys. B 724, 3 (2005),

[hep-ph/0504242].






Martin:2004xw 
A. D. Martin, M. G. Ryskin and G. Watt,
“A QCD analysis of diffractive deep-inelastic scattering data,”
Eur. Phys. J. C 37, 285 (2004),

[hep-ph/0406224].








Martin:2005hd 
A. D. Martin, M. G. Ryskin and G. Watt,
“Diffractive parton distributions from perturbative QCD,”
Eur. Phys. J. C 44, 69 (2005),

[hep-ph/0504132].








Monfared:2011xf 
S. T. Monfared, A. N. Khorramian and S. A. Tehrani,
“A Global Analysis of Diffractive Events at HERA,”
J. Phys. G 39, 085009 (2012),

[arXiv:1109.0912 [hep-ph]].








Ceccopieri:2007th 
F. A. Ceccopieri and L. Trentadue,
“A New Fracture Function Approach to QCD Initial State Radiation,”
Phys. Lett. B 655, 15 (2007),

[arXiv:0705.2326 [hep-ph]].









Camici:1998bg 
G. Camici, M. Grazzini and L. Trentadue,
“Fracture functions and jet calculus,”
Phys. Lett. B 439, 382 (1998),

[hep-ph/9802438].










Vogt:2004mw 
A. Vogt, S. Moch and J. A. M. Vermaseren,
“The Three-loop splitting functions in QCD: The Singlet case,”
Nucl. Phys. B 691, 129 (2004),

[hep-ph/0404111].







Moch:2004pa 
S. Moch, J. A. M. Vermaseren and A. Vogt,
“The Three loop splitting functions in QCD: The Nonsinglet case,”
Nucl. Phys. B 688, 101 (2004),

[hep-ph/0403192].







Alekhin:2017kpj 
S. Alekhin, J. Blmlein, S. Moch and R. Placakyte,
“Parton distribution functions, α_s, and heavy-quark masses for LHC Run II,”
Phys. Rev. D 96, no. 1, 014011 (2017),

[arXiv:1701.05838 [hep-ph]].









Martin:2009iq 
A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt,
“Parton distributions for the LHC,”
Eur. Phys. J. C 63, 189 (2009),

[arXiv:0901.0002 [hep-ph]].







Bertone:2017djs 
V. Bertone, A. Glazov, A. Mitov, A. Papanastasiou and M. Ubiali,
“Heavy-flavor parton distributions without heavy-flavor matching prescriptions,”
JHEP 1804, 046 (2018),

[arXiv:1711.03355 [hep-ph]].








Thorne:2012az 
R. S. Thorne,
“Effect of changes of variable flavor number scheme on parton distribution functions and predicted cross sections,”
Phys. Rev. D 86, 074017 (2012),

[arXiv:1201.6180 [hep-ph]].









Ball:2013gsa 
R. D. Ball et al. [NNPDF Collaboration],
“Theoretical issues in PDF determination and associated uncertainties,”
Phys. Lett. B 723, 330 (2013),

[arXiv:1303.1189 [hep-ph]].









Thorne:2014toa 
R. S. Thorne,
“The effect on PDFs and α _S(M_Z^2) due to changes in flavour scheme and higher twist contributions,”
Eur. Phys. J. C 74, no. 7, 2958 (2014),

[arXiv:1402.3536 [hep-ph]].







Aivazis:1993kh 
M. A. G. Aivazis, F. I. Olness and W. K. Tung,
“Leptoproduction of heavy quarks. 1. General formalism and kinematics of charged current and neutral current production processes,”
Phys. Rev. D 50, 3085 (1994),

[hep-ph/9312318].








Thorne:1997ga 
R. S. Thorne and R. G. Roberts,
“An Ordered analysis of heavy flavor production in deep inelastic scattering,”
Phys. Rev. D 57, 6871 (1998),

[hep-ph/9709442].








Forte:2010ta 
S. Forte, E. Laenen, P. Nason and J. Rojo,
“Heavy quarks in deep-inelastic scattering,”
Nucl. Phys. B 834, 116 (2010),

[arXiv:1001.2312 [hep-ph]].









Thorne:2006qt 
R. S. Thorne,
“A Variable-flavor number scheme for NNLO,”
Phys. Rev. D 73, 054019 (2006),

[hep-ph/0601245].








Bertone:2013vaa 
V. Bertone, S. Carrazza and J. Rojo,
“APFEL: A PDF Evolution Library with QED corrections,”
Comput. Phys. Commun.  185, 1647 (2014),

[arXiv:1310.1394 [hep-ph]].










Bertone:2017tyb 
V. Bertone et al. [NNPDF Collaboration],
“A determination of the fragmentation functions of pions, kaons, and protons with faithful uncertainties,”
Eur. Phys. J. C 77, no. 8, 516 (2017),

[arXiv:1706.07049 [hep-ph]].








Bertone:2018ecm 
V. Bertone et al. [NNPDF Collaboration],
“Charged hadron fragmentation functions from collider data,”
Eur. Phys. J. C 78, no. 8, 651 (2018),

[arXiv:1807.03310 [hep-ph]].






Tanabashi:2018oca 
M. Tanabashi et al. [Particle Data Group],
“Review of Particle Physics,”,
Phys. Rev. D 98, no. 3, 030001 (2018).










Ball:2018iqk 
R. D. Ball et al. [NNPDF Collaboration],
“Precision determination of the strong coupling constant within a global PDF analysis,”
Eur. Phys. J. C 78, no. 5, 408 (2018),

[arXiv:1802.03398 [hep-ph]].








Chekanov:2004hy 
S. Chekanov et al. [ZEUS Collaboration],
“Dissociation of virtual photons in events with a leading proton at HERA,”
Eur. Phys. J. C 38, 43 (2004),

[hep-ex/0408009].







Chekanov:2008fh 
S. Chekanov et al. [ZEUS Collaboration],
“Deep inelastic scattering with leading protons or large rapidity gaps at HERA,”
Nucl. Phys. B 816, 1 (2009),

[arXiv:0812.2003 [hep-ex]].










Aktas:2006hx 
A. Aktas et al. [H1 Collaboration],
“Diffractive deep-inelastic scattering with a leading proton at HERA,”
Eur. Phys. J. C 48, 749 (2006),

[hep-ex/0606003].









Aaron:2010aa 
F. D. Aaron et al.,
“Measurement of the cross section for diffractive deep-inelastic scattering with a leading proton at HERA,”
Eur. Phys. J. C 71, 1578 (2011),

[arXiv:1010.1476 [hep-ex]].








Chekanov:2009ac 
S. Chekanov et al. [ZEUS Collaboration],
“Measurement of dijet photoproduction for events with a leading neutron at HERA,”
Nucl. Phys. B 827, 1 (2010),

[arXiv:0909.3032 [hep-ex]].








Chekanov:2007aa 
S. Chekanov et al. [ZEUS Collaboration],
“Dijet production in diffractive deep inelastic scattering at HERA,”
Eur. Phys. J. C 52, 813 (2007),

[arXiv:0708.1415 [hep-ex]].









Aaron:2011mp 
F. D. Aaron et al. [H1 Collaboration],
“Measurement of Dijet Production in Diffractive Deep-Inelastic Scattering with a Leading Proton at HERA,”
Eur. Phys. J. C 72, 1970 (2012),

[arXiv:1111.0584 [hep-ex]].









Chekanov:2007rh 
S. Chekanov et al. [ZEUS Collaboration],
“Diffractive photoproduction of dijets in ep collisions at HERA,”
Eur. Phys. J. C 55, 177 (2008),

[arXiv:0710.1498 [hep-ex]].






Aaron:2010su 
F. D. Aaron et al. [H1 Collaboration],
“Diffractive Dijet Photoproduction in ep Collisions at HERA,”
Eur. Phys. J. C 70, 15 (2010),

[arXiv:1006.0946 [hep-ex]].








Andreev:2014yra 
V. Andreev et al. [H1 Collaboration],
“Measurement of Dijet Production in Diffractive Deep-Inelastic ep Scattering at HERA,”
JHEP 1503, 092 (2015),

[arXiv:1412.0928 [hep-ex]].








Andreev:2015cwa 
V. Andreev et al. [H1 Collaboration],
“Diffractive Dijet Production with a Leading Proton in ep Collisions at HERA,”
JHEP 1505, 056 (2015),

[arXiv:1502.01683 [hep-ex]].








Schmidt:2018hvu 
C. Schmidt, J. Pumplin, C. P. Yuan and P. Yuan,
“Updating and Optimizing Error PDFs in the Hessian Approach,”
arXiv:1806.07950 [hep-ph].







Khalek:2018mdn 
R. A. Khalek, S. Bailey, J. Gao, L. Harland-Lang and J. Rojo,
“Towards Ultimate Parton Distributions at the High-Luminosity LHC,”
arXiv:1810.03639 [hep-ph].






Willis:2018yln 
C. Willis, R. Brock, D. Hayden, T. J. Hou, J. Isaacson, C. Schmidt and C. P. Yuan,
“Reduction of the PDF Uncertainty in the High-Mass Drell-Yan Spectrum with ePump,”
arXiv:1809.09481 [hep-ex].








Azzi:2019yne 
P. Azzi et al. [HL-LHC Collaboration and HE-LHC Working Group],
“Standard Model Physics at the HL-LHC and HE-LHC,”
arXiv:1902.04070 [hep-ph].








CidVidal:2018eel 
X. Cid Vidal  et al.,
“Beyond the Standard Model Physics at the HL-LHC and HE-LHC,”
arXiv:1812.07831 [hep-ph].









Giele:1998gw 
W. T. Giele and S. Keller,
“Implications of hadron collider observables on parton distribution function uncertainties,”
Phys. Rev. D 58, 094023 (1998),

[hep-ph/9803393].









Giele:2001mr 
W. T. Giele, S. A. Keller and D. A. Kosower,
“Parton distribution function uncertainties,”
hep-ph/0104052.









Pumplin:2001ct 
J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung,
“Uncertainties of predictions from parton distribution functions. 2. The Hessian method,”
Phys. Rev. D 65, 014013 (2001),

[hep-ph/0101032].











Ball:2017nwa 
R. D. Ball et al. [NNPDF Collaboration],
“Parton distributions from high-precision collider data,”
Eur. Phys. J. C 77, no. 10, 663 (2017),

[arXiv:1706.00428 [hep-ph]].












Ball:2012cx 
R. D. Ball et al.,
“Parton distributions with LHC data,”
Nucl. Phys. B 867, 244 (2013),

[arXiv:1207.1303 [hep-ph]].








Ball:2016neh 
R. D. Ball et al. [NNPDF Collaboration],

Eur. Phys. J. C 76, no. 11, 647 (2016),

[arXiv:1605.06515 [hep-ph]].










Khanpour:2016pph 
H. Khanpour and S. Atashbar Tehrani,
“Global Analysis of Nuclear Parton Distribution Functions and Their Uncertainties at Next-to-Next-to-Leading Order,”
Phys. Rev. D 93, no. 1, 014026 (2016),

[arXiv:1601.00939 [hep-ph]].








Hou:2016sho 
T. J. Hou et al.,
“Reconstruction of Monte Carlo replicas from Hessian parton distributions,”
JHEP 1703, 099 (2017),

[arXiv:1607.06066 [hep-ph]].








Martin:2002aw 
A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne,
“Uncertainties of predictions from parton distributions. 1: Experimental errors,”
Eur. Phys. J. C 28, 455 (2003),

[hep-ph/0211080].








Alekhin:2014irh 
S. Alekhin et al.,
“HERAFitter,”
Eur. Phys. J. C 75, no. 7, 304 (2015),

[arXiv:1410.4412 [hep-ph]].











Stump:2001gu 
D. Stump, J. Pumplin, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai and W. K. Tung,
“Uncertainties of predictions from parton distribution functions. 1. The Lagrange multiplier method,”
Phys. Rev. D 65, 014012 (2001),

[hep-ph/0101051].







Blumlein:2006be 
J. Blumlein, H. Bottcher and A. Guffanti,
“Non-singlet QCD analysis of deep inelastic world data at O(alpha(s)**3),”
Nucl. Phys. B 774, 182 (2007),

[hep-ph/0607200].






James:1975dr 
F. James and M. Roos,
“Minuit: A System for Function Minimization and Analysis of the Parameter Errors and Correlations,”
https://doi.org/10.1016/0010-4655(75)90039-9 Comput. Phys. Commun.  10, 343 (1975);





F. James,
“MINUIT Function Minimization and Error Analysis:  Reference Manual Version 94.1,”
https://inspirehep.net/record/1258343/files/minuit.pdf CERN-D-506, CERN-D506.








Gluck:1991ng 
M. Gluck, E. Reya and A. Vogt,
“Parton distributions for high-energy collisions,”
Z. Phys. C 53, 127 (1992).









Aaron:2010ab 
F. D. Aaron et al. [H1 Collaboration],
“Measurement of Leading Neutron Production in Deep-Inelastic Scattering at HERA,”
Eur. Phys. J. C 68, 381 (2010),

[arXiv:1001.0532 [hep-ex]].








Chekanov:2002pf 
S. Chekanov et al. [ZEUS Collaboration],
“Leading neutron production in e+ p collisions at HERA,”
Nucl. Phys. B 637, 3 (2002),

[hep-ex/0205076].








Chekanov:2008tn 
S. Chekanov et al. [ZEUS Collaboration],
“Leading proton production in deep inelastic scattering at HERA,”
JHEP 0906, 074 (2009),

[arXiv:0812.2416 [hep-ex]].










Buckley:2014ana 
A. Buckley, J. Ferrando, S. Lloyd, K. Nordstrm, B. Page, M. Rfenacht, M. Schnherr and G. Watt,
“LHAPDF6: parton density access in the LHC precision era,”
Eur. Phys. J. C 75, 132 (2015),

[arXiv:1412.7420 [hep-ph]].









