Saturation QCD predictions with heavy quarks at HERA

The measurement of the proton structure function at HERA is often seen as a hint for the observation of saturation in high-energy QCD e.g. through the observation of geometric scaling. Accordingly, the dipole picture provides a powerful framework in which the QCD-based saturation models can be confronted to the data. In this paper, we give a parametrisation of proton structure function which is directly constrained by the dynamics of QCD in its high-energy limit and fully includes the heavy quark effects. We obtain a good agreement with the available data. Furthermore, to the contrary of various models in the literature, we do not observe a significant decrease of the saturation momentum due to the heavy quark inclusion.


I. INTRODUCTION
The study of the high-energy limit of perturbative QCD has lead to strong predictions for the scattering amplitudes. One of the most important result is the property of geometric scaling which is a consequence of saturation [1,2,3,4] which extends into the dilute regime where the amplitude is far from the unitarity limit. We shall recall later in this paper how this general property can be proven in perturbative QCD from the Balitsky-Kovchegov (BK) equation [5] or from the Colour Glass Condensate (CGC) formalism [6].
At small x, the confrontation of those predictions with the experimental measurements of the proton structure function at HERA can be achieved within the framework of the dipole model. In the dipole frame, the virtual photon fluctuates into a qq pair of flavour f and size r which then interacts with the proton: The photon wavefunction can safely be computed in QED and is found to be is thus expressed in terms of the dipole-target cross-section σ dip = 2πR 2 p T (r, x) where T is the dipole-target scattering amplitude as entering the QCD evolution equations. R p is often referred to as the radius of the proton and is to be taken as a free normalisation parameter of our model. We are therefore left with the parametrisation of T (r; Y ), with Y = log(1/x) called the rapidity.
This kind of approach is not new in itself. Different approaches to parametrise the dipole-proton scattering amplitude has already been proven successful. One can cite e.g. the pioneering work of Golec-Biernat and Wusthoff [7] which can be improved by adding the collinear DGLAP effects [8,9]. Those approaches, directly based on the gluon distribution function xg(x, Q 2 ), give good description of the data, even including the contributions from heavy quarks [10]. From a different point of view, some other successful approaches [11] use Regge parametrisations for the dipole-proton scattering amplitude.
Beside those various approaches, perturbative QCD provides definite predictions for T , including its approach to saturation. This is perfectly suited for this kind of problem. Indeed, the factorisation formula (1) underlying the dipole picture is valid only in the high-energy (small-x) limit. Hence, using prediction from QCD at high energy to parametrise the dipole-proton scattering amplitude appears as a natural way to proceed. Those predictions that we shall recall in is paper have been successfully gathered into a parametrisation for σ dip (r, x) and tested against the HERA data [12], so far including the contributions from light quarks only. All those approach suggest a saturation scale Q s , the energy-dependent momentum scale below which the amplitude is saturated, of order 1 GeV 2 for x ∼ 10 −4 − 10 −5 at HERA. However, it is rather well-known that, once including the heavy quarks (mainly the charm), all approaches observe a decrease of the saturation scale. Given those a priori important effects of the heavy quarks on the saturation, it appears important to reconsider the predictions of the QCD at high energy to include those contributions. In this paper, we will obtain two results: first that it is possible to accommodate the predictions from [12] with heavy-quark contributions and, second, that this does not lead to a decrease of the saturation scale.
In the next section, we shall recall how to build the dipole-target scattering amplitude T from the equations of QCD at saturation. We also discuss the idea allowing for the inclusion of heavy quarks. We shall then present the fit to the F 2 HERA data in itself, including data selection and parameter adjustment. As a conclusion, we shall finally discuss our results w.r.t. other models as well as with predictions from NLO BFKL [13,14,15].

II. QCD PREDICTIONS FOR THE AMPLITUDE AND HEAVY QUARKS
Since our approach is mainly based on the QCD fit introduced by Iancu, Itakura and Munier (IIM) [12], we start by a presentation of that model. The formula we use to parametrise the dipole-proton scattering amplitude is obtained from our knowledge of the solutions of the BK equation which captures the main ingredients of the high-energy physics with saturation effects. The exact solution to that equation is not known but its asymptotic behaviours, for large and small dipole sizes has been studied in details.
In the last years, it has been shown [4] that the BK equation lies in the same class of universality than the Fisher-Kolmogorov-Petrovsky-Piscunov (F-KPP) equation. The latter has been extensively studied in statistical physics over the past seventy years and it is well-known that its solutions can be written in terms of travelling waves. In the language of the QCD variables we have used so far this means that if one looks at the rapidity evolution of the amplitude T (r; Y ) (seen here as a function of r), the amplitude "front" moves towards smaller values of r without changing its shape. The "position" of the wavefront is then naturally associated with Q −1 s , the inverse of the saturation scale. At asymptotic rapidities, the amplitude T , initially a function of r and Y independently becomes a function of the single variable rQ s (Y ). This very important consequence of saturation is the geometric scaling [17,18,19] property which physically means that the physics remains unchanged when one moves parallel to the saturation line.
This travelling-wave analysis provides two fundamental pieces of information: first, the saturation scale increases exponentially with the rapidity 1 : Q 2 s (Y ) ∝ exp(λY ). Then, the amplitude is known in the small-r region: with ρ = log(4/r 2 ) and ρ s = log(Q 2 s ). The critical slope γ c as well as the parameter λ =ᾱχ ′ c and χ ′′ c are determined from the linear BFKL kernel only. This is an important property: though the scattering amplitude fully satisfies the unitarity constraints and is sensitive to saturation effects, the parameters which describes it do not depend upon the details of how saturation is encoded. We will discuss the value for those parameters later in this section. The fact that they do not depend on the details of the saturation mechanism is another interesting feature of this approach. As a related comment, one can also obtain [1,2] the result (2) by looking at the BFKL equation with a boundary condition at the saturation line: performing a saddle-point approximation gives (2) as well as the exponential behaviour of the saturation momentum. Note finally that the Gaussian part of the exponential in (2) violates geometric scaling as it introduces an explicit dependence in Y . This term however becomes less and less important as rapidity increases. It controls how geometric scaling is approached and allows to deduce that geometric scaling is valid within a window ρ − ρ s 2ᾱχ ′′ c Y . This is an important point that saturation effects are relevant up to large scales above the saturation momentum i.e. in the dilute domain.
The amplitude in the saturated domain is also obtained from the BK equation [20] (it can also be obtained from the Colour Glass Condensate formalism [21]). Putting it together with equation (2) we reach the final expression for our dipole-proton scattering amplitude: where the parameters a and b are fixed so as to ensure that T and its derivative are continuous at rQ s = 2. Going from (2) to (3), we have usedᾱχ ′′ Let us now discuss the parameters in this QCD-saturation-based model. First, the saturation scale can be written This leaves two free parameters: x 0 which is related to the value of the saturation scale at zero rapidity and corresponds to the value of x at which Q s = 1 GeV, and λ which controls the rapidity evolution of the saturation momentum. While leading-order (LO) BFKL predicts λ =ᾱχ ′ c ≈ 0.9, an analysis of the next-to-leading order (NLO) BFKL [3] gives λ ∼ 0.3, a value which is in much better agreement with the phenomenological analysis.
Concerning the parameters in the amplitude itself, we shall fix the matching amplitude T 0 . As in [12], a default value of 0.7 gives good results and variations around that value leads only to small differences, so we will adopt T 0 = 0.7. The value of κ will be set from the LO BFKL kernel which gives κ ≈ 9.9. Note that the NLO BFKL predictions, though a little bit smaller, remain of the same order.
The value of the critical slope γ c is a fundamental issue of this paper. In the original work [12], where only the light quarks were considered, the authors has fixed it to the value obtained from the LO BFKL kernel (γ c ≈ 0.6275). However, as we shall see in section III, when including the heavy quarks, keeping that value leads to a dramatic decrease of the saturation momentum together as well as to a poor χ 2 for the fit. A key issue of the present work is to show that, allowing that parameter to vary, we recover a similar saturation scale and a good fit. In addition, we shall see that the value for γ c coming out of the fit is rather close to what we expect from NLO BFKL (γ c 0.7).
The last parameter entering the dipole-proton cross-section σ dip (r, x) is the radius of the proton R p which fixes the normalisation w.r.t. the dipole-proton amplitude T (r, x). We are thus left with 4 parameters: R p , x 0 and λ which were already present in [12], and the new one: γ c . In the massless case, it was found, for T 0 = 0.7, λ = 0.253, x 0 = 2.67 10 −5 and R p = 0.641 fm for γ c fixed to 0.6275.
Last but not least, we have to specify our heavy-quark prescriptions. While in [12] the sum in (1)  Now that we have fully described the saturation-based model we are using to describe the DIS structure function, we will test its validity by fitting its free parameters to the experimental measurements of the proton structure function F 2 . In this section we present the details of the fit and discuss the results.
We first have to specify which dataset we are working with. Following the recent analysis, we shall use the last HERA data i.e. the last ZEUS [23] and H1 [22] measurements of F 2 . We include a 5% renormalisation uncertainty on the H1 data to account for a normalisation discrepancy between ZEUS and H1. Note that the analysis in [12] only takes into account the ZEUS data. We shall come back on this point later, when we turn to the discussion of our results.
Our approach is focused on the high-energy behaviour of DIS: the dipole-model factorisation 1 is only valid at sufficiently small x and, accordingly, the dipole-proton amplitude is build from the high-energy QCD equations. Hence, we shall limit ourselves to x ≤ 0.01, the usual cut in those approaches. Furthermore, our approach, based on the QCD equations describing saturation, does not takes into account the DGLAP corrections beyond the doublelogarithmic approximation. We shall restrict our analysis to Q 2 ≤ 150 GeV 2 . Note that the former analysis with massless quarks [12] is more conservative and uses Q 2 ≤ 45 GeV 2 . However, in order to check the validity of the results discussed hereafter, we have tested both cuts and observed that the χ 2 and the parameters were not significantly changing when increasing the Q 2 cut.
We have fitted the free parameters of the model discussed in Section II to the 281 data contained in our dataset. In order to grasp the effect of the correct treatment of the heavy-quark masses, we performed the fit with and without   including the charm and bottom contribution to F 2 . In addition, for both cases we give the result of the fit for the critical slope γ c fixed to its LO value or considered as a free parameter. The resulting parameters and χ 2 are presented in table I where we have also added the initial parameters from [12] (with ZEUS data only and Q 2 ≤ 45 GeV 2 ) for better comparison.
Those results deserve some comments: • Concerning the re-analysis of the fit without heavy quarks (second line of table I), we see that with the addition of the H1 data and the extension of the Q 2 domain, the parameters remain similar and the χ 2 is good.
• If one allows γ c to freely vary in the massless case (third line of table I), again, the fit naturally converges to a minimum which is close to the initial one without improving significantly the χ 2 . The LO choice γ c = 0.6275 is even compatible with the error bars of what we obtain when we fit it.
• Once the heavy quarks are taken into account, the situation changes drastically. If the critical slope is fixed to its LO value, the situation becomes dramatic. Indeed, not only the quality of the fit is getting worse, but also the saturation scale is going down by two orders of magnitude (the exponent λ is also decreasing significantly).  [12] (light flavours only, dashed lines), we plot two information: the saturation scale Q 2 s (x) and the extension of the geometric scaling window log(Q 2 /Q 2 s ) ≤ √ 2κλY . We see that the inclusion of the heavy quarks does not lead to a strong decrease both those scales. The points on the plot represent the (x, Q 2 ) position of the HERA measurements.
• If as anticipated in Section II we allow the critical slope γ c to vary, the fit (last line of table I) converges back to a good description (even a better χ 2 than the corresponding massless fit). The parameters describing the saturation scale are also rather close to those obtained in the massless case.
• The value of the critical slope we obtain from the fit seems much larger than the LO result used previously. However, if one extract that value from various renormalisation-group-improved NLO BFKL kernels [15] one get a value of γ c slightly larger than 0.7. This is again in good agreement with the value to which the new fit naturally converges.
• To test the dependence upon the fixed parameters of our model, we have performed various fits varying those parameters around their default value. We have thus changed the Q 2 cut from 150 to 45 GeV 2 , included only the heavy charm (without bottom), varied the masses of the heavy quarks and varied the matching point T 0 . For all those variations, both the quality of the fit and the values of the parameters remained similar, which enforces the robustness of the present parametrisation.
The last line of table I, giving a description of the saturation effects in DIS including the heavy quark effects, has to be considered as the main result of this paper. The corresponding description of the proton structure function is plotted in figure 1.
One of the most interesting point of this parametrisation is that the saturation scale obtained with heavy quarks included is very similar to the one obtained with light quarks only. This is better seen in figure 2 where we have plot the saturation scale (lower curves) as well as the limit of the geometric scaling window both with and without heavy quarks contributions. One clearly see that the addition of the heavy quark contribution only slightly reduces the effect of saturation. This is an important result as previous models including heavy quark effects all report a decrease of the saturation momentum by (roughly) a factor of 2 once those heavy quarks contributions are included. From figure 2 it also appears that a large number of data (all data from small Q 2 up to the limit of the geometric scaling window) are sensitive to saturation effects.
Finally, we can compare the predictions of our parametrisation with the HERA measurements [24,25]    the charm and bottom structure functions respectively. In both cases, we observe a good agreement with the data. Similarly, by taking the contribution coming from the longitudinal part of the wavefunction in (1), we can obtain predictions for the longitudinal structure function. Our result is shown in figure 5 together with the H1 measurements [22]. Again, the present parametrisation gives a good description of the data.

IV. DISCUSSION AND CONCLUSIONS
In this paper, we have shown that it was possible to accommodate the saturation model introduced by Iancu, Itakura and Munier [12] to take into account the heavy quark contribution to the proton structure function. The resulting fit provides a very good description of the HERA measurements of F p 2 and the predictions for the heavy-quark structure functions as well as for the longitudinal structure function are in good agreement with the existing data. The present work, using the dipole picture formalism, has two important features that distinguish it from previous studies in the literature: firstly, following the idea in [12], we have directly use the predictions from high-energy QCD to parametrise the dipole-proton amplitude. The saturation properties, e.g. geometric scaling and its window of validity, are thus parametrised as they are predicted from perturbative QCD. Secondly, the saturation scale that result from our fit is not significantly reduced compared to the saturation scale obtained with light quarks only. This contrasts with previous studies in the literature which report a decrease of the saturation momentum. In the present approach saturation effects cannot be neglected. This new analysis gives an additional argument in favour of saturation.
We also stress a recent analysis [26] also including heavy-quark effects and based on the high-energy QCD properties, where the amplitude is parametrised in momentum space rather than in coordinate space. Compared to that study, the parameters for the saturation scale and geometric scaling window extension, both slightly too small in the momentumspace parametrisation, are more reasonable in the present approach. The analysis in momentum space was however carried with fixed γ c and further studies are requested to check whether a higher critical slope also improves the model.
A last comment concerns the relation between the parameters obtained in our fit and the predictions from the NLO BFKL kernel with a saturation boundary. Indeed, though a bit smaller, the exponent of the saturation scale remains close to the predictions from the NLO BFKL [3]. In addition, we have seen that the critical slope γ c kept fixed to its leading-order value in the massless case, is no longer in agreement with that value once the heavy quarks are taken into account. The value we obtain in that case is rather in agreement with the slope one would obtain from the renormalisation-group-improved BFKL kernels at NLO. This finding is another welcome outcome of our parametrisation.
The last point naturally suggests to directly compare the predictions of NLO BFKL with saturation effects to the proton structure function. However, to rigorously address that question, one probably also needs to correctly introduce the running-coupling corrections to the present formalism. Though they are not expected to lead to significant modifications, they come with a few changes that we leave for future work.
Finally, now that the coordinate-space dipole-proton scattering amplitude is available as directly predicted from the QCD saturation framework, applications to other observables are to be done. This includes the diffractive structure function often considered as an excellent candidate for the observation of saturation and well-described in the dipole picture [7,10,11,27,28]. In addition, extending to non-zero momentum transfer (also predicted by high-energy QCD [29]) or parametrising the impact-parameter dependence, it can also be tested against the diffractive vector-meson production and DVCS cross-sections [30,31,32]. Those complementary analysis are also left for future work.