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

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 {\tt 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.

In the past decades, there has been an increasing interest in attempting to understand the structure of hadron [1,2]. 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) [3][4][5][6] 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) [7] and Future Circular Collider in electron-hadron mode (FCCeh) [8,9]. 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. [10] 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 [11] and the validity of proton vertex factorization [12]. 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. [10,13,14] 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 [15], ZEUS-2010-DPDF [16], and the most recent analysis by GKG18-DPDF [17]. All these of diffractive PDFs determinations were performed at nextto-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. [18]. 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 [19,20] including the most recent H1 and ZEUS combined measurement for the inclusive diffractive DIS cross sections [21]. We do so by using a methodology which has been suggested in Refs. [22,23], and recently used to study leading neutron production [24][25][26] as well as inclusive diffractive DIS [27]: 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 2 0 . 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 [24][25][26].
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. II 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, hardscattering factorization and heavy quark contributions. Then in Sec. III we present our input for the diffractive PDFs at a given initial scale Q 2 = Q 2 0 . In Sec. IV 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. V 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. VI. 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. VII gives a brief summary and critique of the findings.

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

A.
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. 1 and 2 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), 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: 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. (2) 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. 2, the longitudinal momentum fractions, x I P 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, ( 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 I P interpreted as the fractional momentum of interacting partons in the protons with respect to the pomeron IP fractional momentum x I P . 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 leptonproton ( p) cross section, "neutral current diffractive reduced cross section" σ , respectively. In the one-photon exchange approximation, it is given by: In the measurements of inclusive diffractive DIS at HERA, the data were presented in terms of σ D(4) r (β, Q 2 ; x I P ). In this analysis, we analyze all available and up-to-date diffractive DIS data sets [19,20] including the most recent data for combined H1 and ZEUS diffractive DIS cross sections measurements [21].

B.
Diffractive parton distributions and hard-scattering factorization According to the factorization theorem for diffractive DIS in pQCD [11,28,29], 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 D i/p (β, µ 2 F ; x I P , 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. (4) are given by The index i in Eq. (5) runs on the flavors of the interacting partons in the proton [27]. 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 [30]. As one can find in literature, for example [15-17, 19-21, 31-33] 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 I P and the invariant momentum transfer t. In the general factorization theorem [11,28,29] for the diffractive DIS events in the form of above equation holds at fixed values of x I P and invariant momentum transfer t [27]. Hence, the parton content of the colour singlet exchange described by diffractive PDFs, F D i/p , is uniquely controlled by the kinematics of the outgoing protons in which carry the fractional momentum of 1 − x I P .
The diffractive PDFs, F D i/p (β, µ 2 F ; x I P , t), appearing in Eq. (5) are known as "proton-to-proton fracture functions" [22] 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 I P and invariant momentum transfer t [27]. 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 [22,34], therefore the diffractive PDFs also obey the standard DGLAP evolution equations [24-26, 34, 35]. The t-unintegrated diffractive PDFs presented in Eq. (5) also obey the standard DGLAP evolution equations [24,27].
By integrating Eq. (5) 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 [22]. 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 fourmomentum 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. (5) are integrated over t in a limited range, they can obey the standard DGLAP evolution equation [27,34]. Hence, we considered such approximation in our analysis and therefore one can use the standard DGLAP evolution when F D i/p (β, µ 2 F ; x I P , t) are integrated over t in a limited range [27,34] Hence, like for the case of PDFs, the evolution equations of diffractive PDFs in diffractive DIS are easily obtained by the DGLAP evolution equations [22,34] as , correspond to the singlet and gluon diffractive distributions, respectively. These non-perturbative distributions need to be parametrize at an input scale Q 2 0 and their evolution to higher scale, Q 2 > Q 2 0 , can be described by using the evolution equation given above. P Σ and P g in Eq. (7) are the common NLO and NNLO contributions to the splitting functions governing the evolution of unpolarized singlet (S) and non-singlet (N S) 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. (7) are the same as in fully inclusive DIS [36,37].
We have now, in principle, all the essential ingredients to write down the diffractive DIS cross sections of Eq. (4) 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. 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 2 h . It has been shown in literature that the GM-VFNS is the most appropriate scheme to extract PDFs from a global QCD analysis [41][42][43]. 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 [44], the Thorne-Roberts (TR GM-VFNS) scheme [45] and the so-called FONLL GM-VNFS scheme [46]. In our previous diffractive PDFs analysis, GKG18-DPDF [17] we adopted the more recent "TR" prescription [47].
In the present analysis, we prefer to use the FONLL GM-VFNS [46] which provides the charm structure functions by using the exact x-space O(α 2 s ) 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 I P ) is performed within the FONLL GM-VFNS by using the publicly available APFEL program [48]. The QCD parameters are the ones quoted in the analysis by NNPDF Collaboration [49,50]. 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 α n F =5 s,NLO (M Z ) = 0.1185 [51]. This selection of α s (M Z ) is consistent with the very recent determination of the strong coupling constant reported by the NNPDF3.1 [52], in which for the first time uses the jet production and tt differential distributions at NNLO accuracy.

III. INPUT DISTRIBUTIONS
In this section, we present the input distributions for the diffractive PDFs at the reference scale of Q 2 = Q 2 0 . We also glance ahead to mention some of the main feature of this selection. As is clear from the discussion in Sec. II B, 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. [25,26], we take for most PDFs a parameterization of the form where Q 2 0 = 2 GeV 2 is the the input scale, and W(x I P ) is the flux factor which we assume that it depends only to the x I P variable. For the diffractive flux factor, we use the standard functional form which can be written as, 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 I P depends very much on the parameterization that one chooses for the flux factor. In Eq. (8), β f q (β, Q 2 0 ) and β f g (β, Q 2 0 ) are the quark and gluon densities at the reference scale Q 2 0 . The quark and gluon densities depend on the scaled fractional momentum variable β and photon virtuality Q 2 .
For the quark density β f q (β, Q 2 0 ), 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. (10). Like to the case of quark density, the poorly determined gluon density β f g (β, Q 2 0 ) is also taken to has a simpler input functional form as presented in Eq. (10). 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. (9) as well as the quark and gluon densities presented in Eq. (10), one can express the diffractive PDFs in terms of the following set of basis functions for the quark and gluon diffractive PDFs, 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. (11) states that the variables x I P , 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. II B, for a fixed value of x I P , the Q 2 evolution of diffractive PDFs β F(β, Q 2 ; x I P ) is given by the DGLAP equations.

IV. 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. I, where we indicate, for each data set: the name of experiments and the corresponding references, observables, the kinematic range of β, x I P 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 [21], the H1-LRG-12 [20] data and finally the H1-LRG-11 [19] 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 [21]. 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. [17] 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 I P < 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 Col- 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 [20]. 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 I P < 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 [15] 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, σ D(4) r (β, Q 2 ; x I P , t). Finally, we have used H1 measurements of the diffractive DIS reduced cross section [19] 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 I P of 5 × 10 −4 < x I P < 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 [15] 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 I P domain. One can not impose any cut at large x I P , although, at present, there are essentially no diffractive DIS data available probing the x I P > 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 [21,56]. 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 [57][58][59][60][61][62][63], 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 [10]. 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.

V. 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 [64][65][66]), 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 [67], as well as to searches for new physics (NP) beyond the standard model (BSM) [68], especially for the future high-luminosity LHC and highenergy 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 [69,70] and the standard Hessian method [39,71]. 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 [4,52,[72][73][74].
The details of "Hessian approach" is fully addressed in Refs. [39,71,[75][76][77]. In the standard "Hessian approach??, smaller number of' "error sets?? are considered to obtain an estimate of the PDFs errors [71] 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 [39,71]. A version of this method, "Hessian profiling", has been included in the xFitter package [78]. Our previous analysis on diffractive PDFs, GKG18-DPDF [17], 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 χ 2 global ({η 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 2 0 = 2 GeV 2 . This function is given by, where n labels a particular data set and w n is a weight factor for the n th experiment with default value of 1 [79,80]. χ 2 n ({η i }) can be written as The minimization of χ 2 global ({η i }) function presented above is done applying the CERN program library MINUIT [81]. In Eq. (13), 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 χ 2 global ({η i }) value with the final 7 fit parameters. The obtained normalization factors N are presented in Table. I Table III: The values of χ 2 /n data for the H1-LRG-12 [20] included in our QCD analysis. See the caption of Table. II for further details.
In order to show the effects arising from the use of the different diffractive DIS data sets, in Tables. II, III, IV,    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. [39,71,76,77]. The diffractive PDFs, βF(β, Q 2 ; x I P ), defined at the initial scale Q 2 0 , are parametrized by η parameters. As we mentioned, the determination of the diffractive PDFs is obtained using a χ 2 n ({η i }) function given in Eq. (13), 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 χ 2 n ({η i }) function. We entitled this as diffractive PDFs set s 0 . The parameters values of s 0 , i.e. {η 0 1 . . . η 0 n }, in which extracted from QCD fit to diffrac-  tive DIS data, will be presented and discussed in details in Sec. VI. By moving away the parameters from their best fitted values, χ 2 increases by the amount of ∆χ 2 where H ij is the Hessian or error matrix which is given by These covariance or Hessian matrix can be obtained by running the CERN program library MINUIT [81]. Having ay hand the derivative of the observable O with respect to each parameter {η}, one can use the following equation 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. (16) ij 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 The displacement of the parameter {η i } from its obtained minimum values η 0 i can be expressed in terms of the rescaled eigenvectors e ik = √ λ k v ik . It reads Considering the orthogonality of eigenvectors υ ik and putting Eq. (18) in (14), one can obtain 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 Finally the neighborhood parameters can be written as 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 [39,71,77] 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. (21). 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 [39,77] P =ˆ∆ where Γ is the Gamma function. For the case of the onefree-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. (23). Considering the parameter number in our diffractive PDFs analysis, Eq. (23) 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.

VI. 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 [17] analysis.
In Table. VII, we present the fit parameters obtained in this work for our NLO and NNLO QCD analyses at the initial scale of Q 2 0 = 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 I P ) 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 [49,50,72]. Specifically, we use α s (M 2 Z ) = 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 re- gion 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 I P dependence of quark and gluon diffractive PDFs presented in Eq. (11). 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 I P , 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.  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 [20] 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 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. 5, the NLO theory predictions for the diffractive reduced cross sections x I P σ D(3) r (β, Q 2 ; x I P ) have been presented as a function of Q 2 for some selected values of β and for two representative bins of x I P = 0.05 and 0.075. The H1/ZEUS combined [21] 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 [17]. 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 I P ) obtained from this analysis with other results in literature, and discuss how their central values and uncertainties vary. As we presented in Sec. V, 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. 6 detailed comparisons of the NLO and NNLO diffractive quark β F q (β, Q 2 ; x I P ) and gluon β F g (β, Q 2 ; x I P ) PDFs have been presented at the input scale Q 2 0 = 2 GeV 2 for two different x I P 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. 6 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 I P = 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. 6. 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. II to VI, 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. 7, we present the quark β F q (β, Q 2 ; x I P ) and gluon β F g (β, Q 2 ; x I P ) diffractive PDFs at the scale of Q 2 = 100 GeV 2 and for the x I P = 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 [17] 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. 8, we present the results for the diffractive gluon and quark PDFs at Q 2 = 10 GeV 2 for two selected bin of x I P = 0.01 and 0.003. We see that the overall effects of the new methodology are comparable to those reported by [17], but that they act in different kinematical regions of β and x I P and for different diffractive PDFs. From Fig. 8, 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 [82] 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 I P . 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 II C, for the calculation of heavyquark structure functions is performed in the FONLL GM-VFNS [46]. In Fig. 9, the charm βc(β, Q 2 ; x I P ) and bottom βb(β, Q 2 ; x I P ) quark diffractive PDFs obtained from our NLO QCD fits have been shown at Q 2 = 100 GeV 2 and for x I P = 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 [17] 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.

VII. 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 I P . 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 [62,63]. 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 [83][84][85]. The LO, NLO and NNLO diffractive PDFs sets β F(β, Q 2 ; x I P ) presented in this work are available in the LHAPDF format [86] from the author upon request.