First global next-to-leading order determination of diffractive parton distribution functions and their uncertainties within the {\tt xFitter} framework

We present {\tt GKG18-DPDFs}, a next-to-leading order (NLO) QCD analysis of diffractive parton distribution functions (diffractive PDFs) and their uncertainties. This is the first global set of diffractive PDFs determined within the {\tt xFitter} framework. This analysis is motivated by all available and most up-to-date data on inclusive diffractive deep inelastic scattering (diffractive DIS). Heavy quark contributions are considered within the framework of the Thorne-Roberts (TR) general mass variable flavor number scheme (GM-VFNS). We form a mutually consistent set of diffractive PDFs due to the inclusion of high-precision data from H1/ZEUS combined inclusive diffractive cross sections measurements. We study the impact of the H1/ZEUS combined data by producing a variety of determinations based on reduced data sets. We find that these data sets have a significant impact on the diffractive PDFs with some substantial reductions in uncertainties. The predictions based on the extracted diffractive PDFs are compared to the analyzed diffractive DIS data and with other determinations of the diffractive PDFs.


Introduction
High precision calculations of hard scattering cross sections in lepton-hadron deep inelastic scattering (DIS) and hadron-hadron collider experiments can be done within the framework of perturbative quantum chromodynamics (pQCD). The computations of cross sections can be performed using the so-called factorization theorem that allows for a systematic separation of perturbative and nonperturbative physics [1,2]. Some examples for describing the latter in various processes are the well-known parton distribution functions (PDFs) [3][4][5][6][7], nuclear PDFs [8][9][10][11], and polarized PDFs [12][13][14][15][16][17][18], which are rather tightly constrained by global QCD fits to DIS and hadron collider data. In fact, they are crucial assets in all scattering processes involving hadrons (nucleons and nuclei) in the initial state. In this respect, phe-nomenological and experimental studies over the past three decades have provided important information on the structure of hadrons. A significant amount of PDF sets has been determined considering the most precise data from LHC Run I and II [3,5,7,[19][20][21][22][23][24]. In the literature, the relative importance of LHC data has been subject to considerable discussion. These new and up-to-date sets of PDFs have played an important role in the search for new physics, for example in the top quark and Higgs boson sectors [3,25].
Diffractive processes, ep → ep X, where X represents hadronic final state separated from the recoiled proton by a rapidity gap and the proton in the final state carries most of the beam momentum (see Fig. 1), have been studied extensively in the H1 and ZEUS experiments at the electron-proton (ep) collider HERA [2,[26][27][28][29][30][31]. At HERA, a substantial fraction of up to 10% of all ep DIS interactions proceeds via the diffractive scattering process initiated by a highly virtual photon. In the framework of the collinear factorization theorem, the theoretical calculation of diffractive cross sections requires a special type of nonperturbative functions as input, so that the universal diffractive PDFs may be defined. To be more precise, the factorization theorem predicts that the cross section can be expressed as the convolution of nonperturbative diffractive PDFs and partonic cross sections of the hard subprocess calculable within the framework of pQCD. Consequently, the dynamics of the diffractive processes can be formulated in terms of quark and gluon densities. The diffractive PDFs have properties similar to the PDFs of the free nucleon, but with the constraint of a leading proton or its low mass excitations being present in the final state. Like PDFs, it is well established that the diffractive PDFs are universal quantities, which can be extracted from diffractive DIS data through global QCD analyses. The knowledge of diffractive PDFs for different hadron species as well as the estimation of their uncertainties is therefore vital for precise theoretical Fig. 1 Representative Feynman diagram for the neutral current diffractive DIS process ep → ep X and experimental calculations and, hence, has received quite some interests in the past (see, for example, Ref. [32] for a recent review).
The main sources to constrain the diffractive PDFs are the inclusive diffractive DIS data measured at HERA. Given the diffractive PDFs, perturbative QCD calculations are expected to be applicable to other processes such as the jet and heavy quark production in diffractive DIS at HERA [29][30][31][33][34][35]. A full discussion of diffractive dijet production at HERA will be the main subject of our future work. Indeed, the next-toleading order (NLO) QCD predictions using diffractive PDFs describe these measurements rather well. There are several studies in which the diffractive PDFs have been determined from the QCD analyses of diffractive DIS data [27,28,[36][37][38][39][40][41]. In this paper, we present a new set of diffractive PDFs, referred to as GKG18-DPDFs, through a comprehensive NLO QCD analysis. The GKG18-DPDFs diffractive PDFs are determined using all available and up-to-date data from diffractive DIS cross section [42][43][44], including, for the first time, the H1 and ZEUS combined inclusive diffractive cross section measurements [45].
The outline of this paper is as follows: In Sect. 2.1, we briefly present the theoretical formalism adopted for describing the diffractive DIS at HERA. After reviewing the QCD factorization theorem in Sect. 2.2, we explain the heavy flavor contributions to the diffractive DIS structure function in Sect. 2.3. The phenomenological framework used in GKG18-DPDFs global QCD analysis is presented in Sect. 3. This section includes our parametrizations of the diffractive PDFs (Sect. 3.1), a detailed discussion of the description of different data sets included in GKG18-DPDFs global fit (Sect. 3.2), and the method of minimization and diffractive PDF uncertainties (Sect. 3.3). In Sect. 4, we present GKG18-DPDFs results for diffractive PDFs obtained from global fits to H1 diffractive DIS cross sections [42][43][44], and H1 and ZEUS combined inclusive diffractive data [45]. In Sect. 4.1, we compare the diffractive PDFs obtained in this work to the previously determined by other groups. Section 4.2 is also devoted to comparing the theoretical predictions based on the extracted diffractive PDFs with the analyzed diffractive DIS data. Finally, in Sect. 5, we present our summary and conclusions.

Theoretical framework and assumptions
In the following we describe the standard theoretical framework adopted for the diffractive DIS. Although, there are different theoretical approaches to describe the diffractive processes in literature [46], it is well known now that the approach, where the diffractive DIS is mediated by the exchange of the hard Pomeron and a secondary Reggeon can be remarkably successful for the description of most of diffractive DIS data.

Cross section for diffractive DIS
In order to discuss the cross section for diffractive DIS, one needs to introduce the kinematic variables first. The common variables in any DIS process are as follows: the photon virtuality Q 2 = −q 2 , where q = k − k is the difference of the four-momenta of the incoming (k) and outgoing (k ) leptons; the longitudinal momentum fraction x = −q 2 2P.q , where P is the four-momentum of the incoming proton; and the inelasticity y = P.q P.k . The representative Feynman diagram for the neutral current diffractive DIS process ep → ep X, proceeding via a virtual photon exchange, is depicted in Fig. 1. In the case of diffractive DIS, as illustrated in Fig. 1, the additional variables are the squared four-momentum transferred t = (P − P ) 2 , where P is the four-momentum of the outgoing proton, and the mass M X of the diffractive final state, which is produced by diffractive dissociation of the exchanged virtual photon. This mass is much smaller than the invariant photon-proton energy and should be considered as a further degree of freedom. It is usually replaced by the light-cone momentum fraction of the diffractive exchange β, The t-integrated differential cross section for the diffractive process, ep → ep X, is presented in the form of a diffractive reduced cross section σ where x IP = (P−P ).q P.q refers to the longitudinal momentum fraction lost by the incoming proton, which is carried away by the diffractive exchange; and t is the four-momentum transfer squared at the proton vertex. Note that the longitudinal momentum fraction β of the struck parton with respect to the colourless exchange can be also expressed as β = x/x IP . The diffractive reduced cross section is given by

QCD factorization theorem
It has been shown that the diffractive DIS cross sections at HERA [27,28,30] are well interpreted assuming the "proton vertex factorization" approach which provides a good description of diffractive DIS data in terms of a resolved Pomeron (IP) [47,48]. Within the Regge phenomenology [49], the cross sections of diffractive processes at high energies are described by the exchange of so-called Regge trajectories. The diffractive cross section is dominated by a trajectory usually called the Pomeron, while the subleading Reggeon (IR) contribution is significant only for x IP > 0.01. It has been shown that the QCD factorization theorem and the well-known DGLAP parton evolution equations can be applied to describe the dependence of the cross section on β and Q 2 , while a Regge inspired approach is used to express the dependence on x IP and t.
In the QCD factorization approach, the diffractive structure functions can be written as a convolution of hard scattering coefficient functions with the diffractive PDFs, where the sum runs over quarks and gluons. Considering QCD factorization theorem, various hard scattering diffractive processes are calculable by means of diffractive PDFs, such as the diffractive jet production in DIS. The concept of QCD hard factorization of the diffractive PDFs as well as the validity of the assumption of QCD hard factorization have been theoretically predicted to hold in diffractive DIS processes [1]. We should mentioned here that the hard QCD factorization has been tested at HERA in various diffractive processes. In recent H1 analyses the validity of the hard factorization has been successfully examined for open charm production in photoproduction and DIS with D mesons [29,50] and in diffractive production of dijets in DIS [30,34,35,51]. These studies support the validity of QCD hard scattering factorization in diffractive DIS.
We should notice here that in DGLAP NLO QCD global fits, NLO contributions to the splitting functions governing the evolution of unpolarized nonsinglet and singlet combinations of quark densities are the same as in fully inclusive DIS. Hence, the diffractive parton densities satisfy the same (DGLAP) evolution equations as the usual parton distributions in inclusive DIS [52][53][54]. The Wilson coefficient functions C 2 and C L in Eq. (4) are also the same as in inclusive DIS and calculable in perturbative QCD [55]. The diffractive PDFs f D i (β, Q 2 ; x IP , t) are universal and non-perturbative quantities, which can be obtained from the QCD fit to the inclusive diffractive data. Note that diffractive PDFs can be defined in terms of matrix elements of quark and gluon operators; the renormalization of divergencies at next-to-leading order is carried out similarly to the inclusive case and leads to the DGLAP evolution equations.
In GKG18-DPDFs analysis, the proton vertex factorization [47] is assumed, where the x IP and t dependencies of the diffractive PDFs factorize from the dependencies on β and Q 2 . In this framework, the diffractive PDFs can be written as, where f i/IP (β, Q 2 ) and f IR i/IR (β, Q 2 ) are the partonic structures of Pomeron and Reggeon, respectively. The emission of Pomeron and Reggeon from the proton can be described by the flux-factors of f IP/ p (x IP , t) and f IR/ p (x IP , t). The detail discussion on the parametrization of the diffractive PDFs in Eq. (5) will be presented in a separate section.

Heavy flavour contributions to the diffractive DIS structure function
In this section, we discuss a general framework for the inclusion of heavy quark contributions to diffractive DIS structure functions. The correct treatment of heavy quark flavours in an analysis of diffractive PDFs is essential for precision measurements at DIS colliders as well as for the LHC phenomenology. As an example, the cross section for the Wboson production at the LHC depends crucially on precise knowledge of the charm quark distribution. A detailed discussion on the impact of the heavy quark mass treatments in the parton distributions as well as the determination of the their uncertainty due to uncertainty in the heavy quark masses can be found in Ref. [56]. Like to the case of inclusive DIS, the treatment of heavy flavours has an important impact on the diffractive PDFs extracted from the global analysis of diffractive DIS, due to the heavy flavour contribution to the total structure function at small values of z. Recall that there are various choices that can be used to consider the heavy quark contributions. These are the so-called variable flavour number scheme (VFNS), fixed flavour number scheme (FFNS) and general-mass variableflavor-number scheme (GM-VFNS).
In the case of FFNS, Q 2 m 2 c , m 2 b , the massive quark may be regarded as being only produced in the final state and not as partons within the nucleon. Hence, the light up-, down-and strange-quarks are active partons and the number of flavours is fixed to n f = 3. However one can also con-sider charm or bottom quark as light quark at high scales. It has been shown that the accuracy of the FFNS becomes increasingly uncertain as Q 2 increases above the heavy quark mass threshold m 2 H [57]. In the zero-mass VFNS, the massive quarks behave like massless partons for Q 2 m 2 c , m 2 b . The ZM-VFNS misses out O(m 2 H /Q 2 ) contributions completely in the perturbative expansion, and hence, this scheme is not accurate enough to be used in a QCD analysis. One can also see a discontinuity in the parton distributions and total structure function at Q 2 = m 2 H in ZM-VFNS [57]. The GM-VFNS is the appropriate scheme to interpolate between these two regions and could correct FFNS at low Q 2 and ZM-VFNS at high Q 2 → ∞, and hence, could improve the smoothness of the transition region where the number of active flavours is changed by one [57]. Therefore, for a precise analysis of structure functions and other inclusive DIS or hadron colliders data, one can use the GM-VFNS, which smoothly connects the two well-defined scheme of VFNS and FFNS [57]. This scheme is that most commonly approach in variety of global fits. In H1-DPDFs-2006 [27] and ZEUS-DPDFs-2010 [28] diffractive PDFs analyses, the heavy quark structure functions have been computed using the FFNS and general-mass variable-flavor-number scheme of Thorne and Roberts (TR GM-VFNS), respectively. Our approach is based on the TR GM-VFNS [5,58,59] which extrapolates smoothly from the FFNS at low Q 2 to the ZM-VFNS at high Q 2 and produces a good description of the effect of heavy quarks on structure functions over the whole range of Q 2 .
In our analysis, we follow the MMHT14 PDFs analysis and adopt their default values for the heavy quark masses as m c = 1.40 and m b = 4.75 GeV [60]. In Ref. [60], the variation in the MMHT14 PDFs when the heavy quark masses m c and m b were varied away from their default values of m c = 1.40 and m b = 4.75 GeV has been investigated. The dependence of the MMHT14 PDFs and the quality of the comparison to analyzed data, under variations of the heavy quark masses away from their default values has been studied. It has been shown that the effects of varying m c and m b in the predictions of cross sections for standard processes at the LHC are small and the uncertainties on PDFs due to the variation of quark masses are not hugely important [60].

The method of diffractive PDFs global QCD analysis
In the following, we present the method of GKG18-DPDFs global QCD analysis. This section also includes our parametrizations of the diffractive PDFs, the detailed discussion of the description of different data sets included in our global fit, and the method of minimization and uncertainties of our resulting diffractive PDFs.

GKG18-DPDFs parametrizations of the diffractive PDFs
As we already mentioned, the scale dependence of the distributions f i=q,g (β, Q 2 ) of the quarks and gluons can be obtained by the DGLAP evolution equations, provided the diffractive PDFs are parametrized as functions of β at some starting scale Q 2 0 . In our analysis, the diffractive PDFs are modelled at the starting scale Q 2 0 = 1.8 GeV 2 (below the charm threshold) in terms of quark z f q (z, Q 2 0 ), and gluon z f g (z, Q 2 0 ) distributions. Here, z is the longitudinal momentum fraction of the struck parton, which enters the hard subprocess, with respect to the diffractive exchange. Considering the lowest-order quark-parton model process, we have z = β, while the inclusion of higher-order processes leads to 0 < β < z. For the quark distributions we assume that all light-quarks and their antiquarks distributions are equal, The heavy quark distributions f q (=c,b) are generated dynamically at the scale Q 2 > m 2 c,b above the corresponding mass threshold in the TR GM-VFN scheme.
Due to the significantly smaller amount of data for inclusive diffractive DIS data than for the total DIS cross section, we adopt a slightly less flexible, more economical functional form to parametrize the nonperturbative diffractive PDFs at the initial scale Q 2 0 = 1.8 GeV 2 . Our standard parametrizations for the quarks and gluon diffractive PDFs are as follows: An additional factor of e − 0.001 1−z is included to ensure that the distributions vanish for z → 1. Therefore, the parameters γ q and γ g have the freedom to take negative as well as positive values in the fit. We have tested that Eqs. (6) and (7) nevertheless yield a very satisfactory description of the analyzed diffractive DIS data. We found that the two parameters η q and η g had to be fixed to zero since the data do not constrain them well enough. These simple functional forms with significantly fewer parameters have the additional benefit of greatly facilitating the fitting procedure.
The (5) is parametrized by the Pomeron and Reggeon flux factors where the trajectories are assumed to be linear, α IP,IR (t) = α IP,IR (0) + α IP,IR t. The Pomeron and Reggeon intercepts, α IP (0) and α IR (0), and the normalization of the Reggeon term, A IR , are free parameters and should be extracted from the fit to data. Note that the value of the normalization parameter A IP is absorbed in α q and α g .
The Reggeon parton densities f IR i/IR (z, Q 2 ) presented in Eq. (5) are obtained from the GRV parametrization derived from a fit to pion structure function data [61]. The values of the parameters, which are fixed in GKG18-DPDFs fit, are the following: These values are taken from the following experimental measurements [26,62], In total, 9 free parameters are left in GKG18-DPDFs QCD analysis, which are

Diffractive DIS data sets used in GKG18-DPDFs fits
In this section, we present the new experimental data and their treatment in GKG18-DPDFs diffractive PDFs analysis. After reviewing the analyzed data sets, which include the recent H1 and ZEUS combined data, we discuss each of the new data sets in turn. We finally review the way in which the total diffractive DIS data sets are constructed and, in particular, which data and which cuts are included.
A list of all diffractive DIS data points used in GKG18-DPDFs global analysis is presented in Tables 1 and  2. These tables correspond to our two different scenarios for including inclusive diffractive DIS data in GKG18-DPDFs global analyses, namely Fit A and Fit B.
For each data set presented in these tables, we have provided the corresponding references, the kinematical coverage of β, x IP , and Q 2 and the number of data points. We strive to include as much of the available diffractive DIS experimental data as possible in our diffractive PDF analysis. However, some cuts have to be applied in order to ensure that only proper data are included in the analysis.
The first data set we have used in our QCD analysis is the inclusive diffractive DIS data from H1-LRG-11, which were taken with the H1 detector in the years 2006 and 2007. These data correspond to three different center-of-mass energies of √ s = 225, 252 and 319 GeV [42,43]. In this measurement, the reduced cross sections have been measured in the range In addition to the H1-LRG-11 data set, we have used for the first time the H1-LRG-12 data, where the diffractive process ep → eXY with M Y < 1.6 GeV and |t| < 1 GeV 2 has been studied with the H1 experiment at HERA [44]. This high statistics measurement covering the data taking periods 1999-2000 and 2004-2007, has been combined with previously published results [27] and covers the range of 3.5 < Q 2 < 1600 GeV 2 , 0.0017 ≤ β ≤ 0.8, and 0.0003 ≤ x IP ≤ 0.03.
Finally, for the first time, we have used the recent and upto-date H1/ZEUS combined data set for the reduced diffractive cross sections, σ [45]. This measurement used samples of diffractive DIS ep scattering data at a centre-of-mass energy of √ s = 318 GeV and combined the previous the H1 FPS HERA I [63], H1 FPS HERA II [64], ZEUS LPS 1 [65] and ZEUS LPS 2 [26] data sets. This combined data cover the photon virtuality range of 2.5 < Q 2 < 200 GeV 2 , 3.5 × 10 −4 < x IP < 0.09 in proton fractional momentum loss, 0.09 < |t| < 0.55 GeV 2 in squared four-momentum transfer at the proton vertex, and 1.8 × 10 −3 < β < 0.816.
While all H1-LRG data are given for the range |t| < 1 GeV 2 , the combined H1/ZEUS diffractive DIS, which is based upon proton-tagged samples, are restricted to the range 0.09 < |t| < 0.55 GeV 2 , so one needs to use a global normalization factor between those two measurement regions.
Assuming an exponential t dependence of the inclusive diffractive cross section, the extrapolation from 0.09 < |t| < 0.55 GeV 2 to |t| < 1 GeV 2 has been done using the H1 value of exponential slope parameter b 6 GeV −2 [45,64]. The slope parameter can be extracted from fits to the reduced cross section x IP σ D(4) r . With the above choice of constant slope parameter, a good description of the data over the full x IP , Q 2 and β range is obtained [63,64].
In addition to the extrapolation discussed above, distinct methods have been employed by the H1 and ZEUS experiments, and hence, cross sections are not always given with the corrections for proton dissociation background. The different contributions from proton dissociation in the different data sets should be considered by application of different global factors. Proton dissociation is simulated using an approximate dσ [27,41]. The combined H1/ZEUS diffractive DIS are corrected by a global factor of 1.21 to account for such contributions.
It should be noted that the two data normalization factors, which we described above, bring a small systematic uncertainty to the fitted data. However, since the extrapolation in |t| is rather modest and the slope parameter b is experimentally determined with better than 10% accuracy [63] and the factor due to proton dissociation is rather well-constrained phenomenologically and experimentally, this uncertainty is at the level of a few percent. Hence, it can be safely neglected compared to the total experimental error of the H1/ZEUS combined data [45]. To determine our diffractive PDFs, we apply β ≤ 0.80 over the data sets. The data with M X > 2 GeV are included in the fit and the data with Q 2 < Q 2 min are excluded to avoid regions, which are most likely to be influenced by higher twist (HT) corrections or other problems with the chosen theoretical framework.
To ensure the validity of the DGLAP evolution equations, we have to impose certain cuts on the above mentioned data sets. In order to finalize the cut on Q 2 , the sensitivity of χ 2 to variations in Q 2 > Q 2 min is investigated for data used in the analysis. Considering these χ 2 scans, our full diffractive PDFs fits are repeated for each different Q 2 > Q 2 min cut. In Fig. 2, the dependence of χ 2 per number of degrees of freedom, χ 2 /dof, on the minimum cut value of Q 2 has been presented as a function of Q 2 min for all inclusive diffractive DIS data sets used in GKG18-DPDFs (see Table 1). The Q 2 min dependence is reflected from this plot and no further improvement on χ 2 /dof can be expected for larger value of Q 2 > Q 2 min = 9 GeV 2 . Therefore, the lowest Q 2 data are omitted from our QCD fit and Q 2 min ≥ 9 GeV 2 is applied to the diffractive DIS data sets. We refer this fit to Fit A.
However, this choice is somewhat different from the cut used in Refs. [27,28] (Q 2 min > 8.5 GeV 2 ). Since this issue can be related to the possible tension between the H1-LRG-11 and H1-LRG-12 data sets with the H1/ZEUS combined data in low-Q 2 bins, some further investigations are required. To resolve this issue, we also present similar plots for the H1/ZEUS combined data as well as for all H1 LRG data sets. As one can see from the upper panel of Fig. 3, an improvement on χ 2 per number of data points, χ 2 /N pts , can be expected for larger value of Q 2 > Q 2 min = 16 GeV 2 for the H1/ZEUS combined data. In Fig. 3, we have also shown the same plot for the H1 LRG data sets. This plot clearly shows that the appropriate choice for the case of H1 LRG data sets is Q 2 min > 9 GeV 2 . This fact indicates that the choice of Q 2 min > 9 GeV 2 is still suitable for all data sets excluding the H1/ZEUS combined data. Hence, we repeated our analysis by applying an additional cuts on Q 2 min ≥ 16 GeV 2 for the H1/ZEUS combined and keeping Q 2 min ≥ 9 GeV 2 for other H1-LRG-11 and H1-LRG-12 data sets. We refer this fit to Fit B. The number of data points after all cuts for both Fit A and Fit B are summarized in Tables 1 and 2, respectively. Note that since higher twist (HT) can be potentially large is inclusive diffractive DIS [66], the choice of larger Q 2 min also tends to reduce the HT influence.

The method of minimization and diffractive PDF uncertainties
As we already discussed, GKG18-DPDFs diffractive PDFs are provided at NLO in perturbative QCD and the data used in our fits cover a wide range of β, x IP and Q 2 kinematics.
In order to achieve an accurate theoretical descriptions of both the diffractive PDFs evolution and the hard scattering cross sections, a well-tested software package is necessary. In GKG18-DPDFs analysis, we have used the xFitter [67] which is a standard package for performing the global QCD analysis of PDFs. Fortunately, the necessary tools for making theoretical predictions of the diffractive DIS observables have been implemented in the xFitter, allowing one to perform also a global analysis of diffractive PDFs. For the minimization, χ 2 definition and treatment of experimental uncertainties, we used the methodology implemented in xFitter to determine the unknown parameters of diffractive PDFs. The QCD fit strategy follows closely the one adopted for the determination of the PDFs in the HERAPDF methodology [68,69]. The QCD predictions for the inclusive diffractive cross section are obtained by solving the DGLAP evolution equations at NLO. As we mentioned, the heavy quark coefficient functions are calculated in the TR GM-VFNS [5,58] and the heavy quark masses for charm and beauty are chosen as m c = 1.40 GeV and m b = 4.75 GeV [60]. The strong coupling constant is fixed to the α s (M 2 Z ) = 0.1176 [70] which is close to the best-fit value of NNLO MMHT2014 global PDF analysis, α s (M 2 Z ) = 0.1172 + ±0.0013 [71]. The χ 2 function is minimized using the CERN MINUIT package [72]. The form of the χ 2 minimized during our QCD fits is expressed as follows [69], where μ i is the measured value of inclusive diffractive cross section at point i, and T i is the corresponding theoretical predictions. The parameters δ i,stat , δ i,unc , and γ i j are the relative statistical, uncorrelated systematic, and correlated systematic uncertainties. The nuisance parameters b j are associated to the correlated systematics which are determined simultaneously with the unknown parameters {ξ k } of our functional forms of Eq. (6) and (7). We minimize the above χ 2 value with the k = 9 unknown fit parameters {ξ k } of our diffractive PDFs. Table 3 contains the final results of χ 2 /N pts for our global fits. For each data set, the value of χ 2 /N pts has been presented for both Fit A and Fit B. In the last row of the table, the values of χ 2 /dof have also been presented as well. These table illustrates the quality of our QCD fits to inclusive diffractive cross section at NLO accuracy in terms of the individual χ 2 values obtained for each experiment. For Fit A and Fit B, we obtain χ 2 of 322 and 280 with the total 289 and 263 data points, respectively. As one can see from this Table, a Q 2 min ≥ 16 GeV 2 cut on the H1/ZEUS combined data set significantly reduces the χ 2 /N pts from 128/96 to 85/70. Note also that the values of χ 2 /N pts for H1-LRG-11 data sets at √ s = 225 and 252 GeV do not change from Fit A to Fit B and just a very small reduction is observed for the H1-LRG-11 ( √ s = 319 GeV) and H1-LRG-12 data sets. In conclusion, the quality of Fit B is slightly better than that of Fit A, indicating a better description of the inclusive diffractive DIS data. A substantial part of the improvement in the description is driven by the H1/ZEUS combined data.
In order to obtain the uncertainties on the diffractive PDFs, we use the xFitter framework, which includes both the experimental statistical and systematic errors on the data points and their correlations in the definition of the χ 2 function. The uncertainties on the diffractive PDFs as well as the corresponding observables throughout our analysis are computed using the standard "Hessian" error propagation [57,73,74].

Results and discussions
Key results of the current NLO diffractive PDFs fit compared to all previous analyses are the inclusion of all new and up-to-date experimental diffractive DIS data, in particular, the H1/ZEUS combined data set [45], and the error analysis of the extracted diffractive PDFs. Since these new data sets may have the potential to provide more information on the extracted diffractive PDFs, it is important to precisely study their impact on the diffractive PDFs as well as on their uncertainty bands. The second significant addition is the first determination of the diffractive PDFs in the framework of xFitter [67]. The diffractive PDFs in our fits are parameterized at the input scale Q 2 0 = 1.8 GeV 2 according to Eqs. (6) and (7), which provide considerable flexibility. As we mentioned, the available diffractive DIS experimental data are not sufficient enough to constrain all parameters of such a flexible parameterization. However, due to more precise data from H1/ZEUS combined experiments, an enhanced flexibility is maybe allowed for the quark and gluon parameterizations compared to the H1-2006 and ZEUS-2010 fits. We investigated Eqs. (6) and (7) in our analysis and found that relaxing η g and η q does not cause significant changes to the fit results. Therefore, in our Fit A and Fit B QCD analyses, we set these parameters to zero. The details of the fits are summarized in Table 4, which shows our best fit values of the free parameters. In this table, the values of the fixed parameters of α s (M 2 Z ), m c and m b for our Fit A and Fit B QCD analyses are also listed.
The total quark singlet z (z, Q 2 0 ) = q=u,d,s z[q(z, Q 2 0 ) +q(z, Q 2 0 )] and gluon densities zg(z, Q 2 0 ), obtained from our  ever, in the case of the gluon distribution (right panel), the differences between the two analyses are noticeable almost for all kinematic ranges of z. This result can be considered as a evidence for the existence of a possible tension between the low Q 2 data points of the H1/ZEUS combined data. Note that in our Fit A there are more lower-Q 2 data points of the H1/ZEUS combined data than in our Fit B. Overall, it seems that Fit B can be considered as a more conservative analysis because the tension between these data sets has been decreased as much as possible by imposing a more restrictive cut on the H1/ZEUS combined data.
As a last point, we have shown the rations of z Fit B (z, Fig. 4. As illustrated in this figure, in view of the uncertainties of the obtained diffractive PDFs, there is no significant difference between Fit A and Fit B. Consequently, imposing a more restrictive cut on the H1/ZEUS combined data has a slight impact on the central values of the diffractive PDFs, though they do not reduce the uncertainty of the diffractive PDFs. However, from obtained χ 2 /ndf, one can conclude that the GKG18 predictions describe these data very well, particularly for Fit B.
In summary, despite slightly different central values, Fits A and Fit B have overlapping uncertainty bands and, hence, are compatible. The difference comes from the inclusion of the lower-Q 2 region of the combined H1/ZEUS data and thus reflects the overall compatibility of the used data sets. It is in turn related to a few-percent systematic uncertainty in the relative normalization of the data sets, see our discussion above.
The uncertainties on diffractive PDFs need to be improved in the future for very high precision predictions at present and future hadron colliders. Like the total DIS cross section, the diffractive DIS cross section is directly sensitive to the diffractive quark density, whilst the gluon density is only indirectly constrained through scaling violations. Since the gluons directly contribute to the jet production through the boson-gluon fusion process [34,35,50,51,75], one can use quadrature. The combined H1/ZEUS diffractive DIS data are corrected by a global factor of 1.21 to consider the contributions of proton dissociation processes and also corrected by a global normalization factor to extrapolate from 0.09 < |t| < 0.55 GeV 2 to |t| < 1 GeV 2 as described in the text the measurements of dijet production in diffractive DIS to further constrain the diffractive gluon PDF. As an example of the inclusion of dijet production data in the QCD analysis of the diffractive PDFs, one can refer to the ZEUS analysis [28].

Q 2 evolution and comparison to other diffractive PDFs
Having the optimised values of the free parameters, we study next the shape and behaviour of GKG18-DPDFs diffractive PDFs extracted from Fit A and Fit B analyses with an increase of Q 2 and also compare our results with those of other collaborations, in particular with the ZEUS-2010 Fit SJ and H1-2006 Fit B parton sets.
In order to study the scale dependence of diffractive PDFs, in Fig. 5 we show the obtained total quark singlet z (z, Q 2 ) and gluon zg(z, Q 2 ) densities with their uncertainties at some selected Q 2 values of Q 2 = 6, 20 and 200 GeV 2 . These plots also contain the related results of two previous analyses of diffractive PDFs from H1 [27] and ZEUS [28] Collabora-tions. Note that for the H1 analysis we have used the result of their H1-2006 Fit B, while for the ZEUS analysis, their standard analysis of ZEUS-2010 Fit SJ has been considered for comparison.
As can be seen from Fig. 5, due to the evolution effects, both the quark singlet and gluon distributions are undergone an enhancement at low values of z. For large value of z, one can see a reduction of the diffractive PDFs with an increase of Q 2 . For the gluon distributions (left panels), the results of our Fit A and Fit B are in good agreements with the ZEUS-2010 Fit SJ analysis. However, there are some deviations between our results and the H1 ones, especially at smaller and larger values of z. To summarize, the agreement between our results for the gluon diffractive PDFs and the ZEUS-2010 Fit SJ is somewhat better than for H1-2006 Fit B. The discrepancy between our results and H1 fit can be directly attributed to the inclusion of the H1-LRG-12 and H1/ZEUS combined data sets which is not used in the H1 analysis. Fit SJ are inside the error bands of the two Fit A and Fit B total quark singlet distributions. Overall, we have obtained comparable singlet distribution in comparison to the other groups. According to the obtained results, one can conclude that the preliminary impact of these new data sets on the extracted diffractive PDFs is mostly on the behavior of the quark diffractive PDFs.
We conclude this section by presenting the heavy quark diffractive PDFs determined in this analysis in the TR GM-VFNS. In Fig. 6, the charm z(c +c)(z, Q 2 ) (left) and bottom z(b +b)(z, Q 2 ) (right) quark diffractive PDFs obtained from our NLO QCD fits have been shown at selected Q 2 value of Q 2 = 60 and 200 GeV 2 . The error bands correspond to the fit uncertainties derived only from the experimental input. The results from ZEUS-2010 Fit SJ also presented for comparison. As one can see from these plots, only insignificant differences between our results and ZEUS-2010 Fit SJ can be found for all heavy quark diffractive PDFs at low values of z; z < 0.01.

Comparison to the diffractive DIS data
This section presents a detailed comparison of the theoretical predictions based on our diffractive PDFs extracted from the analyses Fit A and Fit B with the experimental data used in these analyses. Note that for all figures, the error bars shown on the experimental data points correspond to the statistical and systematic errors added in quadrature. It should be noted here that the data points excluded from the analysis with Q 2 ≤ Q 2 min = 9 GeV 2 , due to the requirement cuts mentioned in Sect. 3.2, are not shown in the figures in this section. In addition, note that the HERA combined data are corrected by a global factor of 1.21 to consider the contributions of proton dissociation processes as described in Sect. 3.2. As we discussed in Sect. 3.2, while all H1-LRG data sets have been given for the range of |t| < 1 GeV 2 , the combined H1/ZEUS diffractive DIS data are restricted to  [44]. See the caption of Fig. 8 for further details the 0.09 < |t| < 0.55 GeV 2 range. Hence all the combined H1/ZEUS diffractive DIS data sets are corrected by a global normalization factor to extrapolate from 0.09 < |t| < 0.55 GeV 2 to |t| < 1 GeV 2 .
In the following, using the results of Fit A and Fit B, we compare the reduced diffractive cross section In the case of H1-LRG-2011 data [42,43], we present in Fig. 14 Fig. 11 The results of our NLO pQCD fit based on Fit B for the reduced diffractive cross section x IP σ D(3) r as a function of β for x IP = 0.03 in comparison with H1-LRG-2012 data [44]. See the caption of Fig. 8 for further details function of β for x IP = 0.003 and Q 2 = 11.5 GeV 2 in comparison with H1-LRG-2011 data at √ s = 225 GeV (left) and 319 GeV (right). The error bars on the data points and the yellow bands represent the uncorrelated uncertainties and the total uncorrelated and correlated uncertainties, respectively. As can be seen, in the kinematics considered, the theory is again in good agreement with the experiment. From the results presented in this section, one can conclude that our NLO QCD predictions based on the DGLAP approach and using diffractive PDFs extracted from our QCD analysis of inclusive diffraction DIS data describe all analyzed data well.

Summary and conclusions
In this paper, we have presented GKG18-DPDFs, the first global QCD analysis of diffractive PDFs that makes use of the H1/ZEUS combined and the most recent H1 data sets on the reduced cross section of inclusive diffractive DIS. Previous determinations of non-perturbative diffractive PDFs in the parton model of QCD [27,28,41] were based on the older diffractive inclusive DIS data from H1 and ZEUS collaboration. The advent of precise data from the H1 [42][43][44] and H1/ZEUS combined [45] data sets as well as the widely used xFitter package offer us the opportunity to obtain a new set of diffractive PDFs. The TR GM-VFNS provides a rigorous theoretical framework for considering the heavy-quarks We study the impact of the new inclusive diffractive DIS data sets by producing two diffractive PDFs using two different scenarios. Firstly, by considering simultaneously the Q 2 min = 9 GeV 2 cut on all analyzed diffractive DIS data sets, and secondly by removing H1/ZEUS combined data with Q 2 min < 16 GeV 2 in order to investigate possible tension between these data sets at small values of Q 2 . In order to validate the efficiency and emphasize the phenomenological impact of this selection, the differences between these two diffractive PDFs sets are presented and discussed. We find that both of our diffractive PDFs determinations are in very good agreement with the results in the literature for the total quark singlet densities.
We also find differences between our results and the H1-2006 DPDFs fit for the gluon density. There is much better agreement between GKG18 and ZUES-2010 for the gluon density. For the charm and bottom quark densities, there are insignificant discrepancies between GKG18-DPDFs results and ZEUS-2010 for the small values of z; z < 0.01. Our theory predictions based on the determined diffractive PDFs for  [42,43] at √ s = 225 (left) and 319 (right). The error bars on the data points represent the uncorrelated uncertainties and the yellow bands represent the total uncorrelated and correlated uncertainties the reduced diffractive cross section are also in satisfactory agreements with the data sets analyzed as well as with the previous set of H1 data sets. The most significant changes are seen for the heavy quark densities at small values of z and in the increased precision in the determination of the gluon diffractive PDF due to the inclusion of new precise data. For the future, our main aim is to include the very recent diffractive dijet production data, which could provide an additional constraint on the determination of the diffractive gluon density.
A FORTRAN subroutine, which evaluates the leading order (LO) and NLO diffractive PDFs presented here for given values of β, x IP and Q 2 , can be obtained from the authors upon request via electronic mail.