Saturation model of DIS : an update

We present the results of new fits to the recently extracted data on $F_2$ at low $x$ with the GBW saturation model and its modification to cover high values of $Q^2$. We find that the model stands the test of time and gives a good description of the data with slightly modified parameters. All the essential elements of the model, especially the saturation scale, are retained.


Introduction
The recently extracted data [1] on the proton structure function F 2 from the HERA collider summarizes experimental effort [2,3] done by the DESY collaborations H1 and ZEUS to maximally extend our knowledge about the quark-gluon structure of the proton. The new kinematic region of small values of the Bjorken variable x, revealed by HERA, corresponds to the high energy (or Regge) limit of QCD. In this limit, QCD enters the semihard perturbative domain where virtual photon "mass", Q, is much bigger than Λ QCD and much smaller than invariant energy of the virtual photon-proton system, W , i.e. when x ≈ Q 2 /W 2 ≪ 1. The first condition makes the running strong coupling, α s (Q), small enough for perturbative calculations to be applicable, leading to either the DGLAP [4][5][6] or the BFKL [7-10] evolution equations. Both equations predict that the small x structure of the proton is dominated by a strongly rising gluon density when x → 0, which drives a similar rise of the sea quark densities.
The problem of taming the growth of such densities, related to the question about unitarity, was pioneered in [11,12] and developed in [13] in terms of nonlinear modifications of the known evolution equations, resulting from gluon recombination. A further refinement of the problem was proposed in [14][15][16] as an effective field theory approach, which culminated in the formulation of the JIMWLK evolution equations [17][18][19][20][21] describing the over-populated gluonic state, called Color Glass Condensate (CGC). An equivalent formulation was proposed in [22] in terms of the hierarchy of equations for Wilson line operators. This hierarchy simplifies in the limit of large number of colors N c , which leads to the formulation of the Balitsky-Kovchegov (BK) equation for a dipole scattering amplitude in [23,24]. The dipole amplitude is a key element in the computation of the small-x DIS structure functions in the dipole picture in which the virtual photon splits into a qq dipole interacting with the proton through gluonic fields. The dipole scattering amplitude obtained from the BK equation unitarizes the DIS structure functions by taming their power-like growth with x → 0, which is a reflection of a similar growth of the gluon density, see [25] for numerical studies of this effect.
The solution to the BK equation provides theoretical justification of the qq dipole scattering amplitude (or a dipole cross section) in the GBW model proposed in [26] for phenomenological studies of the transition of F 2 to small values of Q 2 and large values of W 2 . With only three fitted parameters, a good description of the DIS data was found. The dipole scattering amplitude in the GBW model vanishes for transverse dipole sizes r → 0 and saturates to a constant value for large r. These features make a connection to QCD for small r (color transparency) and to the idea of gluon saturation for large r in which the number of gluons (or qq dipoles in large N c limit) is tamed due to their recombination.
The key element in the description is an x-dependent saturation scale, Q s (x), which sets the scale for dipole sizes. With the proposed saturation scale, the dipole cross section saturates for smaller and smaller dipole sizes with decreasing x. The GBW model was updated in [27,28] to improve the large Q 2 description of F 2 by a modification of the small r behaviour of the dipole cross section to include the DGLAP evolved gluon distribution. A similar in spirit parameterization of the dipole scattering amplitude, based on the BK equation solution, was proposed in [29].
Since the publication [26], numerous fits with the dipole picture of DIS were performed. Let us mention the fits with the impact parameter dependent cross section [30][31][32] and the fits based on the solution to the BK equation [33][34][35][36]. Such studies are particularly important in view of the analyses [37,38], which show limitations of the linear DGLAP evolution for small x and Q 2 values. However, we should also mention the analyses [39,40] on the limitations of the dipole models.
The main purpose of this presentation is to answer the question whether we still obtain a good description of the new HERA data [1] with the saturation model [26] and its modification [27,28]. To this end, we preform new fits which update the original values of the fit parameters. We include both charm and bottom quark contributions to F 2 (generated radiatively from gluons), which allows us to make predictions for the comparison with data. We also compute the longitudinal structure function F L to be compared with the data.
The paper is organized as follows. In Section 2 we present the results for the fits with the GBW and DGLAP improved models and make a detailed comparison of them. In Section 3 we presents the comparison of the model results with the data on the proton structure functions, F 2 , F c,b 2 and F L . Finally, we summarize our findings in Conclusions.

Fit results
We fit the proton structure function F 2 computed from the HERA data [1] on the γ * p cross section from the relation valid in the small x approximation. In all our fits we take the data with x ≤ 10 −2 and minimal value of Q 2 = 0.045 GeV 2 . We consider both the light and heavy quarks (charm and bottom) in our calculations of the cross sections. The basic theoretical formula for F 2 corresponds to the dipole picture of DIS at small x in which the virtual photon dissociates into a quark-antiquark pair (a qq dipole) and subsequently interacts with the proton. Thus where T, L refer to virtual photon polarization, transverse and longitudinal, and where the sum over quark flavours f is performed. The photon wave function, Ψ T,L , is known [41], but the dipole cross sections is modelled with a few parameters which are fitted to the data. Since the photon wave function depends on mass of the quarks in the qq dipole, we can consider contributions to the structure functions from the individual quark flavour pairs where F l T,L is the sum of the contributions from the light quark pairs, uū, dd and ss, while F c T,L and F b T,L are the contributions from the cc and bb pairs, respectively. In such a case, we also modify the Bjorken variable x in the dipole cross section where W 2 is an invariant energy squared of the γ * p system. The conditionx f ≤ 1 accounts for the quark-pair production threshold. However, to be in accord with the small x approximation, we set the upper limit for the allowed values ofx f substituted in the dipole cross section tox f ≤ 0.1. Above this value, the charm and bottom structure functions are equal to zero. Nevertheless, almost all experimental points from HERA, which we use in our analysis, fulfill the boundx f ≤ 0.1.

Fits with the GBW model
The dipole cross section of the GBW model is given by [26] where the saturation scale Q s is defined as  Table 1. Fit results to the HERA data for Q 2 ≤ 10 GeV 2 with N p = 222 points, using the dipole cross section (2.6). Quark masses are given in units of GeV. χ 2 /N dof computed with the parameters from [26] are given in the first two rows. The parameters of Fit 2 are used for further analysis.
with Q 2 0 = 1 GeV 2 . The above cross section has a property of geometric scaling [42], i.e. it becomes a function of a single variable, rQ s , for all values of r and x. The three parameters of the fits with the GBW model are: σ 0 , λ and x 0 . The GBW model does not incorporate QCD evolution in Q 2 in the small-r part of the dipole cross section. As a consequence, the power λ, which governs change of the saturation scale with x, is independent of Q 2 . Therefore, one does not expect the model to fare well for high photon virtualities. In view of the above, we restrict data points used in the fits to those with Q 2 ≤ 10 GeV 2 and x ≤ 10 −2 , which leads to N p = 222 points.  Table 1 as a function r (left plot) and rQ s (right plot) for x = 10 −6 , . . . , 10 −2 (curves from left to right, respectively). All the curves in the left plot merge into one line in the right plot due to geometric scaling (2.8).  Table 2. Sensitivity of the quality of Fit 2 from Table 1 to the choice of Q 2 max in the data. N p is the number of experimental points selected by Q 2 max .
The results of our fits with the new HERA data are summarized in Table 1. In the first two rows we show the value of , computed with the original parameters of the GBW model [26] (no fits were performed). Rather large values of χ 2 /N dof suggest necessity of new fits.
In the next three rows, Fits 0-2, we present the results of the fits which update the original parameters of the model. The parabolic errors of the fit parameters are given by MINOS from the MINUIT package [43,44], which we use in this paper to minimize χ 2 . We see a good fit quality, taking into account the data precision and a minimal number of fitted parameters. The new parameters are close to the original ones with smaller values of the parameter λ by around 10%. In addition, unlike the original fit results, adding charm into the analysis improves the data description. The bottom quark contribution is very small for Q 2 ≤ 10 GeV 2 , nevertheless, we take the parameters of Fit 2 for further analysis to have a full handle on heavy quarks.
In Figure 1, we show the dipole cross section (2.6) with the parameters from Fit 2 computed for x = 10 −2 , . . . , 10 −6 . These are the curves from right to left in the left plot where the dipole cross section is plotted as a function of the dipole size r. All the curves merge into one solid line in the right plot where the dipole cross section is plotted as a function of the scaling variable rQ s . This is a reflection of geometric scaling (2.8) in the GBW model.
Finally, we study the sensitivity of the fit quality to the choice of the maximal value of the photon virtuality, Q 2 max , in the data. The results of the fits with five flavours for the indicated values of Q 2 max are shown in Table 2. We refrain from quoting the fit parameter errors, whose numerical values are very close to those from Fit 2 in Table 1. We see that the choice Q 2 max = 10 GeV 2 is indeed optimal, although extending the applicability of the GBW model up to Q 2 max = 20 GeV 2 is still acceptable. Summarising, the GBW model gives a good description of the HERA data for small and moderate values of Q 2 , accounting for the transition of the DIS structure functions to small values of Q 2 ≤ 1 GeV 2 This is achieved with the idea of parton saturation and only three fitted parameters. A graphical illustration of the agreement of the model with the data will be presented in Section 3.  Table 3. The results of the fits to the HERA data for Q 2 ≤ 650 GeV 2 with N p = 387 points, using the dipole cross section (2.9) with the scale (2.17). The quark masses m l = 0 and m c = 1.3 GeV. χ 2 /N dof computed with the parameters from [28] and the scale (2.10) are given in the first two rows. The parameters of Fit 2 are used for further analysis.

Fits with the DGLAP improved saturation model
The DGLAP improved saturation model [27,28] implements the dipole cross section given by where g(x, µ 2 ) is the gluon distributions taken at the scale The gluon distribution is evolved with the DGLAP evolution equations truncated to the gluonic sector, from the initial condition taken at the scale Q 0 = 1 GeV. The choice of the power 5.6, which regulates the large-x behaviour, is motivated by global fits to DIS data with the LO DGLAP equations, see [27,28] for more details. The splitting function P gg contains real and virtual terms with the number of active quark flavours n f in the latter one with C A = N c = 3 and T R = 1/2. In the leading order strong coupling constant we set Λ QCD = 300 MeV. Thus, the model has five parameters to fit: σ 0 , A g , λ g , C and µ 2 0 . The dipole cross section (2.9) has the property of colour transparency and tends to the perturbative QCD result in the limit r → 0. Indeed, for small dipoles, the scale µ 2 ≈ C/r 2 and the dipole cross section is proportional to r 2 with the logarithmic modifications due to the scale dependence of the gluon distribution [45], σ dip ≈ π 2 3 r 2 α s (C/r 2 ) xg(x, C/r 2 ) .  Table 3 (dashed lines) as a function r (left plot) and rQ s (right plot) for x = 10 −6 , . . . , 10 −2 (curves from left to right, respectively). The solid lines, corresponding to Eq. (2.15), merge into one solid line in the right plot due to geometric scaling with the saturation scale (2.16). For rQ s ≥ 1, also the dashed curves merge due to geometric scaling in the dipole cross section (2.9) in this region.
These additional logarithms allow better fits to the data for large values of Q 2 . In the limit of large dipoles µ 2 ≈ µ 2 0 , which leads to Thus, at large r, we find the GBW form of the dipole cross section with the saturation scale with the x dependence given by the gluon distribution taken at the scale µ 2 0 . Geometric scaling is strictly valid for (2.15), which is not the case for small dipoles when an additional r dependence is introduced in the dipole cross section (2.9) through the scale (2.10).
The above features of the dipole cross section can also be obtained for a slightly different choice of the scale µ, which interpolates smoothly between the C/r 2 behaviour for small r and the constant behaviour, µ 2 = µ 2 0 , for r → ∞. The fit quality is better for such a choice, thus in the forthcoming, we present the results of the fits with the above scale.
The results of the fits are presented in Table 3. The parabolic errors of the fit parameters are given by MINOS from the MINUIT package. We no longer restrict the data to  Figure 3. The normalized-to-one dipole cross sections (2.6) (solid lines) and (2.9) (dashed lines) with the parameters from Fit 2 in Table 1 and 3, respectively, for x = 10 −6 , . . . , 10 −2 (curves from left to right, respectively) Geometric scaling for the variable rQ s is clearly visible in the right plot.
low Q 2 values since the DGLAP evolution is incorporated in this model. Now, the data set contains N p = 387 points with Q 2 ≤ 650 GeV 2 and x ≤ 10 −2 . We also no longer consider the case with light quarks only since the heavy quark contribution to F 2 (especially from charm) becomes important for large values of Q 2 . In addition, we use the PDG world average value of the charm quark mass, m c = 1.275 ± 0.025 GeV ≈ 1.3 GeV [46].
In the first two rows of Table 3, we show χ 2 /N dof computed with the original parameters of the DGLAP improved model [28] and the prescription (2.10) for the scale µ (no fits were performed). In the next two rows, Fits 1-2, we present the fit results to the new HERA data with the scale (2.17). We find much better description of the data after refitting the parameters. The nonzero light quark mass (equal to 140 MeV in the GBW model) deteriorates the fit quality, thus we set it to zero, m l = 0. For the rest of the presented analysis, we use the parameters from Fit 2 in Table 2 with the charm and bottom contributions included.
The dipole cross section (2.9) from Fit 2 is shown in Figure 2 (left) as the dashed red lines corresponding to different values of x = 10 −2 , . . . , 10 −6 (from right to left). The blue solid lines show to the dipole cross section (2.15), plotted in the whole range of r. Geometric scaling in this dipole cross section is clearly visible as the single blue line in the right plot. Since Eq. (2.15) is a limiting form the dipole cross section (2.9), it is not surprising that geometric scaling is also visible for the dashed curves in the region rQ s ≥ 1, together with its violation for rQ s < 1.
Summarising, the DGLAP improved model allows to extend the GBW saturation model to large values of Q 2 , giving a good description of the HERA data up to Q 2 = 650 GeV 2 . This is achieved by the modification of the dipole cross section for small dipole  Table 1 and 3, respectively. The saturation scale from the original GBW model with charm [26] is shown as the dash-dotted line.
sizes to match the perturbative QCD result in this region. The essential elements of the GBW model, the saturation scale and geometric scaling, are retained in the DGLAP improved dipole cross section in the region rQ s ≥ 1, which is mostly responsible for the transition of F 2 to small Q 2 values.

Fit comparison
In Figure 3 we compare the normalised-to-one dipole cross sections from the GBW (blue solid lines) and DGLAP improved (red dashed lines) models with the parameters from Fit 2 in Table 1 and 3, respectively. We also use the appropriate saturation scales, given by Eqs. (2.7) and (2.16), for the scaling variable,r ≡ rQ s . We see in the left plot that for large values of r the two functions overlap while they differ in the small-r region where the running of the gluon distribution starts to play a significant role. This is clearly seen in the right plot where geometric scaling holds for the DGLAP improved model curve for the scaling variabler ≥ 1 and in the whole region for the GBW model curve. The two model functions overlap forr ≥ 1 due to their universal form, σ dip ∼ 1 − exp(−r 2 /4), in this region.
Notice also that the leftmost dashed curve of the DGLAP model in the left plot in Figure 3, corresponding to x = 10 −2 , lies below the analogous GBW curve for small r. This is due to the suppression term (1 − x) α present in the gluon distribution. Such a term is missing in the GBW dipole cross section, where the saturation scale is always proportional to x −λ .
It is also interesting to compare the saturation scales from the fits. They are shown in Figure 4 in the two analysed models as the blue solid (GBW model) and red dashed (DGLAP improved model) lines. For the reference, we also plot the saturation scale from the original GBW model [26] with charm as the blue dot-dashed line. We see that all lines lie close to each other, although their slope is slightly different. As a result, at x = 10 −6 the saturation scale Q 2 s ≈ 2 − 3 GeV 2 . The structure function F 2 computed with the two dipole cross sections as a function of Q 2 for fixed values of x is shown in the left plot in Figure 5. We see that for Q 2 > 10 GeV 2 , F 2 computed in the DGLAP improved model rises stronger for small values of x. This is an effect of the DGLAP evolution of the gluon distribution in the dipole cross section. In the right plot, we show the logarithmic slope of F 2 in x as function of Q 2 , found from the approximate relation    Table 1 and 3, respectively. The GBW model lines are shown only in the region where the model was fitted.

Comparison to data
In Figures 6-8, we present the comparison of the predictions from the two models discussed above with the newest HERA data. Finally, in Figure 6, we show the structure function F 2 as a function of x for the indicated values of Q 2 . We show separately the comparison with the data in the region up to Q 2 ≤ 10 GeV 2 (left plot) and above this value (right plot). This is to indicate that the GBW model fits were only performed in the first region. We see that good values of χ 2 /N dof of the fits are reflected in a very good agreement with the data. In the region of small and moderate values of Q 2 (left plot), the two models give similar results. Only for the lowest values of Q 2 , the results slightly differ because of nonzero light quark mass in the GBW model. For lager values of Q 2 (right plot), the DGLAP improved model curves also agree very well with the data.
A similar comparison, now for the charm and bottom contributions, F c,b 2 respectively, to the structure function F 2 , is shown in Figure 7 for the HERA data [47] with Q 2 ≥ 5 GeV 2 . Note that the combined data from ZEUS and H1 for the reduced cross section for charm production, σ cc red , was published in [48]. We stick, however, to the comparison with the structure function data. The model values of  Table 1 and 3, respectively. The GBW model lines are only shown in the region where the model was fitted and slightly above.
determined from the fits to the data on the total structure function F 2 . We find good agreement with the data for both the DGLAP improved model (red dashed curves) and the GBW model in the region of Q 2 ≤ 10 GeV 2 (blue solid lines). Finally, in Figure 8 we show the comparison with the HERA data on the longitudinal structure function F L [49]. Both models give predictions which are in the right ballpark. We have to remember that the GBW model was fitted to the data with Q 2 max = 10 GeV 2 , and the curves above this value are only extrapolations. The difference between the two set of curves can be attributed to the lack of the suppression term (1 − x) α in the GBW model structure functions. This term is relevant for the values of x ∼ 10 −2 and the three data points with the largest Q 2 in Figure 8 have such values. The experimental precision of the data prevents us from drawing more precise conclusions.

Conclusions
The new data [1] on the proton structure function F 2 , extracted from the HERA measurements [2,3], prompted us to address the question how the saturation model of DIS [26][27][28] describes this data for small values of the Bjorken variable, x ≤ 10 −2 , and Q 2 in the range between 0.065 GeV 2 and 650 GeV 2 .  Table 1 and 3, respectively, against the HERA data from [49].
In the first part of our analysis, we only considered the data on F 2 for low and moderate values of Q 2 ≤ 10 GeV 2 to fit three parameters of the GBW model [26], which was originally devised to describe the transition of F 2 to small values of Q 2 < 1 GeV 2 . We found a good fit quality with the charm and bottom quark contributions to F 2 included in the fits. The refitted parameters are close to the original values from [26] with around 10 % decrease of the power λ in the saturation scale (2.7), see Table 1. However, as shown in Figure 4, the position of the saturation line on the (x, Q 2 )-plane has not changed significantly. Therefore, we conclude that the GBW model still provides an economical description of the data on the DIS structure functions, F 2 , F L and F c,b 2 , for the values of Q 2 ≤ 10 GeV 2 . In the second part of our study, we used the DGLAP improved model of the dipole cross section with saturation, which allows one to perform a good quality fits to the new data with Q 2 ≤ 650 GeV 2 . In this model, the small dipole size part of the dipole cross section is modified by the presence of the gluon distribution evolved with the DGLAP evolution equation. We proposed a new prescription (2.17) for the scale of the gluon distribution which gives better results than the original scale in [27,28]. As a result, we obtain a dipole cross section which for large dipoles retains the features of the GBW model with the saturation scale given effectively by Eq. (2.16), see Figure 3. This scale is close to the scale from the GBW model, see the red dashed line in Figure 4. Thus, we showed the robustness of the analysed saturation models, verified by the agreement of with the data.
In conclusion, the saturation model has stood the test of time and is still a valuable contribution to our understanding of DIS at small x 1