QCD traveling waves phenomenology revisited

In this paper we review and update the Amaral-Gay Ducati-Betemps-Soyez saturation model, by testing it against the recent H1-ZEUS combined data on deep inelastic scattering, including heavy quarks in the dipole amplitude. This model, based on traveling wave solutions of the Balitsky-Kovchegov equation and built in the momentum space framework, yields very accurate descriptions of the reduced cross section, $\sigma_{r}(x,y,Q^{2})$, measured at HERA as well as DIS structure functions such as $F_{2}(x,Q^{2})$ and $F_{L}(x,Q^{2})$. Additionally, it provides reasonable descriptions of heavy quark structure functions, $F_{2}^{cc}$ and $F_{2}^{bb}$ at small-$x$ and $Q^{2}\lesssim 60$ GeV$^{2}$. Having obtained excellent agreement with HERA data we use the improved model to make predictions for structure functions to be measured in the near future at LHeC.


I. INTRODUCTION
Most of the QCD nonlinear evolution equations at very large rapidities, Y = ln(1/x), possess asymptotic solutions which fall in universality classes. Since a long time [1][2][3], it is known that the nonlinear Balitsky-Kovchegov (BK) equation [4][5][6] can be mapped onto the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation [7,8]. It has universal and analytical traveling wave solutions which do not depend either on initial conditions or on the definite form of the nonlinear correction terms. Specifically, the BK equation derived for the unintegrated gluon distribution (UGD) in momentum space presents traveling-wave solutions in the transition region near the saturation domain despite the precise form of the nonlinear terms. Namely, the corresponding solution is controlled but the linear (dilute) region. Interestingly enough, the geometric scaling property observed in inclusive and exclusive processes at DESY-HERA data at small-x is connected to a traveling-wave structure of the scattering amplitude for a QCD color dipole off nucleons. The underlying quantity is the momentum saturation scale, Q s (Y ), which has its rapidity dependence driven by the velocity of the wave front, v c . The evolution time is t =ᾱ s Y and the position coordinate is ρ ∼ ln(k 2 /k 2 0 ) (where k 0 ∼ Λ QCD is a fixed infrared scale and we use ρ to denote the coordinate to avoid confusion with the Bjorken x), and the function obeying the universal class of equations is u(ρ, t). The position of the wave front is measured by the quantity ρ s = ln[Q 2 s (Y )/k 2 0 ] = v c Y and in the mean field approximation the solution to u presents the form u(ρ, t) = u(ρ − v c t) [1][2][3].
In the large-N c limit and in the mean-field approximation, the small-x behaviour of the forward QCD dipole scattering amplitude, N (r, Y ), follows the BK equation in coordinate space. This equation can be obtained also in momentum space, where it evolves the amplitude N (k, Y ), which is directly related to the UGD, F(k, Y ), through where R p is the proton radius. It can be easily show that the celebrated Golec-Biernat-Wusthoff (GBW) [9] form for the UGD, i.e. F GBW (k, Y ) = F 0 (k 2 /Q 2 s ) exp[−(k 2 /Q 2 s )] (with Q 2 s = κ 2 0 exp(λY ) and F 0 = N c R 2 p /2πα s ) gives N GBW (k, Y ) = 1 2 Γ(0, k 2 /Q 2 s ). Here, Γ(0, x) is the incomplete Gamma function and the amplitude presents clear scaling on τ s = k 2 /Q 2 s (Y ). On the other hand, the complete behavior of the amplitude at fixed QCD coupling (leading logarithmic order, LL) has been extensively investigated and presents a k-diffusion term typical of the Balitsky-Fadin-Kuraev-Lipatov (BFKL) solution [10][11][12] in the region k Q s . However, for the BK equation the saturation scale plays the role of a natural infrared cut-off and the fast broadening of the UGD at small-x is properly controlled. Geometric scaling behavior on the scaling variable τ s is restored in the region where the diffusive factor is negligible and solution is closer to the GBW form. The BK solution at LL accuracy will be revisited in next section, where it will be used in order to describe the recent results for the proton structure function at small-x.
Going beyond the LL approximation, the solutions of the BK equation at next-to-leading logarithmic (NLL) order have been also investigated [13][14][15][16][17]. In Ref. [15] three versions of the NLL BK equation were considered, namely the 1-dimensional BK equation with running coupling and two versions using quark-loop contributions. Moreover, modified BK equations including the renormalization-group corrections to the NLL BFKL kernels were studied. It was shown that there is a unified asymptotic prediction to observables and predictions for the behavior of exact solutions fall upon a large universality class of solutions [15]. That theoretical analysis led to phenomenological models presenting geometric scaling in √ Y rahter than in Y as in the fixed coupling case [18,19]. The role played by the fluctuations effects (Pomeron loops) in the NLL BK solution was analyzed in Ref. [16]. The starting point is a Langevin equation for the forward dipole-target scattering amplitude, N (r, Y ), with a Gaussian white noise. It was verified that a diffusive scaling for large rapidities, Y > Y form , takes place, where Y form is the rapidity interval needed for the solution to form a wave front down to the low density domain, where the noise term is relevant. The semi-analytical solution is somewhat consistent with numerical solutions of the (1+1)-dimensional reaction-diffusion toy model for high energy QCD presented in Ref. [20]. Afterwards, this numerical solution was used to describe inclusive and diffractive deep inelastic scattering (DDIS) in [21]. There, it was found that in DDIS the diffusive scaling is present for fixed coupling and, on the other hand, in the running coupling case geometric scaling takes place and it is reached at smaller values of rapidity than in the case without fluctuations [21]. Furthermore, in Ref. [17] the connection between the BK equation (with non-running and running couplings) in the diffusive approximation with noise and the extension of the stochastic FKKP (sFKPP) to the radial wave propagation in an absorptive medium is done. An important result is that a new geometric scaling domain forward to usual traveling wave front is found. The corresponding extended scaling presents a new scaling variable, with the wave front at position In this work, we revisit the phenomenological model proposed by Amaral, Gay Ducati, Betemps and Soyez (AGBS) [22], based on the analytical solutions of BK equation at leading logarithmic accuracy in the momentum space. An updated AGBS model is provided through fits to the recently extracted combined HERA DIS data on the reduced cross section. Both charm and bottom quark contributions to the proton structure function, F 2 (x, Q 2 ), are included. As a byproduct, charm, bottom and longitudinal structure functions are computed to be compared with the data. The plan of the paper is as follows. In Sec. II we describe the DIS cross section in terms of the AGBS model [22] for the dipole scattering amplitude in momentum space. In Sec. III fitting methods to HERA data on the reduced cross section, σ r (x, y, Q 2 ), are presented along with the fit-tuned parameters to F QQ (Q = c, b) and F L structure functions. In the last section, we discuss the main results of this study and give prospects of possible future studies.

II. DIS CROSS SECTION IN THE MOMENTUM SPACE FRAMEWORK
A. Asymptotic behaviors of N (k, Y ) and the AGBS model The BK equation at LL order in momentum space, whose whose solution is the dipole scattering amplitude N (k, Y ), is written as follows: whereᾱ s = α s N c /π and L = log(k 2 /k 2 0 ) with k 0 being an infrared cut-off scale. The quantity χ(γ) = 2ψ(1) − ψ(γ) − ψ(1 − γ) is the characteristic function of the BFKL kernel [10][11][12]. After an appropriate change of variables, the BK equation reduces to the FKPP equation [7,8] for u(ρ, t) ∝ N (k, Y ) when its kernel is approximated in the saddle point approximation, i.e., to second order in the derivative ∂ L , the so-called diffusive approximation. In this case the equation takes the form, ∂ t u(ρ, t) = ∂ 2 ρ u(ρ, t) + u(ρ, t) − u 2 (ρ, t), with t ∼ Y and ρ ∼ L corresponding to the time and space variables, respectively.
The FKKP equation presents asymptotic solutions described by travelling waves, meaning that at large rapidities (very small Bjorken longitudinal momentum fraction, x) they take the form u(ρ, t) = u(ρ − v c t) of a front travelling to large values of ρ at a speed v c without deformation. This is translated into the geometric scaling property, where the amplitude depends only on the quantity k 2 /Q 2 s , i.e. N (k, Y ) = N (τ s = k 2 /Q 2 s ). At non-asymptotic rapidities geometric scaling is violated, and the forward amplitude takes the following form for k Q s [1][2][3], where χ denotes the second derivative of the BFKL kernel with respect to the anomalous dimension γ. The parameters γ c and v c are obtained uniquely from the BFKL kernel and correspond to the selection of the slowest possible wave, v c =ᾱ s χ (γ c ). For the leading-order (LO) BFKL kernel, one obtains γ c = 0.6275 . . ., v c = 4.88ᾱ s , χ (γ c ) = 4.883 . . . and χ (γ c ) = 48.518 . . .. The rapidity dependence of the saturation scale is explicitly obtained and the leading contribution on the rapidity Y reads Looking at the amplitude in Eq. (3) it can be verified that the geometric scaling is obtained for a kinematic range where k 2 Q 2 s (Y )e β √ Y (the so-called geometric scaling window), with β = 2χ (γ c )ᾱ s . Equation (3) provides the behavior of the amplitude in the dilute region, where k Q s . The other region of interest is the vicinity of the saturation bound, k ≈ Q s , and deep inside the saturation domain, where k Q s . In a rough approximation, by using a Heaviside function for the amplitude in coordinate space, N (r, Y ) = Θ(rQ s − 1), in Ref. [22] it has been shown that the amplitude in momentum space, in the region k Q s , behaves like with a being a constant to be determined by the boundary conditions. The same behavior can be also obtained from the explicit solution of BK evolution equation inside the saturation region. For instance, for small momentum Q s k Λ QCD a similar expression for N (k, Y )can be derived from the Levin-Tuchin (LT) formula [23,24] for the S-matrix valid for larger dipoles, r 1/Q s . Starting from the LT solution, the corresponding unintegrated gluon distribution has been determined in Ref. [25]. Here, one has τ = (1 + 2iν 0 )/4χ(0, ν 0 ) ≈ 0.2. The function χ is an eigenvalue of the BFKL equation and is defined as 1275i. The determination of the UGD is non-trivial, however it can be approximated at leading logarithmic approximation for k 2 Q 2 s as [25] F(k, Y ) where γ E is the Euler-Mascheroni constant. It should be noticed that in the region k Q s , the UGD above goes to zero as k 2 → 0. By making use of the relation in Eq. (1), the dipole amplitude in momentum space is given by, k Qs which has exactly the same parametric form as the simple asymptotic expression in Eq. (5). The AGBS model [22] is an approach which explores the implications of the traveling wave solutions of BK evolution equation to the γ * p scattering. It provides a phenomenological expression for N (k, Y ) which analytically interpolates between the behaviors of the amplitude at the saturated and dilute domains, Eqs. (3) and (5). The final expression for N (k, Y ) [22] reads where and the saturation scale is given by Eq.(4).

B. DIS cross section with dipoles in momentum space
In electron-proton DIS the ep interaction is dominated by the exchange of a virtual photon γ * with virtuality Q 2 . In the dipole model this interaction can be seen in the following way: the virtual photon has enough energy to split into a quark-antiquark pair, a dipole, which then interacts with the target proton via gluon exchanges. This dipole has fixed transverse size r, the quark carrying a fraction z, and the antiquark carrying a fraction 1 − z, of the photon longitudinal momentum. The total γ * p cross section is then given by where Y is the total rapidity interval of the γ * p system and Ψ T,L (r, z; Q 2 ) 2 (well known from QED [26]) give the probabilities for the photon, with transverse (T ) and longitudinal (L) polarization, to split into the dipole. The quantity σ dip (r, Y ) is the total dipole-proton cross section which, according to the optical theorem, is given by where N (r, b, Y ) is the imaginary part of the dipole-proton scattering amplitude in coordinate space. In the general case, the amplitude depends not only on the dipole transverse size, but also on the impact parameter vector, b, of the dipole-proton interaction. If one neglects the b-dependence (which means considering the proton an homogeneous disk), the integration over the impact parameter is simplified and the dipole-proton cross section reads: where R p is the proton radius and N (r, Y ) is the solution of the simplified (b-independent) BK equation in coordinate space. The corresponding equation in momentum space, Eq.(2) (for N (k, Y )), is obtained through the modified Fourier transform [5,6]: This relation allows us to rewrite the total inclusive DIS cross section, Eq.(14), in terms of N (k, Y ) [22] and calculate, for example, the F 2 proton structure function as follows: where α em is the electromagnetic coupling constant. |Ψ T,L (k, z; Q 2 )| 2 now refer to the photon-dipole splitting probabilities expressed in momentum space, whose explicit forms can be straightforwardly obtained through and are given by [22] |Ψ T (k, z; and where 2 q = z(1 − z)Q 2 + m 2 f and m f denotes the mass of the quark with flavor f . Thus, with a model for N (k, Y ) at hand -and the AGBS parametrization is such a model -it is possible to calculate, in a momentum space framework, not only the F 2 structure function, but other physical quantities related to inclusive DIS, for example, the contributions of different flavors (masses) of quarks to the F 2 , as well as the longitudinal structure function, which can be evaluated in this momentum space approach by Its is important to point out that, besides being useful in the description of DIS data, the AGBS model also provides the fundamental tools to study inclusive observables at RHIC and LHC energies. In the paper where the model was proposed [22], it was used to fit measurements of the F 2 proton structure function from H1 [27] and ZEUS Collaborations [28,29] taking heavy-quark (charm) effects into account. Afterwards, in [30] another fit to F 2 has been performed, considering only the contribution of light quarks, but using (more recent) H1 and ZEUS combined HERA data [31]. The model has also been also used to investigate possible pomeron loop effects at HERA [32] and to describe inclusive hadron and photon production at the LHC [33]. In all these phenomenological applications the AGBS parametrization has been shown to be successful in data description. This, together with the fundamental properties underlying the construction of the model, makes its improvement an interesting issue.

III. DIS DATA AND FITTING PROCEDURE
In this paper we make an improvement of AGBS model by updating its parameters with a fitting procedure to recent HERA data including heavy -charm and bottom -quarks. In particular, we include in our fit the recently released data for the reduced cross section [31], which reads: where y = Q 2 /(sx) is the inelasticity variable, √ s denotes the center of mass energy of the ep collision and F L (x, Q 2 ) is the longitudinal structure function.
In fitting σ r a kinematic cut to HERA data is applied to the Bjorken-x variable, namely x ≤ 0.01, since this approach is concieved to describe high-energy amplitudes (the small-x behavior). Two bins of the photon virtuality are considered: Both bins prevent us from the need to include Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) corrections, which must be properly accounted for at too high values of Q 2 . The choice of fitting data in the Bin 1 range can be regarded as a conservative one, with respect to traditional approaches such as, e.g. GBW [34], for which an even lower Q 2 max (=10 GeV 2 , as long as DGLAP corrections are not included) is probed. Moreover, as we take into account heavy quark contributions and since the experimental range considered includes very small values of Q 2 , we perform the usual kinematic shift in the definition of Bjorken-x [9] x for charm and bottom, when the cutx f 0.1 is satisfied. Otherwise the contribution of heavy quarks is switched-off. Fits have been performed using the ROOT framework [35,36], through the members of the TMINUIT class 1 . In specific, we use the MIGRAD algorithm throughout, setting the Confidence Level (CL) to 95% 2 . Goodness-of-fit is evaluated using the standard Chi-squared (χ 2 ) per degrees of freedom (dof) criterion, with s i representing the reduced cross section data (N p = 271 for Bin 1 and N p = 309 for Bin 2), σ i the total uncertainty with respect to central values, s i , and σ r (x, y, Q 2 ) our model, according to Eqs. (18), (22) and (23). We also provide the integrated probability, P (χ 2 ; dof), the well-known p−value, also as goodness-of -fit estimator, with due care, namely limiting to interpret its results in the light of an overall agreement with data sets for the various models tested, specially when comparing fits to Bin 1 and Bin 2, and not in the traditional sense, that is, as a test of hypothesis used to discriminate good from bad models. Concerning the model parameters, the one kept fixed in this analysis isᾱ s = 0.2. For the value of the critical slope γ c , two different values were tested: γ c = 0.6275, which as mentioned before comes from the LO BFKL kernel, and γ c being free parameter. These values can be further compared to the one obtained from the fit performed in [37], γ c = 0.7376, using the Iancu-Itakura Munier (IIM) saturation model for N (r, Y ) including the heavy quarks. This value is also in agreement with what is expected from NLO BFKL (γ c 0.7). Thus, as in the previous studies using AGBS model, we are left with at least four free parameters, v c , k 2 0 , R p and χ (γ c ).  [22,30]. Only the results which provided the best fit quality are presented.
fit quality a Ref. [22] 0 For the quark masses we consider two different situations (i) with only light quarks and (ii) with light and heavy (charm and bottom) quarks. In both situations we use two different values for the light quark masses: m q = m u,d,s = 0.14 and 0.05 GeV. The first value is the most used in DIS phenomenology in the dipole framework, while the second is the one which provided the best fit to previous (not combined) HERA data in the original AGBS model [22]. In the case where heavy quarks are taken into account, charm and bottom quark masses are assumed to be m c = 1.3 GeV and m b = 4.6 GeV, respectively. For the sake of comparison, we show the values obtained in Refs. [22] (light quarks and charm quarks) and [30] (only light quarks) in Table I (we present only the main results).
The main results described above are presented in Tables II and III, where we introduce labels for different fit variants: V i B j , with i = 1, 2, 3, 4 standing for the different values of light quark masses in both situations described above, while j = 1, 2 indicates which bin has been used to tune the model parameters. In Table II we show the best results (Bin 2) of our fits to the DESY-HERA data for the reduced cross section when only light quarks considered. Although the two different choices of quark masses lead to significant differences with respect to the resulting values of the parameters of the AGBS model, they provide fits with similar qualities. In order to perform a cross check, we can compare our results with those obtained from the fits of [30] (see Table I). We see that the results are quite similar (concerning both the parameter values and the quality of the fit) to those obtained in the present work, with the same value, m q = 0.14 GeV, for the light quark masses.
In Table III we summarize the results of the fits to the data on the reduced cross section data with light and heavy quarks, for the two Q 2 -bins. A suited fit quality is found (see variants V 3 B 2 and V 4 B 2 ), given the data precision and a minimal number of fitted parameters. In the Table we present only the results with the parameter γ c fixed at the value γ c = 0.6275, since it provided the best fits to the data. Besides, we have tested the case where γ c is left free and verified a good stability for this parameter, with the fitted one being very close to the that coming from BFKL dynamics. The updated parameters are close to the original ones (Ref. [22], with only charm effects taken into account, see Table I) with v c having lower values by around 15%. We clearly see that the inclusion of heavy quarks still provide good fits to HERA data (see Tables II and III). The bottom quark contribution plays a small role in the Bin 1, whereas in the Bin 2 it is significant, although most of the results present p-values larger than the confidence level considered, α = 0.05, which demonstrate good statistical significance of the analysis. The variability in the fit quality estimators, χ 2 /dof and p− values, between the fits performed in Bins 1 and 2 can be noticeable, even though that does not compromise the goodness of fits by all means. In fact, as we shall see in the following, fits and predictions of models obtained by tuning our model parameters with Bin 1 essentially overlap with the ones from Bin 2. Such behaviour, seems to evidence not only an important effect of high−Q 2 and high−x in our dipole amplitude, but also that fitting a larger Q 2 bin may not be required in order to obtain reasonable predictions for heavy quarks structure functions. In essence, these results demonstrate that the AGBS model remains doing a good job even at large virtualities and small-x, mimicking part of the typical DGLAP evolution (driven by the extended geometric scaling behavior present in large k tale of the dipole amplitude). Parameter v c = λ 0.15 − 0.17 is compatible with λ values found in recent analyses using dipole models with extended geometric scaling in position space. For instance, IIM/CGC model [37,38] gives λ 0.23 whereas b-CGC model [38] found λ = 0.2063. The value of the parameter R p 4.62−5.3 GeV −1 , which is related to the black disc limit of γ * p cross section, σ 0 = 2πR 2 p 52−67 mb, produces larger values compared to corresponding models in configuration space where σ 0 ∼ 30 mb [34,38,39].
In Fig. 1, a comparison between the variants V 1 B 2 (solid lines, only light quarks) and V 3 B 2 (dashed lines, including charm and bottom) is shown against the H1-ZEUS combined F 2 data at Q 2 ∈ [0.1, 150] GeV 2 and x ≤ 10 −2 . A very good agreement with data can be observed and the curves are practically the same at very low-Q 2 . Small deviations appear only at large Q 2 and very small x. The results for light quarks are steeper (v c 0.17) than for those including heavy quarks (v c 0.15). The resulting dipole amplitude in momentum space obtained from present fits can be used for the prediction of LHC cross sections along the lines presented in Ref. [30]. In addition, by using of Eq. 1, the unintegrated gluon distribution in proton can be easily obtained. This is important for physics based calculations in the scope of TMD/k ⊥ -factorization formalisms.
As previously stated, the ABGS model nicely describes HERA data for small and moderate photon vitualities including the transition of the DIS structure functions to small values of Q 2 . It is known that this is achieved by the parton saturation corrections to the BFKL formalism embedded in the approach. This should be more evident in the longitudinal structure function, which is strongly affected by the screening corrections.
With the parameters given in Table III we are able to compute and make predictions for the charm and bottom structure functions. The results are presented in Figs. 2 and 3, respectively, where both contributions for the F 2 structure function are considered in the range 2.5 GeV 2 Q 2 120 GeV 2 [40], and we use the variants V 3 B 1 and V 3 B 2 . In the case of charm structure function, F cc 2 , one sees that AGBS model provides a good description of the data within a wide range of the photon virtuality, up to 60 GeV 2 , and a reasonable description at Q 2 = 120 GeV 2 . We have also made predictions for the bottom structure function, F bb 2 , finding a good agreement with the data, in particular for Q 2 7 GeV 2 , where there is a larger number of experimental points. Finally, we present our predictions for the longitudinal structure function, F L (x, Q 2 ), which in the present analysis can be evaluated using Eqs. (21) and (22). The results are presented in Fig. 4, where we show the behavior of F L as a function of x in the range 1.5 GeV 2 Q 2 120 GeV 2 of the photon virtuality, considering, as before, variants V 3 B 1 and V 3 B 2 . In all the range considered, we see that AGBS model provides a good agreement with the data. Besides, as already mentioned, Figs. 2 and 3 reveal that Bin1 data is sufficient to furnish accurate descriptions of both, F cc 2 and F bb 2 , even at virtualities as large as Q 2 ∼ 100 GeV 2 . Notwithstanding, this effect is even more drastic for F L (x, Q 2 ), as one can see from   [31]. F2 uncertainties are estimated, considering δF2 ≈ δσr, for y < 0.6. Curves: black solid and blue dashed curves are the predictions of variants V1B2 and V3B2, following from fits to σr(x, y, Q 2 ) including heavy quarks and only with light ones. Fit parameters of these curves are given in Table III. As a last analysis, based on V 3 B 2 variant, we give predictions for the Large Hadron Electron Collider (LHeC) [42], which extends the kinematical range of ep DIS to very low-x. It is proposed as a configuration with electrons of 50 − 100 GeV colliding with 7 TeV protons in the LHC accelerator. It is also planned high energy/luminosities configuration in a long term period [43,44] (HE-LHeC, √ s ep 1.7 TeV, and FCC-ep with √ s ep 3.5 TeV). This allows to explore Bjorken-x in DIS down to ∼ 10 −6 with high luminosity. Specifically, here we consider the LHeC scenario with E e = 50 GeV on E p = 7 TeV, √ s ep 1.3 TeV, with a luminosity of 50 fb −1 . This provides access to a kinematic region of 2 × 10 −6 < x < 0.8 and 2 < Q 2 < 10 5 GeV 2 . Our predictions to F 2 and F L are shown in 5 (fit including c and b quarks) for F 2 and F L compared to the simulated LHeC pseudodata electron-proton collisions at Q 2 = 10 GeV 2 and for 10 −6 ≤ x ≤ 10 −2 [42]. The extension of present model to very low-x is reasonably consistent with simulated LHeC data and it is expected that the real measurements can be able to discriminate between the models including saturation physics and constraints on the small x QCD dynamics. Predictions are also shown for the charm and bottom structure functions in Figs. 6. They are compared to the pseudodata generated by RAPGAP Monte Carlo for an LHeC scenario with electrons with E e = 100 GeV and protons with E p = 7 TeV for an integrated luminosity of L int = 10 fb −1 . We present the pseudodata for the configuration where the detector acceptance covers the whole polar angle range as well as events where at least one heavy quark (Q = c, b) is found with polar angles θ Q > 2 (10) degrees. The overall trend of simulated data indicates a possible enhancement of the charm and bottom within the proton at very low-x. For the time being, bearing in mind the recent HERA results on F cc 2 and F bb 2 , it seems premature (if not speculative) to take for granted such a behaviour at LHeC, reason for what our predictions are shown. Future, actual data, shall shed light on this matter.

IV. CONCLUSIONS
In this work we revisited and updated the AGBS color dipole model in the momentum space framework. The amplitude contains the BFKL dynamics at large k (diffusion) and transition to saturation regime using the traveling wave solutions of BK equation at leading order. The parameters have been fitted to the reduced cross section σ r [31] measured at DESY-HERA, taking into account the heavy quark contributions in the theoretical prediction for proton structure function F 2 . The investigation covered data in the region x ≤ 10 −2 and Q 2 ≤ 150 GeV 2 . An excellent quality of fit was found with χ 2 /dof ≈ 1 and good statistical significance with p-value either large. Using a confidence level of 95% (α = 0.05) most of analyzed cases obey p α. The fit quality of the original results for the AGBS model with light (+ charm) quarks remains preserved with heavy flavours included. The parameters have not changed significantly in comparion with previous versions of the model, with and without heavy quarks with exception to the χ (γ c ). Interestingly, the model considering only light quarks still describes the low-x/low-Q 2 data in a nice way. The saturation scale, Q 2 s (x) = k 2 0 x −vc , presents a weaker growth on x for heavy quarks than for only light ones. By using the parameters of the dipole amplitude in momentum space, N (Y, k), determined from the fit to the F 2 data, we predicted other inclusive structure functions. New predictions include the longitudinal, charm and beauty structure functions (F L , F cc 2 , F bb 2 ). It is found remarkable agreement with updated HERA data in all Q 2 bins. This means the model is able to emulate the DGLAP evolution at very large Q 2 and the correct parton saturation effects at low Q 2 . Predictions for the LHeC kinematic range were provided and compared to available pseudodata for that TeV scale ep machine.
Here, we have only considered the simplest scenario of LO expression for dipole amplitude and an extension addressing its NLO correction could certainly be done. Moreover, we envisage as future possibility to further explore the impact parameter dependence of the amplitude at LO and NLO. Finally, the present approach can be regarded as a starting point to study diffractive DIS (DDIS) and exclusive particle production as the Deeply Virtual Compton Scattering (DVCS) and exclusive vector mesons production, which we intend to investigate in future work.