Color Glass Condensate at next-to-leading order meets HERA data

We perform the first dipole picture fit to HERA inclusive cross section data using the full next-to-leading order (NLO) impact factor combined with an improved Balitsky-Kovchegov evolution including the dominant effects beyond leading logarithmic accuracy at low $x$. We find that three different formulations of the evolution equation that have been proposed in the recent literature result in a very similar description of HERA data, and robust predictions for future deep inelastic scattering experiments. We find evidence pointing towards a significant nonperturbative contribution to the structure function for light quarks, which stresses the need to extend the NLO impact factor calculation to massive quarks.


I. INTRODUCTION
The inner structure of protons and nuclei can be accurately determined in deep inelastic scattering (DIS) experiments, where the target stucture is probed by a simple pointlike electron via the exchange of a virtual photon. For proton targets, the combined structure function data from the H1 and ZEUS experiments at HERA [1][2][3][4] have made it possible to extract the parton densities with an excellent precision.
At small momentum fraction x the gluon densities rise rapidly, and one eventually expects non-linear highoccupancy effects to be important and become visible in the weak coupling regime. At high gluon densities, these non-linear effects tame the growth of the gluon density, and a dynamical scale known as the saturation scale Q 2 s is generated. This scale characterizes the region of phase space where the non-linear saturation effects dominate. To describe QCD in this high energy regime an effective theory known as the Color Glass Condensate has been developed, see Refs. [5,6] for a review.
The precise DIS data can provide a crucial test for the saturation picture. Theoretically the inclusive DIS cross section is a relatively simple observable, as the probe has no internal structure and one does not need to consider e.g. fragmentation effects. As the proton structure is not perturbatively calculable, some input from experimental data is needed. In the CGC framework, one can calculate the energy dependence of various observables, e.g. the total photon-proton cross section, perturbatively by resumming contributions enhanced by a large logarithms of energy or ln 1/x. The non-perturbative input in this case is the proton structure at an initial (and smallish) Bjorken-x, which is a parametrized input fitted to the data. The leading order CGC calculations have been able to obtain a good description of the precise HERA data by fitting the initial condition with only a few free parameters [7][8][9]. However, in all these fits one needs to introduce an additional fit parameter to slow down the x-evolution to be compatible with the HERA measurements.
To precisely test the saturation picture of CGC, it is crucial to move beyond leading order accuracy. In recent years the theory has been rapidly developing towards full NLO accuracy. The impact factors, describing the photon-proton interaction, have been calculated at this order in case of massless quarks [10][11][12][13][14][15], and the first numerical results were reported in Ref. [16]. The impact factors need to be combined with evolution equations that describe the Bjorken-x dependence and resum contributions enhanced by large logarithms of energy, (α s ln 1/x) n at leading order and α s (α s ln 1/x) n at nextto-leading order. The Balitsky-Kovchegov (BK) equation describing the evolution of the dipole-target interaction [17,18] is available at NLO accuracy [19] with the higher order contributions enhanced by large transverse logarithms resummed in Refs. [20][21][22][23] and numerical solutions reported in Refs. [24,25].
An additional complication in the small-x evolution is that the Coulomb tails obtained from a perturbative calculation result in the proton size growing much faster than seen in the data, and faster than suggested by the Froissart bound [26] for hadronic collisions. It has been argued [27,28] on the theoretical level that including some non-perturbative damping of the gluon emission at large transverse distance is necessary and sufficient to recover a Froissart behavior for the virtual photonproton cross section. This idea has been studied [29,30] in fits to the HERA data using the impact parameter dependent BK equation supplemented by either a nonperturbative cut-off or collinear resummations. In addition to the BK equation, one can solve the more general JIMWLK evolution equation [31][32][33][34][35][36][37] (available at NLO accuracy [38,39], but no numerical solution exists for the NLO equation). The JIMWLK evolved proton structure was compared with the HERA data in Ref. [40] (see also Refs. [41,42]), where again large non-perturbative contri-butions were needed to describe the system with a finite proton geometry. Due to these additional complications, we only study an impact parameter independent evolution here, and assume that the transverse area of the proton can be factorized in the cross section calculations.
This paper is structured as follows. First, in Sec. II, we will briefly introduce the dipole picture of DIS at leading and next-to-leading order. Then, in Sec. III, we will review the necessary details of the different variants of the BK equation used in this work. Section IV reviews the data sets used in the fits, and Sec. V discusses the results of the fits.

II. DEEP INELASTIC SCATTERING IN THE DIPOLE PICTURE AT NLO
The photon-proton cross section is parametrized in terms of the structure functions F 2 and F L , that are related to the virtual photon-proton cross sections σ γ * p as and Here the subscripts T and L refer to the transverse and longitudinal polarizations of the virtual photon. The experimental data is often reported as a reduced cross section σ r (x Bj , y, Q 2 ) = F 2 (x Bj , Q 2 ) − y 2 1 + (1 − y) 2 F L (x Bj , Q 2 ). (3) Here −Q 2 is the photon virtuality, x Bj is the Bjorken variable and y is the inelasticity.
The focus in this paper is on the next-to-leading order corrections to the total DIS cross section in the dipole picture. As an introduction, let us first briefly describe the process in the leading order dipole picture.
At leading order, the virtual photon-proton scattering in the dipole picture is understood in the following way (see e.g. [75]). First, the incoming photon fluctuates into a quark-antiquark pair. This splitting is described by the photon light cone wave function ψ γ * →qq . Subsequently, the produced dipole interacts with the target. At high energy, the quark-target interaction is eikonal, and the transverse position of the quark does not change during the scattering. Instead, the quark goes through a color rotation in the target color field and picks up a Wilson line V (x 0 ) in the fundamental representation, where x 0 is the transverse coordinate of the quark. Similarly, the antiquark at point x 1 picks up a conjugate Wilson line To calculate the total cross section, one applies the optical theorem and calculates the imaginary part of the forward elastic scattering amplitude for the process γ * p → γ * p. The resulting cross section reads Here, z is the light-cone momentum fraction of the photon carried by the quark. The dipole size is r = x 0 − x 1 and its impact parameter is b = (x 0 + x 1 )/2. In the following, S(r, b, x) is assumed to depend only slowly on b. Thus we will drop this dependence on b and replace the integration over b by a constant d 2 b → σ 0 /2. The dipole scattering matrix S is defined as a two point function of the Wilson lines that the quarks pick up in the scattering process: The brackets refer to the average over the target color charge configurations. Here the momentum fraction x in the subscript stands for the fact that the Wilson lines are evaluated at some energy or rapidity scale corresponding to the kinematics of the process. This dependence is given by the Balitsky-Kovchegov equation which, at leading order, is usually used to evolve the Wilson lines up to an evolution rapidity Y = log 1/x Bj . At NLO the question of the evolution rapidity becomes more complicated, as discussed in more detail in Sec. III. At Next-to-Leading Order (NLO) the virtual photonproton scattering involves Fock states of the photon that contain a gluon in addition to the quark and antiquark, which all scatter off the target. There are also other NLO contributions with only a quark-antiquark Fock state scattering off the target, which include a gluon loop correction to the photon splitting. These NLO qqg and qq contributions have been calculated independently using the conventional dimensional regularization [13,14] and four dimensional helicity schemes [15]. The individual diagrams contain UV divergences that cancel each other in the sum. On top of these, there remains a divergence related to low x gluons, which must be resummed into the evolution of the target. Subtraction schemes for this low x gluon divergence in DIS were devised and tested in Refs. [14,16], and in our present paper we continue to refine the 'unsub' scheme to enable a comparison between the theory and experimental data.
In Refs. [14,16] the low x gluon divergence factorization from the NLO DIS cross sections (for a more detailed discussion in the context of single inclusive particle production see Refs. [67,76,77]) were written in two distinct but equivalent forms: a form where the factorization is implicit, and another where it was made explicit, named 'unsubtracted' and 'subtracted' schemes respectively. In this work we use the unsubtracted form for the cross sections, which can be expressed as Here the first term is the leading order cross section (4) where the dipole scattering amplitude is evaluated at the chosen fixed initial rapidity scale of the target, corresponding to the initial condition of BK evolution. The other terms can be interpreted as arising from the NLO qq diagrams (σ dip L,T ) and from the NLO qqg diagrams (σ qg,unsub.

L,T
), up to subtraction terms used to make the cancellation of UV divergences between these diagrams explicit. In our scheme, the 'unsubtracted' qg term is: and the 'dipole' term is: with the shorthand xi := d 2 xi 2π . The integrand kernels and dipole operators for the leading order and 'dipole' terms are where Here, the rapidity scale which the dipole operator (5) is evaluated at is left implicit. It will be discussed together with the associated small-x evolution in Sec. III. However, we note already now that in the qg-term this rapidity scale must be taken to depend on the gluon momentum fraction z 2 , not just the external kinematical scales x Bj and Q 2 . This is essential for the stability of the factorization scheme, as discussed in great detail e.g. in Refs. [16,20,67,76,77].
For the qg-terms that only appear at NLO, the kernels and Wilson line operators are: Here z 0 , z 1 , z 2 are the longitudinal momentum fractions of the quark, antiquark and gluon, which satisfy i z i = 1. The parton configuration factor QX 3 is interpreted as the ratio of the qqg state formation time to the γ * lifetime [12]. It is defined as X 2 3 := z 0 z 1 x 2 01 + z 0 z 2 x 2 02 + z 2 z 1 x 2 21 . We have also defined a shorthand P (z) := 1 + (1 − z) 2 . The qqg state-target scattering Wilson line operator is In Eq. (7) the lower limit z 2,min in the gluon longitudinal momentum fraction integral is yet undefined, and its proper value will be discussed in the next section.

III. HIGH ENERGY EVOLUTION
A. Balitsky-Kovchegov equation In the calculation of the photon-proton cross section at NLO, as discussed above, the dipole-target scattering amplitude depends on the energy, or equivalently, on Bjorken-x. In the large N c limit, the evolution is given by the Balitsky-Kovchegov (BK) equation [17,18]. At leading order, the BK equation reads The kernel K BK is proportional to the probability density to emit a gluon with transverse coordinate x 2 from the dipole of size x 01 = x 0 − x 1 . The evolution rapidity Y is discussed in detail later. When running coupling corrections following the Balitsky prescription [78] are included, it reads In principle we should use the next-to-leading order BK equation when using the impact factors calculated to the order α s . The required numerical solution of the NLO BK equation exists [24,25,79]. However, the equation is numerically burdensome due to the high dimensional transverse integration (in the NLO BK equation one integrates over the transverse coordinates of the two emitted gluons, instead of just one gluon in the leading order equation). Instead of the full equation, in this work we use prescriptions of BK evolution that capture an important subset of beyond leading order effects. The difference between the studied evolutions reflects some of the uncertainty due to the missing full NLO evolution.
In practice we have chosen three related formulations of the BK equation that resum some or all of the large transverse momentum logarithms in the NLO equation. Firstly we consider the nonlocal evolution equation in terms of the projectile momentum fraction introduced in Ref. [20], where collinear double logarithms are resummed via the inclusion of a Kinematical Constraint: we denote this the KCBK equation. Secondly, we consider the local equation in the projectile momentum fraction of Ref. [21], where the same double logarithms, together with DGLAP-like single logarithms, are explicitly resummed into a kernel that is a nontrivial function of α s : we call this the ResumBK equation. Thirdly we study a nonlocal equation in the target momentum fraction, recently formulated in Ref. [23] and denoted here as the TBK equation. The first two are formulated in terms of the projectile momentum fraction, so that the rapiditylike evolution variable in the BK equation is defined by the momentum fraction of the probe, or equivalently by the plus component of the 4-momentum 1 . Since the projectile momentum fraction is the variable appearing explicitly in the NLO DIS impact factors, using these evolution equations is fairly straightforward. However the fact that the TBK equation is written in terms of the target momentum fraction, i.e. the minus component of 4-momentum, means that the evolution and the perturbative impact factors can only be matched approximatively, and the procedure requires more care.
As we will discuss in detail in the following, we solve the BK equation by taking a parametrized dipole amplitude as an initial condition at an initial rapidity scale. Starting from this initial condition the evolution predicts the behavior of the dipole at higher rapidities, i.e. energies, or correspondingly at smaller Bjorken-x. The phase space available for the emission of the gluon grows with energy and determines the amount of BK evolution. Thus the evolution range is controlled by the lower limit of the gluon momentum fraction z 2,min in Eq. (7), with a smaller lower limit corresponding to longer evolution.

B. Evolution in projectile momentum fraction
Let us first consider the evolution written in projectile momentum fraction, which is the case for the KCBK and ResumBK equations. In this case the high energy evolution is parametrized by the rapidity variable Y , which is defined using the plus components of the gluon momentum k + and a plus momentum scale P + associated with the target as Since the incoming photon energy q + (which is the maximal k + ) in the target rest frame is proportional to the photon-target c.m.s. energy W 2 , one should think of evolution in the rapidity variable Y as evolution in ln W 2 , as we will see more explicitly below.
In the impact factor, Eq. (7), the gluon momentum is parametrized by the momentum fraction z 2 as k + = z 2 q + . Both the probe momentum fraction evolution equations and the NLO impact factor are derived in terms of the same z 2 . Thus it is straighforward to see that the dipole operators in the evaluation of the qgterm in the cross section (7) are always evaluated at the projectile rapidity depending on the integration variable z 2 . We will specify the value of P + below. First, we have to determine the lower limit z 2,min for the z 2 integral in the NLO impact factor, Eq. (7), which controls the amount of evolution. This limit is set by the overall kinematics of the process. One way to understand the existence of this limit is to note that in the limit z 2 → 0 the invariant mass of the qqg-system interacting with the target grows as M 2 qqg ∼ 1/z 2 . The fact that this invariant mass cannot be larger than the c.m.s. collision energy results in a lower limit for kinematically allowed values of z 2 . Since the validity of the eikonal approximation used to derive the dipole picture cross section requires in principle M 2 qqg W 2 , one could require a more strict limit on z 2 than resulting from purely kinematics. Thus there is a choice in how close to the kinematical limit one allows the integral to go, which we quantify by the parameter e Y 0,if 1. In terms of this parameter we have the limit Here we have introduced a non-perturbative target transverse momentum scale Q 2 0 , for which in this work we use the value 2 Q 2 0 = 1 GeV 2 . This allows us to write P + = Q 2 0 /(2P − ), and we used the fact that x Bj = Q 2 /(2P · q) = Q 2 /(2P − q + ). This limit is already derived e.g. in Refs. [14,16,20]. In Ref. [16] the authors for simplicity set Q 2 0 /Q 2 = 1 in practical evaluations of the NLO impact factors.
In principle also the limits z 1 → 0, z 1 → 1 in Eqs. (7) and (8) correspond to the invariant mass of the scattering state becoming infinite, similarly to the limit z 2 → 0. Thus, as discussed in Ref. [14], one could also take the energy or rapidity scale at which the dipoles are evaluated to depend on the (anti) quark momentum fractions z 0 , z 1 . This part of phase space does not, however, generate a parametrically large logarithmic contribution to the cross section. Thus we leave it to further work to explore improved treatments of these "aligned jet" contributions.
When the Bjorken-x of the process is such that the smallest momentum fraction z 2,min is close to 1, i.e. when , the possible phase space for real gluon emission allowed in the expression for the cross section vanishes. Thus the qg contribution to the NLO cross section goes to zero at x Bj ∼ e −Y 0,if by construction. The NLO calculation does not fix any exact value for Y 0,if . A possible choice to consider for Y 0,if would be to take Y 0,if ≈ ln 1/0.01, corresponding to the limit were the dipole picture is usually considered applicable. This was the choice used in Ref. [16]. However, this choice leads to a transient effect in the NLO cross sections at the upper end of the x Bj range x Bj ∼ e −Y 0,if since the positive virtual correction remains large, while the the negative qg contribution vanishes, as demonstrated in Ref. [16].
To avoid this unphysical transient effect, we adopt here instead the maximal (or minimal depending on the point of view) choice Y 0,if = 0. This means that the integral over z 2 in the cross section extends all the way to the kinematical limit, outside of the validity of the eikonal approximation. The contribution from this region is, however, only a parametrically small part of the cross section for small x Bj , which is where we are comparing the cross section to experimental data. Also, since there is a cancellation between the real and virtual contributions to the cross section, and the latter includes a z 2 integral over the full range 0 < z 2 < 1, one could in fact argue that this choice minimizes the net effect of very large invariant mass states in the photon on the cross section.
The above discussion only applies to the qg term (7) in the cross section. The virtual correction in Eq. (8) is already integrated over z 2 and cannot be evaluated at a z 2 -dependent rapidity. Thus for this term dipole operators are taken to be independent of z 2 and evaluated at rapidity Y = ln 1/x Bj . Using a z 2 independent dipole is justified, as the region z 2 1 gives only a negligible contribution to the virtual correction. Including these formally subleading effects, namely the z 2 dependent dipole operator in the virtual term and improving the approximation Y ≈ ln 1/x Bj is left for future work.
The choice Y 0,if = 0 removes the unphysical transient effect, but it forces us to confront another problem that the earlier formulation of Ref. [16] wanted to avoid by choosing a larger Y 0,if . Namely, at the lower end of the z 2 range we are forced to evaluate also the (BK-evolved) dipoles at a rapidity scale that is lower (or x Bj -scale that is higher) than where the BK-equation is normally used. Now we again have different options regarding the rapidity where we start the the BK evolution. We parametrize this choice by another constant Y 0,BK , whose value can also be chosen in different ways.
One way is to take Y 0,BK = Y 0,if = 0, in which case we simply start the BK evolution much earlier (much higher x Bj ) than where we are actually calculating the cross section. Here the contribution of the unphysical small rapidity or large x phase space to the cross section (7) is suppressed, because target gets more and more dilute following the evolution backwards to smaller rapidities. This procedure changes the way the parametrization of the initial condition for the BK-evolution should be interpreted. In this approach, the quantity that can meaningfully be compared to the initial dipole amplitude at x = x 0 ∼ 0.01 in LO fits is not the actual initial condition at Y = Y 0,BK = 0, but the result obtained after Y = ln 1/0.01 units of rapidity evolution.
Another option is to take a more typical initial energy scale for the BK equation, which we here take as Y 0,BK = ln 1/0.01. In this latter case one has to model the dipole amplitude in the region Y 0,if < Y < Y 0,BK . In this case, we simply assume that the dipole operator is independent of Y in this region, which is "before the initial condition" in Y . Assuming an energy independent dipole amplitude in this region is consistent within the accuracy of the framework.
To summarize, we have two parameters that we must choose, Y 0,if and Y 0,BK . In this work we always take Y 0,if = 0 to avoid the large transient effect at x Bj ∼ e −Y 0,if in the data region. We then apply two approaches for the parameter Y 0,BK . The first option is to start the BK evolution at rapidity Y 0,BK = ln 1/0.01 and freeze the dipole amplitude at Y < Y 0,BK . The second option is to also start the BK evolution at rapidity We recall that the dipole amplitudes in the cross section Eq. (7) are evaluated at a rapidity Note that the maximum rapidity Y max encountered is obtained at the z 2 → 1 limit corresponding to values of Y probed by the the z 2integral in the cross section ranging from Y 0,if to Y max , i.e. over a rapidity interval ∆Y = ln 1/z 2,min . When Y 0,BK > Y 0,if , the actual range of BK evolution is smaller by Y 0,BK , and for Y 0,if < Y < Y 0,BK the dipole does not change. We emphasize that the evolution rapidity depends only on the total center-of-mass energy W 2 , and not explicitly on x Bj or Q 2 . This is natural, as the scattering amplitude for a dipole with a fixed transverse size can only be sensitive to the total center-of-mass energy, and strictly speaking the dipole does not exactly know about the photon virtuality or the Bjorken x.

C. Kinematically constrained BK
As discussed previously, in order to keep the computational cost of the fit procedure manageable, we do not use the full NLO BK equation to obtain the rapidity dependence of the dipole scattering matrix S. Instead, we use modified versions of the leading order evolution equation that resums the most important higher order corrections, in particular the collinear double logarithms.
First, we use the KCBK equation [20] that is nonlocal in the projectile momentum fraction (i.e. the evolution variable Y ): This equation explicitly forces time ordering between subsequent gluon emissions. The theta function ensures that only dipoles in the range Y > Y 0,if are included.

D. Rapidity local resummed BK
The most important higher order corrections to the BK equation that are enhanced by double large transverse logarithms can be resummed alternatively into a kernel that is local in the evolution rapidity Y by a method introduced in Ref. [21] 3 . This procedure resums exactly the same contributions that are included in Ref. [20] to derive the kinematically constrained BK equation shown above in Eq. (21). A practical advantage of the approach taken in Ref. [21] is that the resulting equation is local in evolution (projectile) rapidity, and as such numerically easier to solve using standard Runge-Kutta methods. In addition to the double transverse logarithms resummation, the contribution of some of the single transverse logarithms present in the NLO BK equation can be included following Ref. [22], keeping the equation local in rapidity. In Ref. [80] it was shown that this resummed BK equation is in practice close to the kinematically constrained BK equation discussed previously. As the resulting resummed evolution equation is written in terms of the projectile rapidity Y , it can be used with the impact factors exactly as the kinematically constrained BK equation.
The resummed equation is obtained by multiplying the BK kernel (15) by K DLA K STL , where K DLA is a resummation of double and K STL single transverse logarithms. The kernel resumming the double transverse logarithms reads with x = ln x 2 02 /x 2 01 ln x 2 12 /x 2 01 andᾱ s = α s N c /π. If ln x 2 02 /x 2 01 ln x 2 12 /x 2 01 < 0, an absolute value of the argument is used and the Bessel function is changed to J 1 → I 1 , see Ref. [21]. The single transverse logarithms ∼ α s ln 1/(x 2 ij Q 2 s ) are included multiplying the kernel multiplied by In Ref. [25] it was shown that the resummation of single transverse logarithms can be done such that the resummed equation is a good approximation to the full NLO BK evolution by adjusting the constant C sub whose numerical value is not fixed by the resummation procedure. This renders the O(α 2 s ) contributions in the NLO BK equation that are not enhanced by large (single) transverse logarithms minimal. With this procedure, one obtains a rapidity local projectile momentum fraction resummed BK equation which we use as an approximation to the full NLO BK equation (with a resummation of large transverse logarithmic corrections), with C sub = 0.65 determined in Ref. [25].
The resummation of the single transverse logarithms is completely independent of the resummation of the double transverse logarithms, and thus could be included in the same way also in the other studied evolution equations. In this work, however, we only include this contribution in the ResumBK evolution, as we prefer to work with the established versions of the BK evolution. We will discuss the effect of the single transverse logarithm resummation on our fits in Sec. V.

E. Target momentum fraction evolution
As discussed in detail in Ref. [23], it is possible to formulate the evolution in terms of the target rapidity η defined as a logarithm of the minus component of the momentum. This corresponds to a fraction of the total longitudinal momentum of the target, which is the variable used in its DGLAP evolution, and also the usual physical interpretation of x Bj . In order to translate a plus momentum to a minus, one needs to have access to the correct transverse momentum scale. In the case of the whole DIS process this would naturally be Q 2 . Thus we would want to define things in such a way that the largest evolution rapidity reached in the process is with Y max from Eq. (20). Here we see that the target rapidity η is directly related to the Bjorken-x in DIS, as expected. Thus, similarly as one can think of evolution in Y as evolution in ln W 2 , evolution in η corresponds to evolution in ln 1/x Bj . The complication in using the target momentum fraction is that both the evolution equation and the impact factors are written in transverse coordinate space, which is natural for the eikonal interaction with the target. Thus the gluon transverse momentum is not very explicit in either. The usual procedure is to use an uncertainty principle argument and estimate the transverse momentum as the inverse of the corresponding transverse distance. In both the BK equation and the impact factor one integrates over transverse distances up to infinity, which would correspond to zero transverse momentum and infinite η (for a fixed Y ). Distances longer than some nonperturbative scale should, however, not have a significant effect on the physics. Thus we do not want large dipoles with sizes above a (soft) target transverse momentum scale 1/Q 2 0 (the same Q 2 0 that we have already used) to appear in the relation between the rapidities Y and η. In practice we are thus led to consider a dipole of size r at a projectile momentum fraction corresponding to Y , to have a target evolution rapidity η given by We see that with this definition we always have η < Y , which corresponds to the fact that for perturbative size dipoles r 2 < 1/Q 2 0 we always have less evolution in η than in Y (see Eq. (25)).
The evolution equation for the dipole amplitude in terms of the target rapidity η was derived in Ref. [23] as 4 (27) whereS refers to the dipole scattering matrix depending on the target rapidity η, instead of projectile rapidity Y . This evolution equation then needs to be provided with an initial condition at the initial rapidity η 0,BK . We include running coupling corrections and use the kernel K BK (x 0 , x 1 , x 2 ) from Eq. (15). The rapidity shift reads The step function with δ ≡ max{δ 02 , δ 21 } ensures that the equation is a well defined initial value problem and no information about the dipole amplitude for η < η 0,BK affects the evolution. When calculating the cross section, we useS(r, η) =S(r, η 0,BK ) for η < η 0,BK . Equation (27) is the "canonical" BK equation from [23], which contains an all-order resummation of the double collinear logarithm enhanced corrections, and thus is perturbatively correct up to an error of O(ᾱ 2 s ). The full NLO BK evolution in target rapidity has not been solved numerically so it is not known in practice how well this resummation captures the NLO effects. In Ref. [23] a comparison is made between two formulations of the equation with double logarithm resummations, and the differences are minor, and mostly in the early evolution. The resummed evolution is also compared to the LO BK evolution formulated in target rapidity, and the resummed evolution is found to be notably slower.
To use the target rapidity dependent dipole amplitudes in the NLO impact factors, we simply need to replace the dipoles in the impact factor with the η-dependent dipoles. The rapidity argument is determined by using z 2 to obtain Y with Eq. (19), which is then transformed into η using Eq. (26), i.e.
with Y defined in Eq. (19). The regulator ensures that the rapidity shift is always negative, consistent with the definition (26) and with the rapidity shift in the TBK evolution equation (28).
Let us finally discuss the kinematical limits in the z 2 integral and their connection to the target momentum fraction probed by the process. The lower limit z 2 > z 2,min of the z 2 -integral (18) corresponds to the lower limit of the η values probed by the impact factor With our choice Y 0,if = 0 the values of η needed in the cross section always extend down to evolution rapidities before the initial condition that is imposed at η 0,BK . In this region η < η 0,BK the dipole operators are, as in Ref. [81], just frozen to the initial value: S(r, η < η 0,BK ) ≡S(r, η 0,BK ). At the upper limit z 2 = 1, on the other hand, the largest values of η are reached for r > 1/Q 0 , with the range in η extending up to which is the same as the maximum Y reached in projectile momentum fraction evolution. Here let us note two things. Firstly, the min function used in defining the target momentum fraction rapidity (26) prohibits η from getting infinitely large (which would require an infinite amount of evolution) for very large dipoles r > 1/Q 0 that are not expected to contribute significantly to the cross section. Secondly, the amount of evolution, or largest rapidity reached, in target momentum fraction given by x Bj , as in Eq. (25), strictly speaking applies only to typical dipole sizes r ∼ 1/Q. For larger dipoles 1/Q r 1/Q 0 one actually evolves further in the target momentum fraction.

F. Running coupling
For the strong coupling constant in coordinate space we use the expression with β = (11N c − 2N F )/3 and N F = 3, Λ QCD = 0.241 GeV. The parameters µ 0 and c control how the coupling is frozen in the infrared, and we choose µ 0 /Λ QCD = 2.5 and c = 0.2. With this choice, the coupling freezes to α s = 0.762 in the infrared.
We have performed fits with two different running coupling prescriptions. The first one is denoted Balitsky + smallest dipole (Bal+SD) scheme below. In this scheme, we use the Balitsky prescription from Ref. [78] in the BK evolution as in Eq. (15). In the NLO impact factor, Eq. (7), and in the terms resumming large transverse logarithms, Eqs. (23) and (24) Note that the Balitsky prescription reduces to the smallest dipole one when one of the dipoles is much smaller than the others. For comparison we also use another scheme denoted as parent dipole. Here, the scale is always set by the size of the parent dipole, both in the evolution equation and in the impact factor. In the LO-like σ dip L,T term of the impact factor, Eq. (8), there are no daughter dipoles in the scattering state. For this term the smallest dipole scheme is equivalent with the parent dipole scheme.

G. Initial conditions
The initial condition for the (projectile momentum fraction) BK evolution is parametrized at rapidity Y = Y 0,BK . We use the MV γ parametrization, and write the initial condition as The fit parameters in the initial condition are Q 2 s0 which controls the saturation scale at the initial x, and the anomalous dimension γ which determines the shape of the dipole amplitude at small |x ij |. We note that this parametrization results in both a negative unintegrated gluon distribution and negative particle production cross sections in proton-nucleus collisions at high transverse momenta if γ > 1. As the inclusive DIS measurements are not sensitive to asymptotically small dipoles, we do not consider our dipole amplitude to be valid in that region and as such, having an anomalous dimension γ > 1 is acceptable. The practical interpretation of γ in our fit is that it controls the shape of the dipole amplitude in the transient region r ∼ 1/Q s . The leading order BK fits to HERA data generally prefer γ ∼ 1.1 [8,9], and similar results were found in recent fits where the BK equation with some higher order corrections resummed [81] was used. For a detailed discussion related to the Fourier positivity of the dipole amplitude, the reader is referred to Ref. [82].
For the local resummed projectile momentum fraction (ResumBK) evolution, the resummation should also in principle affect the initial condition [21]. However, as the initial condition is in any case a non-perturbative input, we will use the same parametrization of Eq. (34) also for solving the ResumBK equation.
For target momentum fraction evolution, the initial condition for the evolution (27) corresponds to the scattering amplitudeS(r, η = η 0,BK ) at some rapidity η 0,BK . We use the same parametrization, Eq. (34), as in the case of projectile rapidity evolution.

IV. AVAILABLE DIS DATA
The HERA experiments H1 and ZEUS have published their combined measurements for the reduced cross section σ r in Refs. [1,2]. Additionally, the charm and bottom quark contributions to the fully inclusive data are available [3,4]. As the impact factors at next to leading order accuracy in the massive quark case are not available, we only calculate the light quark contribution to the photon-proton cross section. In the leading order fits [8,9] it has been possible to obtain a good description of the fully inclusive data with only light quarks, even though the charm contribution is significant (parametrically up to ∼ 40% at Q 2 m 2 c ). On the other hand, the leading order fits aiming to simultaneously describe the total and charm structure function data require separate parameters (e.g. different transverse areas) for the light and charm quarks [8], or an additional effective soft and non-perturbative contribution [29,40].
In this work we consider two different setups. First, we follow the strategy that has been successfully used at leading order and calculate the light quark contribution to the structure functions, and compare with the inclusive HERA data from Ref. [1]. We note that the newer combined dataset containing data from the HERA-II run is also available [2], but at low x and moderate Q 2 the two datasets result in very similar fits (see e.g. Ref. [83]).
As a second approach, we construct an interpolated dataset that only contains the light quark contribution. Since the charm and bottom data are not measured in the same kinematical x, Q 2 bins as the inclusive data, it is not possible to just subtract the heavy quark contribution from the fully inclusive cross section. Instead, we use a leading order dipole model fit from Ref. [83], where the Bjorken-x and dipole size r dependence is described using the so called IPsat parametrization [84]. This parametrization includes a smooth matching to the DGLAP evolution [85][86][87][88] in the dilute region, and at large dipoles or densities the scattering amplitude saturates to unity. The advantage of this parametrization is that it results in an excellent description of both inclusive and heavy quark datasets. Consequently, it can be used to interpolate the charm and bottom contributions to the structure functions. We use this parametrization to subtract the heavy quark contributions from the measured reduced cross section. We then use this interpolated light quark-only data in the NLO fits. In our procedure we do not modify the uncertainties of the inclusive data in the subtraction (another possible approach would be to reduce the uncertainties proportionally). This is not really a consistent treatment for the errors; ultimately only the experimental collaborations would be in a position to correctly take into account the correlation between errors in the total and heavy quark data. Thus the errors and consequently χ 2 values in the light-quark fits are not correct statistically. However, we expect the magnitude of the uncertainties to only affect the final fit qualities, and to have only a limited effect on the extracted best fit parameter values and the interpretation in terms of physics. This detail must be kept in mind for the interpretation of the χ 2 values from the light quark fits. The total reduced cross section for some Q 2 bins from HERA [1] is shown in Fig. 1, and compared with the result obtained by the IPsat fit mentioned above. The description of the data is excellent. The interpolated light quark data in the same kinematics is also shown, and compared to the light quark reduced cross section computed using the same IPsat fit.
When fitting the initial condition for the BK evolu-tion, we consider datapoints in the region 0.75 < Q 2 < 50 GeV 2 at x < 0.01. This results in N = 187 datapoints to be included in the fit. Although the correlation matrix for the experimental uncertainties is available [1], we do not take these correlations into account as we expect it to have only a negligible effect in our fits.

V. FIT RESULTS
In this section we will look at our fit results. The discussion is divided first by the data that is being fitted, followed by a comparison of the evolution prescriptions in the kinematical domain accessible in future DIS experiments, which lies outside the HERA region included in the fits.
Let us first recall the essential details of our fit schemes. The choice of a fit scheme consists of the version of the BK evolution equation (discussed in Secs. III C, III D and III E), the running coupling scheme (see Sec. III F), and the starting point of the BK evolution, parametrized in terms of Y 0,BK or η 0,BK . The fit results in values for the free parameters characterizing the initial condition as discussed in Sec. III G: Q 2 s0 , σ 0 and γ, and in a value for the parameter C 2 in the scale of the running coupling, see Sec. III F.
Our main fit results are presented in Tables I, II and III classified by the BK equation used, with secondary and tertiary grouping keys being the running coupling scheme and Y 0,BK (or η 0,BK ) controlling the rapidity scale of the BK initial condition used in the fits. The saturation scale Q 2 s defined as N (r 2 = 2/Q 2 s ) = 1 − e −1/2 is also shown at fixed projectile rapidity Y = ln 1 0.01 . We will first discuss in the next subsection the fits to the full HERA reduced cross section data, and in the following subsection the fits to the interpolated light quark pseudodata presented in Sec. IV and labeled as light-q in the tables where the fit results are shown. The two datasets differ enough to to warrant their own discussion.

A. Fitting the HERA reduced cross section
Before we discuss the results and their systematic features in more detail we show in Fig. 2 that all three BK evolutions combined with next to leading order impact factors are capable of describing the HERA data equally well. The results shown are obtained using the Bal+SD running coupling, and Y 0,BK = η 0,BK = ln 1/0.01, but excellent fit results are obtained with other scheme choices, too. Even though the resulting parametrizations for the dipole at initial rapidity can differ significantly, the resulting reduced cross sections are mostly indistinguishable.
We first present in Table I the fit results obtained using the kinematically constrained BK equation as discussed in Sec. III C. We find a very good description (χ 2 /N = 1.49) of the HERA data using our main setup with the  Bal+SD prescription and Y 0,BK = 0. We consider this as our preferred HERA data fit, with a BK equation derived in the same framework as the impact factor, a theoretically preferred running coupling scheme, and only one starting scale Y 0,if = Y 0,BK = 0. We note that starting the BK evolution at Y 0,BK = ln 1 0.01 (and freezing the dipole at smaller rapidities) results in an equally good fit. This suggests that we are only weakly sensitive to the details of extrapolation scheme used to describe the dipole amplitude in the region Y 0,if < Y < Y 0,BK . The parameter C 2 controlling the evolution speed is not required to be large as it is in the case of leading order fits, where one generally finds C 2 ∼ 10 [8,9]. Instead, we find C 2 ≈ 0.85, which is of the same ballpark as the general estimate C 2 = e −2γ E ≈ 0.3 [89,90].
As seen in Table. I, larger values of C 2 are required in the parent dipole scheme fits. This is expected, as C 2 maps the coordinate space scale x 2 ij to momentum space C 2 /x 2 ij , and in the parent dipole scheme the coordinate space scale is generically larger. Consequently a larger C 2 is needed to render the strong coupling values, and the resulting evolution speeds, comparable between the coupling constant scheme choices.
We generically find γ > 1 at the initial condition, with the exception γ ≈ 1 found in the case where the evolution starts at Y 0,BK = ln 1 0.01 and the parent dipole prescription for the running coupling is used. We note that γ > 1 is also required in the leading order fits to obtain a Q 2 dependence at the initial condition compatible with the HERA data [8,9]. The disadvantage of an initial condition with γ > 1 is that, as discussed in Sec. III G, it results in the unintegrated gluon distribution not being positive definite at large transverse momenta.
To understand why different running coupling prescriptions result in different initial anomalous dimen-    sions, we study the slope of the dipole defined as For the KCBK fits this is shown in Fig. 3 as a function of dimensionless dipole size rQ s . The kinematically constrained BK equation is found to keep the anomalous dimension (slope at small r) approximatively constant at very small r, unlike the leading order BK equation. A similar effect was found in case of the ResumBK equation in Ref. [25]. At intermediate r ∼ 1/Q s which dominates the cross section, there is clear evolution towards an asymptotic shape. Let us first focus on results where the smallest dipole Bal+SD coupling is used. Here, the anomalous dimension is large at the initial condition and the evolution decreases the slope at intermediate rQ s , which results in the cross section growing more rapidly with Q 2 . If the BK evolution is started at the rapidity scale from which there is a long evolution before entering the data region (i.e. Y 0,BK = 0), a larger initial anomalous dimension is required in order to obtain the shape dictated by the Q 2 dependence of the HERA structure function data around r ∼ 1/Q s . Let us then consider the evolution with the parent dipole prescription. In this case, we start from a relatively small γ = 0.98, and the evolution increases the slope at small (but not asymptotically small) r. This can be seen to stem from the fact that in the parent dipole prescription the coupling, and consequently the evolution speed of the dipole amplitude N (r), grows more as a function of parent dipole size r in comparison to other running coupling prescriptions. At larger r, the slope evolves only slightly. After a few units of rapidity evolution, the dipole amplitudes have the same shape in the r ∼ 1/Q s region independently of the running coupling prescription. This is expected, as r ∼ 1/Q s size dipoles dominate when calculating the structure functions in HERA kinematics.
Next we move to the local projectile momentum fraction (ResumBK) fits, the results of which are shown in Table II. In general, the results are close to the ones previously discussed in case of the kinematically constrained BK equation. This is not surprising, as both equations are designed to include the same subset of higher order corrections enhanced by large double transverse log- arithms. Similarly to the preferred fit with KCBK, in the ResumBK fit to HERA data with Bal+SD and Y 0,BK = 0 the obtained C 2 is quite small at C 2 ≈ 0.49 ∼ e −2γ E and the anomalous dimension is large, γ = 1.56. The obtained anomalous dimension values behave similarly as in the case of KCBK, and there similarly seems to be a systematic preference for smaller σ 0 /2 with the Balitsky + smallest dipole coupling. The ResumBK equation evolves generically more slowly than the KCBK equation, which is reflected in the required C 2 values being smaller (except in the case Y 0,BK = 0 with parent dipole prescription when the C 2 values are comparable). This is a consequence of the Re-sumBK equation including an additional resummation of some single transverse logarithms. The main effect of this resummation is that it results in a slower evolution, see discussion in Sec. III D. We have confirmed numerically that if the resummation of single transverse logarithms is not included, our fit results are almost intact, except that a larger value for the parameter C 2 is obtained.
In both ResumBK and KCBK fits with Bal+SD running coupling, the obtained values for the proton transverse area σ 0 /2 are generally smaller than what is found in leading order fits with similar running coupling schemes, with or without a resummation of large transverse logarithms [8,9,81]. The obtained saturation scales at Y = ln 1 0.01 , on the other hand, are comparable to the leading order fit results. In the LO fits, one typically obtains σ 0 /2 ∼ 16 mb (proton sizes comparable to our results were found in the leading order fit presented in Ref. [80] where double logarithmic corrections were resummed in the BK equation similarly as in our setup).
We note that the proton transverse area can in principle be obtained by studying the squared momentum transfer t dependence of exclusive vector meson production. If the cross section is written as e −B D |t| at small |t|, the HERA measurements on J/ψ production [91,92] give B D ≈ 4 GeV −2 . Depending on the assumed proton density profile, this corresponds to σ 0 /2 ≈ 9.8 . . . 19.6 mb (using Gaussian or a step function profile). As the vector meson t spectra are not measured precisely enough especially at large |t|, the exact form of the proton density profile can not be deduced. Consequently, we find that all obtained values for the proton transverse size σ 0 /2 in our fits to HERA reduced cross section data are compatible with the J/ψ spectra. However, we also note that the step function profile is not really favored by the HERA data [93]. Thus one would prefer values that are in the lower part of the range σ 0 /2 ≈ 9.8 . . . 19.6 mb. Indeed, especially with the Balitsky + smallest dipole running coupling, our fit results for the proton size also favor such smaller target sizes for the proton.
As we are neglecting the impact parameter dependence, we can not compute the evolution of the proton transverse area and consequently use a fixed σ 0 /2 at all x Bj . We note that the HERA vector meson production data [92,94] suggest that the transverse area depends logarithmically on the center-of-mass energy. This growth is effectively included in the energy dependence of the proton saturation scale in our framework.
Let us finally discuss the results obtained with the third evolution equation considered in this work, the BK equation formulated in terms of the target momentum fraction (TBK). The fit results in this case are shown in Tab. III. While the fit qualities overall are quite similar to the projectile momentum fraction setups, we find some departures from the shared qualitative features of the KCBK and ResumBK fits. With the Balitsky + smallest dipole running coupling the TBK evolution needs to be slowed down more with a larger values of C 2 compared to KCBK and ResumBK equations. The TBK fit with parent dipole coupling is more mixed in this respect: with η 0,BK = 0 setups the C 2 values are quite a bit larger but then with η 0,BK = ln 1 0.01 the HERA data fit is found to require only a small C 2 .
Comparing the initial conditions with η 0,BK = 0 and η 0,BK = ln 1 0.01 we see that every evolution starts from a significantly larger anomalous dimension when η 0,BK = 0. This difference is more pronounced compared to the previously studied KCBK and ResumBK equations. This is because the TBK evolution drives the dipole towards the asymptotic shape with a small anomalous dimension γ ∼ 0.6 [23]. This behavior is similar to the leading order BK equation, in which case it is already known that the asymptotic shape can not be used to parametrize the initial condition [8] 5 . Indeed the development of the geometric scaling regime independently of the initial condition is a theoretically attractive feature of the TBK formulation. However, HERA data seems to prefer to lie in the preasymptotic regime in the fits. Thus, especially in the fits with more evolution before the data region (smaller η 0,BK ), one needs to slow down the evolution more and start with a significantly larger anomalous dimension in order to still have a transient form of the dipole amplitude in the data regime.
The evolution of the dipole slope in TBK evolution is shown in Fig. 4 at different evolution rapidities η. Unlike in the case of KCBK equation discussed earlier and shown in Fig. 3, the slope of the dipole from the TBK evolution is decreasing with both running couplings in the r ∼ 1/Q s regime. As shown in Ref. [23], the asymptotic anomalous dimension γ ∼ 0.6 is obtained only at very large rapidities, and at least the Bal+SD coupling case can be seen to be evolving towards this asymptotic value. In the rapidity range relevant in HERA or even LHeC kinematics the asymptotic anomalous dimension is not reached. With the parent dipole coupling, a significantly longer evolution is needed before evolution towards the asymptotic shape at small r becomes visible in the best fit case with a small anomalous dimension γ = 0.9 in the initial condition.
The dipole amplitudes at different evolution rapidities as a function of dipole size are shown in Fig. 5. Here, results obtained using the all three considered evolution equations are shown at fixed projectile rapidity Y = Y 0,BK + ∆Y . The solution to the TBK evolution is shifted from the target rapidity η to the projectile rapidity Y by performing the shift (29). The shifted TBK solutions are shown in the region where η > η 0,BK . In the region where the dipole amplitude is not small, all evolution equations result in comparable dipole amplitudes. This is expected, as all the shown dipoles result in a compatible description of the HERA structure function data. At small dipole sizes that do not significantly contribute to the structure functions some differences appear. Despite the fact that KCBK and ResumBK equations have very similar initial conditions the resulting amplitudes differ significantly for small dipoles. This is mostly driven by the resummation of the single transverse logarithms not included in the kinematically constrained BK equation, as this resummation is more important at small parent dipole size r. At very small dipoles the TBK evolved dipole also differs significantly from the other dipoles when the shift from target rapidity, Eq. (28), results in the dipole being evaluated close to the initial condition. If the parent dipole scheme for the running coupling were used, the differences between the dipoles obtained from the different evolution equations would be significantly reduced, as in that scheme the coupling constant is generically smaller at small r and differences between the evolution equations are suppressed by the small α s .

B. Fitting the interpolated light quark reduced cross section
Next we consider fits to our interpolated light quark data set. The fit results are also shown in Tables II and III. Figure 6 shows a comparison between the HERA and interpolated light quark data with one of the fits, obtained with the KCBK equation with the Balitsky + smallest dipole running coupling and initial condition parametrized at Y 0,BK = ln 1/0.01. The light quark only fits have quite distinct systematics in comparison to the actual HERA data fits. Every single fit setup used needs a substantially larger C 2 and to a varying degree larger anomalous dimensions. Lastly, and importantly, light quark fits need larger values of σ 0 compared to the corresponding total HERA cross section fit.
The slow evolution speed (visible as a large C 2 especially when using the parent dipole prescription) and a large σ 0 in the light quark pseudodata fits can be understood to result from an effective description of nonperturbative effects. We expect that there is a nonperturbative hadronic contribution in the light quark production cross section which is large (resulting in a large σ 0 ) and evolves more slowly as a function of Bjorken-x Bj than the fully perturbative cross sections, like charm production. In our framework, these nonperturbative effects correspond to large dipoles, with sizes larger than roughly the inverse pion mass. In this case, quark-antiquark dipoles are not the right degrees of freedom, and one should in principle use an another effective desription for the non-perturbative physics, e.g. the vector meson dominance [96][97][98][99] model.
The same non-perturbative effects are there also in the total reduced cross section, and consequently in our fits to full HERA data. However, the full reduced cross section also includes the more reliably perturbative charm production contribution (and a small b quark one), with a much faster x evolution and a smaller magnitude (σ 0 ). Consequently, when performing our (massless) NLO fits to the full HERA data more weight is given to perturbative contributions compared to light quark fits, and there is less need for the fit parameters to adjust to nonperturbative effects with unnatural values.
These observations are compatible with some of the previous analyses. In the study by the AAMQS collaboration [8] it was found that a combined fit to both charm and total reduced cross section requires one to introduce separate fit parameters for the charm quarks, especially the charm quarks require a smaller σ 0 . A slowly evolving non-perturbative contribution to the light quark production was also found to be necessary in Refs. [29,40]. In the dipole picture applied here, one finds that very large dipoles up to a few femtometers contribute significantly to the light quark structure function [83]. In reality, non-perturbative confinement scale effects not included in our perturbative calculation are expected to dominate in these cases as discussed above.
To arrive at one of our central points of this article, we make the observation that even though the HERA DIS data has been described well with leading order dipole picture fits with the BK equation in the past, simultaneous fits to the full data and charm quark data have not been successful with a single BK-evolved amplitude (note however the existence of fits [93,100,101] using parametrizations that mimic BK evolution). Similar results are found in the recent study with the target rapidity BK prescription as well [81]: fits to the full data are excellent but the fit parametrizations do not describe the heavy quark data. Our next-to-leading order analysis, where we separately consider the light quark production only, results in similar conclusions. This indicates that the description of the light quark contribution has a large theoretical uncertainty as well in any such fit to the full DIS data.
Thus we find that it would be preferable to fit the charm quark structure function F 2,c separately (or inclusive F L data, as the longitudinal photon splits generally to smaller dipoles, resulting in smaller non-perturbative contributions). The F L measurements from HERA [102] are however not precise enough for our purposes (see the next section). Very precise F L data (among with inclusive and charm structure functions) can be expected from the future Electron Ion Collider [103,104] or from the LHeC [105].

C. Beyond HERA
Given the equality in the capabilities of the different versions of the BK equation in describing the HERA and light quark data, a question arises if it is possible to distinguish the different fit schemes and find the preferred form of the BK equation. In general, one might expect to see differences in the Q 2 dependence of the structure functions at small x (in the HERA kinematics, the fit procedure ensures a compatible evolution). This is because the Q 2 dependence is controlled by the anomalous At asymptotically small x both approaches can result in the same Q 2 dependence of the cross section in spite of the different anomalous dimensions. This can be seen as follows. Let us first consider the BK equation formulated in the target rapidity, and write the dipole amplitude as N ∼ (Q 2 s r 2 ) γ . The TBK equation results in the saturation scale scaling as Q 2 s ∼ x −λ , as the evolution range is ln 1/x Bj . This gives N ∼ x Bj −λγ (Q 2 ) −γ , and consequently the structure functions behave as where we have scaled out the Q 2 dependence originating from the virtual photon wave function ψ γ * →qq . Substituting an asymptotic anomalous dimension γ ∼ 0.7 we get On the other hand, when applying the KCBK or ResumBK equations formulated in terms of the projectile momentum fraction, the evolution range is controlled by ln W 2 = ln(Q 2 /x Bj ). Consequently, we get In general, in the case of ResumBK and KCBK we expect γ ∼ 1 as the evolution does not change the asymptotic anomalous dimension. Using λ ∼ 0.3 for the generic evolution speed we get which is the same Q 2 scaling as obtained in case of TBK equation, see Eq. (37).
In practice, however, in HERA or even LHeC kinematics the TBK evolution has not reached its asymptotic form, and the anomalous dimension is still close to unity as shown in Fig. 4. Consequently, the Q 2 -dependence is expected to be slower in the TBK evolution in realistic kinematics. We note that the structure functions are not actually sensitive to the slope of the dipole at asymptotically small r but in the region r ∼ 1/Q s or r ∼ 1/Q, which makes it in practice difficult to compare Q 2 dependences analytically. We also note that when computing the structure function at low x Bj , also dipole amplitudes at higher x Bj are probed when performing the z 2 integral.
The numerically calculated Q 2 dependence of the structure functions F 2 and F L is shown in Fig. 7. The results are shown at small x Bj = 5.6 · 10 −5 corresponding to the LHeC kinematics using each of the BK equations, employing the fit to the full HERA data with the Bal+SD running coupling prescription and Y 0,BK = η 0,BK = ln 1 0.01 . For comparison, the leading order result based on Ref. [9] is shown. Compared to the leading order fit, the Q 2 dependence is weaker at next-to-leading order, due to the different asymptotic shape of the dipole amplitude (the leading order BK equation develops a small anomalous dimension γ which results in faster Q 2 dependence).
The different fit schemes that result in an equally good description of the HERA data start to differ slightly at large Q 2 when considering the Bjorken-x Bj region not included in the fits. The longitudinal structure function F L is more sensitive to small dipole sizes, and as such it can be expected to be more sensitive on the details of the evolution. This is especially visible when the ResumBK evolution is compared to other approaches: the Q 2 dependence is much weaker at large Q 2 . This is due to the resummation of single transverse logarithms not included  in other evolution schemes, which has the largest effect at small parent dipole sizes probed at large Q 2 . However, in the realistic kinematical range considered here, the difference between the fits is moderate. This suggests that our next-to-leading order predictions for the structure functions in the future collider experiments are robust. Future high energy DIS data from e.g. LHeC will be extremely precise, with the expected uncertainties in the structure function measurements being even at the per mill level [105]. As such one could be sensitive to details in NLO BK evolution, even though the effects are not large. Ultimately more differential measurements in addition to the reduced cross section will be needed. The most precise measurement of the proton longitudinal structure function F L up to date is performed by the H1 collaboration at HERA [102] (with compatible results obtained by the ZEUS collaboration [106]). In Fig. 8 we compare the F L computed from our fits to the H1 data. Due to the limited statistics, the most precise results are not reported as a function of both x and Q 2 , but at fixed x, Q 2 combinations. Consequently, it is crucial to note that the higher Q 2 points are measured at higher x. All three fit setups result in almost identical F L , as expected as the F L is measured in the kinematical domain mostly included in our reduced cross section fits. Even though the future Electron Ion Collider [103] will not reach as small Bjorken-x values as the LHeC, the F L measurements it can perform will be very useful as the HERA measurement has large uncertainties and it only covers a small fraction of the phase space where the details of the evolution can not be accessed.

VI. CONCLUSIONS
We have performed, for the first time, a fit to the HERA structure function data in the Color Glass Con-densate framework at next to leading order accuracy in the case of massless quarks. As the full next to leading order BK equation is computationally demanding, we approximate it by employing evolution equations that resum higher order corrections enhanced by large transverse logarithms. As a result of the fits, we obtain the initial condition for the perturbative BK evolution. The resulting dipole-target scattering amplitude can be used in other phenomenological applications, for example when calculating particle production in proton-nucleus collisions at next to leading order in α s .
Similarly as in the leading order fits previously studied in the literature, we find that it is possible to obtain an excellent description of the precise combined HERA structure function data. Equally good fits are obtained when using both the BK equation formulated in terms of the projectile momentum fraction, and the recently proposed BK equation where the evolution rapidity is dictated by the fraction of the target longitudinal momentum. When extrapolated to LHeC energies, the different BK evolution prescriptions are found to result in moderate differences in the Q 2 dependence of the structure functions. This suggests that the NLO calculation presented here is robust, and has a strong predictive power for future DIS measurements in new experimental facilities such as the EIC or LHeC.
As next-to-leading order impact factors for massive quarks are not yet available, it is not possible to compute charm and bottom contribution to the structure functions. To perform consistent fits, we also generated an interpolated light quark data set by subtracting the interpolated charm and bottom contribution from the HERA reduced cross section data. Fits to this light quark data require a much slower Bjorken-x evolution than we naturally get from the perturbative evolution equations applied. Additionally, the apparent proton transverse size obtained is significantly larger than seen when fitting the full HERA data. These features we interpret to result from a non-perturbative hadronic component in the light quark production cross section. This component is large (resulting in a large proton transverse area), and evolves more slowly as a function of Bjorken-x, as expected for a hadronic component.
Our results demonstrate the need for massive quark impact factors at next-to-leading order accuracy in the CGC framework, which would allow fits to fully perturbative charm cross section separately. Precise measurements of the charm structure function over a wide range of x and Q 2 , in addition to the longitudinal structrue function, from future experiments will also be useful. The fits to the generated light quark data should in principle be considered our principal preferred fits as there the agreement between the data and the massless theory should be on the most solid footing. However, if used for QCD phenomenology in other observables where the presumed non-perturbative contribution is smaller, the best one can do is use the full HERA data fits.
In addition to inclusion of the quark masses and the usage of the full NLO BK, the NLO DIS calculation can be improved by relaxing some of the kinematical assumptions. First, in addition to the gluon momentum fraction z 2 , the quark momentum fraction z 1 should not be allowed to get arbitrary close to endpoints z 1 → 0, 1 in order to avoid production of qq pairs with invariant mass larger than the center-of-mass energy. Additionally, in the virtual correction one should also perform the integral over the gluon momentum fraction and evaluate the dipole operator at the same rapidity as in the real term. This would make it possible to also consistently include a Q 2 dependent evolution range in the virtual contribution. Finally, when the Balitsky prescription for the running coupling is used in the BK evolution, there is a mismatch in the running coupling schemes between the impact factor and the evolution equation which could be improved. We plan to address these issues in future work.