Collinearly-improved BK evolution meets the HERA data

In a previous publication, we have established a collinearly-improved version of the Balitsky-Kovchegov (BK) equation, which resums to all orders the radiative corrections enhanced by large double transverse logarithms. Here, we study the relevance of this equation as a tool for phenomenology, by confronting it to the HERA data. To that aim, we first improve the perturbative accuracy of our resummation, by including two classes of single-logarithmic corrections: those generated by the first non-singular terms in the DGLAP splitting functions and those expressing the one-loop running of the QCD coupling. The equation thus obtained includes all the next-to-leading order corrections to the BK equation which are enhanced by (single or double) collinear logarithms. We then use numerical solutions to this equation to fit the HERA data for the electron-proton reduced cross-section at small Bjorken x. We obtain good quality fits for physically acceptable initial conditions. Our best fit, which shows a good stability up to virtualities as large as Q^2=400 GeV^2 for the exchanged photon, uses as an initial condition the running-coupling version of the McLerran-Venugopalan model, with the QCD coupling running according to the smallest dipole prescription.


Introduction
The wealth of data on electron-proton deep inelastic scattering collected by the experiments at HERA over 15 years of operation has allowed for stringent tests of our understanding of high-energy scattering from first principles. This refers in particular to the 'small-x' regime where perturbative QCD predicts a rapid growth of the gluon density with increasing energy (or decreasing Bjorken x), leading to non-linear phenomena like multiple scattering and gluon saturation [1,2]. The simplicity of the dipole factorization for deep inelastic scattering at high energy [3,4] has favored the emergence of relatively simple 'dipole models', in which the high-density effects are efficiently implemented as unitarity corrections to the crosssection for the scattering between a quark-antiquark dipole and the proton. Such models allowed for rather successful fits to the small-x HERA data at a time where the theory of the non-linear evolution in QCD was insufficiently developed and the pertinence of gluon saturation for the phenomenology was far from being widely accepted. The first such model -the "GBW saturation model" [5,6] -provided a rather good description of the early HERA data for the inclusive and diffractive structure functions at x ≤ 10 −2 with only 3 free parameters. This success inspired new ways to look at the HERA data, which in particular led to the identification of geometric scaling [7]. The subsequent understanding [8][9][10] of this scaling from the non-linear evolution equations in QCD -the Balitsky-JIMWLK hierarchy [11][12][13][14][15][16][17] and its mean field approximation known as the Balitsky-Kovchegov (BK) equation [18] -has greatly increased our confidence in the validity of the pQCD approach to gluon saturation as a valuable tool for phenomenology. BK equation. Besides, our equation is local in 'rapidity' (the logarithm of the energy, which plays the role of the evolution variable), a property which in this context is rather remarkable since the physics behind the double collinear logarithms is the time-ordering of subsequent, soft, gluon emissions, which is genuinely non-local. 1 The first numerical studies of this collinearly-improved BK equation demonstrate the essential role played by the resummation in both stabilizing and slowing down the evolution [50,57].
In this paper, we shall provide the first phenomenological test of our resummation scheme, by using it in fits to the inclusive HERA data. To that aim, it will be important to first extend this scheme to also include the single transverse logarithms which appear in the NLO correction to the BK equation -that is, the NLO terms expressing the first correction to the DGLAP splitting kernel beyond the small-x approximation and those associated with the one-loop running of the coupling. Indeed, such single-log effects must be kept under control to ensure a good convergence of the perturbative expansion. Besides, the inclusion of running coupling effects is essential for the description of the data, as well known.
The resummation of the DGLAP logarithms to the order of interest turns out to be rather straightforward: it amounts to adding an anomalous dimension (a piece of the leading-order DGLAP anomalous dimension) to both the resummed kernel and the resummed initial condition. For the running coupling corrections, the situation turns out to be more subtle since, strictly speaking, they cannot be encoded into an equation which is local in rapidity. This being said, and following the standard strategy in the literature, we shall propose various schemes for introducing a running coupling directly in the local evolution equation and test these schemes via fits to the HERA data.
After these additional resummations, we are led to a new, more refined, version for the 'collinearly improved BK equation', namely Eq. (9) below, which will be our main tool for phenomenology. By construction, this equation resums the double-logarithmic corrections completely -meaning to all orders inᾱ s ≡ α s N c /π (α s is the QCD coupling and N c is the number of colors) and with the right symmetry factors -, whereas the single-logarithmic terms are resummed only partially (but in such a way to include the respective terms to NLO). It is rather straightforward to extend our resummed equation to full NLO accuracy, by adding the remaining corrections of O(ᾱ 2 s ), as computed in [33]. But the ensuing equation would be very cumbersome to use in practice, due to the intricate, non-local and non-linear, structure of the pureᾱ 2 s corrections. In this first analysis, we shall adopt the viewpoint that the most important higher-order contributions (say, in view of phenomenology) are those enhanced by collinear logs, as explicitly resummed in Eq. (9), and that the pureᾱ 2 s effects are truly small and can be effectively taken care of via the fitting procedure. A similar viewpoint has been advocated in previous fits based on rcBK, but given the importance of the collinear logarithms, this assumption was not so well motivated and led indeed to some tensions in the respective fits, as we shall later explain.
Using numerical solutions to this collinearly improved BK equation together with suitable forms for the initial condition, we have performed fits to the HERA data for the ep reduced cross-section [38] at x ≤ 10 −2 and Q 2 ≤ Q 2 max , where the upper limit Q 2 max on the virtuality Q 2 of the exchanged photon is varied between 50 GeV 2 (a common choice in small-x fits) and 400 GeV 2 . These fits show several remarkable characteristics.
(i) The fits are indeed successful: for Q 2 max = 50 GeV 2 and two types of initial conditions -GBW-like [5] and the running-coupling version of the McLerran-Venugopalan (rcMV) model [58] -, we obtain a χ 2 per number of data points around 1.2 with only 4 free parameters.
(ii) The fits are also very discriminatory: they clearly favor some initial conditions over some others, and some prescriptions for the running of the coupling over the others. For instance, the standard MV initial condition, which truly corresponds to a fixed coupling, appears to be disfavored, whereas a more physical version of it, including a running coupling, works quite well. The latter works also better than the GBW initial condition, in the sense that it provides a fit which remains stable up to Q 2 = 400 GeV 2 .
(iii) Our fits alleviate some tensions (in terms of physical interpretation) which were visible in previous first based on rcBK [21][22][23]36] and could be attributed to the choice to replace all the NLO corrections with the running of the coupling alone (see also the related discussion in [23]). Notably, our fits prefer prescriptions where the QCD coupling α s (µ 2 ) is running according to the smallest dipole size, they do not require any artificial 'anomalous dimension' in the initial condition, and treat the heavy quarks on the same footing as the light ones, in agreement with general expectations from the dipole factorization.

The NLO BK equation and large transverse logarithms
To motivate the resummations that we shall later perform, let us first explicitly exhibit the large transverse logarithms which appear when computing the NLO corrections to the BK equation [31][32][33]. We recall that the BK equation describes the rapidity evolution of the S-matrix S xy = 1 − T xy for the scattering of a color dipole with transverse coordinates (x, y) off a hadronic target. The dipole scattering amplitude T xy is small in the regime where the target is dilute, but it approaches the unitarity (or 'black disk') limit T xy = 1 when the target is dense. The separation between these two regimes is controlled by the saturation momentum Q s (Y ), which increases with the rapidity difference Y between the projectile and the target.
Neglecting the terms suppressed in the limit of a large number of colors N c 1, one finds a closed equation for S xy , whose NLO version reads as follows [33] dS xy dY =ᾱ where N f is the number of flavors,b = (11N c − 2N f )/12N c , andᾱ s = α s N c /π, with the QCD coupling α s evaluated at the renormalization scale µ.
There are two main changes in the structure of the evolution equation as we go from LO to NLO. First, the term with a single integration (SI) over the transverse coordinate z only receives a correction of order O(ᾱ 2 s ) to the kernel, which in particular contains the running coupling corrections proportional tob. Second, there are new terms, of order O(ᾱ 2 s ), which involve a double integration (DI) over the transverse coordinates u and z and which refer to partonic fluctuations involving two additional partons (besides the original quark and antiquark) at the time of scattering. The first such a term, which is independent of N f , represents fluctuations where both daughter partons are gluons. The S-matrix structure therein, that is, S xu S uz S zy − S xu S uy , corresponds to the following sequence of emissions: the original dipole (x, y) emits a gluon at u, thus effectively splitting into two dipoles (x, u) and (u, y); then, the dipole (u, y) emits a gluon at z, thus giving rise to the dipoles (u, z) and (z, y). The 'real' term S xu S uz S zy describes the situation where both daughter gluons interact with the target. The 'virtual' term −S xu S uy describes the case where the gluon at z has been emitted and reabsorbed either before, or after, the scattering. This negative 'virtual' term subtracts the equal-point contribution (z = u) from the 'real' piece, ensuring that the potential 'ultraviolet' singularity associated with the factor 1/(u − z) 4 in the kernel is truly harmless. A similar discussion applies to the second DI term, proportional to N f , except for the fact that the additional partons at the time of scattering are a quark and an antiquark.
In principle, one should be able to undertake the task of solving the NLO BK equation. The hope would be that the solution would only add a relatively small correction to the LO result. However, this is not the case since there are terms in the kernels of the NLO equation which can become large in certain kinematic regimes and thus invalidate the strictᾱ s -expansion. One obvious class of such terms contains the corrections proportional tob in the SI term in Eq. (1), which by themselves bring no serious difficulties: as well known, these corrections can be absorbed into a redefinition of the scale for the running of the coupling, which thus becomes a dynamical scale (see Sect. 4 below for details). Here, we would like to focus on the corrections enhanced by 'collinear logarithms', that is, logarithms associated with the large separation in transverse sizes (or momenta) between successive emissions. These corrections become large only in the weak-scattering regime where all the dipoles are small compared to the saturation scale 1/Q s (Y ) and the equation can be linearized w.r.t. to the (small) scattering amplitude T . This in particular means that one can ignore the last term, proportional to N f /N c , in Eq. (1) since this term vanishes after linearization, as one can easily check (by also using the symmetry of the kernel under the interchange u ↔ z).
To be more precise, let us consider the strongly ordered regime that is, the parent dipole is the smallest one, a gluon is emitted far away at u, a second one even further at z, but with all possible dipole sizes remaining smaller than the inverse saturation momentum. Whenever appropriate, we will denote by r,ū andz the size of the parent dipole, the size of the dipoles involving u and the size of the dipoles involving z, respectively, with r 2 ū 2 z 2 . By inspection of the SI piece in the NLO BK equation, it is quite obvious that the dominant NLO term is the one involving a double transverse logarithm (DTL), that is, the last term within the square brackets. Still within this regime (2), we can approximate the scattering matrices in the SI term as follows: where the second approximate equality follows since the dipole amplitude for a small dipole is roughly proportional to the dipole size squared. Notice that the net result in the approximation of interest fully comes from the 'real' term, which involves the large daughter dipoles.
What is not immediately obvious is the presence of a single transverse logarithm (STL) coming from the DI term. Let us isolate here the relevant part of the kernel, To implement the limit in Eq. (2) we can successively write the expression in Eq. (3) as with φ the angle between r and any of the two dipoles involving u. To obtain (4), we have first set all dipole sizes which include z equal to each other, since any subleading term would be suppressed by inverse powers ofz. Then the only z dependence left is the one explicit in the prefactor. We have subsequently taken the limit r ū (by expanding the logarithm to cubic order) and we have finally averaged over the angle φ between the parent dipole and those involving u. Notice that the would-be leading term, of order 1/z 4 , has cancelled out in these manipulations. The first non-vanishing term, as visible in the r.h.s. of Eq. (4), is suppressed by r 2 /ū 2 , thus creating the conditions for a logarithmic integration overū. To explicitly see this, recall that we consider the weak-scattering regime, where the product of S-matrices multiplying M STL can be linearized. This allows us write S xu S uz S zy − S xu S uy −T uz − T zy + T uy −2T (z). (Once again, the dominant contribution has been generated by the 'real' term.) We see that the net scattering amplitude in this approximation is independent of the intermediate dipole sizeū. Accordingly, when integrating over u, within the range limited by r andz, we find a STL, as anticipated. After also including the LO term and the NLO one enhanced by the DTL, one finds that the NLO BK equation in the strongly ordered region (2) reduces to It is now clear that, if the daughter dipoles are allowed to become sufficiently large, the NLO contributions enhanced by large transverse logarithms become comparable to, or larger than, the LO one. In that case, the present perturbative expansion cannot be trusted anymore. To be more explicit, consider a single step ∆Y in the evolution with the following, simple but physically meaningful, initial condition Thez-integration in Eq. (5) becomes logarithmic and gives This shows that, for sufficiently small rQ s , such thatᾱ s ln 2 (1/r 2 Q 2 s ) 1, the NLO correction becomes larger than the LO term and the perturbation series is unreliable. In particular, the NLO correction is negative rendering the solution unstable, as indeed observed in numerical solutions [48][49][50].

Resumming the large collinear logarithms
The large NLO corrections that we have singled out in the previous section are the lowest-order examples of collinearly-enhanced radiative corrections, which occur to all orders and spoil the convergence of the perturbation theory. When the separation between the transverse scales of the projectile and the target is large enough, as is actually the case in the DIS kinematics at HERA, the higher-order terms of this type become more important than the pureᾱ 2 s NLO terms (i.e. the contributions of O(ᾱ 2 s ) which are not amplified by any transverse logarithms). From now on, we shall focus on this situation, discarding the purē α 2 s NLO corrections, but focusing on the resummation of the large transverse logarithms to all orders. In this section, we consider both the single and double collinear logarithms, thus following and expanding our recent results in [50]. In the next section, we shall explain how the running coupling corrections can be included in this scheme.
In Ref. [50] we have devised a strategy for resumming double-logarithmic corrections to the BK equation to all orders. Our main observation was that these corrections are generated by the diagrams common to the BFKL and DGLAP evolutions -i.e. the Feynman graphs of light-cone perturbation theory in which the successive gluon emissions are strongly ordered in both longitudinal momenta and transverse momenta (or 'dipole sizes') -after enforcing the additional constraint that the emissions must also be ordered in lifetimes (or, equivalently, in light-cone energies [56]). Concerning the single collinear logarithms, it is intuitively clear that they must represent DGLAP-like corrections to BFKL, 'small-x', emissions. For instance, the effect of orderᾱ 2 s ∆Y ρ 2 visible in Eq. (7) is the result of a sequence of two emissions: one small-x emission (in the double logarithmic regime) yielding a contribution ∝ᾱ s ∆Y ρ, and a DGLAP-like emission, characterized by strong ordering in dipole sizes (see eq. Eq. (2)) and which gives an effect of orderᾱ s ρ. Here, ρ ≡ ln(1/Q 2 s r 2 ) is the large logarithm generated by the transverse phase-space. This scenario is corroborated by the following observation: the numerical coefficient A 1 ≡ 11/12 in front of the STL in Eq. (5) can be recognized as the second-order term in the small ω expansion of the relevant linear combination of DGLAP anomalous dimensions: Recalling that one needs one factor of 1/ω in order to generate a small-x logarithm ∆Y = ln(1/x), one sees that the NLO effect ∼ᾱ 2 s ∆Y ρ 2 is indeed produced by combining the singular (1/ω) piece of one emission with the first non-singular piece (A 1 ) of another one. This discussion also instructs us about the strategy to follow in order to resum such STLs to all orders: it suffices to include this piece A 1 as an anomalous dimension, i.e. as an extra power-law suppression, in the evolution kernel previously obtained in Ref. [50]. The "fac" scheme as given in Eq. (12). Blue (dotted): The Balitsky scheme [32]. In all cases the parent dipole size is |x − y| = 1, the coupling is smoothly frozen at the value 0.7 and Λ QCD = 0.2.
We are thus led to the following, collinearly-improved, version of the BK equation, where the overall kernel is written as a product of three factors: the familiar dipole kernel which appears already at leading order, the 'DLA kernel', resuming the double collinear logs to all orders [50] K DLA (ρ) = J 1 2 ᾱ s ρ 2 evaluated at ρ = L xzr L yzr , with L xzr ≡ ln[(x − z) 2 /r 2 ], and a new factor, which features the exponent ±ᾱ s A 1 (the positive sign in the exponent is taken when |x−y| < min{|x−z|, |y−z|} and the negative sign otherwise), which expresses the contribution of the single collinear logarithms. From the above discussion, is should also be clear that the present resummation of STLs is only partial: it refers to the particular class of such corrections which are generated by the first non-singular piece in the expansion in Eq. (8). The higher terms in this ω-expansion will produce single collinear logarithms too, but only starting at higher orders in perturbation theory (NNLO or higher). At the level of the BFKL equation, more complete resummations of the single logarithms have been devised in [52][53][54], but so far it is not clear how to extend these resummation schemes to a non-linear evolution equation like BK.
Returning to Eq. (9), the tilde symbol inT xy is intended to remind that this is truly a suitable analytic continuation of the dipole amplitude which coincides with the physical quantity T xy only for ρ < Y . For ρ > Y , the physical amplitude can be obtained by either solving an equation non-local in Y , or by matching onto the solution to the DGLAP equation [50]. However, explicit numerical studies at DLA level have shown that the solutionT xy to Eq. (9) remains very close to the actual physical amplitude, including for ρ > Y . For this reason, we shall ignore this subtlety (and the related issue of the resummation in the initial condition) for the purpose of the fits to be constructed in the next section. We shall return to a more detailed study of these issues in a forthcoming publication [59].

Prescriptions for the running of the coupling
The last source of potentially large NLO corrections to the BK equation are the running coupling corrections, i.e. the logarithmic terms proportional tob in the SI term in Eq. (1). Such terms can grow large when the scales in their arguments are very disparate. More precisely, the first logarithm can be problematic when r is much smaller or much larger than 1/µ, while the second when the soft gluon at z is collinear to the quark or the antiquark composing the parent dipole. We need to choose µ in such a way to cancel these potentially large logarithms, which could otherwise spoil the convergence of the perturbative expansion 2 . It is clear that there is not a unique choice, but in QCD one usually expects the hardest scale to determine the running of the coupling. Indeed, a quick inspection shows that the smallest dipole prescription α min =ᾱ s (r min ) with r min = min{|x−y|, |x−z|, |y−z|} (11) cancels the large logarithms in all kinematic regions. Another possibility is to choose µ so that all the terms with coefficientᾱ 2 sb vanish. Given that in the current work we neglect all finite (i.e. not enhanced by a large logarithm)ᾱ 2 s terms, this looks like what is called the "fastest apparent convergence" (fac) scheme [60][61][62]. It is convenient in the sense that one is left with just the leading term inᾱ s . We find that and it is an easy exercise to show that it reduces to the minimal dipole choiceᾱ s (r min ) in all limits where one of the three dipoles is much smaller than the other two.
In this work, we shall use both above schemes. Let us add that the most popular prescription, widely used so far in phenomenological applications, is the one due to Balitsky [32], and reads but it will not be adopted here for a number of reasons. First, it is based on an extrapolation to all orders of a coordinate space kernel which includes theᾱ 2 sb terms above as well asᾱ 3 sb 2 corrections. At this order, we would also expect corrections proportional to the two-loop beta function. Second, even though it also reduces toᾱ s (r min ) in the extreme kinematical limits, it does that very slowly for large daughter dipoles (in certain configurations) and this leads to an unphysically small coupling in a large region of phase space, as can be seen in the respective plots in Fig. 1. Finally, and perhaps as a result of the above drawbacks, when used in fitting the DIS data, it gives a much worse fit than Eqs. (11) and (12) and with fit parameters which take somewhat unnatural values.

Fits to the HERA data
We now turn to the description of the HERA reduced cross-section measurements using the resummed BK equation. To this aim several ingredients first have to be specified.
Initial condition. We must fix the initial condition for the dipole amplitude at some Y 0 , which afterwards will be evolved towards higher rapidities using Eq. (9). We consider two choices: the simple parametrisation of the Golec-Biernat and Wüsthoff (GBW) [5] type and the running-coupling version of the McLerran-Venugopalan (rcMV) model [58] T  Table 1: χ 2 and values of the fitted parameters entering the description of the HERA data. The fit includes the 252 σ red data points. The quoted χ 2 for σ cc red and F L are obtained a posteriori.
It is worth noticing that, as dictated by collinear physics, there is no anomalous dimension in the above initial conditions. The extra parameter p determines the shape of the amplitude close to saturation and its approach towards unitarity.
Running coupling. We consider the two prescriptions given by Eqs. (11) and (12). For the explicit expression of the strong coupling in coordinate space in terms of r we introduce a fudge factor as in [22], namely with . This fudge factor is also included in the rcMV type initial condition in (15). The N f -dependent Landau pole is obtained by imposing α s (M 2 Z ) = 0.1185 at the scale of the Z mass [63] and continuity of α s at the flavour thresholds, using m c = 1.3 GeV and m b = 4.5 GeV. To regularise the infrared behaviour, we have decided to freeze α s at a value α sat = 1 and we have checked explicitly that reducing this down to, for example, 0.7 does not affect the fit in any significant manner.
Note that we do not include any form of resummation or matching for ln 1/r 2 > Y , as introduced in [50], in these initial conditions. One of the reasons for not doing so is that the extra factor in the initial condition can always be reabsorbed in a re-parametrisation. Furthermore, a proper matching at small dipole sizes, suited for phenomenological studies, would require a careful treatment of the small-dipole region. In that respect, the resummed BK evolution is expected to perform a better job than a fixed matching with a fixed asymptotic behaviour. We leave a better treatment, e.g. a genuine matching to DGLAP evolution, for future work.
Rapidity evolution. Of course this is determined by the resummed BK equation given in (9). Here, we again consider two separate cases, one in which the evolution resums only the leading double logarithms and one in which it also includes the single ones.
From the dipole amplitude to observables. Once we have the dipole amplitude for all rapidities and dipole sizes, we use the standard dipole formalism to obtain the physical observables: where the transverse and longitudinal virtual photon wavefunctions read  Table 2: Evolution of the fit quality when including data at larger Q 2 (in GeV 2 ).
In the above we have introduced the customary notationQ 2 , and we have assumed a uniform distribution over a disk of radius R p in impact parameter space. The sum in (17) runs over all quark flavours and we will include the contributions from light quarks with m u,d,s = 100 MeV as well as from the charm quark with m c = 1.3 GeV. From the longitudinal and transverse cross-sections, we can deduce the reduced cross-section and the longitudinal structure function as When the quark masses, the value of the strong coupling at the Z mass and its frozen value in the infrared have been fixed, we are left with 4 free parameters according to our choice of initial condition: R p the "proton radius", Q 0 the scale separating the dilute and dense regimes, C α the fudge factor in the running coupling in coordinate space, and p which controls the approach to saturation in the initial condition.
We have fitted these parameters to the combined HERA measurements of the reduced photon-proton cross-section [38]. Since the BK equation is applicable only at small-x, we have limited ourselves to the region x ≤ 0.01. We note that since Eq. (17) probes dipoles at the rapidity ln 1/x f , the exact cut we impose isx c ≤ 0.01 since the most constraining cut comes from the charm, the most massive quark we include in our model. Accordingly, our initial condition for the BK evolution corresponds tox = 0.01. Furthermore, since we do not expect the BK equation to capture the full collinear physics, we impose the upper bound Q 2 < Q 2 max . By default we will use Q 2 max = 50 GeV 2 but we will also give results for extensions to larger Q 2 . In the default case we have a total of 252 points included in the fit. We have added the statistical and systematic uncertainties in quadrature. 3 The results of our fits for the 2 3 = 8 cases, depending on the initial condition, the running coupling prescription and the inclusion or not of single logarithms in the kernel, are presented in Table. 1. The table includes the parameter values obtained from fitting the σ red data and, besides the fit χ 2 , it also indicates the χ 2 obtained a posteriori for the latest σ cc red [64] and F L [65] measurements. These results deserve a few important comments.   (iii) The two initial conditions give similar results, with a slight advantage for the rcMV option, which is likely due to the extra parameter. Note that for a standard MV-type of initial condition T (Y 0 , r) = {1 − exp[−(r 2 Q 2 0 /4 [c + ln(1 + 1/rΛ)]) p ]} 1/p , we have not been able to obtain a χ 2 per point below 1.3 and the parameters, typically c or p, tend to take unnatural values. (iv) As far as the running-coupling prescription is concerned, the smallest dipole prescription given in Eq. (11) tends to give somewhat better fits than the "fac" prescription given in Eq. (12). This can be seen as an estimate of subleading corrections (including the pureᾱ 2 s NLO terms) that we neglect in the present fit. Note also that we have not been able to reach a fit of equivalent quality and robustness with the Balitsky prescription.
(v) The resummation of the single logarithms tends to yield slightly larger values for χ 2 , but the difference is too small to be significative (at least, without performing a full NLO analysis). Perhaps more significantly, this resummation leads to more physical values for some of the parameters, especially C α for the smallest dipole prescription and p for the GBW initial condition. These findings are consistent with the expectation that, once properly resummed, single logarithms should have only a modest impact. Recall however that their resummation is a crucial step towards a full NLO fit -failing to do so could lead to instabilities similar to those observed when double logarithms are not resummed. (vi) The fit remains stable when varying the parameters we have imposed by hand. For example, using α sat = 0.7 instead of 1 has no significant effect on the fit. Varying the light quark masses within the rather wide range 0 ≤ m u,d,s ≤ 140 MeV only slightly changes the quality of the fit. For instance, taking one of our best fits (rcMV initial condition, the smallest dipole prescription for the running of the coupling, and resummation of the single logarithms), we have found χ 2 = {1.180, 1.153, 1.126, 1.159} when choosing m u,d,s = {0, 50, 100, 140} MeV, respectively. This lack of sensitivity to the light quark masses is likely a consequence of saturation, which reduces the dependence of the DIS cross-section to very large dipole fluctuations. (The corresponding amplitudes reach the unitarity, or 'black disk', limit T = 1, so they are independent of the size r of the dipoles fluctuations, as regulated at low Q 2 by the quark masses; see also the discussion of Fig. 3 below.) Also, we have obtained an equally good fit with the slightly larger value m c = 1.4 GeV for the mass of the charm quark, although the quality started deteriorating for significantly larger values m c ≥ 1.6 GeV. Very similar findings have been reported for the saturation fits in [26]. (vii) Trying to extend the fit to larger Q 2 shows an interesting behaviour as seen from Table 2. While the χ 2 obtained using the GBW initial condition increase when including higher-Q 2 data, the fits using the rcMV initial condition remain stable. We suspect that this is due to the fact that this choice of initial condition stays closer to the expected physics at high Q 2 .
In Fig. 2 one can see the quality of our fit and the extracted values of the evolution parameter λ s = d ln Q 2 s /dY . In Fig. 3 we show the value of the saturation momentum in the (x, Q 2 )-plane on top of the data points as well as a few selected initial conditions for the fit. Note that amplitudes which a priori have different functional forms, cf. Eqs. (14) and (15), look nevertheless quite similar in shape (at least in double-logarithmic scale) when plotted for the particular values of the parameters that are selected by the fits.
To conclude, this work can be seen as the first description of small-x DIS data which includes mandatory perturbative QCD ingredients in that region: leading-order small-x evolution, the resummation of large transverse logarithms, and saturation corrections 4 . The dipole amplitude obtained from our fits to inclusive DIS can in principle be used to compute several other observables, like particle multiplicity in hadronic collisions, the diffractive structure functions, the elastic production of vector mesons, or the forward particle production in heavy-ion collisions. This is certainly not the end of the story: beyond what we have included here, it would be interesting to add the pureᾱ 2 s NLO corrections to the BK evolution kernel, thus obtaining a genuine resummed NLO-BK fit, and to perform a proper matching between this small-x evolution and a DGLAP-like evolution at large Q 2 and large x. These steps go beyond the scope of the present paper and are left for future studies.