Evidence of quasi-partonic higher-twist effects in deep inelastic scattering at HERA at moderate $Q^2$

The combined HERA data for the inclusive deep inelastic scattering (DIS) cross sections for the momentum transfer $Q^2>1$ GeV$^2$ are fitted within the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) framework at next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) accuracy, complemented by a QCD-inspired parameterisation of twist 4 corrections. A modified form of the input parton density functions is also included, motivated by parton saturation mechanism at small Bjorken $x$ and at a low scale. These modifications lead to a significant improvement of the data description in the region of low $Q^2$. For the whole data sample, the new benchmark NNLO DGLAP fit yields $\chi^2$ / d.o.f. $\simeq 1.2$ to be compared to 1.5 resulting from the standard NNLO DGLAP fit. We discuss the results in the context of the parton saturation picture and describe the impact of the higher-twist corrections on the derived parton density functions. The resulting description of the longitudinal proton structure function $F_L$ is consistent with the HERA data. Our estimates of higher-twist contributions to the proton structure functions are comparable to the leading-twist contributions at low $Q^2 \simeq 2$ GeV$^2$ and $x \simeq 10^{-5}$. The $x$-dependence of the twist 4 corrections obtained from the best fit is consistent with the leading twist 4 quasi-partonic operators, corresponding to an exchange of four interacting gluons in the $t$-channel.


Introduction and conclusions
Good understanding of the proton structure has been one of the fundamental goals of particle physics over recent decades. Measurements of the deep inelastic e ± p scattering (DIS) performed by H1 and ZEUS collaborations at the HERA collider contributed invaluable experimental input into this task. The combined data of H1 and ZEUS [1] that include all the measurements of the proton structure functions provide the most accurate information on the proton structure over wide range of Bjorken x and the momentum transfer Q 2 , in particular at smaller x and Q 2 . Hence it is crucial to fully use these data to extract the precise information on the parton density functions (PDFs) in the proton.
The standard description of the proton structure function in QCD relies on the operator product expansion (OPE) in which only the leading-twist 2 operators-are retained. The twist 2 contributions to proton structure functions obey the hard factorisation theorem that allows to isolate the universal twist 2 parton density functions. The PDFs drive the proton scattering cross sections and the accuracy of the PDFs determination is crucial for the precision of measurements at proton colliders. It follows from the OPE however, that the twist 2 description of proton scattering is subject of higher-twist (HT) corrections that enter with suppression of inverse powers of the hard process scale squared, Q 2 . Those corrections, although quickly decreasing with Q 2 , may affect the determination of the PDFs from the cross sections. In order to avoid determination error of the PDFs it is necessary to include the higher-twist terms in the analysis. Currently not much is known about higher-twist components of the proton structure. The operator content is increasingly complicated with the increasing twist and the available data are not sufficient to perform a clean and straightforward measurement of the higher-twist terms. Fortunately, the model independent characteristics of higher-twist terms given by their Q 2 scaling provides opportunity to obtain some information on the higher-twist corrections from fits to DIS data extended to low Q 2 . From the theory side, properties of the leading twist 4 singularity at small x were investigated [2], corresponding to a quasi-partonic [3] four-gluon exchange. It was found [2] that the energy dependence of the leading exchange is up to 1/N 2 c corrections, given by a double gluonic ladder exchange in the t-channel. Hence, although the overall magnitude of the twist 4 contributions is currently undetermined, the Q 2 and x-dependencies of these contributions are known from theory and may be used as the higher-twist signatures. Estimates of the higher-twist corrections to DIS at small x [4,5] and fits to the DIS data with higher-twist corrections [6] have been performed since many years. Only recently however, with the most precise set of the combined HERA data these ideas implemented in several fits to diffractive DIS [7] and inclusive DIS [8,9] have lead to accumulating an evidence for higher-twist corrections in diffractive and inclusive DIS.
The recent studies of higher-twist effects in inclusive DIS [8,9] are based on DGLAP fits of the leading-twist contribution complemented by a simple model of twist 4 correction. In the central models elaborated in these analyses a multiplicative twist 4 correction was assumed for the longitudinal structure function F L of the form of AF L /Q 2 and the twist 4 correction to the structure function F 2 was set to vanish. 1 The fit quality increased significantly for the combined HERA data for Q 2 > Q 2 min = 2 GeV 2 for both NLO and NNLO DGLAP approximations of the leading-twist evolution. This simple model provides a surprisingly good description of the DIS data, except of the predicted steep rise of F L towards the low Q 2 for the NNLO DGLAP fit with the twist 4 correction found in Ref. [9]. Such a rise is not inline with the F L data [10][11][12].
In this paper we adopt a more flexible model of the twist 4 contribution motivated by an extraction of the twist 4 corrections to structure functions [13,14] from the Golec-Biernat-Wüsthoff saturation model [15]. The model is inspired by a resummation of multiple scattering in QCD in the eikonal approximation and it is capable to provide more information on the details of the higher-twist corrections and physics insight into their origin. In this approach, the twist 4 corrections to F 2 and F L structure functions have non-trivial Q 2 and x dependence. In addition, we modify the standard form of the DGLAP input for the gluon and sea density, so that they are consistent with general features of parton saturation in QCD at small x. With this model complemented by the NLO or NNLO DGLAP evolution of PDFs we analyse the combined HERA data on the reduced cross sections using the xFitter package [16] with suitable extensions of the code to incorporate the new features of the model. The sensitivity to the higher-twist corrections is enhanced by performing independent fits of the data sets with the momentum transfer constrained by Q 2 > Q 2 min and varying the limit Q 2 min . When Q 2 min is larger than 20 GeV 2 , the χ 2 /d.o.f. measures for the DGLAP fits at NLO and NNLO accuracy, with and without the higher-twist terms are close to 1. 15  Hence the improvement of the fit quality by adding the higher-twist corrections for Q 2 min = 1 GeV 2 is large and particularly pronounced for the fits using the NNLO DGLAP leading-twist part. Also a good description of the HERA F L data [10][11][12] is obtained down to the lowest measured values of Q 2 .
The evidence for sizeable contributions of higher-twist terms is further strengthened by an explicit analysis of the twist composition of the structure functions at small x and moderate and low scales. Consistently we find the growing higher-twist effects when x and Q 2 decrease. The relative importance of the higher-twist corrections is found to be larger in the NNLO fit than in the NLO one. In particular, in the NNLO fit at Q 2 = 1.2 GeV 2 , the twist 4 correction to the reduced cross section is found to be larger than the leading-twist contribution for x < 2 · 10 −4 , and the relative higher-twist correction further grows towards small x, and at the lowest available x 2 · 10 −5 it reaches about 200% of the leading-twist term. The higher-twist effects quickly decrease with increasing Q 2 and reach ∼ 10% level at Q 2 = 6.5 GeV 2 . The higher-twist effects are found to be much stronger in the longitudinal structure function F L . In particular, the twist 4 contribution to F L is larger than the leading-twist contribution for Q 2 < 5 GeV 2 (Q 2 < 6 GeV 2 ) for the NLO DGLAP fit (NNLO fit). Of course, also in F L the effects of higher-twist corrections decrease quickly with Q 2 , but in the NNLO fit, the higher-twist contribution is still visible at 10% level up to a sizeable scale Q 2 20 GeV 2 .
The inclusion of higher-twist corrections is found to affect significantly the fitted gluon and sea density functions at small x < 0.01 and moderate factorisation scales, µ 2 , while the sensitivity of the valence quark distribution to the higher-twist effects is minor and may be neglected. The largest difference in PDFs coming from the higher-twist effects is found in the gluon PDF-the difference is large and much larger than the corresponding uncertainty at small scales, 1 GeV 2 < µ 2 < 3.5 GeV 2 , then decreasing with µ 2 to a few percent level at µ 2 = 50 GeV 2 . The sea distribution at small x is also affected but it exhibits lower sensitivity to the presence of higher-twist corrections.
In conclusion, we have found a consistent evidence of the sizeable twist 4 corrections to proton structure functions. The evidence comes primarily from the χ 2 quality measure of fits to the combined HERA data on the inclusive DIS with the leading-twist component described by the NLO / NNLO DGLAP evolution. This evidence is further strengthened by the strong effects of the higher-twist corrections in the reduced cross section and the structure function F L at small x and moderate / low Q 2 . The fitted twist 4 contributions have the x-dependence that is consistent with the double exchange of hard gluonic ladder at small x, as expected from the QCD analysis of the evolution of leading quasi-partonic operators [2,3].

The model of higher-twist corrections
The Golec-Biernat-Wüsthoff (GBW) saturation model [15] offers a simple and effective description of DIS, DDIS structure functions down to very low Q 2 at small x, and also of the exclusive vector meson production [17]. In particular, with this model one is able to describe reasonably well even the transition from DIS at large Q 2 to the photoproduction limit. This transition may be viewed as a transition from the twist 2 regime to the region that all twist contributions are relevant. From the point of view of perturbative QCD, the GBW model corresponds to multiple independent high-energy scatterings of photon hadronic fluctuations, that is to the eikonal iteration of a single gluon ladder exchange. In particular, the leading behaviour of twist 4 contributions of the GBW cross section is ∼ (x λ Q 2 ) −2 (modulo logarithms) compared to the leading twist cross section ∼ (x λ Q 2 ) −1 (modulo logarithms) [13,14]. Such behaviour of twist 2 and twist 4 amplitudes is in a qualitative agreement with results of the evolution of twist 4 contributions in the Bukhvostov, Frolov, Lipatov and Kuraev framework [2,3,18], where dominant contributions at small x are driven by quasi-partonic operators.
From these studies [13,14] it follows that the twist 4 contributions to the transverse and longitudinal cross sections take the form where A, B are positive constants. The saturation scale depends on x variable as Q 2 sat = Q 2 0 (x 0 /x) λ where Q 2 0 = 1 GeV 2 and x 0 = 3.04 · 10 −4 (the GBW fit without charm) are model parameters. It is important to notice the opposite signs of the corrections in (1), positive for the transverse and negative for the longitudinal part. Additionally, the twist 4 contribution to the longitudinal cross section is logarithmically enhanced. This structure of the corrections can be also deduced from the general QCD analysis which was discussed in [13]. Following the above considerations we shape a saturation model inspired ansatz for twist 4 corrections in the following form in which coefficientsc L/T are left as free and independent parameters. In practical implementation of the fits we rewrite the above parameterisation to the following convenient form of twist 4 corrections to the structure functions: where F L/T = Q 2 σ L/T /4π 2 α em .

Equation (2) contains thec
(log) T parameter which determines the magnitude of the logarithmically enhanced term in the transverse cross section. Due to the properties of the transverse photon impact factor, this parameter vanishes at the leading order in the strong coupling constant expansion [13,14,19]. Therefore, one expects that this parameter is much smaller than parameterc  2 We confirmed that this parameter is indeed small in independent numerical fits.
When one attempts to extend the proton structure function analysis to the region of small Q 2 and very small x it is necessary to consider possible effects of parton rescattering and/or recombination in the dense parton regime. Those effects may lead not only to the higher-twist corrections but also they are expected to influence the form of the input for QCD evolution of the matrix elements in the OPE. In particular, it is natural to require that the input functions for the DGLAP evolution of parton densities are consistent with unitarity constraints relevant for high-energy scattering at very low x.
Hence, a precise analysis of the higher-twist effects which pronounce at low Q 2 and small x requires careful treatment of the gluon and sea input distributions at small x.
Currently, one of the most successful tools for an analysis of γ * -nucleon (or γ * -nucleus) scattering with the unitarity corrections is the Balitsky-Kovchegov (BK) equation [20,21]. In constructing the model for initial conditions of the PDFs we take into account the outcome of this equation analysis. The BK equation, that resums multiple scattering effects in the extended generalised leading logarithmic 1/x approximation [22,23] and the large N c limit, may be described in a natural way in terms of the colour dipole language [21,24,25]. In this approach the high-energy scattering is described in terms of the imaginary part of the BK dipole forward scattering amplitude, N (x, r), where r is the dipole extension vector in the transverse plane. On the other hand the same dynamics may be covered by the BK equation represented in terms of the unintegrated gluon density [26][27][28][29]. It can be shown that at the leading logarithmic 1/x accuracy the imaginary part of the BK dipole forward scattering amplitude, N (x, r), and the BK unintegrated gluon density 3 F(x, k 2 ) (where k is the gluon transverse momentum) are in one-to-one correspondence. Hence the unintegrated gluon density at small x can be recovered from N (x, r), whereÑ is the Fourier transform of the dipole scattering amplitude, and R p is an effective radius of the proton. After employing the leading logarithmic relation of the collinear gluon density f g (x, µ 2 ) and the unintegrated gluon density F(x, k 2 ) one obtains Explicit numerical solutions of the BK equation [29,30] show that at small x and for k below the saturation scale, Q sat (x), generated by the BK evolution, the solution of the BK equation tends to Exactly the same asymptotic behaviour of the unintegrated gluon density is found in the GBW model where one approximates the scattering amplitude with the saturation formula N (x, r) ∼ 1 − exp(−r 2 Q 2 sat (x)). In fact, the small k asymptotics of the saturated unintegrated gluon density may be traced back to the unitarity constraint on the dipole cross section in the position space, N (x, r) ≤ 1. For such a form of F(x, k 2 ) at small k, it is straightforward to show that in the limit µ 2 Q 2 sat (x) the gluon density small x asymptotics is with λ > 0. Therefore the gluon density at a low scale is expected to decrease toward zero with decreasing x. Following this argument our parameterisation of the input for gluon distribution fulfils the condition xf g (x, µ 2 0 ) ∼ x Bg , at x → 0, where B g is a positive fit parameter. Furthermore in our model construction we consider the input for the sea distribution at a small scale and small x. Assuming that the sea quarks at small x are generated predominantly from the gluon DGLAP splitting to quark/antiquark, one may approximate the sea singlet distribution by the LO DGLAP expression describing the feed-down from gluons, Obviously, since the scales probed are low and the impact of non-perturbative effects unknown, the above expression should be treated only as a QCD hint on the actual shape of the sea distribution at low scales. An explicit evaluation of the model expression in Eq. (7) leads to the sea-quark distribution asymptotic behaviour at small x following the gluon asymptotics, xf sea (x, µ 2 ) ∼ x λ . Hence in the fits that include parton saturation effects in the input distributions, we impose the same asymptotic behaviour of the sea and gluon distributions at the initial scale, xf sea (x,

DGLAP framework
The leading twist 2 contributions to the F 2 and F L structure functions are given in terms of PDFs, f k (x, µ 2 ), determined within the DGLAP framework. Our approach follows closely the scheme adopted in the HERAPDF2.0 [1] study, in order to clearly see the effects of higher-twist contributions. The PDFs are parameterised at the starting scale µ F0 and then determined at all scales µ F by solving the DGLAP evolution equations. The factorisation and renormalisation scales are chosen to be equal and in the following we denote them by µ, while the evolution starting scale is denoted by µ 0 .

Scheme description
The light quarks, u, d, s, are taken to be massless. The heavy quarks, c, b, t, are generated radiatively and appear only at transition scales, taken to be equal to the corresponding quark masses, m h . The PDFs of heavy quarks start from 0 once µ goes above m h . In other words, there are no intrinsic heavy flavours. For a simple realisation of this scenario we take the starting scale, µ 0 , below the charm mass.
The coefficient functions of heavy quarks are calculated in the Thorne-Roberts general-mass variableflavour-number scheme called RT OPT [31][32][33]. This scheme is adopted in accordance with the HERAPDF2.0 fit [1].

The input parameterisation
The distributions parameterised at the starting scale include the gluon g, d and u valence quarks, and up-and down-type sea quarks,Ū =ū,D =d +s.

The generic form of the input parton
The relatives contribution to the down-type sea at the starting scale is assumed to be a fixed (x-independent) fraction β ofD, i.e. fs = βfD.
Thus, in general, we have 26 fit parameters to start with. Several assumptions are made in order to make this number smaller.
First, the quark-counting and momentum sum rules are used to fix the valence and gluon PDFs normalisation parameters, A u val , A d val and A g . Next, the u val and d val distributions are assumed to have the same shape at small x i.e. B u val = B d val 4 . Also, to ensure a uniform sea behaviour at low x, (ū d ), the following constraints are imposed: AŪ = (1 − β)AD and BŪ = BD ≡ B sea . Finally, the strange sea fraction, β, is set to 0.4. Based on the HERAPDF2.0 experience and our numerous fit results we set to zero all D k and E k , except of E u val which is left free. With this setup we have 10 free parameters of the input PDFs to be compared with 14 free parameters of the HERAPDF2.0 fit. With this restricted parameterisation at the leading-twist level, we gain on the stability of the fits with the saturation and higher-twist effects included.
The first step towards improved description at low x and low µ 2 is a modification of the basic parameterisations (8), aimed at improving description of parton saturation effects in the gluon and sea input distributions in the small x domain.
In the HERAPDF2.0 fits the input gluon parameterisation is augmented by a negative term −A g x B g (1− x) 25 [1,34]. In the current study we do not include this subtraction. Instead, we consider enhancing the basic parameterisation (8) with saturation-inspired, damping factors for the parton k: applied to the gluon and sea components, with d k + B k > 0. The application of such factors ensures a smooth decrease to zero of xf g (x, µ 2 0 ) and xf sea (x, µ 2 0 ) when x → 0, consistent with the known results from analyses of parton saturation at small x (see Sec. 2 for a more detailed discussion). The damping factors describing the parton saturation effects turn on for x below a specific scalex k , which can be therefore interpreted as the saturation x at Q 0 : Q sat (x k ) = Q 0 . In general,x k and the saturation powers d g , d sea are arbitrary parameters, with the already mentioned constraint d k + B k > 0.
Hence we consider the following input parameterisations of the gluon and sea PDFs: where k = g,Ū ,D. With these parameterisations in the x → 0 limit the input PDFs scale as After a preliminary analysis of the data we found that for the gluon input distribution the saturation damping factor is irrelevant, as the fits yielded B g ∼ 1 which already guarantees power-like approach to zero of xf g (x, Q 2 0 ) for x → 0. Thus we retain the damping factor for the quark sea only, witĥ x D =x U ≡x being a free fit parameter. For the saturation damping exponent for the quark sea, d sea , we impose an additional constraint following from the assumption that the sea input distribution at small x follows the power-like behaviour of the gluon input distribution, see Sec. 2. As a result, a relation of the exponents is obtained: B sea + d sea = B g , resulting in d sea = B g − B sea . In fact, we have checked that leaving d sea as a free parameter does not improve the fit quality (the difference in the χ 2 /d.o.f. is smaller than 0.002). Hence the phenomenological inclusion of the saturation effects in the input of the PDFs is reduced to taking the positive definite gluon input and imposing the sea input damping at small x.

Results
In the current analysis we use the combined HERA data on neutral and charged current e + p and e − p inclusive cross sections, measured at centre-of-mass energies ranging from 225 GeV to 318 GeV [1].
In the fits we use only data points for which Q 2 > 1 GeV 2 . Their kinematic range spans four orders of magnitude in x and Q 2 with lower bounds at x = 1.76 · 10 −5 and Q 2 = 1.2 GeV 2 . The inelasticity y values are between 0.001 and 0.95. The whole data set comprises 1213 data points. A subset of this data set with Q 2 ≥ 3.5 GeV 2 was used to extract the HERAPDF2.0 PDFs [1].
The measured cross sections are presented 5 in terms of the reduced cross section: where In order to fit the data down to Q 2 = 1.2 GeV 2 , we take the starting scale for the DGLAP evolution µ 2 0 = 1 GeV 2 . For the sake of comparison to the HERAPDF2.0 fits we present also some results for the PDFs parameterised at µ 2 0 = 1.9 GeV 2 . The fits are performed using the xFitter package [16] supplemented by us with necessary code extensions including input parameterisation with saturation damping effects, and the higher-twist contributions to F T and F L , as given by Eq. 3. The DGLAP evolution is performed using the QCDNUM program [35].
In the analysis below we use the following fit names: • HTS -Higher Twist + Saturation, µ 2 0 = 1.0 GeV 2 , Q 2 min = 1.2 GeV 2 -Our main fit. It corresponds to the standard DGLAP evolution of the leading twist terms, in which the saturation damping effects are assumed in the input PDFs, complemented by the additive twist 4 corrections. Note that when performing the χ 2 scan (see Sec. 4.1) we use the same name for this model fitted to data with a variable lower cutoff Q 2 min on Q 2 ; • HT -Higher Twist without Saturation, µ 2 0 = 1.0 GeV 2 , Q 2 min = 1.2 GeV 2 -Like HTS but without the saturation damping effects in the initial PDFs; • LT-STD -Standard Leading Twist without Saturation, µ 2 0 = 1.9 GeV 2 , Q 2 min = 3.5 GeV 2 -HERAPDF2.0-like fit -pure DGLAP approach, no saturation damping in the input parameterisation, the initial scale and the data selection as in HERAPDF2.0; • LT-1.0 -Leading Twist, without Saturation, µ 2 0 = 1.0 GeV 2 -Like LT-STD, but with lower initial scale and a variable lower cutoff Q 2 min on Q 2 ; • LT-1.9 -Leading Twist, without Saturation, µ 2 0 = 1.9 GeV 2 -Like LT-STD, but with a variable lower cutoff Q 2 min on Q 2 .
The analysis of the data is carried out in the following stages. First, we perform the scan of the χ 2 /d.o.f. for the data set with Q 2 > Q 2 min as a function of Q 2 min for all chosen setups, that is for NLO and NNLO DGLAP evolution of the leading twist, with and without the higher-twist corrections, with and without saturation damping effects in the input for the PDFs. We have checked that the saturation damping modification of the input distributions becomes important only upon the inclusion NLO LT-1.0, μ 2 0 = 1.0 GeV 2 LT-1.9, μ 2 0 = 1.9 GeV 2 HT, μ 2 0 = 1.0 GeV 2 HTS, μ 2 0 = 1.0 GeV 2  of the twist 4 contributions. Without the higher-twist corrections the damping gives practically no improvement of the fits, so we do not present results for the DGLAP fits without the higher-twist corrections with the saturation damping of the input. Next we describe the features of our best fits with the higher-twist corrections. Furthermore we explicitly study the effects of the higher twists for σ red and F L , as given by the best fits. Finally we show the impact of the inclusion of the higher-twist and saturation effects on the obtained PDFs.

The χ 2 scans
In Fig. 1 we show the χ 2 /d.o.f. for the fits to the reduced cross sections as a function of the lower cutoff Q 2 min imposed on the photon virtuality Q 2 for the data sample taken into account in the fits. The initial scale of the DGLAP evolution is set to µ 0 . The higher-twist parameterisation provides the best description of the data from Q 2 min 16 GeV 2 below. The question of the parton saturation effects in the input sea distribution is more subtle. From the χ 2 scan it follows that within the NLO approximation this is not an important effect. However, in the fits assuming the NNLO DGLAP evolution of the leading twist contribution, both the higher-twist corrections and the parton saturation effects in the input are key ingredients for the best description of the data. The inclusion of both effects improves the data description significantly for Q 2 min < 20 GeV 2 . In Fig. 2 we show a comparison of the χ 2 /d.o.f. between the data and two fits, standard fit LT-STD and the best fit with higher-twist corrections HTS, as a function of the lower cutoff Q 2 cut imposed on the data. Here Q 2 min is kept fixed to 2.0 GeV 2 for the LT-STD fit and to 1.2 GeV 2 for the HTS fit. A systematic improvement of the data description is clearly seen for the HTS fit with respect to the LT-STD fit for Q 2 cut < 5 GeV 2 at NLO and for Q 2 cut < 20 GeV 2 at NNLO. A short comparison of our reference fit LT-STD to the HERAPDF2.0 fit is in order. The resulting values of the χ 2 /d.o.f. for the LT-STD fit are 1.213 and 1.228 at NLO and NNLO, correspondingly, while for the HERAPDF2.0 fit they are 1.200 and 1.205 at NLO and NNLO, correspondingly. Recall however, that the input PDFs parameterisation of the LT-STD fit has 10 free parameters, to be compared with 14 parameters of the HERAPDF2.0 fit, and thus the slight difference in the χ 2 /d.o.f. is acceptable.

Features of the best fits
The parameters of the model obtained from the best fits (HTS) of all the data with Q 2 > 1 GeV 2 with the twist 4 corrections included and the sea input with saturation damping, are given in Tab. 1. The results are displayed both for NLO and NNLO DGLAP evolution of the leading twist terms.
It is interesting to analyse the obtained parameters describing the higher-twist corrections and compare the results to the expectations from the GBW model. First of all, the value obtained from the fit of the saturation scale exponent is λ = 0.351 ± 0.008 for the NLO HTS fit and λ = 0.257 ± 0.016 in the NNLO HTS fit. These values are rather close to the saturation exponents of the GBW model λ = 0.288 (without charm) and λ = 0.277 (with charm). The obtained values of λ ∼ 0.26 − 0.35 are also consistent with the picture of double hard pomeron exchange as the leading contribution to twist 4 corrections at small x. Also the value of the saturation x parameter at Q 0 = 1 GeV obtained from the NNLO fit,x = (2.0 ± 0.4) · 10 −4 compares well to the corresponding GBW saturation x parameters, x 0 = 3.04 · 10 −4 (without charm) and x 0 = 0.41 · 10 −4 with charm. In the NLO fit, however, the obtainedx = 0.09 ± 0.2 · 10 −5 is consistent with zero. Recall that in our approach parameterx is the characteristic x for emergence of the saturation damping effects in the sea distribution. Hence the conclusion implied by the χ 2 scan is confirmed: that the sea saturation damping is important for the NNLO DGLAP fit with twist 4 corrections, while the DGLAP NLO fit with higher twists does not require the saturation input damping in the x range of the fitted data.
Interestingly, the pattern of the twist 4 multiplicative coefficients is found to differ significantly from the predictions of the GBW model. At small x the model yields a sizeable negative twist 4 correction to F L and a positive correction to F T . The performed fits exhibit a different pattern -both at NLO and NNLO we find a small positive twist 4 correction to F T and a larger positive correction to F L . The difference of the sign of the higher twist correction to F L at small x between the GBW model prediction and the fit results, occurs both in the leading logarithmic term ∝ c (log) L log(Q 2 /Q 2 sat (x)) for Q > Q sat (x) and in the constant term c

Comparison with the data for σ red and F L
In the approach presented here the relative importance of the higher-twist corrections to the proton structure functions may be estimated for different x and Q 2 . Such an estimate provides a measure of both the expected accuracy of the leading twist description and the sensitivity to the higher-twist contribution. It also permits to determine the kinematic region in which the higher-twist corrections are most important and, in this way, the evidence of the higher-twist contribution to the structure functions is strengthened. Indeed, we find that the most important effect in the structure function deviations from the DGLAP leading twist description comes from the region of small x and Q 2 , where the higher-twist effects are strongest.
In Fig. 3 the data for the reduced cross sections for 1.2 GeV 2 < Q 2 < 15 GeV 2 are compared to the HTS fit results. In order to illustrate the twist content of the proton structure function we also show the twist decomposition of the reduced cross section at the NLO level. The twist 4 contribution makes up to 75% of the twist 2 contribution at Q 2 = 1.2 GeV 2 and x = 3 · 10 −5 , where the higher-twist effects are largest. In this kinematic region the twist 4 effects are estimated to provide about 40% of the reduced cross section. As expected, the higher-twist contribution is suppressed with increasing x, and at x = 3 · 10 −4 the twist 4 correction is reduced to about 20% of the total value. The relative importance of the higher-twist correction decreases also with growing Q 2 . Indeed, at x = 3 · 10 −5 the relative twist 4 contribution is about 30% at Q 2 = 2 GeV 2 and about 20% at Q 2 = 3.5 GeV 2 . In the HTS fit at NLO, the relative higher-twist effect is below 10% at Q 2 = 6.5 GeV 2 .
The results of a similar investigation at NNLO are displayed in Fig. 4 for Q 2 up to 15 GeV 2 . In the NNLO DGLAP fits the higher-twist effects are found to be significantly stronger than in the NLO fits over the whole probed range of Q 2 . In particular, in the NNLO fit at Q 2 = 1.2 GeV 2 , the twist 4 correction is found to be larger than the leading-twist contribution for x < 2 · 10 −4 , and the relative correction further grows towards small x, to reach about 200% of the leading-twist term at the lowest available x 2 · 10 −5 . At Q 2 = 2 GeV 2 and the smallest x, the twist 4 correction reaches about 80% of the twist 2 contribution, and at Q 2 = 3.5 GeV 2 the higher-twist term is still around 25% of the leading-twist term. Finally, the higher-twist correction reaches ∼ 10% level at Q 2 = 6.5 GeV 2 and quickly decreases for larger Q 2 .
The characteristic behaviour of the data at moderate Q 2 is a turn-over at small x. This feature is not reproduced by the DGLAP fits without higher-twist corrections [1] and the inclusion of higher-twist effects is necessary to provide a good description of this behaviour [8,9]. Hence the turn-over may be considered to be a signature of the higher-twist contributions. The HTS fits reproduce well this shape both at NLO and NNLO.
In the existing analyses of the combined HERA data with higher-twist corrections [9] a satisfactory description of the F L data at smaller Q 2 has not been achieved within the NNLO framework [9]. The predictions for F L obtained in our approach are shown in Fig. 5 in comparison to the experimental data from H1 [11]. The F L data are well described down to Q 2 = 1 GeV 2 . Note that the plotted F L data were not directly fitted, the F L contribution was treated in the fits only as a part of the reduced cross section σ red . The higher-twist contributions are found to be important in F L at small and moderate Q 2 . In particular, the twist 4 term dominates for Q 2 < 5 GeV 2 (Q 2 < 6 GeV 2 ) for the NLO fit (NNLO fit). Remarkably, in the NNLO fit, the higher-twist contribution is still visible at 10% level up to a sizeable scale of Q 2 20 GeV 2 . This shows that the longitudinal structure function is particularly sensitive to the higher-twist effects and it may be used as their effective probe. F L data HTS twist-2 twist-4 Figure 5: The predictions for F L from the HTS fit compared to the H1 data [11]. Also shown are twist 2 (dashed green line) and twist 4 (dotted purple line) contributions.

Impact of the higher-twist effects on the PDFs
One of the key goals of the present analysis is to understand the impact of the higher-twist corrections on the parton density functions. It is expected that inclusion of the higher-twist corrections alters the resulting PDFs, hence affecting the predictions in which the PDFs are used. Below we present the PDFs obtained from the HTS fits (including the higher-twist effects) compared to the PDFs obtained within the standard framework, which we denote as the LT-STD fit. The determined PDFs are presented with the corresponding experimental uncertainties. The included higher-twist corrections are larger in the small x region, so they affect mostly the gluon and sea PDFs at small x. Hence, in the figures we show only these parton distributions. We have checked that the effect of the higher-twist corrections on the valence quark distribution is much smaller than for the gluon and sea PDFs.
In Fig. 6 we compare the gluon and sea PDFs at small scales, µ 2 = 1.2, 2.0 and 3.5 GeV 2 , with (HTS) and without (LT-STD) the higher-twist corrections, at the NLO and NNLO accuracy. The input condition for the HTS fit is assumed at µ 2 0 = 1 GeV 2 and it includes the saturation damping effect in the gluon and sea PDFs. The standard fit (LT-STD) starting scale is µ 2 0 = 1.9 GeV 2 , so for µ 2 = 1.2 GeV 2 the LT-STD PDFs should be treated as extrapolations only. The dotted vertical lines in the plots mark the minimal kinematically allowed value of x in HERA measurements. Hence the PDFs values to the left of this line are also extrapolations from the region of available data towards smaller x.
It is clear from Fig. 6 that inclusion of the higher-twist and saturation effect leads to sizeable changes of the gluon and sea PDFs at small x < 0.01 and small scales µ 2 < 3.5 GeV 2 both at the NLO and NNLO accuracy. The effects found are larger for the NLO fits, where e.g. for µ 2 = 3.5 GeV 2 and x = 3 · 10 −5 the HTS gluon is larger than the LT-STD gluon by about 60% and the HTS sea PDF is smaller by more than 10% than the LT-STD one. In the NNLO fits the differences between the gluon PDFs with and without the higher-twist effects tend to be smaller than in NLO fits, while the effect is larger in the sea quarks at NNLO. In general, at the smaller scales, the inclusion of the higher-twist effects leads to a larger gluon PDF and the reduced sea PDF. The differences found are not only in the shape and values of the PDFs but also in the relative behaviour of the gluon and sea distributions.  In the standard approach those two are decoupled at the input scale whereas in the HTS approach the sea distribution at small x follows the gluon distribution already at the input scale. In a similar way, Fig. 7 displays the impact of higher-twist effects on the PDFs at larger scales, µ 2 = 20 and 50 GeV 2 . The higher-twist effects are still significant in the gluon distribution at small x. As for the smaller scales, the inclusion of the higher-twist effects leads to a larger gluon PDF. The change of the gluon PDF at the small x values corresponding to the kinematic lower limit at HERA is about 10% for µ 2 = 20 GeV 2 , and slightly below 10% for µ 2 = 50 GeV 2 . At the larger scales the sea distribution is not significantly affected by the higher-twist corrections. Note finally that the higher-twist corrections lead to larger changes in the PDFs for the NLO DGLAP framework than for the NNLO DGLAP one. Fig. 8 shows the independent PDFs at µ 2 = 10 GeV 2 obtained from our best fits with higher-twist corrections and the parton saturation effects. The valence up and down quark, the sea and the gluon PDFs are shown with their experimental uncertainties at the NLO and NNLO accuracy.

Discussion
Let us recall the main results of this paper that substantiate the evidence of the significant highertwist contributions in the inclusive DIS at HERA at small x and Q 2 . The strongest point is the comparison of the fits to the combined HERA data on the reduced cross sections based on the DGLAP evolution with and without twist 4 corrections. For the data sample analysed, with Q 2 > 1 GeV 2 , the inclusion of the twist 4 terms improves the χ 2 /d.o.f. of the NLO fit from about 1.35 (for the pure DGLAP) to about 1.22. With about 1200 degrees of freedom of the fit, the statistical significance of the χ 2 change corresponds to an improvement of the p-value by more than seven orders of magnitude. For the NNLO DGLAP fit we find an improvement of the χ 2 /d.o.f. from about 1.5 to less than 1.2, so the statistical significance of the improvement is greater by many orders of magnitude that in the NLO DGLAP fits. On the other hand, the values of the χ 2 /d.o.f. about 1.2 found in the DGLAP fits with twist 4 corrections are still uncomfortably larger than 1.0. Here, however, being aware of difficulties in treatment of the systematic errors of the measurement, we take as the reference level the value χ 2 /d.o.f. = 1.15 that is obtained in all the fits for Q 2 min > 10 GeV 2 , where the higher-twist corrections should be negligible. Given the crudeness of the higher-twist model applied, we find the change of the χ 2 /d.o.f. from the reference level 1.15 to less than 1.2 in the best NNLO fit with higher twists to be sufficiently small to accept the model.
An issue that should be raised is a possible alternative explanation of the data by a pure DGLAP fit with a more flexible input parameterisation. Taking into account this possibility we studied in detail the effect of variation of the sea and gluon input PDFs at small x due to the inclusion of the fitted saturation damping effects. We found that these variations of the input PDFs do not change the conclusion on the presence of the higher-twist contributions. Moreover, the contribution of the higher twists found is characterised by a strong Q 2 dependence which cannot be mimicked by the features of the input parameterisations that are given at a fixed scale µ 2 0 . The explicit analysis of the Q 2 and x-dependencies of the higher-twist terms leading to the improvement of the fit quality shows full consistency with the assumption of the leading twist 4 exchange as the contribution necessary to well describe the data. In particular the obtained twist 4 correction is characterised by a leading x −2λ behaviour with λ = 0.26 for the NLO fit and λ = 0.35 for the NNLO fit-the values strongly supporting the interpretation of the twist 4 terms as coming from the double hard gluon ladder exchange, as expected from the QCD analyses of the twist 4 evolution at small x.
The higher twist interpretation advocated in this paper should be also confronted with recent interesting results obtained for the vector meson production [36], where power corrections to the coefficient functions in high energy scattering were derived, that do not come from the leading quasipartonic 4-gluon twist 4 exchange 6 . Although we agree that such corrections should also contribute to the DIS cross sections, we stress that the x-dependence of the twist 4 corrections (i.e. the power corrections in 1/Q 2 ) found in the present analysis is consistent with the dominance of the exchange of two hard gluonic ladders.
The results presented here are compatible with the two recent papers [8,9] in which a similar conclusion was reached that the combined HERA data require an inclusion of the twist 4 corrections at small x and Q 2 . These two analyses have been performed using a similar method of the χ 2 scan as a function of Q 2 min for DGLAP fits with and without the higher-twist corrections. In both approaches the data were selected with Q 2 > 2 GeV 2 and no sensitivity was found to the higher-twist corrections in F 2 . Hence, in the central model of both the approaches, only F L receives the twist 4 correction in the form of ∆ (τ =4) F L (x, Q 2 ) = AF L (x, Q 2 )/Q 2 . The improvements and extensions of the present analysis are the following. The model of higher-twist corrections we propose here is based on general features of multiple scattering in QCD at high energies. In particular it includes the known results on the leading evolution of twist 4 quasi-partonic operators [2,3]. Therefore the model offers an explicit physics interpretation of the result and allows to obtain a deeper insight into the x-dependence of higher-twist corrections in both F L and F T (or F 2 ). Furthermore, in order to enhance the sensitivity to the highertwist effects, we extended the Q 2 range of the fitted data down to Q 2 = 1 GeV 2 . Indeed, with this extension the fits are sensitive to twist 4 effects in both F L and F 2 and the x-dependence of the twist 4 corrections is well constrained. Hence, the additional conclusion we find w.r.t. the previous studies is the steep x −2λ dependence of the twist 4 corrections, while the twist 4 correction ∆ (τ =4) F L (x, Q 2 ) of Refs. [8,9] is characterised by a single gluonic ladder exchange ∼ x −λ . It is intriguing why the fits do not distinguish clearly between these two different scenarios. We propose the following interpretation of this apparent discrepancy. As follows from the performed twist decomposition of the reduced cross section in the best fits, the higher-twist corrections become significant in the bins of the lowest Q 2 , below Q 2 = 5 GeV 2 and at the smallest available x. Therefore there is only a moderate region of the kinematic space (and a moderate number of the data points in this space) to pin down the x-dependence of the higher-twist corrections with a high precision, see a discussion of this issue performed in Ref. [8]. In Refs. [8,9] the data points were taken with Q 2 > 2 GeV 2 so their sensitivity to the x-dependence of the higher-twist corrections is smaller than in the present analysis. Furthermore in the present approach we introduce additional freedom of the sea input parameterisation at small x motivated by parton saturation damping effects. This feature leads to a reduced sea PDF at small x and small scales w.r.t. the standard sea PDF and a modified physics picture emerges. In [8,9] the sizeable higher-twist corrections are found in F L and they are negligible in F 2 . In the HTS NNLO fit obtained in the present analysis we find the sizeable positive higher-twist corrections in both F L and F 2 at small x and Q 2 , and certain suppression of the leading twist contribution to F 2 due to the saturation damping of the sea. Finally the model of higher-twist correction of Ref. [9] is simple and efficient, however it leads to unphysical behaviour of the predicted F L at small Q 2 at NNLO. In summary, our results are consistent with the main findings of Refs. [8,9] but we offer a more flexible approach to the higher-twist corrections with a clear physical interpretation. The proposed framework can be successfully extended down to Q 2 = 1 GeV 2 and gives a good description of F L over the whole Q 2 range.
A comment is in order on the fact that the sign of the obtained twist 4 correction to F L is positive and it disagrees with the predictions [13,14] from the Golec-Biernat-Wüsthoff saturation model [15]. The result is also not intuitive in a simplified picture of the multiple scattering, in which the longitudinal virtual photon-proton scattering total cross section receives a positive correction from double scattering, while one could expect a negative shadowing correction. The sign of the twist 4 correction depends however, on the sign of the corresponding coefficient function. In the GBW model the twist 4 coefficient function is implicitly computed assuming eikonal coupling of two gluonic ladders to the virtual quark loop coming from the virtual photon fluctuation. Such a coupling is not what is predicted by QCD at small x. The γ * coupling to four-gluon exchange was analysed in the leading logarithmic 1/x approximation in the framework of extended generalised leading logarithmic approximation [19,22,23,37] and it was found that the coupling occurs via triple pomeron vertex, which is distinctly different from the eikonal coupling. A more detailed analysis of the four-gluon coupling to the virtual photon was performed in Ref. [38]. The definite prediction of the sign of the twist 4 corrections to the proton structure functions in this approach is not available yet, so it is interesting to obtain it and confront with the results of the fit.
Finally, let us discuss shortly the concept of the gluon and sea saturation damping effects in the input PDFs at small x. The motivation to propose this damping for the gluon PDF comes from the properties of the solution of the Balitsky-Kovchegov equation [20,21] and for the sea PDFs we additionally assume that the input sea PDF at small x is driven by the gluon PDF. The damping may be understood as a result of unitarity corrections due to the multiple scattering effects below the scale of the input parameterisation. For the gluon PDF, the DGLAP fits (with or without higher-twist corrections) with the conventional parameterisation of the input gluon PDF, with the small x asymptotics xf g (x, µ 2 0 ) ∼ x Bg , yield a positive B g when µ 2 0 is sufficiently small, inline with the saturation damping property. So, for the gluon PDF the saturation damping is implicit in the standard parameterisation. For the sea PDFs, however, the saturation damping leads to a significant modification of the input shape at µ 2 0 = 1 GeV 2 . The saturation damping effect of the sea distribution is found to be supported by the data in the NNLO DGLAP fit with the higher-twist corrections. The inclusion of the saturation damping effects in this fit improves the χ 2 /d.o.f. significantly for moderate Q 2 . In particular for the complete data sample the χ 2 /d.o.f. goes down from 1.25 to 1.2. In this fit, the characteristic value of x for which the saturation damping effects turn on at µ 2 0 = 1 GeV 2 was found to bex 2 · 10 −4 , in consistence with the earlier analyses of the parton saturation at small x [15,17].