HERA data and collinearly-improved BK dynamics

Within the framework of the dipole factorisation, we use a recent collinearly-improved version of the Balitsky-Kovchegov equation to fit the HERA data for inclusive deep inelastic scattering at small Bjorken $x$. The equation includes an all-order resummation of double and single transverse logarithms and running coupling corrections. Compared to similar equations previously proposed in the literature, this work makes a direct use of Bjorken $x$ as the rapidity scale for the evolution variable. We obtain excellent fits for reasonable values for the four fit parameters. We find that the fit quality improves when including resummation effects and a physically-motivated initial condition. In particular, the resummation of the DGLAP-like single transverse logarithms has a sizeable impact and allows one to extend the fit up to relatively large photon virtuality $Q^2$.


Introduction
Derived via systematic approximations within perturbative QCD, the Colour Glass Condensate (CGC) effective theory [1][2][3][4][5] is a powerful framework for computing high-energy processes in the presence of nonlinear effects associated with high parton densities. There are intense ongoing efforts towards extending this effective theory to next-to-leading order (NLO) accuracy, as required by realistic applications to phenomenology [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. These efforts refer both to the Balitsky-JIMWLK equations [20][21][22][23][24][25][26], which govern the high-energy evolution of the scattering amplitudes, and to the impact factors, which represent cross-sections at relatively low energy. The first NLO results, obtained more than a decade ago [6][7][8], refer to the Balitsky-Kovchegov (BK) equation [20,27]. The latter is a non-linear equation emerging from the B-JIMWLK hierarchy in the limit of a large number of colours (N c → ∞). It describes the evolution of the elastic scattering amplitude between a colour dipole and a dense hadronic target. Via appropriate factorisation schemes, like the "dipole factorisation" or the "hybrid factorisation" [28], the BK equation also governs the high-energy evolution of the cross-sections for processes of phenomenological interest, like deep inelastic scattering (DIS) at small Bjorken x, or forward particle production in proton-nucleus collisions.
A few years after the full NLO BK equation was first presented [8], its numerical study in [29] showed that it is unstable, as anticipated in [30,31]. Similar problems had been identified, and cured [32][33][34][35][36][37], for the NLO version of the BFKL equation [38][39][40] -the linearised version of the BK equation, valid for weak scattering. The origin of this difficulty has been clearly identified, both numerically [29] and conceptually [31,41]: it is associated with large and negative NLO corrections enhanced by a double "anti-collinear" logarithm, generated by integrating out gluon emissions with small transverse momenta (see Sect. 2). Such doublelogarithmic corrections spoil the convergence of the fixed-order perturbative expansion of the high-energy evolution, unless they are properly resummed to all orders. Refs. [31,41] proposed two different strategies for resumming this series of double-logarithmic corrections to all orders.
At a first sight, these strategies seemed to be successful, leading to stable evolution equations [41,42] and allowing for good fits to the small-x HERA data [43,44]. However, a recent study [45] revealed some inconsistencies in the original analyses in [31,41,43,44]. In particular, there was a confusion concerning the meaning of the rapidity variable which plays the role of the evolution time. The variable which a priori enters the perturbative calculations at NLO [8] and in the resummed equations proposed in [31,41] is the rapidity Y ≡ ln(s/Q 2 0 ) of the dipole projectile, with s the centre-of-mass energy squared and Q 0 a typical transverse momentum scale for the target. It is different from the rapidity η ≡ ln(1/x) of the hadronic target, which is the variable used in DIS (see Sect. 2 for details). When translated to this physical variable, the results of the original resummations in [31,41] show a strong scheme dependence, preventing any meaningful phenomenological applications [45]. The success of the corresponding fits to the HERA data [43,44] was rather fortuitous and can be attributed to several factors: (i) these fits blindly assumed the physical rapidity variable η = ln(1/x) (although this was inconsistent with the resummation scheme), (ii) the rapidity interval over which one can probe high-energy evolution is relatively limited making it delicate to critically probe the effects of resummation. In other words, even though the fits can hint at an evolution being better than another (e.g. via a better χ 2 , or more physical fit parameters), it remains difficult to draw a firm conclusion.
To overcome these difficulties, Ref. [45] proposed a reorganisation of perturbation theory in which the evolution time is the physical rapidity η. This program lead to a new version of "collinearly-improved BK equation", shown below in Eq. (6). This equation is non-local in rapidity. In that sense it looks formally similar to the equation proposed in [31] but these two equations are fundamentally different: (i) the evolution variable, η, occurring in Eq. (6) is the physical rapidity ln(1/x), and not the rapidity of the dipole projectile; (ii) the rapidity shift in Eq. (6) accounts for an all-order resummation of double collinear logarithms, and not anti-collinear (see Sect. 2). The collinear logarithms are generated by gluon emissions with relatively large momenta. Such emissions are atypical in the context of DIS. Moreover, they are partially suppressed by non-linear effects, so their resummation is somewhat less critical. As a consequence, a comparison between various resummation prescriptions shows only little scheme dependence [45], at the level of the expected perturbative accuracy of the resummed equation.
The evolution equation we use in this work, Eq. (6), actually goes beyond the original proposal from Ref. [45] by additionally including a class of DGLAP-like single transverse logarithms (both collinear and anti-collinear) which appears in the NLO BK equation. While in principle one can also include [45] the full set of NLO BK corrections in our evolution equation, this is numerically cumbersome (see however [42]) and it goes beyond the scope of this Letter. Instead, we directly confront the relatively simple equation (6) with the most recent HERA data [46] for inclusive DIS at small Bjorken x ≤ 0.01. This dataset includes more data points as compared to the one used in previous fits [43,44].
We use the fits to test the resummation of the aforementioned double logs and single transverse logarithms as well as various prescriptions for the running of the QCD coupling which enters the BK equation. We study two models for the initial condition to the BK equation: the simple Golec-Biernat and Wüsthoff (GBW) model [47,48] and a a running-coupling version of the McLerran-Venugopalan (MV) model [49]. Both models include two free parameters which are fitted to the experimental data. Two additional parameters (leaving aside the quark masses that we keep fixed) are associated with our ansatz for the running coupling and with the overall normalisation of the DIS cross-section.
With this setup, we find a good overall agreement with the HERA data at small x. Furthermore, it appears that the physically-motivated MV model with running coupling is preferred over the GBW model, with the former giving both better χ 2 and more reasonable parameter values than the latter. As with earlier studies, the fits appear to favour an initial condition for the dipole scattering amplitude with a fast and abrupt approach to the unitarity limit. This is the main reason why better fits are obtained with the running coupling version of the MV model as compared to its original fixed-coupling version [49]. This finding is also in agreement with the fact that previous fits [50,51] using the original MV model favoured a modified dependence of the initial amplitude on the dipole size, decreasing like r 2γ0 with γ 0 > 1. Another interesting observation of our fits is that the inclusion of the DGLAP-like single logarithms significantly improves the description of the data at large Q 2 , allowing for good descriptions up to maximum Q 2 of 400 GeV 2 .
This Letter is organised as follows. Sect. 2 provides a qualitative and (hopefully) pedagogical summary of the arguments justifying the replacement of the original collinear-improved BK equation formulated in terms of the rapidity of the dipole projectile [31,41], by a new version formulated in terms of the rapidity of the hadron target, or Bjorken x. (We refer to Ref. [45] for more details.) Sect. 3 discusses the main physical consequences of the various resummations on the solution to the collinearly-improved BK equation. Finally, Sect. 4 presents the main original results of this paper: the setup and the results for the fits together with their physical discussion.

Collinearly-improved BK evolution in the target rapidity
The main difference between our present approach and previous saturation fits to DIS refers to our specific choice of evolution equation used to describe the evolution of the dipole S-matrix with increasing energy. More precisely, we use a collinearly-improved version of the BK equation -recently proposed in [45] -in which the rapidity variable playing the role of the evolution time is the proton rapidity η = ln(1/x), with x the standard Bjorken variable. This contrasts previous studies (see e.g. [8,31,41,43]) where the evolution was formulated in terms of the rapidity Y of the dipole projectile. Beyond leading order [20,27], the choice of η over Y has important consequences. To make this clear, we first summarise the main findings of [45] which are relevant for the fit to DIS data described in the next section.
Basic kinematics, target and projectile rapidity. We use a frame in which the virtual photon, γ * , is an energetic right-mover with (light-cone) 4-momentum q µ ≡ (q + , q − , q ⊥ ) = (q + , − Q 2 2q + , 0 ⊥ ), whereas the proton target is a left-mover with P µ = δ µ− P − . 1 In the high-energy or small Bjorken x regime, the coherence time ∆x + 2q + /Q 2 of the virtual photon, i.e. the typical lifetime of its quark-antiquark (qq) fluctuation, is much larger than the longitudinal extent 1/P − of the target. This justifies the use of the dipole picture in which the γ * fluctuates into a qq colour dipole long before the collision, which then scatters inelastically off the proton. Via the optical theorem, the total dipole-hadron scattering cross-section is related to the S-matrix for the elastic scattering. At high energy, one can work in the eikonal approximation where the transverse coordinates x of the quark and y of the antiquark are not affected by the collision. The elastic S-matrix S xy also depends on the rapidity difference between the dipole and the proton through the high-energy evolution. The physical picture of this evolution and its analytical description depend on how the total energy is divided between the (dipole) projectile and the (proton) target i.e. upon the choice of the "dipole frame" in which one is working. It is useful to consider the two extreme situations: the "target frame", in which most of the total energy (and hence the whole high-energy evolution) is carried by the proton, and the "projectile frame", where the energy is mostly carried by the dipole (and the highenergy evolution is encoded in the dipole wavefunction). Importantly the rapidity variable which represents the "evolution time" for the high-energy evolution, is different in these two situations: For the target rapidity, the typical gluon from the proton which interacts with the dipole has a longitudinal momentum k − = Q 2 /2q + = |q − |, and hence a longitudinal extent ∼ 1/k − of the order of the lifetime ∆x + of the qq pair. For the projectile rapidity, the softest dipole to participate in the collision has a longitudinal momentum q + 0 = Q 2 0 /2P − -i.e. a lifetime 2q + 0 /Q 2 0 equal to the longitudinal extent 1/P − of the proton -, where the scale Q 0 Q is the transverse momentum scale for the onset of unitarity (multiple scattering) effects in the (unevolved) proton. The two rapidities differ by ρ ≡ ln(Q 2 /Q 2 0 ) which is large when Q Q 0 . The non-linear effects associated with the high gluon density are described differently in the two frames. In the target frame, soft gluon emissions occur in the proton wavefunction which is a dense environment. Accordingly, these emissions are modified by non-linear effects like gluon recombinations. The non-linear evolution of the dense hadron wavefunction has been computed only to leading order, yielding the (functional) JIMWLK equation [21][22][23][24][25][26]. Conversely, in the dipole frame one views the evolution as successive, soft, gluon emissions within the dipole wavefunction, a dilute hadronic system. Gluon emissions from the dipole occur like in the vacuum and non-linear effects exclusively refer to multiple scattering. This leads to the Balitsky hierarchy (and the BK equation), currently known to NLO accuracy [8][9][10][11][12]. Since our purpose in this work is to go beyond LO accuracy, we systematically use the dipole frame in what follows.
Time ordering and collinear improvements in Y (dipole frame). Besides being less suited for applications to DIS, the evolution with Y has a more severe conceptual drawback: the typical emissions contributing to this evolution at leading order can violate proper time ordering, that is, the condition that the formation time of a daughter gluon be smaller than the lifetime of its parent. 2 To understand this, we first recall that, when Q 2 Q 2 0 , the typical emissions associated with the high-energy evolution of the dipole wavefunction are strongly ordered both in longitudinal momenta (k + ) and in transverse momenta (k ⊥ ): This corresponds to soft and collinear emissions which yield the dominant, double-logarithmic, contributions proportional to powers ofᾱ s Y ρ. However, an explicit calculation of the relevant Feynman graphs shows [41] that this double-logarithmic enhancement only holds so long as the gluon lifetimes are ordered as well: This condition reduces the rapidity phase-space available for the evolution from Y to Y − ρ ≡ η. The condition Eq. (5) is already violated (due to the emission of daughter gluons with sufficiently soft k ⊥ ) in the LO BK evolution which resums an infinite series inᾱ s Y ρ, instead of the correct series in powers ofᾱ s (Y −ρ)ρ.
The difference between the two corresponds to an alternating series of double "anti-collinear" logarithms proportional toᾱ s ρ 2 which spoil the convergence of the perturbative expansion in Y . In particular, the NLO BK equation includes the first (negative) contribution proportional toᾱ s ρ 2 [8] which makes the evolution unstable [29] and hence unsuitable for physical studies. 3 To overcome this difficulty, it was originally proposed [31,41] to enforce time-ordering directly in the dipole frame evolution. Two "collinearly improved" BK equations have been proposed: in the first [31] the evolution is non-local in rapidity and has the same kernel as at LO, while in the second [41] the evolution is local in Y , but both the kernel and the initial condition receive corrections to all orders inᾱ s ρ 2 . Both methods allow for a faithful resummation of the dominant series in powers ofᾱ s ρ 2 , but the subleading terms (proportional to powers ofᾱ k s ρ 2 with k ≥ 2) are not under control. At a first sight, both strategies appear to be successful: the respective equations are stable [41,42], they can be extended to full NLO accuracy [42], and moreover they allow for good fits to the HERA data for DIS at small x [43,44].
Recasting dipole evolution in terms of η. A more recent study has revealed that these apparent successes were in fact deceiving [45]. The numerical studies in [42][43][44] have been presented in terms of Y instead of the physical rapidity η = Y − ρ = ln(1/x), and in the DIS fits in [43,44], the variable Y has been abusively interpreted as ln(1/x). The correct procedure would require to first transform the results from Y to η by a simple change of variables, before attempting a physical interpretation or a fit to the data. When following this correct procedure, one finds [45] an unacceptably large resummation-scheme dependence. For example when solved with the same initial condition the two "collinearly-improved" BK equations introduced in [31] and [41] yield very different predictions for the evolution in η. 4 More generally, different choices for the "rapidity shift" in the non-local equation in Y , albeit equivalent to DLA, lead to very different predictions for the saturation exponent λ s [45]. This strong scheme dependence forbids any physical interpretation of the results. It demonstrates that the uncontrolled, subleading, double-logarithmic corrections -which generally differ from one resummation scheme to another -are numerically important.
The problem is further complicated by the fact that the resummed BK evolution in Y cannot be formulated as a genuine initial-value problem. The non-local equation introduced in [31] must be solved as a boundary-value problem (on a line of constant Y − ρ) which seriously complicates even numerical studies of the equation. Moreover, even if the local equation with a resummed kernel [41] does admit an initial-value formulation, the corresponding initial condition must itself be resummed to account for double-logarithmic corrections to all orders, a task which appears to be intractable beyond strict DLA.
These difficulties with the evolution in Y can be avoided altogether by using η as an "evolution time" [45]. This choice has some obvious virtues in practice: η = ln(1/x) is the right variable to be used in phenomenological studies of DIS and, clearly, a boundary-value condition formulated at constant Y − ρ becomes an initial condition for the evolution in η. Most importantly, one can show that the evolution in η naturally ensures the proper time ordering of the successive emissions. 5 Instead of computing directly the target evolution in η which would be delicate in the presence of saturation, Ref. [45] proposed to reformulate the dipole evolution in Y , computed in pQCD, as an evolution in η via the change of variables Y ≡ η + ρ. Such a change of variables is unambiguous in perturbation theory and has been used [45] to obtain the NLO BK equation in η from the respective equation in Y [8].
Resummation of atypical collinear double logarithms. When using η as an evolution variable, the BK equation is not affected by the large anti-collinear logarithms that were present in the evolution in Y . This is a consequence of time ordering, Eq. (5), being automatically satisfied for the evolution in η. Furthermore, (5) also guarantees that typical anti-collinear emissions -i.e. those strongly ordered in transverse momentum according to (4) -automatically satisfy the proper ordering in longitudinal momentum k + (cf. (4) again).
However, even with the proper time-ordering, the correct ordering in k + can still be violated by a series of collinear emissions with a strong transverse momentum ordering opposite to that of Eq. (4). These violations yield double-logarithmic corrections to the BK kernel, starting at NLO (cf. [45]). In principle, this problem is as severe as the one of time-ordering violations in the evolution in Y : these collinear logarithms have to be resummed to all orders in the evolution in η, as the anti-collinear were resummed in the evolution in Y . However, these collinear emissions, where the daughter gluon has a much larger transverse momentum than the parent one are atypical in the context of DIS. One can further show that they are also strongly suppressed by saturation which freezes the evolution for emissions with sufficiently soft transverse momentum.
The NLO evolution in η nonetheless has a contribution from collinear logarithms which, albeit suppressed, eventually translates into an instability at large-enough rapidity. In practice one would therefore resum it to all orders (see [45]) using a rapidity shift leading to a non-local evolution in η (see Eq. (6) below). As for the resummations in Y , this resummation scheme is not unique beyond DLA. But unlike what happens with the resummations in Y , the scheme dependence for the resummations in η is reasonably small, in agreement with the expected perturbative accuracy of the resummed equations. For example, choosing different resummation schemes (e.g. by varying the η shift in (6)), one finds a small impact on the saturation exponents, comparable with missing perturbative contributions of O(ᾱ 2 s ). 4 In particular, they predict widely different values for the saturation exponent λs at large η, see the right panel of Fig. 1 in [41], where even the sign of the deviation w.r.t. the LO value appears to be different in the two schemes. 5 In a nutshell, for a gluon of (projectile) rapidity Y k = ln(k + /q + 0 ) and transverse momentum k ⊥ , one has η k ≡ Y k −ρ k = ln τ k τ 0 with ρ k = ln(k 2 ⊥ /Q 2 0 ) and τ k = 2k + /k 2 ⊥ is the gluon lifetime. Ordering in η therefore coincides with ordering in proper time.
The collinearly-improved BK equation in η. We are finally in a position to present the evolution equation in the target rapidity, η which reads where z is the transverse position of the gluon emitted by either the quark or the antiquark. In the large-N c approximation, implicit in (6), this can be viewed as the splitting of the original dipole (x, y) into two daughter dipoles (x, z) and (z, y). The other notations are explained below. Compared to the LO BK equation (in η) a few key differences are worth noting: (i) the use of the one-loop running couplingᾱ s (r min ) with the running scale set by the size r min of the smallest dipole: r min ≡ min{|x − y|, |x − z|, |z − y|}. Alternative schemes are discussed in the next section.
(ii) the rapidity shifts in the arguments of the S-matrices for the daughter dipoles ensure the resummation of the leading double logarithms associated with the k + ordering. They are given by with r ≡ |x − y|, and similarly for δ zy;r . They are non-zero only for emissions in which one of the daughter dipoles is (much) smaller than the parent one, in agreement with our earlier discussion about collinear logarithms. Expanding (6) to first non-trivial order in δ would give the double-logarithmic contribution to the BK kernel at NLO which eventually yields an instability unless properly resummed as in (6).
(iii) Eq. (6) also includes the resummation of the first set of single DGLAP logarithms (either collinear, or anti-collinear), via the factor r 2 /z 2 ±A1 , wherez ≡ min{|x − z|, |z − y|}. The number A 1 = 11/12 is related to the DGLAP splitting function via the following Mellin transform: The singular piece 1/ω generates the η = ln(1/x) logarithm contributing to the LO BK evolution, while the non-singular piece (−A 1 ) contributes at NLO and is enhanced by a single transverse logarithm ln(r 2 /z 2 ). The sign in the exponent, ±A 1 , is taken to be plus for an anti-collinear emission (r 2 <z 2 ) and minus for a collinear one, so this factor is always suppressing the evolution. Eq. (6) has to be solved as an initial value problem, albeit a somehow unusual one due to its non-locality in η. Indeed, since the shift introduces a dependence to rapidities smaller than η, if we want to start the evolution at some rapidity η 0 we should specify the initial condition for η ≤ η 0 . Our prescription is to assume a constant behaviour in η (see Sect. 9 of Ref. [45]) i.e.
With this prescription, Eq. (6) can be solved and used for DIS fits. A final note concerns the perturbative accuracy of Eq. (6). We have argued that it includes all the NLO corrections enhanced by at least one transverse logarithm. Hence, from the viewpoint of a strict weak-coupling expansion, it is accurate up to pure NLO corrections, of O(ᾱ 2 s ) without any logarithmic enhancement. Furthermore, the resummation-scheme dependence of (6) is also coherent with missing O(ᾱ 2 s ) terms. It is possible to extend this equation to full NLO accuracy by adding the missing pureᾱ 2 s corrections. The resulting equation, which can be found in Ref. [45], is substantially more complex and we postpone its applications to DIS to future work.

Illustrating the impact of running coupling and resummation effects
Before turning to the description of inclusive DIS data, it is helpful to briefly illustrate how the various ingredients in the BK equation, namely running-coupling (RC) corrections and the resummation of transverse logarithms, affect the evolution in η. For this purpose, we choose a homogeneous target, i.e. take S xy (η) = S(η, r) with r = |x − y|, together with the simple Golec-Biernat-Wüsthoff (GBW) initial condition [47,48]: S 0 (r) = exp(−r 2 Q 2 0 ) with Q 2 0 = 1 GeV 2 . This Gaussian Ansatz does not capture the proper behaviour for r 2 Q 2 0 1 but is enough for illustrating our points in this section. Let us first discuss running-coupling effects. For the sake of the present discussion, we usē withb 0 = 0.75 (corresponding to n f = 3) and Λ = 0.2 GeV. The Landau pole is avoided by freezing the coupling at the valueᾱ sat = 0.67. There is some freedom in implementing RC corrections in the BK equation and we consider two different prescriptions: the minimal dipole prescription,ᾱ s (r min ), where Eq. (10) is evaluated at r = r min ≡ min{|x − y|, |x − z|, |z − y|} and the BLM prescription (also dubbed as "fast apparent convergence" [43]), defined as Other prescriptions, not studied here, are also possible (see e.g. [6][7][8]). They all have in common that they reduce toᾱ s (r min ) when one of the three dipoles is much smaller than the other two, as one can easily check forᾱ BLM . This minimises the NLO correction to the BK equation proportional to the one-loop β-function.
Besides varying the prescription for the running coupling, we also aim to probe the effect of the resummation of transverse logarithms. Single logarithms can be switched off by removing the factor r 2 /z 2 ±A1 in Eq. (6), while to remove double transverse logarithms we set the η shifts, δ xz;r and δ zy;r , to zero.
In practice, we solve Eq. (6) numerically up to η = 20 and study the saturation exponent defined as with the saturation momentum Q s (η) numerically obtained from the condition that S(η, r) = 1 2 when r = 2/Q s (η). The rapidity range under study is not sufficient to reach the universal asymptotic behaviour but is representative for the actual range covered by the HERA data, η 10.
Our results are presented in Fig. 1 for a fixed-coupling prescription as well as for our two RC prescriptions. In each case we show the saturation exponent for the LO BK equation without any transverse logarithm resummation ("LO", dotted red), when double transverse logarithms are resummed ("DL", dashed green) and when both double and single transverse logarithms are resummed ("DL+SL", solid blue). The first observation is that running coupling effects have a large impact, reducing the saturation exponent by ∼75% for both RC prescriptions at η = 10 and even more at larger rapidities. The effect of resumming large transverse logarithms is smaller but still clearly visible: with RC, we see an additional ∼ 10 − 25% reduction coming from the resummation of double logarithms and a ∼ 10 − 20% reduction coming from single logarithms. (The effect is considerably larger when using a fixed coupling.) The fact that single-log effects are almost as large as double-log effects is most likely due to the fact that while the latter are only relevant for large dipole sizes -where they are reduced by saturation effects -the former have an impact on both large (collinear) and small (anti-collinear) dipoles.
Furthermore, one sees that, without the resummation of the transverse logarithms, λ s is still quite large (e.g. λ s (η = 10) = 0.26 for rcBK with the BLM prescription). This seems to be still too large to optimally accommodate the small-x HERA data. It also likely explains why previous fits based on rcBK [50][51][52][53] appeared to prefer another RC prescription, due to Balitsky [6], which predicts smaller values forᾱ s .
Beside the saturation exponent which characterises the "speed" of the evolution, it is interesting to look at the form of the amplitude itself. Both the initial condition (dash-dotted black) and the amplitude at η = 20 for different degrees of resummation (cf. Fig. 1) are shown in Fig. 2, for the same three configurations as in Fig. 1. Results are plotted as a function of the dimensionless ratio rQ s /2 so as to better highlight the shape of the dipole amplitude. The bottom part of the plot shows the effective slope Two striking features can be observed from these plots. First, the resummations of double and single transverse logarithms have a very small impact on the form of the amplitude. Even though a small effect is visible for the fixed-coupling evolution, almost no effect is observed with either running-coupling prescriptions. Then, one sees that the evolved amplitude has a less sharp transition between the dilute (small-r) and saturation (large-r) regions. This is particularly visible for the effective slope γ s which grows very slowly towards 1 at small dipole sizes.

Fits to the HERA data for inclusive DIS
We now come to our main task which is to describe the HERA data [46] for the inclusive DIS cross-section at x ≤ 0.01.

The dipole factorisation and the fit set-up
The dipole factorisation for DIS at small x (and at leading order in pQCD) expresses the physical picture in which the virtual photon fluctuates into a qq pair with a lifetime much longer than the longitudinal extent of the proton target. The total γ * p cross-section therefore factorises as a wavefunction for the γ * → qq splitting and an interaction between the dipole and the proton: where the squared light-cone wavefunctions (below,Q 2 represent the probability that a (longitudinal or transverse) virtual photon splits into a qq colour dipole with transverse size r and with "plus" longitudinal momentum fractions z and 1 − z for the quark and antiquark respectively. The sum in Eq. (14) runs over the quark flavors and our fit includes contributions from the three light quarks with m u,d,s = 100 MeV and from the charm quark, with m c = 1.3 GeV. 6 Furthermore, σ dipole (η f , r) is the dipole-proton total cross-section, evaluated for a proton (target) rapidity In principle this cross-section should be computed by integrating the dipole scattering amplitude T = 1 − S over all impact parameters, but since the impact-parameter dependence is non-perturbative and not properly encoded in the BK equation, we simply assume, in the spirit of a mean field picture, that the proton is a uniform disk of radius R p (treated as a fit parameter): σ dipole (η f , r) = 2πR 2 p T (η f , r), with T (η, r) obtained from numerical solutions to the homogeneous version of Eq. (6).
One peculiar feature about Eq. (14) is the fact that the dipole cross-section is evaluated at the rapidity scale η f = ln(1/x f ), which refers to the virtual photon, and not to the dipole. This looks natural from the perspective of the target evolution in k − (recall the discussion following Eq. (1)), but it might look less obvious when thinking about the evolution of the dipole. Clearly, the dipole and the virtual photon have different rapidities, due to the splitting fraction z, which can be very asymmetric, i.e. z 1 or 1 − z 1, especially for a transverse photon . We show however in Appendix A that changes in the longitudinal and the transverse phase-spaces associated with the γ * → qq splitting compensate each other and that the use of the photon rapidity η f in (14) is valid.
The quantity we actually fit is the reduced cross-section, related to the γ * p cross-sections (14) by y is the inelasticity parameter defined through Q 2 = xys, with s the squared centre-of-mass energy of the ep collision. In order to solve the evolution equation (6), we still need to specify the running coupling prescription and the initial condition. For the running of the QCD coupling, we use the one-loop expression where the dependence on the number of active flavours N f is made explicit. The value of Λ 5 is determined by imposing α s (M 2 Z ) = 0.1181 [54], and Λ 3,4 are fixed by the continuity of α s at the flavour thresholds. We use m c = 1.3 GeV and m b = 4.5 GeV for the charm and bottom quark masses. In coordinate space we use where we included a fudge factor C α , which will be one of the parameters to be fitted to the data. Finally, we freeze α s at a value α sat = 1 to regularise its infrared behaviour.
For the initial condition of the BK evolution, we use two parametrisations, most conveniently written for the dipole scattering amplitude T (η, r) ≡ 1 − S(η, r). The first choice corresponds to a modified version of the Golec-Biernat and Wüsthoff (GBW) "saturation model" [47] (used already in Fig. 1): where Q 0 and p are parameters to be fitted. The parameter p, not present in the original GBW model, controls the shape of the amplitude close to saturation. Our second choice is a running-coupling version of the McLerran-Venugopalan (MV) model [49], called rcMV in the following, which reads: Once again, Q 0 and p are free parameters (and α sat = 1). The other two parameters of the fit are C α and R p , the proton radius. In practice, we restrict p and C α to values smaller than 4 and 10, respectively. Indeed, the natural value for p in the MV model is p = 1 and it would be reassuring that the fits prefer such a value. Similarly, a natural value for C α should be of order one.

Fit results
In Table 1 we quote the results of the fit to the combined HERA data for the reduced cross section [46]. We include in the fit data points with x < 0.01 and Q 2 < 50 GeV 2 . For both initial conditions (GBW and rcMV) we show the results obtained using the two running coupling prescriptions (smallest dipole or BLM) in three cases corresponding to different resummations of transverse logarithms: pure LO (rcBK) evolution, resumming only double logarithms, or resumming both double and single logarithms. One can see that the resummation of the double and of the single logarithms are both improving the agreement with the data. This is particularly the case for the resummation of single transverse logarithms which not only have a significant impact on the quality (χ 2 ) of the fit, but also lead to more physical values for the fit parameters (in particular C α , for which values much larger than 1 mean that the evolution needs to be artificially slowed down by an unnaturally small value of the coupling in order to be compatible with data). 7 These two resummations allow to reach values of χ 2 per data point of less than 1.2 for all the combinations of initial conditions and running coupling prescriptions considered here. We believe (see e.g. the discussion in section 3) that this good agreement with the data is largely due to the reduction of the saturation exponent, λ s , due to the resummation of transverse logs. A particular consequence is that, contrary to some previous rcBK fits, we do not need to use the peculiar Balitsky prescription for the running coupling [6] to obtain a good agreement with the data. One notices also that the fit shows a slight preference for the BLM running-coupling scheme which predicts slightly smaller saturation exponents compared to the smallest dipole prescription.
In Fig. 3 (left) we show the shape of the initial condition as a function of r for the four fits which take into account the resummation of both single and double logarithms. The results look quite similar despite the different functional forms, which is due to the rather strong constraints from the data. In the right panel of Fig. 3, we show the saturation scale as a function of x for each of these four fits, superimposed over the HERA data points that we use in the fit.
In Table 2 we show how the fit is affected by including data with Q 2 larger than 50 GeV 2 . A priori, one would not expect a very good agreement with data at such high transverse scales, where DGLAP effects are expected to be essential. Nevertheless one can see that, when including the single logarithms resummation, the fit quality remains about the same up to Q 2 max = 400 GeV 2 for a rcMV initial condition. This is likely due to the fact that, as explained in section 2, these logarithms represent an important subset of the DGLAP contributions (at least at small x) and therefore improve the description at large Q 2 . Even when resumming the single transverse logarithms, the fit quality gets worse when going to larger Q 2 with a GBW-type initial condition (21). This is probably related to the fact that the GBW model does not have the correct physical behaviour at large transverse momenta, or small dipole sizes.
Since we take into account the charm contribution to the inclusive reduced cross section σ red in Eq. (18), we could also in principle compute the charm production cross section σ cc red . In [43] it was found that,  after fitting the data for the inclusive cross-section σ red , a very good χ 2 /ndf (< 0.7) was also obtained -without further tuning the fit parameters -for the σ cc red data presented in [55]. We believe that this was mostly a coincidence since a comparison between the best fits in [43] and the more recent, combined, HERA data for charm production [56], with more points and smaller uncertainties, shows a much worse agreement (χ 2 /ndf > 4). The situation is similar with the present fit. This should not be a surprise as the inclusion of heavy flavour data in the saturation fits based on the BK equation is a longstanding issue. The resummations considered here are not expected to bring any concrete improvement on this issue, as they are not adapted to the inclusion of heavy quarks.

Positivity of the dipole amplitude's Fourier transform
In this work we used the initial condition parametrisations (21) and (22) proposed in [43]. As was later shown [57], a drawback of these expressions is that their Fourier transforms are not positive-definite, which can lead to unphysical results in momentum space such as a negative unintegrated gluon distribution. The original GBW [47] and MV [49] parametrisations are not affected by this issue, however we were only able to obtain rather poor fits (with χ 2 /ndf > 1.4) when using these expressions. A similar discussion applies to the more recent MV e form, introduced in [53], which involves an extra parameter e c and reads T (η 0 , r) = 1 − exp − r 2 Q 2 0 4 ln 1 rΛ + e c · e .
That said, it should be pointed out that the high-energy evolution tends to improve the situation. If we concentrate on our fits using the rcMV initial condition, which is physically most appealing, we numerically find that, even though the Fourier transform of the initial condition has negativity issues, the evolved amplitude becomes positive after a few units of rapidity, cf. Fig. 4 (left). (Some small oscillations remain in an intermediate rapidity range before disappearing at larger rapidities.) This happens because the solution develops an effective slope, defined in Eq. (13), which is smaller than 1 even for small r as can be easily seen in the lower row in Fig. 2. This should be contrasted to the small-r behaviour of the rcMV initial condition for T (r), which vanishes faster than r 2 when r → 0. In fact, after evolution, not only S(r) but also T (r)/r 2 (whose Fourier transform is proportional to the gluon occupation number in the proton) satisfy all the conditions given in [57] which are necessary for a function to have a positive Fourier transform. This is confirmed by the numerical results displayed in Fig. 4.
Finding a functional form which has a positive-definite Fourier transform while preserving the good agreement with the HERA data would be extremely interesting, as this could open the way towards a unified description (via the dipole factorisation) of inclusive DIS and of particle production in "dilute-dense" (ep,  Table 1.) The discontinuities in the curves corresponding to the initial condition (η = 0) and the early rapidities reflect the negativity problem mentioned in Sect. 4.3.
eA, pp, pA) collisions. However, this seems to be also very challenging. At this level, it is legitimate to ask whether this difficulty solely reflects our inability to imagine versatile enough parametrisations for the initial condition, or it rather points towards a deeper problem with the dipole picture in the presence of a running coupling (perhaps similar to the problem discussed in [16] in the context of pA collisions). However, addressing such deep issues goes well beyond the scope of the present work. Our main goal here was to show that, for a given and physically-motivated initial condition, the use of the collinearly-improved version of the BK evolution improves over the standard rcBK dynamics when it comes to a description of the HERA data for inclusive DIS at small x.