A global QCD analysis of diffractive parton distribution function considering higher twist corrections within the xFitter framework

We present {\tt SKMHS22}, a new set of diffractive parton distribution functions (PDFs) and their uncertainties at next-to-leading-order accuracy in perturbative QCD within the {\tt xFitter} framework. We describe all available diffractive DIS data sets from HERA and the most recent high-precision H1/ZEUS combined measurements considering three different scenarios. First, we extract the diffractive PDFs considering the standard twist-2 contribution. Then, we include the twist-4 correction from the longitudinal virtual photons. Finally, the contribution of subleading Reggeon exchange to the structure-function $F_2^D$ is also examined. For the contribution of heavy flavors, we utilize the Thorne-Roberts general mass variable number scheme. We show that for those corrections, in particular, the twist-4 contribution allows to include the high-$\beta$ region and leads to a better description of the diffractive DIS data sets. We find that the inclusion of the subleading Reggeon exchange significantly improves the description of the diffractive DIS cross-section measurements. The resulting sets are in good agreement with all diffractive DIS data analyzed, which cover a wider kinematical range than in previous fits. The {\tt SKMHS22} diffractive PDFs sets presented in this work are available via the {\tt LHAPDF} interface. We also make suggestions for future research in this area.

One of the most important experimental findimgs reported by the H1 and ZEUS collaborations at the HERA collider, working at the center of the mass-energy of about √ s = 318 GeV, is the observation of a significant fraction of about 8-10% of large rapidity gap events in a diffractive deep inelastic scattering (DIS) processes [1][2][3][4][5].
Such diffractive DIS processes allow us to define a non-perturbative diffractive parton distribution functions (diffractive PDFs) that can be extracted from a QCD analysis of relevant data [2,3]. According to the factorization theorem, the diffractive cross-section can be expressed as a convolution of the diffractive PDFs and partonic hard scattering cross-sections of the subprocess which is calculable within perturbative QCD.
The diffractive PDFs have properties very similar to the standard PDFs, especially they obey the same standard DGLAP evolution equation [6][7][8][9]. However, they arXiv:2206.13788v1 [hep-ph] 28 Jun 2022 have an additional constraint due to the presence of a leading proton (LP) in the final state of diffractive processes, (k) + p(P ) → (k ) + p(P ) + X(p X ).
Considering the diffractive factorization theorem, the diffractive PDFs can be extracted from the reduced crosssections of inclusive diffractive DIS data by a QCD fit. It should be noted here that if the factorization theorem would be violated in hadron-hadron scattering, then there is no universality for example for the diffractive jet production in hadron-hadron collisions [10,11]. Starting from perturbative QCD, in the first approximation, diffractive DIS is described in the dipole framework and formed by the quark-antiquark (qq) and quark-antiquarkgluon (qq) system. So far all, several groups extracted the diffractive PDFs from QCD analyses of the diffractive DIS data at the next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) accuracy in perturbative QCD [2,3,[12][13][14][15]. In Ref. [14], the authors presented the first NLO determination of the diffractive PDFs and their uncertainties within the xFitter framework [16,17]. Ref. [13] presented the first NNLO determination of the diffractive PDFs, and the framework of fracture functions is used in the QCD analysis [15,18,19]. Some of these studies, such as ZEUS-2010-dPDFs [3], also include the dijet cross-section measurement to determine the wellconstrained gluon PDFs.
In this paper, we report a new QCD fit of diffractive PDFs to the HERA inclusive data in diffractive DIS at the NLO in perturbative QCD within the xFitter framework [16]. The present fit also includes the high-precision H1/ZEUS combined measurements of the diffractive DIS cross-section [20]. The inclusion of the most recent HERA combined data, together with the twist-4 corrections from the longitudinal virtual photons and the contribution of subleading Reggeon exchanges to the structure-function, provide well-established diffractive PDFs sets. We show that these corrections, in particular, the twist-4 contribution allows us to explore the high-β region, and the inclusion of the subleading Reggeons gives the best description of the diffractive DIS data. In addition, by considering such corrections, one could relax the kinematical cuts that one needs to apply to the data. This paper is organized in the following way: In Sec. II we discuss the theoretical framework of the SKMHS22 diffractive PDFs determination, including the computation of the diffractive DIS cross-section, the evolution of diffractive PDFs, and the corresponding factorization theorem. This section also includes our choice of physical parameters and the heavy quark contributions to the diffractive DIS processes. The higher twist contribution considered in SKMHS22 QCD analysis also discussed in detail in this section. In Sec. III we present the details of the SKMHS22 diffractive PDFs global QCD analysis and fitting methodology. Specifically, we focus on the SKMHS22 parametrization, the minimization strategy, and the method of uncertainty estimation. We also present the diffractive data set used in SKMHS22 analy-sis, along with the corresponding observables and kinematic cuts applied to the data samples. In Sec. IV we present in detail the SKMHS22 sets. The perturbative convergence upon inclusion of the higher twist corrections is also discussed in this section. The fit quality and the theory/data comparison are presented and discussed in this section as well. Finally, in Sec. V we summarize our findings and outline possible future developments.

II. THEORETICAL FRAMEWORK
In the following section, we describe the standard theoretical framework which is perturbative QCD (pQCD) for the typical event with a large rapidity gap (LRG) for diffractive DIS processes. We discuss in detail the calculation of diffractive DIS reduced cross-section, the relevant factorization theorem, and our approaches to consider the heavy flavors contributions. We also provide the details of the diffractive structure-function (diffractive SF) taking the twist-4 and Reggeon corrections into account.

A. Diffractive DIS cross section
In Fig. 1, we display the Feynman diagram for diffractive DIS in the single-photon approximation. In the neutral current (NC) diffractive DIS process ep → epX, we have the incoming positron or electron in the initial state in which scatters off an an incoming proton with the fourmomentum k and P , respectively. As one can see from the Feynman diagram, in the final state, the proton with the four-momentum of P remains intact and there is a rapidity gap between the proton in the final state and the diffractive system X and outgoing electron with fourmomentum k .
In order to calculate the reduced diffractive crosssection for such a process, one needs to introduce the standard set of kinematical variables which are the photon virtuality Q 2 , the inelasticity y, and the Bjorken variable x, respectively. In the case of diffractive DIS, one needs to introduce an additional kinematical variable β which is defined to be the momentum fraction carried by the struck parton with respect to the diffractive exchange. The kinematical variable β is given by, where M X is the invariant mass of the diffractive final state, produced by the diffractive dissociation of the exchanged virtual photon, and the variable t = (P − P ) Figure 1: Diagram for diffractive DIS ep → epX. The four-momenta are indicated as well (in teh round brackets). The diffractive scattered proton is distinguished from the diffractive mass X.
is the squared four-momentum transferred at the proton vertex.
The experimental diffractive DIS data sets are provided by the H1 and ZEUS collaborations at HERA in the form of the so-called reduced cross section σ D3 r (β, Q 2 ; x I P , t), where x I P is defined to be the longitudinal momentum fraction lost by the incoming proton. The longitudinal momentum fraction x I P satisfies the relation x = βx I P . The t-integrated differential crosssection for the diffractive DIS processes can be written in terms of the reduced cross-section as, In the one-photon approximation, the reduced diffractive cross-section can be written in terms of two diffractive structure functions [2,3,14]. This reads, It should be noted here that for the y not close to the unity, the contribution of the longitudinal structure function F D(3) L to the reduced cross-sections can be neglected. Since we use the diffractive DIS data sets at HERA for the reduced cross-section, we follow the recent study by GKG18 [14] and consider this contribution in our QCD analysis.

B. Factorization theorem for the diffractive DIS
The main idea of diffractive DIS was proposed for the first time by Ingelman and Schlein [21]. According to the Ingelman-Schlein (IS) model, the diffractive processes in DIS are interpreted in terms of the exchange of the leading Regge trajectory. The diffractive process includes two steps which are the emission of the pomeron from a proton and subsequent hard scattering of a virtual photon on partons in the pomeron. Therefore the pomeron is considered to have a partonic structure as do hadrons. Hence, the diffractive structure functions factorizes into a pomeron flux and a pomeron structure function [22][23][24]. In analogy to the inclusive DIS, the diffractive structure functions can be written as a convolution of the nonperturbative diffractive PDFs which satisfy the standard DGLAP evolution equations [6][7][8][9], and the hard scattering coefficient functions. It is given by, where the sum runs over all parton flavors (gluon, dquark, u-quark, etc.).
Here the long-distance quantity f D i (z, Q 2 ; x I P , t) denotes the non-perturbative part which can be determined by a global QCD analysis of the available diffractive experimental data sets. The Wilson coefficient functions C 2/L,i in the above equation describe the hard scattering of the virtual photon on a parton i and are the same as the coefficient functions known from the inclusive DIS and calculable in perturbative QCD [25]. It has been shown that the description of the experimental data is very good when factorization is assumed [2,3,14]. The vertex factorization states that the diffractive PDFs should be factorized into the product of two terms, one of them depends on the x I P and t, and the another one is a function of β and Q 2 . Hence, the diffractive PDFs f D i/p (β, Q 2 ; x I P , t) is given by, where f I P/p (x I P , t) and f I R/p (x I P , t) are the Pomeron and Reggeon flux-factors, respectively. These describe the emission of the Pomeron and Reggeon from the proton target. The Pomeron and Reggeon partonic structures given by the parton distributions f i/I P (β, Q 2 ) and f I R i/I R (β, Q 2 ). The parametrization and determination of these functions will be discussed in detail in section III.

C. Heavy flavour contributions
The study of the heavy quark flavor contributions to the DIS processes enables us to precisely test QCD and the strong interactions. In this respect, their contributions have an important impact on the PDFs [26][27][28][29], FFs [30][31][32], and diffractive PDFs [12][13][14] extracted from the global QCD analysis.
Generally speaking, there are two regimes for the treatment of heavy quark production. The first region is Q 2 ∼ m 2 h where m h is the heavy quark mass. The massive quarks are produced in the final state and they are not treated as an active parton within the nucleon. This regime is interpreted using the "Fixed Flavour Number Scheme" (FFNS). In this scheme, the light quarks are considered to be the active partons inside the nucleon, and the number of flavors needs to be fixed to n f = 3. The FFNS is not accurate and reliable for scales much greater than the heavy quark mass threshold m 2 h . At higher energy scales, Q 2 m 2 h , the heavy quarks behave as massless partons within the hadron. It that case, logarithmic terms ∼ Q 2 /m 2 h are automatically summed through the solutions of the DGLAP evolution equations for the heavy quark distributions. The simplest approach which describes the Q 2 /m 2 h → ∞ limit is the "Zero Mass Variable Flavor Number Scheme" (ZM-VNS), which ignores all the O(m 2 h /Q 2 ) corrections. In summary, at very large scales, where the resummation of large logarithms are pertinent, the ZM-VFNS would be more precise, and in the limits close to the heavy quark mass threshold m h , the FFNS works well enough. In order to present the correct scheme and to obtain the best description of these two limits of Q 2 ≤ m 2 h and Q 2 m 2 h , one needs to take into account the "General Mass Variable Flavor Number Scheme" (GM-VFNS) [33]. In this scheme, the DIS structure function can be written as follows, where n f is the number of active light quark flavors and m is the number of heavy quarks. Unlike the ZM-VFNS, the hard scattering coefficient functions C FF,n f k depend on the (Q 2 /m 2 h ) but reduce to the zero mass approach as Q 2 /m 2 h → ∞. Considering the transition from n f active quarks to n f + 1 ones, one could write, where the matrix elements A jk (Q 2 /m 2 h ) are calculated and presented in Ref. [34]. The C FF,n f k (Q 2 /m 2 h ) coefficient in the above equation is given by, As we mentioned before, the coefficient functions must behave towards the massless limit as Q 2 /m 2 h → ∞, as given by Eq. (9). The analysis presented in this work is based on the Thorne and Roberts (TR) GM-VFNS which covers smoothly from FFNS scheme at low Q 2 to the ZM-VFNS descritption at high energy scale Q 2 , and in this regard, our choice gives the best description of the heavy flavor effects on the diffractive structure functions.
In this study, following the MMHT14 collaboration, we adopt the values for heavy quarks as m c = 1.40 GeV and m b = 4.75 GeV [35]. The strong coupling constant is fixed to α s (M 2 Z ) = 0.1185 [36].

D. Twist-4 correction
In this section, we describe in detail the twist-4 correction considered in the SKMHS22 QCD analysis. Generally speaking, the higher twist contribution is proportional to the terms depending on 1/Q 2 , and hence, would be highly suppressed at a high-energy scale Q 2 in the DIS process. Nonetheless, in the color dipole framework this term for M → 0 or β → 1 dominates over the twist-2 contribution. Hence, the dipole picture provides a powerful framework in which the QCD-based saturation models can be used to investigate the diffractive DIS data.
For interaction of the color dipole with a proton, one could consider two different models which are the quark dipole (qq system) and gluon dipole (qqg) system. We refer the reader to the Ref. [37] for a clear review.
As we mentioned earlier, for the diffractive crosssection, the virtual photon could be transversely or longitudinally polarized. Therefore, according to the different polarization, one can decompose it into a transverse and a longitudinal part. It turns out for the leading order in Q 2 for β → 1, the qq and qqg from transverse virtual photons vanish proportional to the (1 − β), whereas the longitudinal part which is of higher twist and gives a finite contribution. In diffractive DIS, for M Q 2 , the qq contribution to the final state dominates over qqg of the photon wave function. Thus, the Lqq component is not negligible at higher value of β and has an important contribution in the longitudinal diffractive structure function. The F D Lqq can be written in term of the Bessel functions J 0 and K 0 and dipole cross sectionσ(x I P , r) [37,38], The main idea for the dipole approach is that the photon splits up into a quark-antiquark pair (dipole), which then scatters on the proton target. Following the studies presented in Refs. [39], the dipole cross-section is considered to have the following simple form in our QCD analysis,σ where r indicates the separation between the quark and the antiquark. The saturation momentum Q 2 s = (x I P /x 0 ) −λ GeV 2 is responsible for the transition to the saturation regime. The parameters σ 0 = 29 mb, x 0 = 4 × 10 −5 , and λ = 0.28 are taken from Ref. [39]. This definition of the dipole-proton cross-section presents a good description of the inclusive HERA diffractive DIS data sets. Our QCD analysis which includes the twist-4 correction is called SKMHS22-tw2-tw4. The effect of such corrections on the extracted diffractive PDFs is discussed in section. IV.

E. Reggeon contribution
For the higher value of x I P , these diffractive DIS data sets include the contributions which decrease with energy. In order to truly describe this effect one can include the contributions of the subleading Reggeon which breaks the factorization of the diffractive structure function. In order to take into account such contributions, we add the following Reggeon contribution to the diffractive structure function where the f R (x I P , t) is the Reggeon flux, and the F R 2 (β, Q 2 ) is the Reggeon structure function. In principle, we should consider different Regge pole contributions and sum over them and include interference terms. Approximately, one can neglect the interference terms between Reggeons and the Pomeron and also between different Reggeons. Therefore for the Reggeon flux, we consider the following formulas: and Here, C i (t) = 4 cos 2 [πα i (t)/2] and C i (t) = 4 sin 2 [πα i (t)/2] are the signature factors for the even Reggeon (f 2 ) and the odd Reggeon (ω), respectively. The Reggeon trajectory is given by GeV −2 which denote the Reggeon couplings to the proton. λ 2 i = 0.65 GeV is known from Reggeon phenomenology in hadronic reactions. The analysis for the isoscalar Reggeons f 2 , ω shows that the Reggeon contribution to the diffractive structure-function becomes important for x I P > 0.01 [42]. Our assumption for the Reggeon structure functions is that they are the same for all Reggeons and F R 2 (β, Q 2 ) is related to the parton distributions in the Reggeons in a conventional way and can be determined by the fit to diffractive DIS data from HERA [20,43,44]. For the Reggeon structure-function, we consider the following parametrization form, where w 1 , w 2 , w 3 and w 4 are free fit parameters and they need to be determined from a global QCD analysis. The parameter w 2 controls the shape of F R 2 (β) in the low-β region, while w 3 controls the high-β region. The parameters w 4 and w 5 are considered to be fixed to zero, as the current diffractive DIS data sets do not have enough power to constrain all the shape parameters. Our QCD analysis which includes the Reggeon contribution is called SKMHS22-tw2-tw4-RC, and the effect arising from the inclusion of such correction is discussed in section. IV.

III. DIFFRACTIVE PDFS GLOBAL QCD ANALYSIS DETAILS
The analysis of diffractive PDFs is a QCD optimization problem and can be seen as having four elements to it. The first and most important element is of course the experimental data and their uncertainty. The second element is a theoretical model used to describe these phenomenon. For the analysis of non-perturbative objects such as diffractive PDFs this theoretical framework can have some adjustable parameters as well as a number of constraints (say for ensuring of factorization). Another part of the problem is the objective function which we wish to minimize/maximize. Here, this function is a χ 2 function that is defined by the theoretical framework and includes the data uncertainty. The last part of the problem is the optimization method which is the essential step to finding the solution.
In this section, we present the methodology of our global QCD analysis to obtain the diffractive PDFs and their uncertainties at NLO accuracy in perturbative QCD. We describe the details of the data sets used for the following fit in section III A , the parametrization of diffractive PDFs is discussed in section III B, and the methodology of minimization and uncertainties of our diffractive PDFs are described in section III C.

A. Experimental data sets
In this section, we describe all available inclusive diffractive DIS data sets in detail. However, one needs to apply some kinematical cuts in order to avoid nonperturbative effects and ensure that only the data sets for which the available perturbative QCD treatment is adequate are included in the QCD analysis. We attempt to include more data sets, and hopefully by including the twist-4 and Reggeon contributions, one could relax some kinematical cuts that need to be applied to the data.
The kinematic coverage in the (β, x I P , Q 2 ) plane of the complete SKMHS22 data set is displayed in Fig. 2. The data points are classified by different experiments at HERA. As customary, some kinematic cuts need to be applied to the diffractive DIS cross-section measurements. The list of diffractive DIS data sets and their properties used in the SKMHS22 global analysis are shown Table I: List of all diffractive DIS data points with their properties used in the SKMHS22 global QCD analysis. For each data set we provide the kinematical coverage of β, xIP , and Q 2 . The number of data points is displayed as well. The details of the kinematical cuts applied on these data sets are explained in the text.   in Table I. All the measurements are presented in terms of the t-integrated reduced diffractive DIS cross-section measurements D

D(3) r
(ep → epX). In the SKMHS22 QCD analysis, we also analyze the combined measurement of the inclusive diffractive cross-section presented by the H1 and ZEUS Collaborations at HERA (H1/ZEUS combined) [20]. They used samples of diffractive DIS ep data at the center-of-mass energy of √ s = 318 GeV at the HERA collider where the leading protons are detected by appropriative spectrometers. This high-precision measurement combined all the previous H1 FPS HERA I [42], H1 FPS HERA II [45], ZEUS LPS 1 [46], and ZEUS LPS 2 [1] data sets. These measurements cover the photon virtuality interval 2.5 < Q 2 < 200 GeV 2 , 3.5×10 −4 < x I P < 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. We should highlight here that all H1-LRG data are published for the range of |t| < 1 GeV 2 , while the recent combined H1/ZEUS data sets are restricted to the range of 0.09 < |t| < 0.55 GeV 2 , and hence, one needs to use a global normalization factor to extrapolate from 0.09 < |t| < 0.55 GeV 2 to |t| < 1 GeV 2 [14]. Therefore, the combined H1/ZEUS data are corrected for the region of |t| < 1 GeV 2 .
Another data set that we have used in our QCD analysis is the Large Rapidity Gap (LRG) data from H1-LRG-11, which was measured by the H1 detector in 2006 and 2007. These data are derived for three different center of mass energies of √ s = 225, 252 and 319 GeV [43]. The H1 collaboration measured the reduced cross-section in the photon virtualities range 4 GeV 2 ≤Q 2 ≤44 GeV 2 for the center of mass √ s = 225, 252 GeV, and 11.5 GeV 2 ≤Q 2 ≤44 GeV 2 for the center-of-mass energy of √ s = 319 GeV. The diffractive final state masses and the proton vertex are in the range of 1.25 < M X < 10.84 GeV and |t| < 1 GeV 2 , respectively. The diffractive variables are considered in the range of 5×10 −4 < x I P < 3×10 −3 , 0.033 < β < 0.88 for √ s = 225, 252 GeV, and 0.089 < β < 0.88 for √ s = 319 GeV. Finally the last data set that we have employed is the H1-LRG-12 [44]. These data for the process ep → eXY have been derived by the H1 experiment at HERA. H1-LRG-12 data covers the range of 3.5 < Q 2 < 1600 GeV 2 , 0.0003 ≤ x I P ≤ 0.03, and 0.0017≤β≤0.8.
We applied some kinematical cuts on β, M X , and Q 2 . The kinematical cuts we applied are similar to those used in Refs. [3,14], except for the case of β. By considering the twist-4 and Reggeon corrections, we can relax the cuts to include more data sets than those have been used in Ref. [14]. For the extraction of diffractive PDFs we apply β ≤ 0.90 and M X < 2 GeV over all data sets used in this analysis. The sensitivity of data sets to the Q 2 has been tested in Refs. [3,14]. The authors finalized the cut on Q 2 by making a χ 2 scan. In Ref. [12] the authors used standard higher twist for structure functions and they obtained the best χ 2 by considering the Q 2 min = 6.5 GeV 2 , and in Ref. [14] to extract the diffractive PDFs they found the best cut is Q 2 min = 9 GeV 2 . In this work and after testing the results for the χ 2 we found the best agreement of theory and data will be achieved by taking Q 2 min = 9 GeV 2 . After applying these kinematical cuts the total number of data points is reduced to the 302.

B. SKMHS22 Diffractive PDFs parametrization form
Diffractive PDFs are nonperturbative quantities and cannot be calculated in perturbative QCD. Therefore, for their functional dependence a parametric form with some unknown parameters should be considered at the input scale. Due the lack of diffractive DIS experimental data, for the quark distributions we consider all light quarks and anti-quarks densities to be equal, f u = f d = f s = fū = fd = fs. The scale dependence of the quarks and gluon distributions f q,g (β, Q 2 ) needs to be determined by the standard DGLAP evolution equations. We fit the diffractive quark and gluon distributions at the starting scale Q 2 0 = 1.8 GeV 2 which is below the charm threshold (m 2 c = 1.96 GeV 2 ). Since the diffractive DIS data sets can only constrain the sum of the diffractive PDFs and due to the small amount of the inclusive diffractive DIS data sets, one needs to consider the less flexible parametrization form for diffractive PDFs. We fit the diffractive parton distribution function at the initial scale Q 2 0 = 1.8 GeV 2 with the following pomeron parton distributions which have been used in several analysis [2,12,14]: , (17) where z in above equation is the longitudinal momentum fraction of struck parton with respect to the diffractive exchange. At the lowest order, z = β but by including the higher orders this parameter differs from β and this leads to the 0 < β < z. In order to let the parameters γ q and γ g have enough freedom to achieve negative or positive values in the QCD fit we follow the QCD analyses avalible in the litrature and add an additional term e − 0.001 1−z to ensure that the distributions vanish for z → 1. It turned out that four parameters η q , η g , ξ q and ξ g are not well constrained by the experimental data, and hence, one needs to set them to zero. The heavy flavors are generated through the DGLAP evolution equations at the scale Q 2 > m c,b . As we mentioned, to consider the contributions of heavy flavors, we apply the Thorne-Roberts scheme GM-VFN scheme.
The x I P dependence of the diffractive PDFs f D i (β, Q 2 ; x I P , t) in Eq. 6 is determined by the Pomeron and Reggeon flux with linear trajectories, α I P ,I R (0) + α I P ,I R t [2,12,14] f I P ,I R (x I P , t) = A I P ,I R e B I P ,I R t where the normalization factor of Reggeon A I R and the Pomeron and Reggeon intercepts, α I P (0) and α I R (0) are free parameters and will be determined from fit to the diffractive DIS data. We should note here that, according to Eq (6) the parameter A I P is absorbed into the α q and α g . The remaining parameter appearing in Eq. (18) are taken from [14]. In general we have twelve free fit parameters: six from Eq. (16) and (17), three form Eq. (18), and three from the Reggeon structure function in Eq. (15).

C. Minimization and diffractive PDF uncertainty method
In this section, we present two important parts of the SKMHS22 QCD analysis. First, we discuss the optimization method, and then we show how to take into account the data uncertainty in the final results. In order to estimate the free parameters of the diffractive PDFs given the experimental data of section III A we apply the maximum log-likelihood method. If we assume that the data arise from a Gaussian distribution as is usually done, this method coincides with minimizing the χ 2 estimator. Here we adopt the following form for the χ 2 function [14,47], with E is the measured experimental value, and T is the theoretical prediction based on the fit parameters {ζ k }.
The parameters δ i,stat , δ i,unc , and γ i j denote the relative statistical, uncorrelated systematic, and correlated systematic uncertainties, respectively. The nuisance parameters b k are related to correlated systematic uncertainty and are determined with the {ζ k } parameters simultaneously in the QCD fit.
The above χ 2 function is incorporated in the xFitter framework, in conjunction with other tools required for a perturbative QCD analysis of diffractive PDFs. One then can use this package to perform all the essential operations such as DGLAP evolution up to NNLO accuracy, theoretical calculation of the relevant observables and finding out optimal parameter values and deducing their uncertainty collectively. As we specified above in order to find the optimal fit parameter values one needs to minimize the χ 2 function, this is achieved by utilizing the MINUIT CERN package [48]. This package finds the parameter uncertainties by considering the χ 2 function around its minimum. MINUIT has five minimization algorithms, here we choose to work with MIGRAD which is the most commonly used method of minimization.
In practice, we need to propagate these parameter uncertainties to the observables or other quantities like the diffractive PDFs themselves. For this aim, a set of eigenvector PDF sets along with the the central values of the diffractive PDFs is formed, for our fit that has 9 free parameters, the total number of PDF sets will be 19. Each member of the error set is derived by either increasing or decreasing one of the parameters by its uncertainty. Then, the uncertainty of a quantity O due to its dependence on PDFs is given as in Ref. [49], where O (±) i refer to the values of O which are calculated from PDF sets of the i-th parameter along with the ± directions. In the derivation of this relation, it is assumed that the variation of O can be approximated by linear terms of its Taylor series and then the gradient is approximated, which produces the above result. In the next section, we present our main results and findings of the SKMHS22 diffractive PDFs and some observables in order to show the quality of the analysis.

IV. FIT RESULTS
This section includes the main results of the SKMHS22 diffractive PDFs analysis. As we discussed in detail earlier, in this QCD analysis we present three different QCD fits to determine the diffractive PDFs which are SKMHS22-tw2, SKMHS22-tw2-tw4 and SKMHS22-tw2-tw4-RC. In this section we will present and discuss these results in turn. The similarity and difference of these results over different kinematical ranges will be highlighted, and the stability of the results upon the inclusion of higher-twist corrections will be discussed in detail. This section also includes a detailed comparison with the diffractive DIS data analyzed in this work.
The best fit parameters for three sets of SKMHS22 diffractive PDFs are shown in Table II along with their experimental errors. Considering the numbers presented in this table, some comments are in order. The parameters η g and η q are considered to be fixed at zero as the present diffractive DIS data sets do not have enough power to constrain all shape parameters of the distributions.
The parameter {w i } for the Reggeon structure function for the SKMHS22-tw2-tw4-RC analysis presented in Eq. 15 are determined along the fit parameters and then keep fixed to their best fitted values. The parameters w 4 and w 5 are keep fixed at zero. The values for the strong coupling constant α s (M 2 Z ) and the charm and bottom quark masses also are shown in the table as well.
In the following, we turn our attention to the detailed comparison of the three SKMHS22 diffractive PDFs sets. We display the SKMHS22-tw2, SKMHS22-tw2-tw4 and SKMHS22-tw2-tw4-RC diffractive PDFs parameterized in our QCD fits, Eqs. (16) and (17), along with their uncertainty bands in Fig. 3 at the input scale Q 2 0 = 1.8 GeV 2 . The higher Q 2 value results at 6 and 20 GeV 2 are shown in Figs. 4 and Fig. 5, respectively. The lower panels show the ratio of SKMHS22-tw2-tw4 and SKMHS22-tw2-tw4-RC to the SKMHS22-tw2.
In Fig. 6, we display the perturbatively generated SKMHS22 diffractive PDFs for the charm and bottom quark densities along with their error bands at the scale of Q 2 = 60 GeV 2 and 200 GeV 2 . All three diffractive PDFs sets are shown for comparison.
A remarkable feature of the SKMHS22 diffractive gluon and quark PDFs shown in these figures is the difference both in the shape and error bands, which reflects the effect arising from the inclusion of higher twist corrections. As can be seen, for all cases the inclusion of the twist-4 and Reggeon corrections lead to slightly small error bands. This is consistent with the χ 2 value presented in Table III. The inclusion of the twist-4 corrections and Reggeon contributions leads to an enhancement of the gluon diffractive PDFs at high β, and a reduction of the singlet PDFs as well. Considering such corrections also affects the small regions of β and leads to the reduction of the central values of all distributions. These findings indicate that for the given kinematical range of the  For the gluon PDFs, we see a very large error band for a small region of β, namely β < 0.01, which indicates that the available diffractive DIS data do not have enough power to constrain the low-β gluon density. To better constrain the gluon PDFs, the diffractive dijet productions data need to be taken into account [50,51]. In terms of future work, it would be very interesting to repeat the QCD analysis described here to study the impact of diffractive dijet production data on the diffractive gluon PDF and its uncertainty.
The same discussion also holds for the case of heavy quark diffractive PDFs. As one can see from Fig. 6, the SKMHS22-tw2-tw4-RC charm and bottom quark densities lie below other results for all regions of β.
In Table. III, we present the values for the χ 2 per data point for both the individual and the total diffractive DIS data sets included in the SKMHS22 analysis. The values are shown for all three sets of SKMHS22-tw2, SKMHS22-tw2-tw4 and SKMHS22-tw2-tw4-RC. Concerning the fit quality of the total diffractive DIS data sets, the most noticeable feature is an improvement in the inclusion of the Reggeon contribution. The improvement of the individual χ 2 per data point is particularly pronounced for the H1-LRG 2012 when the Reggeon contribution in considered in the analysis. This finding demonstrates the fact that the inclusion of the higher-twist corrections improves the description of the diffractive DIS data sets.
In order further high-light the SKMHS22 analysis, in the following, we present a detailed comparison of the diffractive DIS data set analyzed in this work to the correspond-ing NLO theoretical predictions obtained using the NLO diffractive PDFs from all three sets. As we will show, in general, overall good agreements between the diffractive data and the theoretical predictions are achieved for all H1 and ZEUS data.
We start with the comparison to the H1-LRG-11 √ s = 225, 252 and 319 GeV data. In Fig. 7, the SKMHS22 theory predictions are compared with some selected data points. The comparisons are shown as a function of β and for selected bins of x I P =0.0005 and 0.003. As one can see, very nice agreement is achieved for all regions of β The same comparisons are also shown in Fig. 8 for the H1-LRG-12 for some selected bins of x I P namely 0.01, 0.03, and 0.003. In these plots, the SKMHS22 theory precision for all three different sets are compared with the H1-LRG-12 diffractive DIS data. The comparisons are shown as a function of β and for different values of Q 2 . As one can see, very good agreements between the data and theory predictions are achieved, consistent with the total and individual χ 2 presented in Table III.
Finally, in Fig. 9, we display a detailed comparison of the theory predictions using the three different sets of SKMHS22 and the H1/ZEUS combined data. The plots are presented as a function of x I P and for different values of β and Q 2 . As one can see, very good agreement is achieved.
As one can see from Figs. 7, 8 and 9, in general, overall good agreement between the diffractive DIS data and the NLO theoretical predictions are achieved for all experiments, which is consistent with the individual and total χ 2 values reported in Table. III. Remarkably, the NLO theoretical predictions and the data are in reasonable agreement in the small-and large-β regions, and the whole range of x I P .  Figure 3: The SKMHS22 diffractive PDFs at the input scale of Q 2 0 = 1.8 GeV 2 . All three diffractive PDFs sets are shown for comparison. The lower panels represent the ratio of SKMHS22-tw2-tw4 and SKMHS22-tw2-tw4-RC to the SKMHS22-tw2.

V. DISCUSSION AND CONCLUSION
In this work, we performed fits of the diffractive PDFs at NLO accuracy in perturbative QCD called SKMHS22 using all available and up-to-date diffractive DIS data sets from the H1 and ZEUS Collaborations at HERA. We present a mutually consistent set of diffractive PDFs by adding the high-precision diffractive DIS data from the H1/ZEUS combined inclusive diffractive cross-sections measurements to the data sample.
In addition to the standard twist-2 contributions, we also considered the twist-4 correction and Reggeon contribution to the diffractive structure function, which dominate in the region of large β. The effect of such contributions on the diffractive PDFs and cross-sections were carefully examined and discussed. The twist-4 correction and Reggeon contribution lead to the gluon distribution which is peaked stronger at high-β than in the case of the standard twist-2 QCD fit. The well-established xFitter fitting methodology, widely used to determine the unpolarised PDFs, are modified to incorporate the higher-twist contributions. This fitting methodology is specifically designed to provide faithful representations of the experimental uncertainties on the PDFs.
The QCD analysis presented in this work represents the first step of a broader program. A number of updates and improvements are foreseen for the future study.
The important limitation of the SKMHS22 QCD analysis is the fact that it is based on the inclusive diffractive DIS cross-section measurements. Despite the fact that the diffractive DIS is the cleanest process for the extraction of diffractive PDF, it is scarcely sensitive to the gluon density. For the near future, our main aim is to include very recent diffractive dijet production data [50,51], which we expect to provide a good constraint on the determination of the gluon diffractive PDFs. This will require the numerical implementation of the corresponding observables at NLO and NNLO accuracy in perturbative QCD in the xFitter package. A further improvement for the future SKMHS22 analyses, as a long-term project, is the inclusion of other observables from hadron colliders which could carry some information on flavor separation. A FORTRAN subroutine, which evaluates the three sets of SKMHS22 NLO diffractive PDFs presented in this work for given values of β, x I P and Q 2 can be obtained from the authors upon request. These diffractive PDFs sets are available in the standard LHAPDF format.