Parton distributions with small-x resummation: evidence for BFKL dynamics in HERA data

We present a determination of the parton distribution functions of the proton in which NLO and NNLO fixed-order calculations are supplemented by NLLx small-x resummation. Deep-inelastic structure functions are computed consistently at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NLO+NLL}x$$\end{document}NLO+NLLx or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {NNLO+NLL}x$$\end{document}NNLO+NLLx, while for hadronic processes small-x resummation is included only in the PDF evolution, with kinematic cuts introduced to ensure the fitted data lie in a region where the fixed-order calculation of the hard cross-sections is reliable. In all other respects, the fits use the same methodology and are based on the same global dataset as the recent NNPDF3.1 analysis. We demonstrate that the inclusion of small-x resummation leads to a quantitative improvement in the perturbative description of the HERA inclusive and charm-production reduced cross-sections in the small x region. The impact of the resummation in our fits is greater at NNLO than at NLO, because fixed-order calculations have a perturbative instability at small x due to large logarithms that can be cured by resummation. We explore the phenomenological implications of PDF sets with small-x resummation for the longitudinal structure function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_L$$\end{document}FL at HERA, for parton luminosities and LHC benchmark cross-sections, for ultra-high-energy neutrino–nucleus cross-sections, and for future high-energy lepton–proton colliders such as the LHeC.


Introduction
The experiments at CERN's large hadron collider (LHC) continue to explore particle physics both at the high-energy and the high-precision frontiers. The outstanding quality of current and forthcoming LHC data challenges the theory community to perform more precise calculations, so that meaningful conclusions can be drawn when comparing these theoretical predictions to experimental measurements. In this respect, the tremendous effort to be made in order to arrive at precision calculations for hard-scattering matrix elements and final-state parton evolution has to be accompanied by a comparable level of understanding of the internal structure of the initial-state hadrons.
Global analyses of PDFs [1][2][3][4][5][6] (see [7][8][9][10][11] for recent overviews) are generally based on fixed-order perturbative calculations, at LO, NLO and NNLO. However, it is well known that further logarithmic enhancements can affect partonic cross-sections and DGLAP evolution kernels order by order in perturbation theory. If we denote by Q the hard scale of the process of interest and by √ s the center-of-mass energy of the colliding protons, we have logarithmic enhancements in two opposite limits, namely Q 2 ∼ s (the threshold region) and Q 2 s (the high-energy region). Introducing the variable x = Q 2 /s, the threshold limit corresponds to large x, while the high-energy limit to small x.
The LHC is exploring a vast kinematic range in x, potentially covering both extreme regions. It is therefore crucially important to consistently assess the role of logarithmic corrections both at large and small x. For instance, searches for new resonances at high mass are sensitive to PDFs in the region between 0.1 x 0.7 [12]. On the other hand, processes such as forward production of Drell-Yan lepton pairs at small di-lepton invariant masses [13] and of D mesons at small p D T [14], both measured by the LHCb collaboration, probe values of x at the other end of the spectrum, down to x ∼ 10 −6 .
Calculations that aim to describe these extreme regions of phase space should in principle include resummation in the calculations of matrix elements and should make use of PDFs that were determined with a consistent theory. Threshold (large-x) resummation has already been included in PDF fits [15] (see also Ref. [16]) and dedicated studies which include threshold resummation in both the coefficient functions and in the PDFs have been performed in the context of heavy supersymmetric particle production [12]. The inclusion of threshold resummation in PDF fits is straightforward because in the widely used MS scheme the DGLAP evolution kernels are not enhanced at large x [17,18], so threshold resummation is only necessary for the coefficient functions, and can thus be included rather easily.
The situation is rather more intricate for small-x resummation, because both coefficient functions and splitting func-tions receive single-logarithmic contributions to all orders in perturbation theory. Small-x resummation is based on the BFKL equation [19][20][21][22][23]. However, the naive application of the fixed coupling leading-log x (LLx) BFKL equation to small-x deep-inelastic scattering (DIS) structure functions predicted a much steeper growth than that actually observed by the first HERA measurements [24,25], which instead were well reproduced by the predictions of LO and NLO running coupling DGLAP [26][27][28][29][30][31]. This paradox was compounded by the computation of next-to-leading logarithmic (NLLx) corrections to the evolution kernels [32][33][34][35][36], which turned out to be large and negative, destabilizing the LLx BFKL result. The correct implementation of small-x resummation turns out to require the simultaneous resummation of collinear and anti-collinear singularities in the small-x evolution kernels, together with a consistent resummation of running coupling effects.
On the other hand, while fixed-order DGLAP theory can provide a reasonable fit to the inclusive HERA data, several groups have found indications that the description of the most precise legacy datasets is not optimal in the small-x and small-Q 2 region, especially at NNLO 1 [64][65][66][67][68][69][70]. Currently, the evidence that this tension is related to lack of small-x resummation is inconclusive. The only way to show that it is due to resummation would be to perform a complete global PDF analysis including small-x resummation. Since the effect of resummation is known to be small, at least in the kinematic region explored at HERA, it is necessary that these fits are free of methodological bias. The NNPDF framework [71][72][73][74][75][76][77][78][79], having been validated by a closure test, is thus ideal in this respect.
With these motivations, the goal of this paper is to present a state-of-the-art PDF determination in which NLO and NNLO fixed-order perturbation theory is matched to NLLx small-x resummation. This will be done by supplementing the recent NNPDF3.1 PDF determination [79] with small-x resummation of DGLAP evolution and DIS coefficient functions using HELL, thereby leading to resummed PDF sets. We will show that the inclusion of small-x resummation significantly improves the quantitative description of the small-x and small-Q 2 HERA data, in particular at NNLO, both for the inclusive and for the charm structure functions. Our results fulfill a program that was initiated more than 20 years ago, when the first measurements of F 2 (x, Q 2 ) at HERA stimulated studies on the inclusion of small-x resummation in perturbative evolution [80][81][82][83].
The outline of this paper is as follows. First, in Sect. 2 we review the implementation of small-x resummation that we will use, and illustrate how resummation affects PDF evolution and DIS structure functions. Then in Sect. 3 we present the settings of our fits, which we dub NNPDF3.1sx, and in particular we discuss the choice of kinematic cuts. The results of the fits with small-x resummation are discussed in Sect. 4. In Sect. 5 we show the comparisons with the HERA experimental data, and provide detailed evidence for the onset of resummation effects in the inclusive and charm-production structure functions. We then perform a first exploration of the phenomenological implications of the NNPDF3.1sx fits at the LHC and beyond in Sect. 6, and finally in Sect. 7 we summarize and outline possible future developments.

Implementation of small-x resummation
Here we briefly review the implementation of small-x resummation which will be adopted in the sequel. First, we summarize the general features of small-x resummation theory, its main ingredients, and available approaches to it. We then discuss separately the implementation and general phenomenology of small-x resummation of perturbative evolution, and of deep-inelastic structure functions.

Basics of small-x resummation
In collinear factorization, the deep-inelastic scattering structure functions can be expressed as where x = Q 2 /s, μ R and μ F are the renormalization and factorization scales, the sum runs over partons, and we have factored out for convenience the Born cross-section σ 0 . Similarly for hadronic processes where M 2 is the invariant mass of the particles produced in the final state, x = M 2 /s, and the parton luminosities 3) The scale dependence of the PDFs f i x, μ 2 is controlled by the DGLAP evolution equations (2.4) and knowledge of the splitting kernels P i j (x, α s ) to (k + 1)loops allows for the resummation of collinear logarithms at N k LO accuracy. The evolution kernels are currently known to NNLO (3 loops) [84,85], and partially even to N 3 LO (4 loops) [86,87]. Single logarithms of x affect higher order corrections to both splitting functions and hard cross-sections. Specifically, the generic all-order behavior of the gluon-gluon splitting function is P gg ∼ 1 x n α n s ln n−1 1 x . Small-x logarithms are mostly relevant for PDFs in the singlet sector, i.e. the gluon and the quark singlet: small-x (double) logarithms in nonsinglet PDFs are suppressed by an extra power of x. Partonic cross-sections (either inclusive, or differential in rapidity or transverse momentum) can also contain small-x logarithms, which depend on the process and the observable. For gluoninduced processes (such as Higgs or top production) resummation affects the leading-order cross-section and it is thus a leading-log x (LLx) effect, while for quark-induced processes (such as Drell-Yan or deep-inelastic scattering) there must be a gluon-to-quark conversion, which makes it a NLLx effect. In either event at small x and low scales the combination α s ln 1 x can become large, spoiling fixed-order perturbation theory. In these circumstances it becomes necessary to resum the large logarithms in both splitting and coefficient functions in order to obtain reliable predictions.
In the ABF approach, one constructs perturbatively stable resummed results by combining three main ingredients: duality, i.e. consistency relations between the DGLAP and BFKL evolution kernels [37,38,97,98], which are used to construct a double-leading evolution kernel that simultaneously resums both collinear and small-x logarithms; symmetrization of the BFKL kernel in order to stabilize its perturbative expansion both in the collinear and anti-collinear regions of phase space [43,47], and thus in the region of asymptotically small x; and resummation of running coupling contributions, which despite being formally subleading are in fact dominant asymptotically, since they change the nature of the small-x singularity [41,42,52,53,58,99]. The resummation of gluon evolution with all the above ingredients consistently combined was originally achieved to NLO+NLLx in Refs. [43,53], while the inclusion of the quark contributions and the rotation to the physical basis of the singlet sector was completed in Refs. [46,57]. The matching to NNLO has been recently achieved in [63] and represents an important new development since it makes it possible to compare NNLO results with and without NLLx small-x resummation included.
Thanks to high-energy factorization [100][101][102][103] (generalized in Ref. [104] to rapidity and in Refs. [105,106] to transverse momentum distributions) it is possible to also perform resummation of the leading small-x logarithms in the coefficient functions both in deep-inelastic cross-sections Eq. (2.1) and hadronic cross-sections Eq. (2.2). The resummation relies on the resummation of the splitting function, which must then be combined with a computation of the hard cross-section with incoming off-shell gluons. Such calculations have been made for a range of processes: heavyquark production [100,101,107,108], DIS structure functions [103,109,110], Drell-Yan production [111,112], direct photon production [113,114] and Higgs production [115][116][117]. The use of these expressions to resum coefficient functions at fixed coupling is straightforward, but becomes more complicated when the coupling runs, due to the presence of anti-collinear singularities. This issue was resolved (both for photoproduction and hadroproduction processes) in Ref. [44], and used in Ref. [46] to compute running coupling coefficient functions for DIS.
In order to discuss NLLx resummation, we have to carefully specify the choice of factorization scheme. The socalled Q 0 MS scheme is often introduced [55,88,102,103], and is preferred to the traditional MS because it gives more stable resummed results. When expanded to fixed order, the scheme-change factor between the two is O(α 3 s ), so NLLx resummation in Q 0 MS can be matched directly to the usual fixed-order NNLO MS scheme calculation.

Resummation of DGLAP evolution
Resummed splitting functions take the generic form where the first contribution is the splitting function computed to fixed-order k (so k = 0, 1, 2 for LO, NLO and NNLO) and the second term is the resummed contribution, computed to either LLx (h = 0) or NLLx (h = 1), minus its expansion to the fixed-order k to avoid double counting. We note that the splitting functions in the gluon sector (P gg and P gq ) contain LLx and NLLx contributions, while in the quark sector (P qg and P qq ) they only start at NLLx. For this reason, there have been attempts to partially extend the resummation to the next logarithmic order (see [118]) which, however, are not considered in this work.
In Fig. 1 we show a comparison of the fixed-order gluongluon x P gg (x, α s ) (left) and the quark-gluon x P qg (x, α s ) (right plot) splitting functions with the resummed counterparts. The comparison is performed in the Q 0 MS factorization scheme, with n f = 4 active quark flavors and at a small scale such that α s = 0.2. We consider LLx resummation matched to LO (for the gluon-gluon case), and NLLx resummation matched to both NLO and NNLO. All calculations are performed using the HELL (version 2.0) implementation of the ABF construction, and thus incorporate a number of technical improvements which makes the numerical implementation more robust, and allow the matching to NNLO fixed order as well as NLO: a detailed discussion and comparison is given in Refs. [62,63]. The resummation of small-x logarithms is more important at NNLO than at NLO, since Fig. 1 Comparison of the fixed-order gluon-gluon x P gg (x, α s ) (left) and the quark-gluon x P qg (x, α s ) (right) splitting functions with the corresponding LO+LLx, NLO+NLLx and NNLO+NLLx results including small-x resummation. The comparison is performed at a scale such that α s = 0.2 and in the Q 0 MS scheme with n f = 4 active quark flavors at NNLO the fixed-order small-x logarithms give rise to perturbative instabilities at small x, as visible from a comparison of the NLO and NNLO curves in Fig. 1. Indeed, from the left hand plot, one can immediately see that for moderately small values of x NLO gluon evolution is closer to the all-order result at small x than NNLO evolution, since for 10 −6 x 10 −3 the NLO splitting kernels are closer to the best prediction, NNLO+NLLx, than the NNLO ones. Additionally, from the right plot, both resummed results for the gluon-to-quark spitting function are closer to NLO than to NNLO for 10 −5 x 10 −1 . N 3 LO evolution, when available [86,87], will lead to even more significant instabilities at small x, due to the appearance of two extra powers of the small-x logarithms (the leading NLO and NNLO logarithms are accidentally zero), and will make the inclusion of small-x resummation even more crucial.
To facilitate the use of small-x resummation, the HELL code has been interfaced to the public code APFEL [119,120]. Thanks to this APFEL+HELL interface, it is straightforward to perform the PDF evolution (and the computation of DIS structure functions) with the inclusion of small-x resummation effects. Note that APFEL+HELL only implements the so-called "exact" solution of DGLAP evolution, rather than the "truncated" solutions used in ABF (for example in Refs. [44][45][46]), and nowadays routinely in NNPDF fits, in which subleading corrections are systematically expanded out [72]. For this reason we will use the exact solution throughout in this paper, to facilitate comparison between fixed-order and resummed results. Since the difference between the two solutions becomes smaller and smaller when increasing the perturbative order, this choice does not affect significantly our NNLO(+NLLx) results, but care should be taken when comparing the NLO PDFs from those of other NNPDF fits.
We now investigate the effects induced by evolving the PDFs with resummed splitting kernels as compared to standard fixed-order DGLAP splitting functions. In order to illustrate these effects, we take a given input PDF set as fixed at a low scale Q 0 , that is, a common boundary condition, and then evolve it upwards using APFEL+HELL with either fixed-order (NLO or NNLO) or resummed (NLO+NLLx or NNLO+NLLx) theory. In this way, we can determine what are the main differences induced at high scales by small-x resummation in the PDF evolution; we stress, however, that the physical meaning of the resulting comparison is limited, as in a PDF fit with small-x resummation the PDFs at low scales, now taken to be equal to their fixed-order counterparts, are likely to change significantly.
The results of this comparison are collected in Fig. 2, where we show the ratio of the gluon (upper plots) and quark singlet (lower plots) as a function of x for the evolution from a fixed boundary condition at Q 0 = 1.65 GeV up to Q = 100 GeV using either (N)NLO fixed-order theory or (N)NLO+NLLx resummed theory for the DGLAP evolution. In this specific case, the input boundary condition has been chosen to be NNPDF3.1 (N)NLO. We observe that the effects of the different PDF evolution settings are negligible at large and medium x, but can reach up to a few percent at the smallest values of x relevant for the description of the data included in a PDF fit, in particular the HERA structure functions. Specifically, we observe that resummation effects change the NLO evolution quite substantially for both the gluon and the quark singlet, an effect which is reduced at NNLO for the gluon, while it remains of the same size (if not larger) for the quark singlet. Although this study is purely illustrative and by no means predictive, it allows us to conclude that the effect of small-x resummation in PDF evolution Fig. 2 The ratio of the gluon (upper plots) and quark singlet (lower plots) for the evolution from a fixed boundary condition at Q 0 = 1.65 GeV up to Q = 100 GeV using either fixed-order theory (NLO left, NNLO right) or resummed theory (NLO+NLLx left, NNLO+NLLx right) for the DGLAP evolution. In this specific case, the input boundary condition has been chosen to be NNPDF3.1 NLO (NNLO) is in general sizable and will certainly impact the determination of PDFs at small x.

Resummation of DIS structure functions
Resummed results for DIS structure functions, including mass effects, have been recently implemented in the public code HELL, version 2.0 [63]. Analogously to Eq. (2.5), resummed and matched results can be written as where the index a denotes the type of structure function, a = 2, L , 3, while the index i refers to the incoming parton i = q, g. Note that in this paper we only consider NLLx resummation of the partonic coefficient functions, since in DIS there are no LLx contributions. Consistently with the choice made for the evolution, we work in the Q 0 MS scheme. A consistent PDF fit which spans several orders of magnitude in Q 2 further requires us to consider a different num-ber of active quark flavors at different energies, to account for potentially large collinear logarithms due to massive quarks. When crossing the threshold of a given heavy quark, matching conditions which relate the PDFs above and below threshold are needed. These matching conditions also contain small-x logarithmic enhancements, which one can consistently resum. As for DIS coefficient functions, the matching conditions are NLLx, and their resummation, as well as the resummation of the massive coefficient functions [63,110] is available in HELL 2.0. These last ingredients make it straightforward to implement a resummation of the FONLL variable flavor number scheme [121] used in the NNPDF fits.
A careful treatment of charm is essential when addressing the impact of small-x resummation on DIS structure functions, since the kinematic region where resummation is expected to be important (small x and low Q 2 ) is rather close to the charm threshold. We thus fit the initial charm distribution, as in Ref. [78]. The FONLL scheme can be readily extended to fitted charm, in the process receiving an extra contribution [122], denoted IC , which is The proton neutral-current (NC) structure function F 2 (x, Q) as a function of Q for two different values of x (left: x = 10 −3 ; right: x = 10 −5 ) and using different calculational schemes. In the top panels we show the structure function computed in fixed-order perturbation theory (NLO and NNLO). In the middle and bottom panels we show the ratio of resummed results (NLO+NLLx and NNLO+NLLx) to their fixed-order counterparts. In particular, in the middle panel the resum-mation is included in the coefficient function but not in the evolution, while in the bottom panel we resum both coefficient functions and parton evolution. The input boundary condition at Q 0 = 1.65 GeV has been chosen to be NNPDF3.1 NLO (NNLO), and all calculations are performed with α s (m Z ) = 0.118, and a (pole) charm mass m c = 1.51 GeV currently known only at O(α s ) [123,124]. When IC is included, the phenomenological damping adopted in the original FONLL formulation to smooth the transition to the regime in which collinear logarithms are resummed does not have any effect [122,124], and is therefore omitted. Since the O(α s ) IC contribution is then a small correction, we expect the NNLO (O(α 2 s )) and small-x resummation corrections to IC to be practically insignificant (see Ref. [78] for a detailed discussion of this issue).
To obtain a first qualitative estimate of the impact of smallx resummation in the DIS structure functions, we can compare theoretical predictions at (N)NLO with predictions that include resummation. To disentangle the effect of resummation on PDF evolution from that in the coefficient functions in the Q 0 MS scheme, we take into account the effect of resummation in two steps. First, we compute structure functions with the same (fixed-order) input PDFs and include small-x resummation in the coefficient functions only. As a second step, we include resummation also in the DGLAP evolution, using a fixed input PDF boundary condition at a small scale Q 0 = 1.65 GeV, as previously done in Fig. 2. Since, as already noticed, the use of a fixed boundary condition at a small scale is not particularly physical, these results should be interpreted with care.
The proton structure function F 2 (x, Q) in neutral-current (NC) DIS is shown in Fig. 3 as a function of Q for two values of x, one moderate (x = 10 −3 , left plot) and one small (x = 10 −5 , right plot). The upper panel of each plot shows the NLO and NNLO results. The middle panel shows the ratio of resummed (N)NLO+NLLx theory over the fixed-order (N)NLO results, including resummation only in coefficient functions. The lower panel, instead, shows the same ratio but with resummation included also in PDF evolution. In all cases, we take the NNPDF3.1 boundary condition at (N)NLO at Q 0 = 1.65 GeV. As mentioned above, heavy-quark mass effects are included using the FONLL-B (C) scheme [121,122,124] for the NLO (NNLO) calculations, supplemented with small-x resummed contribution for the (N)NLO+NLLx as described in Ref. [63].
The comparison in Fig. 3 is interesting from several points of view. First of all, we observe that when resummation is included only in the coefficient functions its effect is rather mild, almost negligible when matched to NNLO, even at rather small x and at low scales. On the other hand, when including resummation in the PDF evolution, the situation changes. In this case, we note that the differences between fixed-order and resummation are larger, thus showing that in F 2 much of the impact of small-x resummation arises from the PDF evolution. Moreover, the effects are always greater at NNLO than at NLO: at NNLO, effects of smallx resummation can reach 10 percent already for x 10 −3 , and 20 percent for x 10 −5 . This discussion suggests that at the level of PDF fits we expect little differences between the fixed-order and resummed cases at NLO, but more significant differences at NNLO. Next, in Fig. 4 we show the same comparison as in Fig. 3 but now for F c 2 (x, Q), the charm component of the proton structure function F 2 (x, Q). By comparing Figs. 3 and 4 we observe that the impact of small-x resummation for inclusive and charm structure functions is similar, except just above the charm threshold where the effects of the resummation in the charm coefficient function can be substantial. From this comparison, we see the importance of a careful treatment of mass effects close to the charm threshold, since these can change the size of the effect of small-x resummation.
Finally, in Fig. 5 we show the corresponding comparison but this time for the longitudinal structure function F L (x, Q) in neutral-current DIS. Here we find that resummation effects in the coefficient functions only are substantially larger than in F 2 , and are now larger when matching resummation to NNLO than to NLO. When resummation is included also in PDF evolution, the overall effect of resummation on F L is somewhat reduced at NLO, thus showing some sort of compensation of the effects in PDF evolution and in partonic coefficient functions, while it is enlarged at NNLO, which now reaches about a 30% deviation at x = 10 −5 at small Q ∼ 5 GeV. The global pattern is similar to F 2 , with differences smaller at NLO and more significant at NNLO, though overall effect is somewhat bigger, consistently with the fact that F L is singlet dominated. Given that F L contributes to the measured reduced cross-sections σ r,NC at high y, which for the HERA kinematics corresponds to small x and Q 2 , this effect should be relevant for PDF fits.

Fitting strategy
In this section we discuss the settings of the NNPDF3.1 fits with small-x resummation, as well as of their fixed-order counterparts, which are used as baseline comparisons. In the following, we will denote these fits as NNPDF3.1sx, each of them consisting of N rep = 100 Monte Carlo replicas. We briefly present the input dataset, and review the theoretical treatment of the deep-inelastic and hadronic data used in the fit. We also discuss the strategy adopted for choosing appropriate kinematic cuts for both DIS and hadronic processes.

Fit settings
The settings of the fits described in this work follow closely those of the recent NNPDF3.1 global analysis [79]. In particular, the same input dataset is used, which includes fixedtarget [125][126][127][128][129][130][131][132] [170][171][172][173][174] at √ s = 7 and 8 TeV. As in the NNPDF3.1 analysis, the charm PDF is fitted alongside the light quark PDFs [78], rather than being generated entirely from perturbative evolution off gluons and light quarks. As usual in NNPDF, we use heavy-quark pole masses [175], and the charm quark pole mass is taken to be m c = 1.51 GeV. In all the results presented here we take α s (m Z ) = 0.118.
The initial scale Q 0 at which PDFs are parametrized is chosen to be Q 0 = 1.64 GeV, i.e. Q 2 0 = 2.69 GeV 2 , which is slightly smaller than the initial scale adopted in the NNPDF3.1 analysis, namely Q 0 = 1.65 GeV. The main motivation for this choice of initial scale is to be able to include the Q 2 = 2.7 GeV 2 bin in the HERA inclusive struc-ture function data [68], which is expected to be particularly sensitive to the effects of small-x resummation, and that was excluded from NNPDF3.1. At the same time, the initial scale cannot be too low, to avoid entering a region in which α s is too large and the numerical reliability of the small-x resummation implemented in the HELL code would be lost. 2 In this work we have produced fits at fixed-order NLO and NNLO accuracy and corresponding resummed fits at NLO+NLLx and NNLO+NLLx accuracy. In the resummed fits, small-x resummation is included both in the solution of the evolution equations and in the deep-inelastic coefficient functions as discussed in Sect. 2. Heavy-quark mass effects are accounted for using the FONLL-B and FONLL-C general-mass scheme [121,122,124] for the NLO and NNLO fits, respectively, modified to include small-x resummation effects when NLO+NLLx and NNLO+NLLx theory is used as previously described.
Theoretical predictions for the Drell-Yan fixed-target and the hadron collider (Tevatron and LHC) cross-sections are obtained using fixed-order or resummed DGLAP evolution for (N)NLO and (N)NLO+NLLx fits, respectively, but with their partonic cross-sections always evaluated at the corresponding fixed order. This approximation is due to the fact that the implementation of hadronic processes in HELL is still work in progress. To account for this limitation, we cut all data in kinematic regions where small-x corrections are expected to be significant, as explained in Sect. 3.2 below.
The settings for the evaluation of the hadronic hardscattering matrix elements are the same as in NNPDF3.1, namely we use fast NLO calculations as generated by APPLgrid [176] and FastNLO [177] tables, which are combined before the fit with the DGLAP evolution kernels by means of the APFELgrid interface [178]. For the NNLO fits, NNLO/NLO point-by-point K -factors are used [79] using specific codes for each process: we use the code of [179,180] for tt differential distributions [181]; for the Z p T distributions we use the calculation of [182,183]; for Drell-Yan production we use FEWZ [184]; while jet cross-sections are treated using NLO matrix elements supplemented by scale variation as additional theory systematics.
For comparison purposes, we have also produced DISonly fits for which small-x resummation is included in both evolution and coefficient functions for all data points included in the fit. That is, in such fit, fully consistent smallx resummed theory is used for the entire dataset. Moreover, while PDF uncertainties are of course much larger due to the lack of hadronic data, the constraints from the HERA structure functions are still the dominant ones in the smallx region. The comparison between the global and DIS-only NNPDF3.1sx fits is discussed in Sect. 4.2.1.

Kinematic cuts
In the NNPDF3.1sx analysis, we apply the same experimental cuts as those of the NNPDF3.1 fit [79] with two main differences. First, as discussed above, the lower Q 2 cut is reduced from Q 2 min = 3.49 GeV 2 in NNPDF3.1 to Q 2 min = 2.69 GeV 2 here. Thanks to this lower cut, we can now include a further bin of the HERA inclusive crosssection data, specifically the one with Q 2 = 2.7 GeV 2 . In turn, this allows us to slightly extend the kinematic coverage of the small-x region, from x min 4.6 × 10 −5 before, down to x min 3 × 10 −5 now. This lower cut also affects a handful of points at low Q 2 (although at larger values of x) of other fixed-target DIS experiments, which are therefore also included in the NNPDF3.1sx fits but not in NNPDF3.1. The cut on W 2 ≥ 12.5 GeV 2 remains the same.
Moreover, no additional cuts are applied to the HERA charm-production cross-sections as compared to the inclusive structure functions. This was not the case in NNPDF3.1, where some points at small-x and Q 2 were excluded in the NNLO fit, specifically those with Q 2 ≤ 8 GeV 2 . We have explicitly verified that the inclusion of these extra points does not affect the resulting PDFs, though the χ 2 of the F c 2 data becomes somewhat worse at NNLO. Taking into account these two differences, from HERA we fit 1162 points for the inclusive structure functions and 47 points for the F c 2 data, to be compared with 1145 (1145) and 47 (37) in NNPDF3.1 NLO (NNLO), respectively. The number of data points N dat for each of the DIS experiments included in NNPDF3.1sx is collected in Table 1.
The second main difference with respect to the NNPDF3.1 kinematic cuts is related to hadronic data. As already discussed, for hadronic processes small-x resummation effects are included only in PDF evolution but not in the partonic cross-sections. Therefore, in order to avoid biasing the fit results, in the NNPDF3.1sx fits we include only those hadronic data for which the effects of small-x resummation on the coefficient function can be assumed to be negligible.
Quantifying the impact of small-x resummation on the partonic coefficient functions would require the knowledge of such resummation. Therefore, in order to estimate the region of sensitivity to small-x logarithms, we resort to a more qualitative argument. The foundation of this argument is the observation that in a generic factorization scheme large logarithms appear both in the partonic coefficient functions and in the partonic evolution factors; in general, resummation corrections are thus expected to have a similar size both in the evolution and in the coefficient functions. This naive expectation is indeed confirmed by explicit calculations of hadronic resummed cross-sections [44,185], where it was found that the most common situation is a partial cancellation between the resummation corrections from evolution and those in the partonic cross-section. It follows that estimates based on the corrections due to resummed evolution alone will probably be conservative, in the sense that they will over-estimate the total resummation correction to the hadronic cross-section.
In order to implement these cuts, we first introduce a parametrization of the resummation region in the (x, Q 2 ) plane. Small-x logarithmic corrections should in principle be resummed when α s (Q 2 ) ln 1/x approaches unity, since the fixed-order perturbative expansion then breaks down. We thus define our kinematic cut to the hadronic data in the NNPDF3.1sx fits such as to removes those data points for which where H cut 1 is a fixed parameter: the smaller H cut , the more data are removed. Assuming one-loop running for the strong coupling constant (which is enough for our purposes), Eq. (3.1) can instead be expressed as where 88 MeV is the QCD Landau pole for n f = 5, and β 0 0.61. Thus the cut is a straight line in the plane of ln 1 x and ln Q 2 2 , with gradient β 0 H cut . Note that the variable x used in the definition of the cut, Eq. (3.1), can in general only be related to the final-state kinematic variables of hadronic observables by assuming leadingorder kinematics. To see how this works in practice, consider for example weak gauge boson production: then Q 2 = M 2 V , and for fixed √ s the cut translates into a maximum rapidity Thus in the case of W boson production at √ s = 7 TeV, a cut of the form of Eq. (3.2) with H cut = 0.5 (0.7) would imply that cross-sections with rapidities above y max 0.3 (1.3) Fig. 6 The ratio of hadronic cross-sections included NNPDF3.1 computed using a fixed input PDF at Q 0 = 1.65 GeV (in this case NNPDF3.1 NNLO) using either NNLO+NLLx or NNLO theory for PDF evolution, always with NNLO partonic cross-sections. We show the results for ATLAS, CMS, LHCb, and the Tevatron, indicating the division of each experiment into families of processes. The empty blue triangles indicate those data points that are excluded from the NNPDF3.1sx fits with the default cut H cut = 0.6, while the filled red ones indicate the points that satisfy the condition Eq. (3.2) would be excluded from the fit. In this case, the first (tighter) cut excludes all the LHC gauge boson production data except for a handful of points from the ATLAS and CMS measurements in the most central rapidity region. The second (looser) cut instead allows one to include most of the ATLAS and CMS gauge boson production data. However, the LHCb measurements are removed altogether for both values of the cut, highlighting the sensitivity of forward W, Z production data to the small-x region. It remains to determine the optimal value of H cut , in a way that minimizes at the same time the amount of information lost from the dataset reduction, but also the possible theoretical bias due to the missing small-x resummed coefficient functions. In this work we will present results with three different values, namely H cut = 0.5, 0.6 and 0.7. In Sect. 4 we will motivate the choice of H cut = 0.6 as our default value, and show explicitly how the main findings on this work are independent of the specific value of H cut adopted.
Here we attempt to provide an a priori argument to justify our choice by estimating the size of the resummation corrections through a comparison of the results obtained with fixedorder and resummed parton evolution. Specifically, we take a fixed input PDF set (NNPDF3.1 NNLO) at Q 0 = 1.65 GeV and evolve it using either NNLO or NNLO+NLLx theory, and then compute the convolution with fixed-order partonic coefficient functions. The comparison is represented in Fig. 6, where we show the ratio of hadronic cross-sections computed using NNLO+NLLx evolution over those computed using NNLO evolution. We show the results for ATLAS, CMS, LHCb, and the Tevatron data points included in NNPDF3.1, indicating the division of each experiment into families of processes. From this comparison, we see that the effects of small-x resummation are likely to be significant only for the W and Z Drell-Yan data, where they could be as large as up to ∼ 5% for ATLAS and CMS, and up to ∼ 8% for the forward LHCb measurements, while they are most likely 2), below which the hadronic data is excluded from the fit. For hadronic processes, the LO kinematics have been used to determine the (x, Q 2 ) values associated to each data bin negligible for all other collider processes, such as jets, the Z p T , and top-quark pair production. Given that the collider DY data have rather small experimental uncertainties, of the order of a few percent or even smaller, we should ensure that we cut data where the effects of small-x resummation could be larger than ∼ 2% (to be conservative). We see from Fig. 6 that this is indeed achieved with the default value of H cut = 0.6: for the included points, differences are always smaller than this threshold.
To summarize this discussion of the kinematic cuts in the NNPDF3.1sx fits, we show in Fig. 7 the kinematic coverage in the (x, Q 2 ) plane of the data included in the present analysis, for the default value H cut = 0.6 of the cut to the hadronic data. As mentioned above, for hadronic processes the LO kinematics have been used to determine the values of x and Q 2 associated to each data bin. The diagonal line indicates the region below which the cut defined in Eq. (3.2) removes hadronic data. As a consequence of the kinematic cuts, the hadronic dataset is restricted to the large-Q 2 and medium-and large-x region.
In Table 2 we show the number of data points for the hadronic data in the NNPDF3.1sx NNLO fits for with H cut = 0.5, 0.6 and 0.7. The number in brackets corresponds to the values for the NLO fits, since the kinematic cuts of the NNPDF3.1 fits [79] are slightly different at NLO and at NNLO. The main effect of the H cut is on the Drell-Yan pre-diction measurements from ATLAS and CMS, which in turn affects the quark and antiquark flavor separation, and the Z p T distributions, which provide information on the gluon. On the other hand, the inclusive jet and top-quark pair production data, which are mostly sensitive to the large-x region, are essentially unaffected by the cut. For completeness, we also provide the values of N dat when no cut is applied at all (H cut = ∞). In the latter case, the fit also includes 85 (93) LHCb experimental points at NNLO (NLO).

Parton distributions with small-x resummation
In this section we present the main results of this work, namely the NNPDF3.1sx fits including the effects of smallx resummation. We will present first the DIS-only fits and then the global fits, based on the dataset described in Sect. 3. Unless otherwise specified, for the global fits we will use the default cut H cut = 0.6 for the hadronic data.
In the following, we will first discuss the DIS-only fits, showing how small-x resummation improves the fit quality and affects the shape of the PDFs. We then move to the global fits, and compare them to the DIS-only ones. We find that the qualitative results are similar, though PDF uncertainties are reduced. We show the impact of resummation on the PDFs, and study the dependence on the cut used to remove the hadronic data potentially sensitive to small-x logarithms and for which we do not yet include resummation. We show how our default choice for H cut does not bias the fit, and still allows us to determine PDFs whose uncertainties are competitive with those of NNPDF3.1. We discuss in detail the role of the additional low-Q 2 HERA bin that we include in this fit for the first time, and how small-x resummed theory is able to fit it satisfactorily. We will further inspect the improved description of the HERA data in Sect. 5, where we will perform a number of diagnostic studies aimed at quantifying the onset of BFKL dynamics in the inclusive HERA structure functions.

DIS-only fits
Let us start our discussion by considering the DIS-only fits, in which we include all the DIS data from fixed-target and collider experiments described in Sect. 3. For all these data, we have a complete theoretical description at resummed level, thus allowing us to perform a fully consistent small- Table 3 The values of χ 2 /N dat for the total and the individual datasets included in the DIS-only NNPDF3.1sx NLO, NLO+NLLx, NNLO and NNLO+NLLx fits. The number of data points N dat for each experiment is indicated in Table 1. In addition, we also indicate the absolute difference χ 2 between the resummed and fixed-order results, Eq. (4.1). We indicate with a dash the case | χ 2 | < 0.5 x resummed fit. First of all, in Table 3 we collect the χ 2 /N dat values for the total and individual datasets computed with the PDFs fitted using NLO, NLO+NLLx, NNLO and NNLO+NLLx theory. The χ 2 values are computed using the experimental definition of the covariance matrix, while the t 0 definition [186] was instead used during the fits, as customary in the NNPDF analyses. In addition, we also show the difference in χ 2 between the resummed and fixed-order results, which is useful to gauge how statistically significant are the differences between the fixed-order and resummed results for each experiment. We immediately observe that the NNLO+NLLx fit has a total χ 2 /N dat that improves markedly with respect to the NNLO result, which instead gives the highest value of χ 2 /N dat . The total χ 2 /N dat is essentially the same in the NLO, NLO+NLLx, and NNLO+NLLx fits. As illustrated by the χ 2 values of Table 3, the bulk of the difference in the fit quality between the NNLO and NNLO+NLLx fits arises from the HERA inclusive neutral-current and charm datasets, which probe the smallest values of x, and whose χ 2 /N dat decrease from 1.16 to 1.12 ( χ 2 = −47) and from 2.45 to 1.24 ( χ 2 = −57), respectively.
We note that the χ 2 /N dat of the charm dataset is rather high at NNLO. In fact, the description of the charm data can be rather sensitive to the details of the heavy quark scheme. For instance, we can set to zero the IC term discussed in Sect. 2.3, thus allowing the inclusion of a phenomenologically induced damping factor which has the role of suppressing formally subleading terms numerically relevant at scales close to the charm threshold (see [121,122,124]). 3 When the damping is included, we find that recomputing the χ 2 /N dat of the charm dataset it becomes 1.10 at NNLO. On the other hand, the quality of resummed theory is very stable with respect to such a variation, and the χ 2 /N dat of the charm data becomes 1.23 ( χ 2 = +6). The rather high value of the charm data χ 2 at NNLO with our default settings is mostly driven by a poor description of the low-x and low-Q 2 bins. Indeed, if we restrict our attention to the region which survives the more conservative cut used in NNPDF3.1 (Q 2 ≥ 8 GeV 2 for the HERA charm data), we obtain χ 2 /N dat = 1.38 at NNLO and 1.35 at NNLO+NLLx ( χ 2 = −1) using our default settings. The low-Q 2 region is somewhat affected by how the subleading terms are treated -ultimately, this choice is driven by phenomenological reasons, and therefore it is possible that by tuning them one may achieve a satisfactory description of the data at NNLO, for instance by mimicking a perturbative behavior 4 ; however, the same choice may be suboptimal at the resummed level. Since at NLO(+NLLx) and with FONLL-B we achieve a satisfactory description of the charm data for all 47 points both at fixed-order and at resummed level, here we shall use the same theory settings of the NNPDF3.1 paper, and interpret the more marked dependence on the subleading terms as a limitation of the fixed-order theory at NNLO.
We further observe that the description of the fixed-target DIS experiments, sensitive to the medium and small-x region, is not significantly affected by the inclusion of small-x resummation, giving us confidence that the resummed and matched predictions reduce to their fixed-order counterpart where they should. The only exception is the slight decrease in fit quality between the NNLO and NNLO+NLLx fits for BCDMS and CHORUS ( χ 2 = +14 and +26, respectively). As we will show in the next section, most of these Another interesting result from Table 3 is that the effect of resummation is instead much less marked at NLO. Indeed, the NLO and NLO+NLLx fits have very similar χ 2 /N dat : in particular the χ 2 change of the HERA inclusive (charm) dataset is rather small, χ 2 = +6 (−5). This is again not surprising, as the whole point of resummation is to cure instabilities in the fixed-order perturbative expansion, by removing the large logarithms causing the instability and replacing them with all-order results. Thus the resummation is more important at NNLO than at NLO, and indeed would probably be yet more important at the next perturbative order (N 3 LO).
We can see this result more clearly by considering the resulting fitted PDFs and their uncertainties. In Fig. 8 we show the ratio between the gluon (left) and the total quark singlet (right) at Q = 100 GeV in the NLO+NLLx fit as compared to the NLO baseline (upper plots) and in the NNLO+NLLx fit as compared to the NNLO baseline (lower plots). In this comparison, as well as in subsequent PDF plots, the bands represent the 68% confidence level PDF uncertainty. Consider first the NLO+NLLx fit. Here the resummation has a moderate effect: the resummed gluon PDF is somewhat enhanced between x = 10 −5 and x = 10 −2 , with the PDF uncertainty bands only partially overlapping, whilst the shift in central values for the singlet is well within the PDF uncertainties. This remains true down to the smallest values of x: even for values as small as x 10 −6 the shifts of the central value of the singlet and the gluon PDF due to the resummation are less than 10%. This is a consequence of the fact that, as discussed in Sect. 2, NLO theory is a reasonably good approximation to the fully resummed result at small x, and any differences are such that can be reabsorbed into small changes in the gluon PDF.
The situation is rather different at NNLO+NLLx. In this case, we see that starting from x 10 −3 the resummed gluons and quarks are systematically higher than in the baseline NNLO fit, by an amount which ranges from 10% for x ∼ 10 −4 up to 20% for x ∼ 10 −5 (though note that Note that we are performing these comparisons at the electroweak scale Q ∼ 100 GeV, where there are no DIS data and where the effect of resummed evolution is combined with the change of the fitted PDFs at low scales. This has the advantage of showing that several observables at the LHC characterized by electroweak scales are likely to be sensitive to small-x resummation through the PDFs, particularly when measurements can be performed at high rapidities. Therefore, for such observables, the use of small-x resummed PDFs (and coefficient functions) is probably going to be necessary in order to obtain reliable theoretical predictions.
In Fig. 8 we observed that including resummation leads to a significantly larger shift in the small-x quark singlet and gluon PDFs at NNLO than at NLO. This is so despite the fact that from the point of view of small-x resummation the information added is the same in both cases, and that the resummed splitting and coefficient functions at small x are quite similar whichever fixed-order calculation they are matched to. The explanation of this paradoxical result is that fixed-order perturbation theory is unstable at small x due to the small-x logarithms, and while this instability is quite small at NLO, due to accidental zeros in some of the coefficients, it is significant at NNLO, and would probably become very substantial at N 3 LO. To better illustrate this effect, and the way it is cured by resummation, in Fig. 9 we compare the NLO, NNLO and NNLO+NLLx results for the gluon and singlet PDFs in the baseline fits at Q = 100 GeV, normalized to the NLO prediction. We find that the NNLO results are systematically below the NLO ones for x ≤ 10 −2 , and that the net effect of adding NLLx resummation to the NNLO fit is to bring it more in line with the NLO (and thus as well with the NLO+NLLx) result. This provides an explanation of our previous observation that NNLO theory fits small-x DIS data worse than NLO, while NNLO+NLLx provides the best description of all.
So far we focused on the gluon and quark singlet, as smallx resummation affects PDFs in the singlet sector. To quantify the effect of resummation on the PDFs in the physical basis it is convenient to use a distance estimator, as defined in Refs. [74,77]. This allows us to represent in a concise way how two PDF fits differ among themselves, both at the level of central values and of PDF uncertainties. In Fig. 10 we show these distances between the central values (left) and the PDF uncertainties (right) of the NNPDF3.1sx NNLO and NNLO+NLLx fits at Q = 100 GeV. Since these fits are based on N rep = 100 replicas each, a distance of d ∼ 10 corresponds to a variation of one sigma of the central values or the PDF uncertainties in units of the corresponding standard deviation.
From the comparison in Fig. 10 we see that the impact of using NNLO+NLLx theory peaks between x 10 −3 and x 10 −5 , where d 30, meaning that the central value shifts by more than three times the corresponding PDF uncertainty. The gluon is the most affected PDF, followed by the charm and then by the light quark PDFs. Note that the differences are not restricted to the region of very small-x, since for gluons d ∼ 10 already at x 5 · 10 −3 , relevant for the production of electroweak scale particles such as W and Z bosons at the LHC. On the other hand, the impact of using NNLO+NLLx theory is as expected small for the PDF uncertainties, since from the experimental point of view very little new information is being added into the fit. However, as we will discuss in greater detail in Sect. 4.2.4, adding small-x resummation has allowed us to lower the minimum value of Q 2 for the HERA data included in the fits -which in turn extends to smaller x the PDF kinematic coverage, thus reducing PDF uncertainties in the very small-x region.  Comparison between the gluon (left) and the total quark singlet (right plots) from the NNLO and NNLO+NLLx DIS-only fits, including the variant of the resummation which differs by subleading terms, as discussed in the text Before moving to the global fits, we want to briefly investigate how our results are sensitive to unknown subleading logarithmic contributions. Indeed, the results of Ref. [63] are provided with an uncertainty band aimed at estimating the impact of subleading (NNLLx) contributions not predicted by NLLx resummation. Ideally, the uncertainty band should be included as a theory uncertainty in the fit procedure; however, at the moment the inclusion of theory uncertainties in PDF fits is still under study. Nevertheless, we can investigate the effects of such uncertainties by performing another fit in which we change the resummation by subleading terms. A simple way to do it in a consistent manner is to vary by subleading terms the anomalous dimension used for the resummation of coefficient functions and of P qg . As the resummed gluon splitting function depends on the resummed P qg , all splitting functions and coefficient functions are affected by this change. More specifically, the so-called LL anomalous dimension used in HELL 2.0 (and hence in this work) is replaced with the full NLLx anomalous dimension, as proposed originally in Ref. [46]. The effect of this variation is contained within the uncertainty bands of Ref. [63].
The result of this fit, based on the same DIS-only dataset considered so far and performed at NNLO+NLLx accuracy, is fully consistent with that obtained with the baseline theory settings. The fit quality is essentially unaffected, and the χ 2 variations with respect to the numbers in Table 3 are compatible with statistical fluctuations. Most PDFs are not sensitive to this variation, except the gluon and the quark singlet, which do change a little, to accommodate the different subleading terms in the splitting functions and coefficient functions. These PDFs are shown in Fig. 11 and compared with the default HELL 2.0 result. In both cases the new PDFs are smaller than our default ones, i.e. closer to the NNLO results. This is mostly due to a harder resummed P qg in the varied resummation, which is therefore closer to its NNLO counterpart, at intermediate values of x, than our default resummation. For the gluon in particular, the new results are not compatible within the uncertainty bands with our default fit, highlighting that the PDF uncertainty does not cover the theory uncertainty from missing higher orders. However, all the qualitative conclusions remain unchanged.

Global fits
We now turn to consider the global fits, based on the complete dataset described in Sect. 3.2. We first show the results of the fits, obtained with the default cut parameter H cut = 0.6, highlighting similarities and differences with respect to the DIS-only fits, and we discuss the impact of resummation on the PDFs. We then study the dependence of our results upon variation of the value of H cut . Finally, we discuss in some detail the description of the low-Q 2 HERA bin which we include in the NNPDF31sx fits.

Fit results and comparison to the DIS-only fits
We start by considering the quality of the global NNPDF3.1sx fits at NLO, NLO+NLLx, NNLO and NNLO+NLLx, using the default value of H cut = 0.6 for the hadronic data cut discussed in Sect. 3.2. The values of the χ 2 /N dat for the total and the individual datasets are shown in Table 4. As in the DIS-only case, in this table we also include the absolute χ 2 difference between the resummed and fixed-order results, χ 2 Eq. (4.1). We observe that the NNPDF3.1sx fit based on NNLO+NLLx theory leads to the best overall fit quality, χ 2 /N dat = 1.100. The NNLO fit, on the other hand, has again the highest χ 2 /N dat = 1.130, so that the overall improvement is χ 2 = −121. Whilst resummation proves particularly beneficial at NNLO, the effect at NLO is very mild; the χ 2 /N dat 1.120 at NLO+NLLx is compatible, within statistical fluctuations, with the 1.117 obtained with fixed-order theory, that is, χ 2 = +11. Note that in the NNPDF3.1 fits the NNLO χ 2 was markedly better than the NLO one [79]: this is no longer the case here, since the high-precision Drell-Yan and Z p T data points, which are poorly described by NLO theory, are now partly removed by the H cut cut.
The improvement of the χ 2 at NNLO+NLLx is essentially due to the HERA charm and neutral-current structure function data. On one hand, as we already noticed in the DIS-only fits, by using NNLO+NLLx theory one achieves an improved description of the precise HERA NC inclusive structure function measurements, whose χ 2 /N dat decreases from 1.17 in the NNLO fit to 1.11 in the NNLO+NLLx fit, A marked improvement is also achieved for the HERA charm cross-sections, whose χ 2 /N dat goes down from 2.33 to 1.14, χ 2 = −56. These two datasets are thus sufficient to explain the overall improvement in the total χ 2 .
We also find that NNLO theory describes better than the corresponding NLO theory the ATLAS and CMS measurements, particularly the recent high-precision data such as the ATLAS W, Z 2011 rapidity distributions, and the ATLAS and CMS 8 TeV Z p T distributions. Specifically, the χ 2 /N dat total values for ATLAS and CMS is 1.18 (1.16) and 0.97 (0.92) in the NLO(+NLLx) fits, respectively, decreasing to 0.99 (0.98) and 0.86 (0.85) when using NNLO (+NLLx) theory. It is interesting that in all cases the resummed fits are slightly better than their fixed-order counterparts.
Despite the improved description of the large-Q 2 collider data with respect to the NLO theory, the NNLO fit turns out to have the highest χ 2 of the four theories, as in the DIS-only case. The main reason is the poor description of the HERA inclusive and charm dataset, which contain almost one third (N dat = 1209) of the number of data points included in the fit (N dat = 3930). Moreover, we observe that the effects of small-x resummation at NNLO are confined to the HERA data; the differences between the χ 2 values of the (N)NLO and (N)NLO+NLLx fits for the other datasets are being all rather small. This is in agreement with the findings of the DIS-only fits, and with the fact that hadronic data potentially sensitive to small-x effects have been cut. Specifically, in the NNLO fits there is no other dataset besides the HERA inclusive and charm data with | χ 2 | ≥ 10.
Comparing the values of the χ 2 /N dat for the DIS experiments in the global and DIS-only fits, we notice that once resummation is accounted for, the global fit is if anything slightly better than the DIS-only fit. In particular for the inclusive HERA data, where χ 2 /N dat is 1. 16 (1.12) at NNLO(+NLLx) in the DIS-only fits, we have χ 2 /N dat = 1.17 (1.11) in the global fits, so that χ 2 decreases from −47 to −62 in the global fit. The other significant difference between the global and DIS-only fits appears in the NuTeV dimuon data, which is fit somewhat less well in the global fit (irrespective of resummation) due to the tension with the LHC data relative to the proton strangeness, especially with the ATLAS W, Z 2011 rapidity distributions [79].
We now move to the impact of small-x resummation on the global dataset PDFs. First, we quantify the differences between the global and the DIS-only fits, taking as a representative the baseline fixed-order NNLO fit. We start by showing the distance estimator in Fig. 12, both for the central value (left) and the PDF uncertainty (right), at Q = 100 GeV. Due to the conservative kinematic cut imposed on the collider observables, the distances between the global and DIS-only fits are moderate and localized to the medium and large-x region, while the small-x region is pretty much unchanged. The PDF flavor which is most affected is the charm PDF, whose distance is about 10 for x ∼ 10 −2 . The decrease in PDF uncertainties in the global dataset at medium and large-x is clearly visible, especially for the gluon PDF which is only constrained in an indirect way by the DIS structure function data.
In Fig. 13 we show a direct comparison between the gluon (left) and the total quark singlet (right) at Q = 100 GeV between the NNPDF3.

Features of the small-x resummed PDFs from the global fit
The comparison done so far demonstrates that the use of the global dataset is very beneficial from the point of view of the PDF uncertainties, while it does not affect the qualitative and quantitative results at small x. Therefore, the global fits will be considered from now on the baseline NNPDF3.1sx fits, and we will focus on these results for subsequent applications and studies. Therefore, before moving forward, it is interesting to analyse the features of these fits in more detail. We focus on the results at NNLO and NNLO+NLLx, as at NLO the impact of resummation is less significant (just as in the DIS-only fits) and also less important from the point of view of applications to the LHC and future high-energy collider physics. In Fig. 14 we show the same distance comparison as in Fig. 10 but now for the NNPDF3.1sx global fits. By comparing this figure with the corresponding DIS-only case, we see that in the global fits the qualitative features are the same. The increased significance of the distances at large x observed in the global fit as compared to the DIS-only is a direct consequence of the reduced PDF uncertainties in the global fit, rather than to a shift in the central values.
To visualize these effects, in Fig. 15 we show the flavor combinations most affected by resummation (as indicated in the distance plot of Fig. 14), namely the gluon, charm, up and down PDFs, at a typical electroweak scale of Q = 100 GeV. The impact of NLLx resummation is very similar for all the quark combinations: the effect is mild for x 10 −3 , whilst it increases at small x, by an amount which is, however, mostly consistent with the one or two sigma PDF uncertainties. The effect is rather more marked for the gluon, where the NNLO+NLLx fit can be up to 30% bigger at x 10 −6 , well outside the uncertainty band. Thus the main impact of high-energy resummation is to strongly enhance the gluon and mildly enhance the quarks at small-x.
To conclude the discussion of the results of the global NNPDF3.1sx fits, we move away from the electroweak scale and consider the PDFs at the input parametrization scale Q 0 . This comparison is interesting because it disentangles the effects of small-x resummation on the fitted PDFs from those due to the evolution from low to high scales. With this motivation, we show in Fig. 16 the gluon and the quark singlet at the fit scale Q 0 = 1.64 GeV. In the case of the total quark singlet, we see that the impact of resummation is moderate, with a one-sigma increase at small x in the NNLO+NLLx fit which helps to improve the fit to the low Q 2 HERA data. The slightly larger effects seen at higher scales are thus mostly driven by the evolution that mixes the singlet with the gluon. On the other hand, the effects of resummation are more marked for the fitted gluon, where we see explicitly a drop in the NNLO gluon at small x driven by perturbative instability, which disappears on resummation in such a way that the NNLO+NLLx gluon is rather flat, and indeed very close to the NLO and NLO+NLLx gluon. Note that the resummation thus extends the perturbative region at small x: even at Q 0 = 1.64 GeV the fitted gluon remains stable, and it seems likely that one would have to go to even lower scales (below the charm threshold) before the kind of instability seen in NNLO fixed-order perturbation theory sets in. Note that we would not expect the same to be true of N 3 LO perturbation theory: the unresummed logarithms at N 3 LO are considerably larger than those at NNLO, and thus the need for resummation at N 3 LO would be even more pressing than at NNLO.

Dependence on the value of H cut
Thus far we have only discussed the results of the global fit obtained using the default cut to the hadronic data, identified as H cut = 0.6. We now discuss the dependence of the fit results with respect to variations of this choice, both from the point of view of the fit quality and of the impact at the PDF level. In doing so, we provide further motivation for the choice of H cut = 0.6 for our default global fits.
To begin with, we study the dependence of the quality of the NNPDF3.1sx fits as a function of the value of the cut parameter H cut applied to the hadronic data. In Table 5 we show a comparison of the NNLO and NNLO+NLLx values of the χ 2 /N dat for the fits with H cut = 0.5, 0.6 and 0.7. In addition, to better appreciate the variations for χ 2 for the fits with different H cut cuts, in Table 6 we also show the differences cut for each experiment.
The main general feature that we note from the comparisons in Tables 5 and 6 is that the χ 2 /N dat values exhibit a rather moderate dependence on the specific value of the kinematic cut to the hadronic data.
Concerning the total dataset, the χ 2 /N dat values slightly increase as H cut is raised and the dataset is enlarged: in particular, for the NNLO (NNLO+NLLx) fits, the values of χ 2 /N dat for the total dataset are 1.120, 1.130, and 1.142 (1.085, 1.100, and 1.112) for H cut = 0.5, 0.6 and 0.7 respectively. The fact that the fit quality of both the fixed-order and the resummed fits is slightly better for H cut = 0.5 is a direct consequence of the more restrictive dataset.
In the case of the NNLO+NLLx fits, the difference between the χ 2 /N dat of the fit with H cut = 0.6 and the fit with H cut = 0.7 is larger than a statistical fluctuation. This might be an indication that the deterioration of the fit with H cut = 0.7 could be related to non-negligible effects of unresummed small-x logarithms in the extra hadronic data that are included in this fit. This conjecture is supported by the fact that, while with H cut = 0.6 the resummation improves the total χ 2 over the fixed order by around 120 points, for H cut = 0.7 the improvement is reduced to less than 100 points. On the other hand, the same trend is also visible in the NNLO fits, and there it can be partly explained by the contributions from some collider points that are in tension between the DIS data, for instance, the ATLAS W, Z 2011 rapidity distributions and the neutrino data. We also find that the more conservative fit with H cut = 0.5 also improves with resummation by even more than the H cut = 0.6 fit (around 140 points), thus suggesting that our default cut value is safe, in the sense that it is not affected by large unresummed logarithms in the hadronic processes.
We further investigate the impact on the PDFs of the various choices of H cut . We show in Fig. 17 the distance estimator to compare the default NNPDF3.1sx NNLO+NLLx global fit with H cut = 0.6 with the corresponding fits with the H cut = 0.7 and H cut = 0.5 fits. In terms of central values, we see that differences are well below PDF uncertainties (which corresponds to d 10) when comparing H cut = 0.7 to H cut = 0.6. On the contrary, the distances between the H cut = 0.5 fit and the H cut = 0.6 fit are larger, especially for charm and strangeness at x 10 −3 . This comparison indicates that there is no real benefit in loosening the cut from H cut = 0.6 to H cut = 0.7 (since differences at the PDF level are small, and the possibility of biasing the fit higher) whilst it is indeed advantageous to use H cut = 0.6 rather than the tighter cut H cut = 0.5, thanks to the increase in PDF constraints provided by the additional data.
Finally, in Fig. 18 we show the relative PDF uncertainties in the NNPDF3.1sx NNLO fits with the three different values of the H cut cut on the hadronic data. For completeness, we also include in this comparison the results of the NNPDF3.1 fit. Specifically, we show the gluon, the quark singlet, the anti-up quark, and the total strangeness, at Q = 100 GeV. From the comparison we see that, as expected, the smaller the value of H cut , the more marked the increase in PDF uncertainties. On the other hand, we see that for H cut = 0.6 the results are already competitive with those of NNPDF3.1. We also find that in the small-x region, PDF uncertainties are smaller in the NNPDF3.1sx fits than in the NNPDF3.1 ones, Summarizing, we have provided here a number of indications that the NNPDF3.1sx fit with H cut = 0.6 is not biased by hadronic data sensitive to small-x resummation, and at the same time we have demonstrated that the resulting PDF uncertainties are competitive, though still larger, with those of NNPDF3.1. These considerations provide further weight for our default choice of the H cut cut to the hadronic data.

The role of the Q 2 = 2.7 GeV 2 bin
We have stressed in Sect. 3 that, as opposed to NNPDF3.1, we include in the NNPDF3.1sx fits an additional low Q 2 bin of the inclusive HERA dataset, specifically the one with Q 2 = 2.7 GeV 2 . This choice has the important advantage of extending the kinematic coverage of the fits from x min 5 × 10 −5 down to x min 3 × 10 −5 . The main reason why this bin was excluded from previous NNPDF fits (as well as in most other global PDF fits) is its low value of Q 2 , which lies at the boundary between perturbative and nonperturbative dynamics, and where fixed-order perturbation theory might not be fully appropriate. Here we show that this failure is not actually due to non-perturbative dynamics, but rather it represents a limitation of the fixed-order expansion in the small-x region enhanced by the larger value of α s . Indeed, we find that once NNLO fixed-order perturbation theory is supplemented by NLLx resummation, this bin can be described with similar quality as the rest of the HERA data.
To illustrate this point, we have computed the values of χ 2 /N dat for the N dat = 17 data points that constitute the Q 2 = 2.7 GeV 2 bin of the inclusive HERA structure function dataset. We find that the values of χ 2 /N dat for this bin are 1.64 and 1.34 in the NNPDF3.1sx NLO+NLLx and NNLO+NLLx fits. These results can be compared with the corresponding values in the NLO and NNLO fits, which turn out to be 2.04 and 3.04, respectively. The trend is the same as that for the total NC HERA inclusive dataset (see Table 4), namely with the NNLO+NLLx (NNLO) fit leading to the best (worst) overall description, and with the NLO and NLO+NLLx values in between. Interestingly, we also see that for this specific fit NNLO+NLLx theory leads to a rather better χ 2 than the NLO+NLLx one, although the small number of data points prevents drawing any strong conclusion from this observation. Fig. 18 The relative PDF uncertainties in the NNPDF3.1sx NNLO fits with the three different values of the H cut cut on the hadronic data, compared with those from NNPDF3.1. We show the gluon, the quark singlet, the anti-up quark, and the total strangeness, at Q = 100 GeV Once we have established that the fit quality to the Q 2 = 2.7 GeV 2 HERA bin is satisfactory when NLLx resummation is included, we can next turn to study the constraints that this bin has on the small-x PDFs. With this motivation, we have performed a global fit at NNLO+NLLx with the same settings as the NNPDF3.1sx baseline but raising the low Q 2 kinematic cut from Q 2 min = 2.69 GeV 2 to Q 2 min = 3.49 GeV 2 , as in NNPDF3.1, so that the HERA bin with Q 2 = 2.7 GeV 2 is excluded. In the latter case, the lowest HERA bin included is the one with Q 2 = 3.5 GeV 2 .
The results are shown in Fig. 19, where the gluon and the quark singlet PDFs obtained in the NNPDF3.1sx NNLO+NLLx fits with and without this additional bin are compared at the input parametrization scale of Q 0 = 1.64 GeV, together with their relative PDF uncertainties. We find that the inclusion of this extra Q 2 bin leads to a significant reduction of small-x uncertainty of the gluon in the region which is constrained by the data (x 10 −5 ), while the quark singlet is essentially unaffected. These results illustrate how the use of an improved theory, NNLO+NLLx in this case, can lead indirectly to a decrease of the PDF uncer-tainties due to the possibility of including more data in the fits from a wider kinematic range.

Small-x resummation and HERA structure functions
The results of the previous section provided two main pieces of information. First of all, the inclusion of small-x resummation improves the description of those datasets which represent the best probe of the small-x region, namely the inclusive and charm HERA structure functions. Second, the impact of resummation at the level of PDFs can be sizable. In this section, we focus on the HERA data in the small-x and small-Q 2 region, in order to further quantify the improvement in its description when fixed-order theory is supplemented by NLLx resummation.
We first compare the HERA structure functions at low x with fixed-order and resummed theoretical predictions, both for the inclusive and charm reduced cross-sections as well as for the longitudinal structure function F L . In all cases, we highlight the improved description that is achieved once NNLO+NLLx theory is used in all cases. To quan- titatively investigate the evidence for the onset of small-x resummation in the HERA data, we introduce several estimators building upon the set of diagnostic tools first presented in Refs. [65,66]. We finally study how removing HERA data at low-x and low-Q 2 affects global NNLO fits, and we discuss how the resulting PDFs are modified at mediumand large-x. This way, it is possible to assess whether the inclusion of data poorly described in a fixed-order analysis might be a source of bias for high-Q 2 phenomenology.

The HERA data in the small-x region
In order to investigate in greater detail how well resummed theory describes the low-Q 2 HERA cross-sections, we first perform a comparison of the theoretical predictions obtained using the results of the NNPDF3.1sx NNLO and NNLO+NLLx global fits to the experimental data. To begin with, in Fig. 20 we show the neutral-current (NC) reduced cross-section, defined as where Y + = 1 + (1 − y) 2 and y = Q 2 sx is the inelasticity. This comparison is performed for the first four bins in Q 2 above our Q 2 min kinematic cut of the √ s = 920 GeV dataset, corresponding to Q 2 = 2.7, 3.5, 4.5 and 6.5 GeV 2 respectively. In the left plots, the uncertainty of the experimental data points is given by the sum in quadrature of the various sources of uncorrelated and correlated uncertainties, whereas the theoretical predictions include the associated PDF uncertainty. In the right plots, instead, only the uncorrelated uncertainties are shown in the data, and the correlations are taken into account via shifts which modify the theoretical prediction [187] and facilitate the graphical comparison. Note that these correlations are included in the χ 2 definition. However, unlike in a Hessian approach, in a Monte Carlo method one does not determine the best-fit systematic shifts. Rather, here we have computed them a posteriori, under the assumption Fig. 20 Comparison between the HERA NC reduced cross-section from the √ s = 920 GeV dataset and the results of the NNLO and NNLO+NLLx fits with the corresponding PDF uncertainties. We show the results for the first four bins in Q 2 above the Q 2 min kinematic cut.
For each bin we also show in the bottom panel the ratio of the theory predictions to the experimental data. The plots on the right show the theoretical prediction including the shifts as discussed in the text that the uncertainties are gaussian, which is not necessarily true in a Monte Carlo fit. Therefore, this comparison must be interpreted with care. From this comparison, we see that for x 5 × 10 −4 the results of the NNLO and NNLO+NLLx fits are essentially identical; in both cases, the theoretical predictions undershoot the data. The trend changes for values of x smaller than 5 × 10 −4 , where the NNLO and the NNLO+NLLx predictions start to differ. Around this value, we observe that the reduced cross-section exhibits a slope change too: the data stop rising and, after a turnover, the reduced cross-section starts decreasing. As a result, the NNLO prediction starts to overshoot the data, whereas the NNLO+NLLx prediction is in reasonable agreement with the data for x 10 −4 . It is worth observing that the differences between the NNLO and NNLO+NLLx predictions are relatively small and concern only a limited number of points. By looking at the bottom panels in Fig. 20, where we show the ratio to the experimental data, we see that the two predictions differ by at most 10% and only for the smallest values of x. Yet the combined HERA dataset is so precise that the improvement in the description provided by small-x resummation is clearly visible at the χ 2 level, as was shown in Table 4, and will be discussed further below in Sect. 5.2.
The improved description of the inclusive reduced crosssection data at small x can be in part traced back to the role of the longitudinal structure function F L (x, Q). As reviewed in Sect. 2, F L is particularly sensitive to the effects of smallx resummation, and in particular to deviations from the DGLAP framework. The reason is that it vanishes at the Born level, and therefore it receives gluon-initiated contributions already at its first non-trivial order. As shown in Fig. 5, the differences between the NNLO and NNLO+NLLx can be as large as ∼ 30% at the lowest x and Q 2 bins for which there are data available. As a consequence, at small x and small Q 2 the contribution of F L to σ r,NC can be significant, see Eq. (5.1), thus partly explaining the differences between the NNLO and NNLO+NLLx predictions observed in Fig. 20. Therefore, it is useful to compare the predictions also for the longitudinal structure function F L in the NNLO and NNLO+NLLx fits.
In Fig. 21 we compare the latest measurements of F L from the H1 collaboration [188] 5 with the predictions from the NNPDF3.1sx NNLO and NNLO+NLLx fits. Note that our fits already include the constraints from F L , not directly but rather via its contribution to the NC reduced cross-section, Eq. (5.1). In this comparison, the experimental uncertainties have been added in quadrature, and each value of Q 2 corresponds to a different x bin as indicated in the plot. The Fig. 21 The longitudinal structure function F L (x, Q 2 ) as a function of Q 2 for different x bins for the most recent H1 measurement [188], comparing the results of the NNLO and NNLO+NLLx fits NNPDF3.1sx results are shown down to the smallest scale for which one can reliably compute a prediction, 6 which is set by the initial parametrization scale Q 2 0 = 2.69 GeV 2 . We see that for Q 2 100 GeV 2 there are significant differences between the NNLO+NLLx and the NNLO predictions, which can be traced back to a combination of the corresponding differences for the input small-x gluon and those in the splitting and coefficient functions (see Fig. 5). The NNLO+NLLx result is larger than the NNLO result by a significant amount: at Q 2 10 GeV 2 , the resummed calculation is more than a factor 2 larger than the NNLO result. Moreover, while at NNLO F L starts becoming negative at small x and Q 2 (below the scale where the positivity constraints are imposed in the NNPDF fits) the NNLO+NLLx result instead exhibits a flat behavior even for the smallest values of Q 2 . The larger value of F L with the NNLO+NLLx theory leads to a lower reduced cross-section at high y, with a more pronounced turnover, thus giving a better description of σ r,NC at small x, as shown in Fig. 20.
Finally, in Fig. 22 we show a similar comparison to that of Fig. 20, this time for the HERA charm-production reduced cross-sections. Here we also show the two Q 2 bins about the lower Q 2 min cut, which in this case correspond to the Q 2 = 5 and 7 GeV 2 bins. We find that especially for the bin with Q 2 = 5 GeV 2 , the NNLO+NNLx prediction agrees well with the HERA data while the NNLO one overshoots it. We remind the reader again that these graphical comparisons do not take into account the correlations between systematic uncertainties. The large difference between the χ 2 at NNLO and at NNLO+NLLx is therefore only partially captured by Fig. 22. As we shall see in greater detail in Sect. 5.2, also in this case the deterioration of the NNLO χ 2 with respect  Table 4 stems mostly from the low-Q 2 , low-x bins.
Note that the HERA charm cross-sections are extracted from the experimentally measured fiducial cross-section [133] by extrapolation to the full phase space using the fixed-order O(α 2 s ) calculation of the HVQDIS program [190], based on the fixed-flavor number scheme. This should be contrasted with the inclusive neutral-current structure function measurements, which are determined from the outgoing lepton kinematics and therefore do not assume any theory input. Given that we have shown that fixed-order and resummed predictions for F c 2 can exhibit important differences at small x, such theory-based extrapolation based on the O(α 2 s ) fixedorder calculation might introduce a bias whose size is difficult to quantify. It is quite possible that a more consistent analysis of the raw data based instead on an extrapolation using resummed theoretical predictions might further improve the already good agreement of the extracted charm cross-section with the NNLO+NLLx fit.

Quantifying the onset of small-x resummation in the HERA data
In this section we resort to a number of statistical estimators to identify more precisely the onset of small-x resummation in the inclusive and charm HERA measurements. First, we perform a detailed χ 2 analysis, which we then complement by a study of the pulls between theory and HERA data.

χ 2 analysis
The χ 2 /N dat values summarized in Table 4 indicate that the fit quality of the inclusive HERA structure functions improves when resummation effects are included: this is particularly true at NNLO, where the total χ 2 drops by χ 2 = −121 units in the NNLO+NLLx fit. We now want to identify the origin of this improvement, and investigate to what extent it arises from a better description of the data in the small-x and small-Q 2 region where the effects of small-x resummation are expected to be most important.
To achieve this goal, we have recomputed the χ 2 /N dat values of the HERA inclusive and charm cross-sections using the NNPDF3.1sx NLO, NNLO, NLO+NLLx, and NNLO+NLLx global fits with the default choice H cut = 0.6, excluding those data points for which The condition Eq. (5.2) is designed to exclude data for which the small-x logarithmic terms are expected to be of the same size at all orders in the coupling α s , thus potentially spoil- ing the perturbative behavior of the theoretical predictions at fixed order. From basic considerations (see also Sect. 3.2), one would expect fixed-order perturbation theory to break down for α s (Q 2 ) ln 1 x of order 1. The parameter D cut should thus be of order 1 as well. By varying the value of D cut , we can vary the number of data points excluded from the computation of the χ 2 /N dat . For sufficiently small values of D cut , all contributions which potentially spoil perturbation theory should be cut away, and we should thus find that small-x resummation does not improve the quality of the fit. Then as we increase D cut , more data points at small x and Q 2 will be included, and the effects of the resummation should become apparent. A kinematic plot showing the HERA structure function data which are cut for various values of D cut is shown in the left panel of Fig. 23. We emphasize that this cut should not be confused with the H cut cut defined in Eq. (3.1), which was used to determine which hadronic data enter in the fit; here the parameter D cut applies only to DIS structure functions and is used as an a posteriori diagnosis tool after the fit has been performed.
In Fig. 24 we display the values of χ 2 /N dat for the HERA neutral-current inclusive (top left) and charm (bottom left) reduced cross-sections as a function of D cut . First of all, we observe that at NNLO the χ 2 /N dat increases sharply for D cut 2, or, equivalently, as more data from the small-x and small-Q 2 region are included, both for the inclusive and the charm data. On the other hand, this trend disappears for the NNLO+NLLx fits: in this case the value of χ 2 /N dat is flat for all D cut values in the studied range.
Another interesting feature of these plots is that the stability with respect to the value of D cut is also present for the NLO and NLO+NLLx fits. Indeed, the χ 2 /N dat values for the NLO, NLO+NLLx, and NNLO+NLLx fits all exhibit a rather similar shape. This is of course a consequence of the fact that, as shown in Sect. 4, the PDFs obtained from the fits using these three theories are rather close to each other, whereas the NNLO PDFs are very different at small x. Remarkably, for the inclusive data especially the NNLO+NLLx fits lead to a better χ 2 /N dat than the NLO and NLO+NLLx ones, presumably due to the additional NNLO corrections included in the NNLO+NLLx matched calculations. This result highlights the importance of the NNLO corrections for the optimal description of the medium and large-x HERA data.
The results of Fig. 24 demonstrate that fixed-order NNLO theory does not provide a satisfactory description of either the inclusive or charm DIS data at small x and small Q 2 . The better description is instead achieved by including NLLx effects, providing direct evidence of the need for small-x resummation at small x. Moreover, we observe that the rise in the χ 2 /N dat values of the NNLO fits becomes very significant for D cut 2. This means that BFKL effects at NNLO approximately start to become important when To study whether the treatment of the hadronic data in the PDF fits can modify this conclusion, in the upper right panel of Fig. 24 we also compare the χ 2 /N dat values as a function of D cut for the NNPDF3.1sx NNLO+NLLx global fits with the three H cut values discussed in Sect. 4.2.3, namely H cut = 0.5, 0.6 and 0.7, as well as with the global H cut = 0.6 NNLO fit and the NNLO+NLLx DIS-only fit. These comparison illustrate that our quantitative conclusions are to a very good approximation independent of the specific cut applied to the hadronic data: very similar NNLO+NLLx results are found in the global fit irrespective of the value of H cut , as well as for the corresponding DIS-only fit. We have also verified that the same conclusion holds for the NLO and NLO+NLLx fits.
In Refs. [65,66], a similar cutting exercise was performed, but in that case the specific form of the cut to the smallx and small-Q 2 data was inspired by saturation arguments. Specifically, the condition used to exclude data points was with λ = 0.3. The value of A cut determines how stringent is the cut: the larger its value, the more data points excluded (so 1/A cut behaves qualitatively in the same way as D cut ). While the inspiration for the cut Eq. (5.4) is different from that of Eq. (5.2) (which is based instead on perturbative considerations), the practical result is the same, with only some differences on the exact shape of the cut in the (x, Q 2 ) plane (see the right panel of Fig. 23). The results for the χ 2 /N dat as a function of 1/A cut are shown in the bottom right panel of Fig. 24, and indeed confirm that the trend is essentially the same, irrespectively of the specific details of how the small-x and Q 2 data are cut.
In summary, the results collected in Fig. 24 clearly demonstrate the onset of BFKL dynamics in the small-x and Q 2 region for both the inclusive and charm HERA data. Specifically, we find that the use of NNLO+NLLx theory gives the best description of the HERA data in the small-x region, while NNLO theory gives a significantly worse description. Moreover, our results also allow us to determine the kinematic region where small-x resummation effects start to become phenomenologically relevant, thus providing useful guidance to estimate their reach at the LHC as well as for future colliders.

Pull analysis
A complementary approach to further investigating the onset of BFKL dynamics in the low-x region, and to make connection with the analysis of Refs. [65,66], is provided by the calculation of the relative pull between experimental data and theory. This relative pull is defined as where the normalization is given by the average of central values. This estimator allows us to quantify the absolute size of the differences between data and theory in units of the crosssection itself. Here we focus on the results computed with NNLO and NNLO+NLLx theory, using the NNPDF3.1sx sets obtained in the respective global fits with the default cut H cut = 0.6.
To visualize the differences between data and theory in the small-x and small-Q 2 region, we can represent the relative pull P rel i (x, Q 2 ), Eq. (5.5), as a function of (x, Q 2 ) in the relevant region of the kinematic plane. In Fig. 25 we show an interpolated representation of P rel (x, Q 2 ) for the HERA neutral-current dataset at √ s = 920 GeV and the NNLO and NNLO+NLLx fits. In the case of the NNLO fit, the relative differences between theory and data can be up to ∼ 20% at small-x and Q 2 , and reduce to less than a few percent at larger x or Q 2 . On the other hand, the agreement between data and theory is markedly improved in the case of the NNLO+NLLx fit: the quality of the data description is essentially the same everywhere in the region considered, and the relative differences between data and theory are everywhere below the 8% level. These plots show that by using NNLO+NLLx theory, one can achieve a satisfactory description of the inclusive HERA measurements in the entire region spanned by the available data.
In order to further quantify differences and similarities between the NNLO and NNLO+NLLx theoretical predictions, in Fig. 26 we show a similar relative pull as in Eq. (5.5), now between the theoretical predictions for the HERA reduced cross-sections obtained with the NNLO and the NNLO+NLLx theory and fits, namely Note that in this comparison both the color code and the (x, Q 2 ) ranges are different from those of Fig. 25. From  Fig. 26 Same as Fig. 25, now for the relative difference in the theoretical predictions of the HERA reduced cross-sections between the NNLO and NNLO+NLLx fits, Eq. (5.6). Note the different color code and x and Q 2 ranges with respect to Fig. 25 the results of Fig. 26 we see that the differences between the cross-sections computed with NNLO and NNLO+NLLx theory are between 5% and 10% for Q 2 10 GeV 2 and x 2 × 10 −4 . Once we move away from this region, differences become smaller. For x 2 × 10 −4 , we find that the differences are always smaller than ∼ 3%, for any value of Q 2 . This comparison provides a detailed snapshot of the region in the (x, Q 2 ) plane where the impact of resummation is phenomenologically more relevant, and is consistent with the results shown in Fig. 20 and the conclusions of the χ 2 profile analysis of Sect. 5.2.1.

Sensitivity to subleading logarithms
Finally, we can study if our conclusions are stable with respect to variations of unknown subleading logarithms. To this end, in the left panel of Fig. 27 we show the values of χ 2 /N dat for the HERA NC inclusive reduced cross-section as a function of D cut , now comparing the DIS-only fit at NNLO and NNLO+NLLx to the NNLO+NLLx fit where subleading logarithms are introduced as described in Sect. 4.1. We observe that the two NNLO+NLLx profiles are very similar, with the χ 2 /N dat of the alternative fit being slightly lower at larger D cut . To better quantify the differences between the two variants, in the right panel of Fig. 27 we also show the relative pull Eq. (5.6) between the theoretical predictions for the two NNLO+NLLx fits for the HERA neutral current √ s = 920 GeV dataset. We observe that the relative difference is at most 2% for the smallest values of x and Q 2 probed by the dataset, and is below 0.5% for all x 3×10 −4 , independently of the value of Q 2 . This analysis shows that our results are stable with respect to variations of subleading logarithms.

Impact of the small-x HERA data on PDFs at medium and large-x
In the last part of this section, we present results of additional NNPDF3.1sx NNLO fits where we have removed a number of HERA structure function data points in the small-x and Q 2 region, in order to study how the resulting PDFs are affected by the use of such reduced dataset. This exercise allows us to understand to what extent existing NNLO global PDF fits might be biased by fitting low-x data while neglecting the  effects of small-x resummation. Since we have just demonstrated that at small x the HERA structure functions prefer NNLO+NLLx theory to the NNLO one, it may be advisable to apply dedicated kinematic cuts in the small-x and Q 2 region in standard NNLO analyses. This would ensure on one hand that the fitting dataset corresponds to a region where a fixed-order perturbative expansion is reliable, and on the other hand that the estimate of the uncertainties at small x is more reliable. For this purpose, we have performed variants of the NNPDF3.1sx NNLO global fit without any cut on the hadronic data (that is, H cut = ∞) but where instead we impose the cut Eq. (5.2) to the DIS structure function data before fitting, thus reducing the number of data points in the small-x and small-Q 2 region. Specifically, we have performed NNLO fits with D cut = 1.7, 2.0 and 2.3, as a well as a fit without cutting any data (D cut = ∞) as a reference. The motivation for this range of D cut values is the observation (see Fig. 24) that D cut 2 indicates the region where the effects of small-x resummation start to become significant.
The comparison between the NNPDF3.1sx NNLO fits with different cuts to the DIS data is shown in the upper plots of Fig. 28. Specifically, we show the comparison between the gluon and quark singlet from the fits with various values of D cut with the corresponding fit without that cut (D cut = ∞) at Q = 100 GeV. From this comparison we can see that -as expected -the higher the value of D cut , the smaller are the PDF uncertainties at small x due the increase in kinematic coverage of the fitted HERA data. However, the central values remain very stable, even at the lowest values of x. Additionally, we also see that for x 5 × 10 −4 the gluon and quark singlet are extremely stable with respect to the D cut variations, both in terms of central value and of PDF uncertainties. Therefore, we conclude that current NNLO fits are not biased in the region relevant for precision LHC phenomenology, even if the fits include points at small x, while neglecting resummation effects.
It is also interesting to compare the NNPDF3.1sx NNLO fit with D cut = 1.7 with our default NNLO+NLLx global fit, namely the one with the cut in hadronic data corresponding to H cut = 0.6, but D cut = ∞. This comparison allows us to understand if the PDF uncertainties of the conservative NNLO fits, where data points at small x and small Q 2 have been removed, account for the PDF shift induced by using the more accurate NNLO+NLLx theory. We show this comparison in the bottom plots of Fig. 28 at the scale Q = 100 GeV. Whereas for the quark singlet the shift between the NNLO+NLLx and NNLO fits is mostly covered by the corresponding PDF uncertainties, for the gluon the shift in central values is bigger and is not covered by the PDF uncertainties, despite the larger PDF errors of the fit with the conservative dataset.
The results of Fig. 28 suggest that at small x the theoretical uncertainties associated with the NNLO gluon are comparable to or larger than the PDF uncertainties. Moreover, we observe that the shift induced by NNLO+NLLx theory is covered by the PDF uncertainties only for values of x larger than x 3 × 10 −3 . Therefore, for processes sensitive to the small-x region, including a number of LHC cross-sections, current NNLO PDF uncertainties do not fully account for the total theoretical uncertainty, suggesting that the use of NNLO+NLLx theory would lead to more reliable theoretical predictions.

Small-x phenomenology at the LHC and beyond
In this section we explore some of the phenomenological implications of the NNPDF3.1sx fits. First of all, we present a first assessment of the possible impact of NLLx resummation at the LHC. We then move to DIS-like processes, for which we can produce fully consistent NNLO+NLLx predictions. In this context, we consider the implications of the NNPDF3.1sx fits for the ultra-high-energy (UHE) neutrino-nucleus cross-sections as well as for future highenergy lepton-proton colliders such as the LHeC [191] and the FCC-eh [192,193], illustrating the key role that small x resummation could play in shaping their physics program.

Small-x resummation at the LHC
In this section we perform a first exploration of the potential effects of small-x resummation on precision LHC phenomenology. We start by considering parton luminosities and then we estimate the effects of small-x resummation for electroweak gauge boson production at the LHC. The latter study will, however, be necessarily only qualitative, since as explained in Sect. 2 the relevant small-x resummed partonic cross-sections are not yet implemented in a format amenable for phenomenological applications. The studies of this subsection will thus be complementary to previous estimates of the effects of small-x resummation on inclusive LHC processes in which both evolution and cross-section were resummed, but the PDFs used were fixed rather than refitted [44,185].

Parton luminosities
In order to provide a first insight on the possible impact of NLLx resummation effects on hadronic cross-sections, it is useful to consider its effects on the parton luminosities. We consider both the integrated parton luminosity Eq. (2.3) and also luminosities differential in rapidity (see e.g. [194]) We assume the production of a hypothetical final state with invariant mass M X , so that x = M 2 X /s with √ s being the LHC center-of-mass energy, and we take the factorization scale to be μ 2 = M 2 X . Despite offering only a qualitative estimate of the effects of small-x resummation, parton luminosities contain the bulk of the information as regards the partonic contributions to a given process. In particular, rapidity-dependent PDF luminosities provide a direct mapping between regions in the (x, Q 2 ) plane (PDF sensitivity) and those in the (M X , y) plane (kinematic coverage for collider production), assuming leading-order production kinematics.
Let us start with the integrated parton luminosities, Eq. (2.3). In Fig. 29 we show the gluon-gluon, quarkgluon, quark-antiquark and quark-quark PDF luminosities at √ s = 13 TeV as a function of the invariant mass M X , comparing the NNPDF3.1sx global fits based on NNLO and NNLO+NLLx theory respectively.
For the two gluon-initiated luminosities, the effects are very large for M X 100 GeV, and smaller above that value.
For the gg luminosity, for instance, the ratio between the NNLO and NNLO+NLLx results can be more than ∼ 30% for M X 10 GeV, a region relevant for instance for open B-meson production.
Even larger effects can be expected for processes at smaller invariant masses, such as D-meson or J/ production.
At M X 100 GeV, a region relevant for e.g. top-quark pair production, the gluon-induced luminosities are instead reduced, albeit by only a few per cent.
In the case of the quark-antiquark and quark-quark luminosities, the differences due to resummation are less significant, with the NNLO and NNLO+NLLx luminosities in agreement within PDF uncertainties for the entire M X range.
However, the effects of small-x resummation are not negligible; for instance, they could still be as large as 10% at M X = 10 GeV, a region probed by the LHC in processes such as low-mass Drell-Yan production. Thus, in this region the effects are small, but nevertheless of the same order as the experimental uncertainties of recent high-precision LHC measurements, such as for instance the ATLAS 2011 W, Z rapidity distributions [152] or the CMS Z p T distributions [169].
It is important to emphasize here that the luminosity comparison in Fig. 29 provides only a rough estimate of the actual differences between the NNLO and the fully resummed NNLO+NLLx cross-sections, since a quantitative assessment requires the resummation of the partonic cross-sections for the relevant processes, and this can be as large as the difference in the luminosities [44,185].
This said, the results of Fig. 29 show that the effects of small-x resummation are potentially significant for LHC cross-sections, in particular for those with large gluoninitiated contributions.
Next, in Fig. 30 we show same comparison but this time between the NLO and NLO+NLLx fits. As discussed in Sect. 4, we expect the differences to be more moderate com-pared to the NNLO fits case. Indeed, the differences are now much smaller, both for the gluon-initiated and for the quarkinitiated luminosities. The most significant effect of resummation can again be seen in the gg luminosity, but now only at the 10% level at M X 10 GeV. The other luminosities all agree within uncertainties. Henceforth, we will focus on the comparison between the NNLO+NLLx and the NNLO fits, as in all cases the corresponding differences between NLO+NLLx and NLO would always be much smaller.
Now we move to compare PDF luminosities which are differential in rapidity, Eq. (6.1). As already mentioned, these luminosities allow for a more direct mapping between the final-state kinematics and the regions of x, Q 2 of the underlying PDFs. For simplicity, we focus here on the gluon-gluon and quark-antiquark luminosities, as the behavior of the gluon-quark and quark-quark is closely related to these two. In Fig. 31 we compare the PDF luminosities of the NNLO and NNLO+NLLx fits, normalized to the central value of the former.
We show the results as a function of y for three different values of M X , namely 10, 30, and 100 GeV.  Fig. 31 we see that the impact of small-x resummation depends on the final-state rapidity y, and it increases close to the kinematic endpoints. This is expected as large (or small) values of the rapidity probe smallx values in one of the two partons that make up the parton luminosity, see Eq. (6.1). For instance, in the case of the gg luminosity, for the production of a final state with invariant mass M X = 10 (30) GeV, the ratio between NNLO+NLLx and NNLO is between 30 and 50% (10 and 20%), depending on the specific value of the rapidity. The differences are smaller in the case of the quark-antiquark luminosities, though we note that they could become more relevant if the PDF uncertainties were reduced by including in the fit data sensitive to the small-x region of the PDFs, such as the LHCb W, Z forward production cross-sections. This would, however, require the inclusion of small-x resummation in the partonic cross-sections for the relevant processes.

Implications for Drell-Yan production
We now present a first exploration of the possible phenomenological consequences of small-x resummation for LHC cross-sections, specifically for the Drell-Yan production process. We do this by providing estimates for some recent Drell-Yan cross-section measurements from the LHC, focusing on those more sensitive to the possible presence of small-x effects, and comparing the results of the predictions from the NNPDF3.1sx NNLO and NNLO+NLLx fits, using in both cases the fixed-order NNLO hard-scattering cross-sections. These differences likely over-estimate the real effect, and in particular might be reduced once the resummation in the partonic cross-sections is taken into account [185]. However, we believe they provide a reliable though conservative estimate of the possible size of the resummation effects that can be expected.
Specifically, we show in Fig. 32 the predictions for the low-mass DY cross-sections from ATLAS at 7 TeV [156], the lowest invariant mass bin for the CMS Drell-Yan crosssections double-differential in y and M ll at 8 TeV [195], as well as for forward W + and Z production at 8 TeV from LHCb [174]. Note that none of these datasets is included in the NNPDF3.1sx fits, since as discussed in Sect. 3 they are removed by the H cut cut, Eq. (3.1). We stress once again that we calculate the NNLO+NLLx and NNLO cross-sections using the corresponding PDFs from the NNPDF3.1sx fits, but using in both cases the fixed-order NNLO coefficient functions, using the same settings described in Sect. 3. We do not show the experimental data points in this comparison, as our aim is to focus on the impact of the resummation rather than a comparison with the measured cross-sections.
From the results shown in Fig. 32, we find that the NNLO and NNLO+NLLx predictions are consistent within uncertainties in almost all cases. The differences are more marked for the kinematic regions directly sensitive to small-x, such as small M ll for the ATLAS data and large rapidities in the case of the LHCb and CMS measurements. In the latter case, Fig. 32 Comparison between the NNPDF3.1sx NNLO and NNLO+NLLx predictions for selected Drell-Yan measurements at the LHC. From left to right and from top to bottom, we show the ATLAS low-mass measurements at 7 TeV, the CMS low-mass measurements at 8 TeV, and the LHCb W + and Z rapidity distributions at 8 TeV. For the NNLO+NLLx predictions, the effects of small-x resummation are included in the PDF evolution but not in the partonic cross-sections the shift due to NNLO+NLLx theory could be as large as ∼ 5% at the largest rapidities, and the two PDF bands do not overlap for y > 2 in the invariant mass bin considered.
Moreover, since the experimental uncertainties for the cross-sections shown in Fig. 32 can be smaller than the corresponding PDF errors (and in some cases also smaller than the shift between the NNLO and NNLO+NLLx curves), we can conclude from this exercise that the inclusion of these data into a fully consistent small-x resummed global PDF fit might provide further evidence for BFKL dynamics, this time from high-precision electroweak LHC cross-sections as opposed to from lepton-proton deep-inelastic scattering.

The ultra-high-energy neutrino-nucleus cross-section
We next briefly explore the implication of the NNPDF3.1sx fits for the calculation of the total neutrino-nucleus crosssections at ultra-high-energies (UHE). The interpretation of available and future UHE data from neutrino telescopes, such as IceCube [196] and KM3NET [197], requires precision predictions for the UHE cross-sections. With this motivation, a number of phenomenological studies of the UHE cross-sections and the associated uncertainties has been presented, both in the framework of collinear DGLAP factorization [198][199][200][201][202][203][204] and beyond it [205][206][207][208][209], the latter including for instance the effects of non-linear evolution or saturation.
Here we focus on the charged-current (CC) neutrinonucleus inclusive cross-sections. Measuring neutrinonucleus interactions at the highest values of E ν accessible at neutrino telescopes explores values of x down to ∼ 10 −9 for Q ∼ M W , thus representing a unique testing ground of smallx QCD dynamics. We have computed the theoretical predictions with APFEL+HELL for NNLO and NNLO+NLLx theory, using the corresponding NNPDF3.1sx fits as input. Heavy-quark mass effects are included using the FONLL scheme, although these mass corrections are negligible at the relevant intermediate and high neutrino energies, so the calculation is effectively a massless one. Fig. 33 The UHE neutrino-nucleus charged-current cross-section σ CC (E ν ) as a function of the neutrino energy E ν , comparing the results obtained using the NNPDF3.1sx NNLO fits with those of the its resummed NNLO+NLLx counterpart In Fig. 33 we show the UHE neutrino-nucleus chargedcurrent cross-section σ CC (E ν ) as a function of the neutrino energy E ν for the fixed-order and for the resummed predictions. We show both the absolute cross-sections and the crosssections normalized to the central value of the NNLO prediction. The error bands indicate the one-sigma PDF uncertainties. The upper limit in E ν corresponds to the foreseeable range of the current generation of neutrino telescopes.
As we can see from the comparison of Fig. 33, the main effect of small-x resummation is to increase the cross-section at the highest energies, by an amount that can be as large as 50% or more. The PDF errors are, however, large, and the NNLO and NNLO+NLLx predictions agree at the onesigma level on the whole range of energy considered. Given the large PDF uncertainties, it appears difficult to tell apart distinctive BFKL signatures in the total UHE inclusive crosssection. However, it is interesting to note that the effect of small-x resummation on the UHE cross-sections is the opposite of that obtained in calculations based on non-linear QCD dynamics, which instead predict a smaller cross-section at high energy (see e.g. [205]).
A promising strategy towards reducing the large PDF errors that affect σ CC (E ν ) in Fig. 33 is provided by the inclusion of charm-production data from LHCb [210][211][212] in the PDF fits. As demonstrated in [14,213,214], the inclusion of LHCb D-meson production cross-sections gives a significant reduction in PDF uncertainties in the small-x region (up to an order of magnitude at x 10 −6 ), which in turns leads to UHE cross-sections with few-percent theory errors up to E ν = 10 12 GeV [213]. In this respect, the combination of NLLx resummation and the additional constraints provided by the LHCb charm data would make possible a calculation of the UHE cross-sections with unprecedented theoretical and experimental uncertainties.

Small-x resummation at future electron-hadron colliders
Since we have demonstrated the onset of BFKL dynamics in the inclusive HERA structure function data, it is natural to expect that the effects of small-x resummation will become even more relevant at the proposed future highenergy electron-hadron colliders: the higher their center-ofmass energy √ s, the smaller the values of x kinematically accessible in the perturbative region of Q 2 .
One such future ep collider is the large hadron-electron collider (LHeC) [191,215]. In its latest design, a proton beam from the LHC with E p = 7 TeV would collide with an electron/positron beam with E e = 60 GeV coming from a new LinAc, thus enabling to access the region down to x min 2 × 10 −6 at Q 2 = 2 GeV 2 . A more extreme incarnation of the same idea corresponds to colliding the same E e = 60 GeV electrons with the E p = 50 TeV beam of the proposed 100 TeV Future Circular Collider (FCC) [192,193]. The resulting collider, dubbed FCC-eh, would be able to reach x min 2 × 10 −7 at Q 2 = 2 GeV 2 . We show in Fig. 34 the kinematic coverage in the (x, Q 2 ) plane of the two machines, compared with that of the HERA structure function data included in the NNPDF3.1sx fits. It is clear that these two machines would probe into the small-x region much deeper than HERA, thus allowing an unprecedented exploration of new QCD dynamics beyond fixed-order collinear DGLAP framework.
In the following, we perform an initial exploration of the potential of the LHeC/FCC-eh for small-x studies. We use APFEL in conjunction with HELL to produce NNLO and NNLO+NLLx predictions for various DIS structure functions, assuming the latest version of the simulated LHeC/FCC-eh kinematics. 7 In Fig. 35 we show these predictions for the F 2 and F L structure functions using the NNPDF3.1sx NNLO and NNLO+NLLx fits at Q 2 = 5 GeV 2 for the kinematics of the LHeC and the FCC-eh. In the case of F 2 , we also show the expected total experimental uncertainties based on the simulated pseudo-data, assuming the NNLO+NLLx curve as central prediction. To compare with the kinematic region within the reach of HERA data, we also show in the inset of the left plot the values of F 2 in the range x > 3 × 10 −5 . The total uncertainties of the simulated pseudo-data are at the few-percent level at most; hence they are much smaller than the PDF uncertainties in most of the kinematic range. No simulated pseudo-data is currently available for F L using the latest scenarios for the two colliders; thus in this case we show only the theoretical predictions.
We now discuss in turn some of the interesting features in Fig. 35. First of all, we clearly see how with the FCC-eh one can probe the small-x region deeper than the LHeC by about an order of magnitude. Second, we find that the differences between NNLO and NNLO+NLLx are moderate for F 2 , especially if we take into account the large PDF uncertainties. The difference between the central values is in fact at the 15% level at x 10 −6 , but the current PDF uncertainties are much larger. However, given the precision that the data could have, measuring F 2 (or alternatively the reduced cross-section σ r,NC ) at the LHeC/FCC-eh would provide discrimination between the two theoretical scenarios of small-x dynamics. Indeed, we see that the differences between the central values of the fixed-order and resummed fits in the restricted kinematic region covered by HERA are already comparable or larger than the size of the simulated pseudodata uncertainties.
This suggests that the inclusion of the LHeC/FCC-eh data for F 2 into a global fit would also provide discrimination power between the two theories, even if restricted to the HERA kinematic range. Finally, we see that differences are more marked for F L , with central values differing by several sigma (in units of the PDF uncertainty) in a good part of the accessible kinematic range. This is yet another illustration of the crucial relevance of measurements of F L to probe QCD in the small-x region (as highlighted also by Fig. 21).
The comparisons of Fig. 35 do not do justice to the immense potential of future high-energy lepton-proton colliders to probe QCD in a new dynamical regime. A more detailed analysis, along the lines of Ref. [216], involves including various combinations of LHeC/FCC-eh pseudodata (σ red NC , F L , F c 2 , etc.) into the PDF global analysis, allowing one to use the pseudo-data to reduce the PDF uncertainties and to quantify more precisely the discriminating power for small-x resummation effects with various statistical estimators, generalizing the analysis of the HERA data presented in Sect. 5. Such a program would illustrate the unique role of the LHeC/FCC-eh in the characterization of small-x QCD dynamics, and would provide an important input to strengthen the physics case of future high-energy lepton-proton colliders.
As a first step in this direction, we have performed variants of the NNPDF3.1sx fits including various combinations of the LHeC and FCC-eh pseudo-data of σ red NC . Specifically, we have used the LHeC (FCC-eh) pseudo-data on E p = 7 (50) TeV + E e = 60 GeV collisions, where the central value of the pseudo-data has been assumed to correspond to the NNLO+NLLx prediction computed with the corresponding resummed PDFs. All experimental uncertainties of the pseudo-data have been added in quadrature. The fits have been performed at the DIS-only level, since we have demonstrated in Sect. 5 that the small-x results are independent of the treatment of the hadronic data. Here we will show results of the fits including both LHeC and FCC-eh pseudodata, other combinations lead to similar qualitative results.
First of all we discuss the fit results at the χ 2 /N dat level. For simplicity, we show only the results of the HERA inclusive cross-sections as well as that of the LHeC and FCC-eh pseudo-data: for all other experiments, the values presented in Table 3 are essentially unchanged. As shown in Table 7, it is not possible to find a satisfactory fit to the LHeC/FCC-eh pseudo-data on inclusive cross-sections using NNLO theory while assuming that NNLO+NLLx theory is the correct underlying theory, as we have done here. As expected, the most marked differences are observed for the FCC-eh  pseudo-data. Note that the last row in Table 7 corresponds to the sum of the three experiments listed on the table. By performing the same analysis as in Fig. 24, we have verified that the significant improvement in χ 2 /N dat between the NNLO and NNLO+NLLx fits arises from the bins in the small-x and small-Q 2 region. Next in Fig. 36 we show the comparison between the gluon and the singlet PDFs at Q = 100 GeV in the NNPDF3.1sx NNLO+NNLx fits without and with the LHeC+FCC-eh pseudo-data on inclusive structure functions. Note that the latter is a DIS-only fit, hence the differences observed at large-x. For completeness, we also show the results of the corresponding NNPDF3.1sx NNLO fit with LHeC+FCC-eh pseudo-data. In the case of the NNLO+NLLx, we see that the central values coincide within uncertainties (as expected by construction) and there is a significant uncertainty reduction both for the gluon and for the singlet. In particular, the LHeC+FCC-eh kinematic coverage ensures that a precision measurement of the small-x gluon, with few-percent errors down to x 10 −7 , would be within reach.
From Fig. 36 we also see that for the gluon case, the NNLO and NNLO+NNLx fits with LHeC+FCC-eh pseudo-data are very different from each other. For instance, at x 10 −5 , where we gluon can be pinned down with 1% errors, the central values of the two fits differ by ∼ 15%. This comparison highlights that the fixed-order description of the small-x region at these future high-energy colliders would be completely unreliable, and that accounting for the effects of resummation at small x is required for any quantitative prediction. Indeed, the LHeC and FCC-eh would be truly unique machines in their potential to unveil the new dynamical regimes of QCD that arise in the deep small-x region.

Summary and outlook
The search for evidence of novel dynamics at small x beyond the linear fixed-order DGLAP framework has been an ongoing enterprise ever since the HERA collider started operations about 25 years ago. While some tantalizing hints have been reported, until now no conclusive evidence had been Fig. 36 Comparison between the gluon (left plot) and the singlet (right plot) PDFs in the NNPDF3.1sx NNLO+NNLx fits without and with the LHeC+FCC-eh pseudo-data on inclusive structure functions. For com-pleteness, we also show the results of the corresponding NNPDF3.1sx NNLO fit with LHeC+FCC-eh pseudo-data found in the HERA inclusive deep-inelastic structure functions. On the contrary, fixed-order perturbative QCD calculations have been remarkably successful, leading to good agreement with experimental data even in kinematic regions where they might naively be expected to fail.
From the theoretical point of view, formalisms for consistently including small-x resummation in DGLAP evolution and partonic coefficient functions were developed more than a decade ago . While these were sufficient to explain the success of fixed-order perturbation theory in describing the data, a state-of-the-art global PDF fit including the effects of small-x resummation was never performed. 8 It was the main goal of this study to bridge this gap, and to present the first genuine attempt at cutting-edge global NLO and NNLO PDF analyses which include the subtle effects of small-x resummation.
This has been made possible thanks to a number of developments both from the theory and from the implementation points of view. These include the consistent matching of NLLx small-x resummation to both the NLO and the NNLO fixed-order results, the resummation of the heavyquark matching conditions and DIS coefficient functions, as well as the implementation of these theoretical developments in the public code HELL and its interface with the APFEL program [62,63]. Also crucial to the success of the enterprise was the development of the NNPDF fitting technology, which is sensitive enough to identify small effects without them being masked by systematic methodological uncertainties from the fit procedure.
The main result of this work is the demonstration that including small-x resummation stabilizes the perturbative expansion of the DIS structure functions at small x and Q 2 , and thus also of the PDFs extracted from them. Specifically, the PDFs obtained with small-x resummation using NLO+NLLx and NNLO+NLLx theory are in much closer agreement with each other at medium and small x than the corresponding fixed-order NLO and NNLO PDFs. This suggests in turn that the theoretical uncertainty due to missing higher order corrections in a NNLO+NLLx resummed calculation is rather less at small x than that of the fixedorder NNLO calculation. This result is reflected in the marked improvement in the quantitative description of the HERA inclusive structure function data at small x when NNLO+NLLx resummed theory is used rather than NNLO: the NNLO+NLLx theory describes the low Q 2 and low x bins of the HERA data just as well as it describes the data at higher Q 2 and larger x. This effect is seen both in the inclusive neutral current and in the charm cross-sections, as expected from small-x resummation. We thus find no need for higher-twist contributions at low x, as proposed e.g. in Refs. [64,69], at least in the region where the resummed perturbative calculation is valid.
We have also presented here a first exploration of the phenomenological implications of our results. It has been understood for some time that the effect of resummation on the evolution of the PDFs can have a significant impact on the shape of parton luminosities and thus of hadronic cross-sections at LHC [44]. We have now shown that the further effect on parton luminosities of including small-x resummation in a PDF fit to low-Q 2 data at small x also remains sizable even at higher scales. We therefore expect that at the LHC smallx resummation might have significant effects, at either low invariant masses or at high rapidities, and thus that the accu-rate description of processes in these kinematic regions will require small-x resummation Conversely, present and future LHC measurements might provide further evidence for the onset of BFKL dynamics, this time in proton-proton collisions.
Small-x resummation also plays a crucial role in shaping the physics case for future high-energy lepton-proton colliders such as the LHeC and the FCC-eh, which would extend the coverage of HERA by up to two orders of magnitude into the small-x region. In this respect, the NNPDF3.1sx fits can be used to improve the accuracy of existing calculations of deep-inelastic scattering processes at these new machines. We have also demonstrated that a clear probe of BFKL dynamics is provided by the UHE neutrino-nucleus crosssections, where differences in event rates could be observed by upcoming measurements with neutrino telescopes such as IceCube and KM3NET.
The main limitation of the present analysis is the need to impose stringent cuts to the fitted hadronic data, in particular for Drell-Yan production, in order to ensure that the contamination from unresummed partonic cross-sections is kept to a minimum. On the one hand, it is well understood how to combine resummation corrections to partonic cross-sections with resummed parton luminosities to obtain fully resummed cross-sections even when the coupling runs [44], and smallx resummed partonic cross-sections have been computed for many of the relevant collider processes [100,101,[104][105][106][107][108][111][112][113][114][115][116][117]. On the other hand, these calculations are still not available in a format amenable to systematic phenomenology, and some effort is still required before they can be used in PDF fits. Future work in this direction will allow us to include a wider range of hadron collider data into a fully consistent small-x resummed global fit by removing the need for such cuts, and therefore allow us to achieve the same experimental precision for the resummed PDFs as is now possible in fixed-order fits. Moreover, an accurate description of processes at high rapidity, such as forward Drell-Yan and D meson production at LHCb, is likely to require the simultaneous resummation of both small-x and large-x logarithms, since at high rapidity while one of the partons is at very small x, the other is at very large x.
Finally, we would like to emphasize that the implications of our results go beyond what is traditionally thought of as "small-x physics". As LHC data become ever more precise, the theoretical challenge is to reduce theoretical uncertainties down to the 1% level, and this will require consistent calculations in perturbative QCD at N 3 LO. Recent progress with four-loop splitting functions [86,87] suggests that this may be possible rather sooner than was previously thought. However, while at NNLO the most singular term in the gluon splitting function is of order α 3 s x ln 1 x (the term with two logarithms being accidentally zero), at N 3 LO the most singu-lar term is or order α 4 s x ln 3 1 x . We thus expect the instability in fixed-order perturbative evolution at small x to be rather worse at N 3 LO than it was at NNLO. Small-x resummation would then be mandatory for improved precision, and this would require N 3 LO + NNLLx calculations to properly resum all the small x logarithms. While there has been some progress in extending the BFKL kernel to NNLLx [88][89][90][91][92][93][94], much work remains to be done.

Delivery
The fits presented in this work are available in the LHAPDF6 format [217] from the webpage of the NNPDF collaboration: http://nnpdf.mi.infn.it/nnpdf3-1sx These sets are based on the global dataset and contain N rep = 100 replicas. Specifically, the following fits are available: • Baseline NLO and NNLO NNPDF3.1sx sets, which are based on the global dataset with the kinematical cut of H cut = 0.6 applied to the hadronic data: NNPDF31sx_nlo_as_0118 NNPDF31sx_nnlo_as_0118 • Resummed NLO+NLLx and NNLO+NLLx NNPDF3.1sx sets, which are the resummed counterparts of the baseline sets above, based on an identical input dataset with the only difference of the theory settings:

NNPDF31sx_nlonllx_as_0118 NNPDF31sx_nnlonllx_as_0118
In addition, the other NNPDF3.1sx fits presented in this work, such as the DIS-only fits, are available upon request from the authors.
The DIS-only fits with various combinations of the LHeC and FCC-eh pseudo-data, discussed in Sect. 6.3, are also available in the same webpage: