Proton structure functions at small x

Proton structure functions are measured in electron-proton collision through inelastic scattering of virtual photons with virtuality Q on protons; x denotes the momentum fraction carried by the struck parton. Proton structure functions are currently described with excellent accuracy in terms of scale dependent parton distribution functions, defined in terms of collinear factorization and DGLAP evolution in Q. With decreasing x however, parton densities increase and are ultimately expected to saturate. In this regime DGLAP evolution will finally break down and non-linear evolution equations w.r.t x are expected to take over. In the first part of the talk we present recent result on an implementation of physical DGLAP evolution. Unlike the conventional description in terms of parton distribution functions, the former describes directly the Q dependence of the measured structure functions. It is therefore physical insensitive to factorization scheme and scale ambiguities. It therefore provides a more stringent test of DGLAP evolution and eases the manifestation of (non-linear) small x effects. It however requires a precise measurement of both structure functions F2 and FL, which will be only possible at future facilities, such as an Electron Ion Collider. In the second part we present a recent analysis of the small x region of the combined HERA data on the structure function F2. We demonstrate that (linear) next-to-leading order BFKL evolution describes the effective Pomeron intercept, determined from the combined HERA data, once a resummation of collinear enhanced terms is included and the renormalization scale is fixed using the BLM optimal scale setting procedure. We also provide a detailed description of the Q and x dependence of the full structure functions F2 in the small x region, as measured at HERA. Predictions for the structure function FL are found to be in agreement with the existing HERA data.


Introduction
The description of the proton in terms its elementary constituents, quarks and gluons, remains one of the big unsovled problems of nuclear-and elementary particle physics. At the typical energy scale of the proton, which is of the order of Λ QCD 200 MeV, Quantum Chromodynamics, the Quantum Field Theory description of strong interactions is strongly coupled, and quarks and gluons are subject to confinement. It is however possible to obtain very valuable information about the structure of the proton from collision processes of protons with leptonic projectiles, such as the electron. Due to the point-like structure of the electron and a very good theoretical understanding of electromagnetic interactions, the electron provides the perfect probe to explore the nucleon. To leading order in Quantum Electrodynamics (QED), scattering of the electron and the proton takes place through the exchange of a virtual photon with virtuality q 2 = −Q 2 , see Fig. 1. If the photon virtuality is large, the proton is destroyed during the scattering and the process is generally referred to as Deep Inelastic Scattering (DIS). The cross-section for neutral-current DIS on unplolarized nucleons can be written in terms of two Lorentz invariant structure functions F 2 and F L in the following way d 2 σ dxdQ 2 = 4πα 2 e.m. xQ 4 1 − y + y 2 2 Here y = (q ·p)/(k ·p) denotes the inelasticity with 0 < y < 1, see also Fig. 1. The structure functions themselves depend on only two Lorentz invariants, the photon virtuality Q 2 and Bjorken x = Q 2 /(2p · q). Within the parton model, to be discussed below, x denotes the momentum fraction of the parton hit by the virtual photon.
If the photon virtuality Q 2 is significantly larger than the non-perturbative energy scale of the proton, asymptotic freedom provides for such processes a weak strong coupling constant α s (Q 2 ) 1 and a description within perturbative theory becomes possible. The conventional theoretical framework for such DIS processes is based on the collinear factorization theorem [2]. At leading order, the essential physics is captured by the parton model [3] of the proton. Within this model, the highly virtual photon interacts not with the entire proton with characteristic size ∼ 1/Λ QCD , but with a single, essentially point-like, parton, i.e., a quark or gluon, with effective size 1/Q. Interference effects with spectator quarks or gluons are on the other hand suppressed by powers of Q 2 . To arrive at the complete cross-section, the "partonic" interaction of virtual photon and quark, needs to be convoluted with parton distribution functions (PDFs), f i (x, Q 2 ), i = q,q, g, which encode the probability to find a parton with a certain proton moment fraction inside the proton. Higher order corrections to such partonic cross-sections, calculated within QCD perturbation theory, possess reveal then a new kind of singularity, apart from the conventional ultra-violet singularity, which is removed through renormalization of the QCD Lagrangian. This new singularity is of infra-red type and can be associated with configu-rations where an additionally emitted parton is collinear to the proton momentum. In physical terms, this initial state singularity reflects interference between the perturbative calculable partonic interactions at the hard scale Q, with spectator quarks and gluons, characterized by the hadronic scale Λ QCD . Collinear factorization provides then a systematic framework to remove such singularities from the perturbative hard cross-section, resulting into finite Wilson coefficients [4][5][6][7], absorbing them into parton distribution functions, which encode the long-distance, non-perturbative physics; for a recent review see, e.g., [8].
To make the separation between long-and short-distance physics manifest, one needs to introduce some arbitrary factorization scale µ f , apart from the scale µ r appearing in the renormalization of the strong coupling α s . The independence of physical observables such as F 2,L on µ f can be used to derive powerful renormalization group equations (RGEs) governing the scale dependence of PDFs in each order of perturbation theory, known as the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation. The corresponding kernels are the anomalous dimensions or splitting functions associated with collinear two-parton configurations [9][10][11]. DGLAP evolution has been impressively confirmed by experiment, in particular through the very accurate DIS data on the structure function F 2 from the DESY-HERA experiment, see and evolved with the DGLAP equations, have been shown to describe F 2 data over orders of magnitude, both in x and Q 2 [13]. Despite of the impressive success of the DGLAP evolution equations, theoretical considerations suggests that at some point in phase space this description is supposed to break down. With decreasing x, logarithms ln 1/x increase and are capable of balancing the strong coupling, α s ln 1/x ∼ 1, leading to a break-down of the naïve perturbative expansion, see Fig. 3. The necessary resummation of enhanced terms (α s ln 1/x) n is the achieved by the Balitsky-Fadin-Kuarev-Lipatov (BFKL) evolution equation [14]. One of the main predictions of BFKL evolution is a power-like rise of the gluon density. If continued to ultra-small x, the 1/Q 2 expansion, on which collinear factorization is based, will eventually break down. A description of proton structure functions in this region of phase space is provided by the Balitsky-Kovchegov (BK) [15] and JIMWLK [16] evolution equations, which provide a non-linear extension of BFKL evolution, resumming corrections due to high gluon densities to all orders. While theoretical arguments suggests the relevance of such corrections already at current collider energies, current data provide no clear evidence for deviations from linear DGLAP evolution, which would provide a signature for the onset of a non-linear kinematic regime dominated by high, or saturated, gluon densities. Definite evidence for such a regime of QCD requires therefore new experiments such as a future high-luminosity electron-ion collider, i.e.., the EIC [1] and the LHeC [17] projects, whose physics case is currently studied. In particular these projects plan to measure both structure functions F 2 and F L and their scaling violations very precisely at small x both in electron-proton and in electron-heavy ion collisions.
From the theory side this requires the development of suitable tools which allows to pin down possible deviations from DGLAP evolution. In particular it is necessary to reduce the large freedom in fitting initial conditions of parton distribution functions. With non-linear saturation effects most likely to manifest themselves at small values of Q 2 , the large number of free parameters used for the description of initial conditions in PDF fits, does not allow to exclude the possibility that saturation effects, while present in reality, are currently hidden in the initial conditions of DGLAP evolution.
In the following we present two approaches which have the potential to restrict this large freedom at low scales. Section 2 is dedicated to the concept of physical evolution kernels, which allows to reduce the number of independent PDFs in DIS fits and eliminates scale-and scheme dependence in their definition. Section 3 contains results of a recently achieved BFKL fit of the combined HERA data. While both theoretically and experimentally less explored than collinear factorization, BFKL evolution has the potential to reveal the emergence of non-linear effects more easily than DGLAP evolution. Unlike DGLAP evolution, BFKL drives the system into the saturated regime, making a detection of high density effects more likely. For details we refer to [18][19][20]

Physical evolution kernels for DIS observables
Since collinear factorization can be carried out in infinitely many different ways, one is left with an additional choice of the factorization scheme for which one usually adopts the MS prescription. Likewise, the RGE governing the running of α s with µ r can be deduced from taking the derivative of F 2,L with respect to µ r . Upon properly combining PDFs and Wilson coefficients in the same factorization scheme, any residual dependence on µ f is suppressed by an additional power of α s , i.e., is formally one order higher in the perturbative expansion but not necessarily numerically small. Alternatively, it is possible to formulate QCD scale evolution equations directly for physical observables without resorting to auxiliary, convention-dependent quantities such as PDFs. This circumvents the introduction of a factorization scheme and µ f and, hence, any dependence of the results on their actual choice. The concept of physical anomalous dimensions is not at all a new idea and has been proposed quite some time ago [4,21,22] but its practical aspects have never been studied in detail. The framework is suited best for theoretical analyses based on DIS data with the scale µ r in the strong coupling being the only theoretical ambiguity. In addition, F 2,L or their scaling violations can be parametrized much more economically than a full set of quark and gluon PDFs, which greatly simplifies any fitting procedure and phenomenological analysis. The determination of α s from fits to DIS structure functions is the most obvious application, as theoretical scheme and scale uncertainties are reduced to a minimum.
Here we largely focus on the practical implementation of physical anomalous dimensions in analyses of DIS data up to next-to-leading order (NLO) accuracy. We shall study in detail potential differences with results obtained in the conventional framework based on scaledependent quark and gluon densities, which could be caused by the way how the perturbative series is truncated at any given order.

Theoretical Framework
This gist of the factorization scheme-invariant framework amounts to combine any two DIS observables {F A , F B } and determine their corresponding 2 × 2 matrix of physical anomalous dimensions instead of the scale-dependent quark singlet, Σ ≡ q (q +q), and gluon distributions appearing in the standard, coupled singlet DGLAP evolution equations. Instead of using measurements of F 2 and F L (actually their flavor singlet parts), one can also utilize their variation with scale for any given value of x, i.e., dF 2,L (x, Q 2 )/d ln Q 2 as an observable. The required sets of physical anomalous dimensions for both {F 2 , F L } and {F 2 , dF 2 /d ln Q 2 } have been derived in [22] up to NLO accuracy. The additionally needed evolution equations for the non-singlet portions of the structure functions F 2,L are simpler and not matrix valued. As we shall see below, the required physical anomalous dimensions comprise the inverse of coefficient and splitting functions and are most conveniently expressed in Mellin n moment space. The Mellin transformation of a function φ given in Bjorken x space, such as PDFs or splitting functions, is defined as where n is complex valued. As an added benefit, convolutions in x space turn into ordinary products upon applying (2), which, in turn, allows for an analytic solution of QCD scale evolution equations for PDFs. The corresponding inverse Mellin transformation is straightforwardly performed numerically along a suitable contour in n space, see, e.g., Ref. [23] for details. The necessary analytic continuations to non-integer n moments are given in [24,25], and an extensive list of Mellin transforms is tabulated in [26]. We will work in Mellin space throughout this review. Assuming factorization, moments of DIS structure functions F I at a scale Q can be expressed where the sum runs over all contributing n f active quark flavors with electric charge squared e 2 q and the gluon g, each represented by a PDF f k . For e 2 g , the averaged quark charge factor e 2 q = (1/n f ) q e 2 q has to be used. µ r and µ f specify renormalization and factorization scale, respectively. The scale Q 0 defines the starting scale for the PDF evolution, where a set of nonperturbative input distributions needs to be specified. For simplicity we identify in the following the renormalization scale with the factorization scale, i.e., µ r = µ f ≡ µ. The coefficient functions C I,k are calculable in pQCD [4][5][6][7] and exhibit the following series in a s ≡ α s /4π where m 0 depends on the first non-vanishing order in a s in the expansion for the observable under consideration, e.g., m 0 = 0 for F 2 and m 0 = 1 for F L .
where the l → k splitting functions have a similar expansion [9][10][11] as the coefficient functions in Eq. (4): The P kl (n) relate to the corresponding anomalous dimensions through γ kl (n) = −2P kl (n) in the normalization conventions we adopt, where we use the leading order (LO) and NLO expressions for γ kl (n) given in App. B of the first reference in [10]. We note that the same normalization is used in the publicly available Pegasus evolution code [23]. In practice one distinguishes a 2 × 2 matrix-valued DGLAP equation evolving the flavor singlet vector comprising Σ(n, µ 2 /Q 2 0 ) and g(n, µ 2 /Q 2 0 ) and a set of n f − 1 RGEs for the relevant non-singlet quark flavor combinations. The scale-dependent strong coupling itself obeys another RGE governed by the QCD beta function with β 0 = 11 − 2n f /3 and β 1 = 102 − 38n f /3 up to NLO accuracy. To compare below with the results for the physical anomalous dimensions in Ref. [22] we also introduce the evolution variable Instead of studying F I (n, Q 2 ) in (3) in terms of scale-dependent PDFs, which are obtained from solving the singlet and non-singlet DGLAP equations (5) in a fit to data [13], one can also derive evolution equations directly in terms of the observables F I (n, Q 2 ). To this end, we consider a pair of DIS observables F A and F B , to be specified below, whose scale dependence is governed by a coupled matrix-valued equation for the remainders.
The required physical anomalous dimensions in Eqs. (9) and (10), obey a similar perturbative expansion in a s as in (6). The singlet kernels in (9) are constructed by substituting into the left-hand side of Eq. (9) and taking the derivatives. Note that we have normalized the quark singlet part of F A,B with the same averaged charge factorē 2 f which appears in the gluonic sector. Upon making use of the RGEs for PDFs and the strong coupling in Eq. (5) and (7), respectively, one arrives at where we have introduced 2 × 2 matrices for the relevant singlet coefficient and splitting functions, respectively. An analogous, albeit much simpler expression holds for the NS kernel K (N S) in (10). As has been demonstrated in [22], the kernels (12) are independent of the chosen factorization scheme and scale but do depend on µ r and the details of the renormalization procedure. We also note that the inverse C −1 in (12), appearing upon re-expressing all PDFs by F A,B , can be straightforwardly computed only in Mellin moment space.
2.2. Example I: F 2 and F L Let us first consider the evolution of the pair of observables {F 2 , F L }. A precise determination of F L in a broad kinematic regime is a key objective at both an EIC [1] and the LHeC [17].
Since the perturbative series for F L only starts at O(a s ), one wants to account for this offset by actually considering the evolution of either L,q )}. Both sets of kernels K AB show a rather different behavior with n, as we shall illustrate below, but without having any impact on the convergence properties of the inverse Mellin transform needed to obtain x dependent structure functions. The kernels K AB at LO and NLO accuracy for L,g )} can be found in [22]. Note that evolution in [22] is expressed in terms of t. Using (8), d/da s = −2/(a s β 0 )d/dt, and (7) to compute to extra terms proportional to β(a s ), we fully agree with their results.   45) in Ref. [22]. For NLO kernels we refer to [18].  L,i and β 0 , respectively, to K (1) AB have been omitted. A global factor of α s /4π has been ignored in the perturbative expansion, i.e., In Fig. 4 we illustrate the n dependence of the LO and NLO singlet kernels K AB for the L,q )} assuming α s = 0.2 and n f = 3. As can be seen, NLO corrections are sizable for all singlet kernels, in particular, when compared with the perturbative expansion of the singlet splitting functions P kl (n) in (6); see Figs. 1 and 2 in Ref. [11]. This is, however, not too surprising given that the known large higher order QCD corrections to the Wilson coefficients C L,g and C L,q [27] are absorbed into the physical anomalous dimensions K AB for the evolution of the DIS structure functions F 2 and F L . The impact of contributions from the NLO coefficients C L,g and C L,q on the results obtained for K AB is illustrated by the dash-dotted lines in Fig. 4. Another source for large corrections are the terms proportional to β 0 in the NLO corrections as can be inferred from the dotted lines; note that K L,q . In Sec. 2.4 we will demonstrate how the differences between the LO and NLO kernels become apparent in the scale evolution of F 2,L (x, Q 2 ).
gq ∼ 1/n, and P (0) gq ∼ ln n. The NLO kernel K L2 exhibits an even stronger rise with n. In the same way one obtains, for instance, that K L,q )} only grows like ln n, see Eq. (14).
Despite this peculiar n dependence and the differences between the singlet kernels shown in Fig. 5, both sets of observables, L,g )} can be used interchangeably in an analysis at LO and NLO accuracy. Results for the QCD scale evolution are identical, and one does not encounter any numerical instabilities related to the inverse Mellin transform, which we perform along a contour as described in Ref. [23]. In fact, it is easy to see that the eigenvalues which appear when solving the matrix valued evolution equation (9), are identical for both sets of kernels and also agree with the corresponding eigenvalues for the matrix of singlet anomalous dimensions P kl ; see also the discussions in Sec. 2.4 2.3. Example II: F 2 and dF 2 /dt Of future phenomenological interest could be also the pair of observables {F 2 , dF 2 /dt}, in particular, in the absence of precise data for F L . Determining experimentally the t or Q 2 slope of F 2 is, of course, also challenging. Defining F D ≡ dF 2 /dt, we obtain the following physical evolution kernels at LO to be used in Eq. (9); for NLO kernels see [18]. The kernels K AB in (16) exhibit more moderate higher order corrections, mainly through terms proportional to β 0,1 , than those listed in Sec. 2.2. This shall become apparent in the next Section when we discuss results for the scale dependence of both {F 2 , dF 2 /dt} and {F 2 , F L }.

Numerical Studies
In this Section we apply the methodology based on physical anomalous dimensions as outlined above and compare with the results obtained in the conventional framework of scale-dependent quark and gluon densities and coefficient functions. Due to the lack of precise enough data for F L or dF 2 /dt we will adopt the following realistic "toy" initial conditions for the standard DGLAP evolution of PDFs at a scale Q 0 = √ 2 GeV [23] xu v (x, for all our numerical studies. The value of the strong coupling α s at Q 0 is taken to be 0.35. For our purposes we can ignore the complications due to heavy flavor thresholds and set n f = 3 throughout. We use this set of PDFs to compute the flavor singlet parts of F 2 , F L , and dF 2 /dt at the input scale Q 0 using Eq. (11). For studies of DIS in the small x region, say x 10 −3 , in which we are mainly interested in, the flavor singlet parts are expected to dominate over NS contributions and, hence, shall be a good proxy for the full DIS structure functions. Results at scales Q > Q 0 are obtained by either solving the RGEs for PDFs or by evolving the input structure functions directly adopting Eq. (9). For the solution in terms of PDFs we adopt from now on the standard choice µ = Q.
For completeness and to facilitate the discussions below, let us quickly review the solution of the matrix-valued RGEs such as Eqs. (5) and (9). While one can truncate the QCD beta function and the anomalous dimensions consistently at any given order in a s , there exists no unique solution to beyond the LO accuracy. The matrix-valued nature of (9) only allows for iterations around the LO solution, which at order a k s can differ in various ways by formally higher-order terms of O(a l>k s ). To this end, we employ the standard truncated solution in Mellin moment space, which can be found, for instance, in Ref. [24], see also [23], and reads where the evolution operator up to NLO is defined as with Here, a 0 = a s (Q 0 ), Γ P = Σ g , and Γ K = F A F B , i.e., the index i = P refers to the coupled RGE for the quark singlet and gluon and i = K to the RGE for the pair {F A , F B } of DIS structure functions in (9). For i = K one has with a corresponding definition for i = P in terms of the 2 × 2 matrices of singlet splitting functions P (0) and P (1) . λ ± denote the eigenvalues given in Eq. (15) and e ± the projection operators onto the corresponding eigenspaces; see Refs. [23,24]. As has been mentioned already at the end of Sec. 2.2, the eigenvalues λ ± (n) are identical when computed for the kernels K AB and P kl . This in turn implies that as long as, say, F 2 and F L are calculated at µ = Q 0 with LO accuracy, their scale evolution based on physical anomalous dimensions reproduces exactly the conventional results obtained with the help of scale-dependent PDFs. Figure 6 shows our results for the scale dependence of the DIS structure functions F 2 and F L . The input functions at Q 0 = √ 2 GeV are shown as dotted lines. While LO results are identical, starting from NLO accuracy the comparison between the two methods of scale evolution becomes more subtle, and results seemingly differ significantly as can be inferred from the middle panels of Fig. 6.
The origin of the differences between F 2,L (x, Q 2 ) computed based on Wilson coefficients and scale-dependent PDFs and physical anomalous dimensions can be readily understood from terms which are formally beyond NLO accuracy. For instance, upon inserting the NLO Wilson coefficients (4) and the truncated NLO solution (18)- (20) into Eq. (11), F 2 at O(a s ) contains spurious terms of both O(a 0 a s ) and O(a 2 s ). Since F L starts one order higher in a s , similar terms are less important here. On the other hand, when we evolve F 2,L with the help of physical anomalous dimensions we first compute, due to the lack of data, the input at a 0 based on Eq. (11), which then enters the RGE solution (18)- (20). Again, this leads to terms beyond NLO. In case of F 2 they are now of the order O(a 0 a s ) and O(a 2 0 ), i.e., even more relevant than in case of PDFs since a 0 > a s .
To test if the entire difference between the two evolution methods shown in Fig. 6 is caused by these formally higher order contributions, one can easily remove all O(a 2 s ), O(a 0 a s ), and O(a 2 0 ) contributions from our results. Indeed, the scale evolution based on physical anomalous dimensions and the calculation of F 2,L from PDFs then yields exactly the same results also at NLO accuracy. We note that this way of computing properly truncated physical observables from scale-dependent PDFs beyond the LO accuracy has been put forward some time ago in Ref. [28,29] but was not pursued any further in practical calculation.
Another interesting aspect to notice from Fig. 6 are the sizable NLO corrections illustrated in lower panels, in particular, for F 2 in the small x region. For this comparison, LO results refer to the same input structure functions F 2,L as used to obtain the NLO results but now evolved at LO accuracy, i.e., by truncating the evolution operator in Eqs. (18)-(20) at L (0) K . At first sight the large corrections appear to be surprising given that global PDF fits in general lead to acceptable fits of DIS data even at LO accuracy [13]. However, this is usually achieved by exploiting the freedom to have different sets of PDFs at LO and, say, NLO accuracy. The framework based on physical anomalous dimensions does not provide this option as the input for the scale evolution is, in principle, fully determined by experimental data, and only the value for the strong coupling can be adjusted at any given order. In this sense it provides a much more stringent test of the underlying framework and perhaps a better sensitivity to, for instance, the possible presence of non-linear effects to the scale evolution in the kinematic regime dominated by small x gluons. . Same as in Fig. 7 but now for the pair of observables F 2 and F D ≡ dF 2 /dt.
In Fig. 7 we show the corresponding results for the scale dependence of the DIS structure function F 2 and its slope F D = dF 2 /dt. Again, any differences between the scale evolution performed with physical anomalous dimensions and based on PDFs are caused by formally higher order terms O(a 2 s ), O(a 0 a s ), and O(a 2 0 ), which can be removed with the same recipe as above. As for {F 2 , F L }, NLO corrections are sizable in the small x region due to numerically large contributions to K DD from the QCD beta function.

Summary
We have presented a phenomenological study of the QCD scale evolution of deep-inelastic structure functions within the framework of physical anomalous dimensions. The method is free of ambiguities from choosing a specific factorization scheme and scale as it does not require the introduction of parton distribution functions. Explicit results for the physical evolution kernels needed to evolve the structure functions F 2 , its Q 2 slope, and F L have been presented up to next-to-leading order accuracy.
It was shown that any differences with results obtained in the conventional framework of scaledependent quark and gluon densities can be attributed to the truncation of the perturbative series at a given order in the strong coupling. At next-to-leading order accuracy the numerical impact of these formally higher order terms is far from being negligible but, if desired, such contributions can be systematically removed. A particular strength of performing the QCD scale evolution based on physical anomalous dimensions rather than auxiliary quantities such as parton densities is that the required initial conditions are completely fixed by data and cannot be tuned freely in each order of perturbation theory. Apart from a possible adjustment of the strong coupling, this leads to easily testable predictions for the scale dependence of structure functions and also clearly exposes the relevance of higher order QCD corrections in describing deep-inelastic scattering data. Next-to-leading order corrections have been demonstrated to be numerically sizable, which is not too surprising given that the physical evolution kernels absorb all known large higher order QCD corrections to the hard scattering Wilson coefficients.
Once high precision deep-inelastic scattering data from future electron-ion colliders become available, an interesting application of our results will be to unambiguously quantify the size and relevance of non-linear saturation effects caused by an abundance of gluons with small momentum fractions. To this end, one needs to observe deviations from the scale evolution governed by the physical anomalous dimensions discussed in this work. The method of physical anomalous dimensions can be also used for a theoretically clean extraction of the strong coupling and is readily generalized to other processes such polarized deep-inelastic scattering or inclusive one-hadron production.
3. F 2 and F L at small x using collinearly-improved BFKL resummation 3.1. Structure functions within the BFKL framework At small x and center-of-mass energy s = Q 2 /x, we can apply high energy factorization and write the structure functions F I , I = 2, L as where q ≡ q 2 ⊥ ) . Φ P is the non-perturbative proton impact factor which we model using where we have introduced two free parameters and a normalization. Φ I is the impact factor associated to the photon which we treat at leading-order (LO), i.e. where , Ω 2 = (11 + 12ν 2 )/8, Ω L = ν 2 + 1/4, and the strong coupling α s is fixed at the renormalization scale µ 2 . In the present work we will also use the kinematically improved impact factors proposed in [30,31], which include part of the higher order corrections by considering exact gluon kinematics. Its implementation requires to replace the functions c I (ν) byc I (γ, ω) wherec ψ(γ) is the logarithmic derivative of the Euler Gamma function and ξ = 1 − 2γ + ω, while ω is the Mellin variable conjugate to x in the definition of the gluon Green function F, see Eq. (28) below. The main difference between these impact factors is that the LO ones roughly double the value of their kinematically improved counterparts in the region with small |ν|, while being very similar for |ν| ≥ 1.
The gluon Green function can be written in the form withᾱ s = α s N c /π. The collinearly improved BFKL kernel as introduced in eq. (28) is an operator consisting of a diagonal (scale invariant) pieceχ(γ) with eigenvalue where , plus a termχ RC (γ) proportional to β 0 which contains the running coupling corrections of the NLO kernel [32]:χ The precise form of the NLO kernel χ 1 can be found in [19,33]. The resummation of collinear logarithms of orderᾱ 3 s and beyond is realized by the term [19,34,35] Our final expression for the structure functions reads where M 2 and M L can be found in [19]. For the kinematical improved version of F I we replace c I (ν) byc I (1/2 + iν, χ(1/2 + iν)). In Eq. (33) the scale of the running of the coupling has been set to µ 2 = QQ 0 . Building on the work of [36] we found in [19] that in order to obtain a good description of the Q 2 dependence of the effective intercept of F 2 , λ, for x < 10 −2 , it is very useful to operate with non-Abelian physical renormalization schemes using the Brodsky-Lepage-Mackenzie (BLM) optimal scale setting [37] with the momentum space (MOM) physical renormalization scheme [38]. For technical details on our precise implementation we refer the reader to [19] (see also [39] for a review on the subject and [40] for a related work). More qualitatively, in these schemes the pieces of the NLO BFKL kernel proportional to β 0 are absorbed in a new definition of the running coupling in order to remove the infrared renormalon ambiguity. Once this is done, the residual scheme dependence in this framework is very small. We also found it convenient [19] to introduce, in order to describe the data with small Q 2 , an analytic parametrization of the running coupling in the infrared proposed in [41].

Comparison to DIS experimental data
In the following we compare our results with the experimental data for F 2 and F L . Let us first compare the result obtained in [19] for the logarithmic derivative d log F 2 /d log(1/x) using Eq. (33) with a LO photon impact factor and our new calculation using the kinematically improved one. In Fig. 8 we present our results with the values of our best fits for both types of impact factors and compare them with the H1-ZEUS combined data [12] for x < 10 −2 . The values of the parameters defining the proton impact factor in (23) and the position of the (regularized) Landau pole (we use n f = 4) for the strong coupling are δ = 8.4, Q 0 = 0.28 GeV, Λ = 0.21 GeV for the LO order case and δ = 6.5, Q 0 = 0.28 GeV, Λ = 0.21 GeV for the kinematically improved (note that the normalization C does not contribute to this quantity).  Figure 8. Fit to λ for F 2 with the LO photon impact factor (solid line) and the kinematically improved one (dashed line). The data set has been extracted from [12].
The LO impact factor generates lower values than the kinematically improved one in the high Q 2 region and slightly higher ones when Q 2 2 GeV 2 . It is interesting to see how the approach presented here allows for a good description of the data in a very wide range of Q 2 , not only for high values, where the experimental uncertainties are larger, but also in the non-perturbative regions due to our treatment of the running of the coupling. Encouraged by these positive results we now turn to investigate more differential distributions. We select data with fixed values of x and compare the Q 2 dependence of our theoretical predictions with them, now fixing the normalization for the LO impact factor to C = 1.50 and 2.39 for the kinematically improved. Our results are presented in Fig. 9. The equivalent  Figure 9. Study of the dependence of F 2 (x, Q 2 ) on Q 2 using the LO photon impact factor (solid lines) and the kinematically improved one (dashed lines). Q 2 runs from 1.2 to 200 GeV 2 .
comparison to data, this time fixing Q 2 and looking into the evolution in the x variable, is shown in Fig. 10. We observe that our predictions give a very accurate description of the data for both types of impact factors. Let us remark that the values for the parameters in this fit are in syntony with the theoretical expectations for the proton impact factor since Q 0 is very similar to the confinement scale and  x Q² 120 GeV² Figure 10. Study of the dependence of F 2 (x, Q 2 ) on x using the LO photon impact factor (solid lines) and the kinematically improved one (dashed lines). Q 2 runs from 1.2 to 120 GeV 2 .
the value of δ sets the maximal contribution from the impact factor also in that region. This is reasonable given that the proton has a large transverse size.
The longitudinal structure function is an interesting observable which is very sensitive to the gluon content of the proton. We will now present our predictions for F L using the best values for the parameters previously obtained in the fit of F 2 . We will see that the agreement with the data is very good. First, Q 2 is fixed and the x dependence is investigated in Fig. 11. The experimental data have been taken from [43]. To present the Q 2 dependence it is convenient to calculate, for each bin in Q 2 , the average value of x, see Fig. 12. In some sense this is a similar plot to the one previously presented for λ in the F 2 analysis and we can see that the effect of using different types of impact factors is to generate a global shift in the normalization. Again we note that we have an accurate description of the transition from high to low Q 2 , which was one of the main targets of our work.    Figure 11. Fit to F L with the LO photon impact factor (solid lines) and the improved one (dashed lines). The experimental data are taken from [43].

Predictions for future colliders
While our predictions for the structure functions are in agreement with the data from the HERA collider experiments H1 and ZEUS, these observables are too inclusive to provide unambiguous evidence for BFKL evolution (for other recent studies in this context see [42]). Comparable in quality fits can be obtained by both DGLAP evolution and saturation models, see e.g. [43,44]. In order to distinguish among different parton evolution pictures new collider experiments are needed, such as the proposed Electron-Ion-Collider (EIC) at BNL/JLab (USA) [1] and the Large Hadron Electron Collider (LHeC) at CERN (Switzerland) [45], which will be able to measure both F 2 and F L at unprecedented small values of Bjorken x. In Fig. 13 we present two studies with our predictions for F 2 and F L down to values of x = 10 −6 .  Figure 12. The proton structure function F L as a function of Q 2 . The average x values for each Q 2 of the H1 data (black) are given in Figure 13 of [43]. ZEUS data are taken from [44]. The solid line represents our calculation with the LO photon impact factor and the dashed line using the kinematically improved one.  Figure 13. Predictions for F 2 (left) and F L (right) for LHeC. On the left plot, the curve with Q 2 = 10GeV 2 can be compared with Figure 4.13 of [45]. Simulated measurements for F L in the kinematic range plotted here (right) can be found in Figure 3.7 of the same reference.

Conclusions
We have presented an application of the BFKL resummation program to the description of the x and Q 2 dependence of structure functions as extracted from Deep Inelastic Scattering data at HERA. We have also provided some predictions for these observables at future colliders. In order to obtain the correct dependence on the virtuality of the photon at high values of the scattering energy, we have included in the BFKL kernel the main collinear contributions to all orders. We have also used optimal renormalization and an analytic running coupling in the infrared in order to accurately describe the regions of low Q 2 .

Summary
DIS scattering experiments allow to explore the proton in terms of its QCD content, i.e., quarks and gluons. At small values of x, the description in terms of linear DGLAP evolution and parton distribution functions is expected to break down and non-linear effects, associated with high gluon densities, are believed to set in. To pin down such effects in DIS, new collider experiments are needed, which allow to measure both structure functions F 2 and F L with high accuracy for both electron-proton and electron-nucleus scattering. DGLAP evolution formulated in terms of physical evolution kernels allows for a direct evolution of structure function doublets. Apart from removing scale-and scheme dependence in the description of structure functions, it further reduces the number of free parameters, used in the parametrization of non-perturbative initial conditions. BFKL evolution describes on the other hand directly evolution in x and hence drives the system into the non-linear regime, promising higher sensitivity to non-linear effects. Being less explored than DGLAP evolution, we took a first step towards such applications by confronting BFKL evolution with the combined HERA data. In particular we demonstrated that NLO BFKL evolution is capable to describe the combined HERA data, if the NLO BFKL kernel is supplemented with collinear resummation and optimal scale setting for the QCD running coupling is used.