In depth analysis of the combined HERA data in the dipole models with and without saturation

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ In depth analysis of the combined HERA data in the dipole models with and without saturation © 2018 the Author(s) Published version Mäntysaari, Heikki; Zurita, P. Mäntysaari, H., & Zurita, P. (2018). In depth analysis of the combined HERA data in the dipole models with and without saturation. Physical Review D, 98(3), Article 036002. https://doi.org/10.1103/PhysRevD.98.036002 2018 In depth analysis of the combined HERA data in the dipole models with and without saturation H. Mäntysaari Department of Physics, University of Jyväskylä, P.O. Box 35, 40014 University of Jyväskylä, Finland P. Zurita Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA (Received 24 April 2018; published 1 August 2018) We present an updated impact parameter dependent saturation model determined through a fit to the combined HERA I and I þII reduced cross section data. The same HERA data are used to fit the linearized version of the applied dipole amplitude, which makes it possible to estimate the magnitude of the saturation effects in various experiments. We find that both parametrizations provide comparable descriptions of the considered data when an effective confinement scale dynamics is incorporated with quark masses. Moreover, it is possible to consistently determine the light and charm quark masses. The role of potentially nonperturbatively large dipoles is examined in detail, with the result that, especially in the case of the structure function F2, their contribution is numerically significant. The potential to discriminate between the two models in future eþpand eþAexperiments is also illustrated. DOI: 10.1103/PhysRevD.98.036002 I. INTRODUCTION Quantum chromodynamics, the theory of the strong interaction, is a vast field with a plethora of diverse phenomena still unexplored. Despite being the object of study of the high energy and nuclear physics communities both theoretically and experimentally, the true nature of a proton,itsconstituents,their interactions,and how theycome together to conform it remain elusive. In the collinear framework at a given resolution scale Q2, the proton can be described as composed of quarks and gluons carrying a fraction xof the proton momentum. Once known at an initial scale Q2 0, the partonic densities can be determined at any scale Q2using the Dokshitzer-Gribov-Lipatov-AltarelliParisi (DGLAP) evolution equations [1–4]. This picture is successfully supported by extensive experimental evidence; however, it cannot be valid for all the kinematic range: as one explores lower xvalues, the DGLAP equations predict an infinite rise of thegluondensity which would break unitarity. It follows that some phenomena, e.g. gluon recombination, have to enter in order to tame this dangerous behavior. This dynamically generates the saturation scale Q2 s, which determines when the transition to the nonlinear regime of QCD takes place. Moreover, if Q2 sis large, perturbative calculations become possible as the strong coupling constant isweak.Asuccessful theoreticalframeworktodescribeQCD in this region is known as the color glass condensate [5]. There are many theoretical models that incorporate saturation to QCD calculations with different approaches and considerations, a popular one being the impact parameter dependent saturation model (IPsat) parametrization [6–10]. The parameters that provide the necessary nonperturbative input to these models are determined through fits to available data, the bulk of them being high precision deep inelastic scattering (DIS) data measured at HERA in electron-proton (eþp) collisions [11–14]. Despite thewide kinematic range covered by this collider there are no spectacular signals of deviation from the DGLAP predictions and some observed discrepancies might be due to reasons other than saturation. Recently a possible hint of a nonlinear regime from HERA data at low xwas shown to be feasibly explainable by the inclusion of resumed logarithmic corrections [15]. Saturation model calculations have also been able to provide a natural description for the nearly flat center-of-mass energy dependence of the diffractive to total cross section ratio at HERA [16,17] (see also Ref. [18]). In general, there is no clear consensus on whether the onset of saturation has been reached, and it will be necessary to perform a thorough and detailed exploration of the kinematic space beyond our current knowledge in order to observe the nonlinear regime of QCD. Future facilities such as the Electron-Ion Collider (EIC) [19,20], the Large electron-Hadron Collider (LHeC) [21], and the Future Circular Collider (FCC-eh) [22] hold the key to this door. Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 98, 036002 (2018) 2470-0010=2018=98(3)=036002(17) 036002-1 Published by the American Physical Society In this study we present a new determination of the IPsat and its linearized version (“IPnonsat”) description of the HERA combined data [11–14] in the framework of the dipole model. What is new here, compared to the previous literature [6–8], is that we also fit the IPnonsat model parametrization to the precise combined HERA data which allows us to explore the expected magnitude of saturation effects in current and future collider experiments. In addition, by simultaneously fitting the total cross section and the charm contribution to it, it becomes possible to determine the quark masses in this framework. For consistency, a variable flavor number scheme is also applied. This work is organized as follows. In Sec. II we describe the inclusive photon-proton interaction in terms of the dipole model for both IPsat and IPnonsat formulations. The analysis of the combined inclusive and charm data from HERA is the subject of Sec. III, while in Sec. IV we present an analysis of the obtained dipole-proton scattering amplitude. The application of our determined parameters to the exclusive production of vector mesons is discussed in Sec. V. In Sec. VI we consider the potential of the EIC and the LHeC to provide a signal of saturation in both eþpand eþAcollisions. Finally Sec. VII summarizes our findings. II. PHOTON-PROTON SCATTERING IN THE DIPOLE PICTURE The most precise study of the proton structure has been performed in eþpDIS experiments at HERA [11–14]. The process is described as the electron emitting a virtual photon with momentum q, which then probes the proton with a resolution scale Q2¼−q2. In the dipole picture, applicable at high energy and not too large Q2, the virtual photon-proton scattering process can be factorized in two parts: the γsplitting into a q¯ qpair, and the dipole-target scattering. The total γpcross section is subsequently obtained as the imaginary part of the forward elastic γp→ γpscattering amplitude using the optical theorem. The photon splitting into a dipole with transverse separation ris described in terms of the photon wave function Ψf L;Tðr; z; Q2Þ, where fis the quark flavor, Land Trefer to transverse and longitudinal polarizations, jrj¼r, and z is the fraction of the longitudinal momentum of the photon carried by the quark. The total photon-proton cross section is then given by σγp L;Tðx;Q2Þ¼Zd2bd2rZ1 0 dz 4πjΨf L;Tðr;z;Q2Þj2dσdip d2b;ð1Þ where dσdip d2bis the proton-dipole cross section with b denoting the impact parameter, and one has to sum over all the quark flavors included in the analysis (u,d,s,c, and bin this work). The photon wave functions for the transverse and longitudinal polarizations summed over spins and helicities read [23] jΨTðr; z; Q2Þj2 ¼2Nc παeme2 qf½z2þð1−zÞ2ε2K2 1ðεrÞþm2 fK2 0ðεrÞg ð2Þ and jΨLðr; z; Q2Þj2¼8Nc παeme2 qQ2z2ð1−zÞ2K2 0ðεrÞ;ð3Þ with ε2¼zð1−zÞQ2þm2 q. Here, eqis the fractional charge of the quark qand mqis the quark mass. The proton structure functions F2and FLcan be written in terms of the total photon-proton cross section as F2¼Q2 4π2αem ðσγp Lþσγp TÞ;ð4Þ FL¼Q2 4π2αem σγp L:ð5Þ The most precise combined HERA data are released in the form of the reduced cross section σrðx; y; Q2Þ¼F2ðx; Q2Þ−y2 1þð1−yÞ2FLðx; Q2Þ;ð6Þ where y¼Q2=ðxsÞis the inelasticity of the eþpscattering with center-of-mass energy ffiffiffi s p. In the IPsat model the saturation scale of the target depends on the impact parameter, and the cross section is written as dσdip d2b¼2½1−exp ð−r2Fðx; rÞTpðbÞÞ:ð7Þ The proton density profile TðbÞ¼e−b2=ð2BpÞ=ð2πBpÞis assumed to be Gaussian, and we use fixed Bp¼4GeV−2 based on HERA exclusive J=Ψproduction data. Thus, the effective transverse area of the proton is not a free parameter in the model, and the root mean square radius of the proton is ffiffiffiffiffiffiffiffi 2Bp p. However, we note that this parametrizationdescribes only part of the observed growth of the proton already at the HERA energies [24,25] due to the Gribov diffusion. Including this effect would require us to either parametrize the proton width Bpand try to fit it simultaneously to the HERA data or solve impact parameter dependent small-x evolution equations such as in Refs. [26–28]. We do not want to include exclusive data in our fit due to additional model uncertainties, and we leave it for a future work. At the lowest order in perturbation theory the function F is proportional to the DGLAP evolved gluon distribution function H. MÄNTYSAARI and P. ZURITA PHYS. REV. D 98, 036002 (2018) 036002-2 Fðx; r2Þ¼ π2 2Nc αsðμ2Þxgðx; μ2Þ;ð8Þ with xbeing the momentum fraction of the proton carried by the gluon, and the scale μ2is a function of r2. This parametrization gives the correct perturbative Quantum Chromodynamics (pQCD) limit for the dipole cross section, σdip ∼r2. At large dipoles, the saturation effects are described by having an eikonalized gluon distribution function, which gives dσdip=d2b→2at large r, corresponding with the interaction probability being unity at large dipoles. The scale at which the gluon density and the strong coupling constant are evaluated is chosen to be μ2¼ μ2 0þC=r2. Here, unlike in previous fits [7,8],wefixμ2 0¼ 1.1GeV2and let Cbe a free parameter. This allows us to force μ2to remain in the perturbative region. In our fit we only include data in the Q2bins that satisfy Q2>μ2 0, which guarantees the applicability of the perturbative calculation. We also consider data in the kinematical region x<0.01 where the dipole picture can be considered to be most reliable. The gluon density at the initial scale μ0is parametrized as xgðx; μ2 0Þ¼Agx−λgð1−xÞ6;ð9Þ where Agand λgare free parameters to be determined by the fit. Unlike previous works, we use a variable flavor number scheme when evaluating the strong coupling constant αsand when solving the DGLAP evolution for the gluon distribution. Neglecting the change in the number of flavors and usinga fixed number of quark flavors as thedata movesin Q2 are not the most adequate strategies from the theoretical point of view but, in practice, are solely reflected in different values of the fitted parameters, without a sizable effect in the description of the data in the currently probed kinematic range. For simplicity we refrain from including the quark singlet contribution to the DGLAP evolution which should also be present. However, we have checked that the fit quality and the resulting dipole amplitude are not significantly affected by its inclusion, and that the fit drives the quark singlet contribution to zero at the initial scale. Furthermore we choose the high-xbehavior to be an integer exponent (6) instead of the standard 5.6 in order to speed up the DGLAP evolution performed in Mellin space. We have checked that this exponent does not have a significant impact on the determination of the parameters. The strong coupling constant is required to satisfy αsðMz¼ 91.1876 GeVÞ¼0.1183 [13]. When evaluating the heavy quark contribution, the Bjorken-xis replaced by xq¼x1þ4m2 q Q2ð10Þ in order to take into account the kinematical shift, where q refers to the quark flavor cor b. As we stay in the perturbative (large Q2) region in this work, the shift (10) would have negligible effect in the case of light quarks. The quark masses mffor the light and charm quarks are kept as free parameters and constrained by the fit. In the fit we only include data that satisfy xc<0.01 for the charm quark. The bquark mass is set to 4.75 GeV, and the bquark contribution to the structure function is included if the corresponding Bjorken-xsatisfies xb<0.1. Effectively in the IPsat model we fit the xdependence of the cross section to the HERA data and extrapolate it down to smaller values of Bjorken-xwhen calculating predictions e.g. for the future DIS experiments. The other approach used in small-xphenomenology is to evolve the dipole amplitude in xusing the perturbatively derived BalitskyKovchegov evolution equation [29,30] and incorporate some of the DGLAP effects in the initial condition of the evolution, which is fitted to the HERA data. Currently these fits are done at the leading logarithmic accuracy with running coupling corrections, and a very good description of the HERA data is obtained if the xevolution speed is also adjusted in the fit process by fitting the coordinate space scale at which the running coupling constant is evaluated [31,32] [note that our fit parameter Chas a similar effect by controlling the scale at which we evaluate αsðμ2Þand xgðx; μ2Þ]. Recently, there has been a lot of progress in developing the theory to next to leading order (NLO) accuracy; see e.g. Refs. [33–40]. In order to quantify the magnitude of the saturation effects, we also study the linearized version of the IPsat parametrization (7), dσdip d2b¼2r2Fðx; rÞTpðbÞ;ð11Þ to which we refer as the IPnonsat model. We emphasize that the rigorous way to look for the saturation effects is to compare saturation model calculations with the perturbative QCD results obtained by applying collinear factorization. In practice, however, comparing IPsat and IPnonsat results can be used as a first estimate for the strength of the saturation effects in the given process. In order to enable such a comparison, we fit the IPnonsat model parameters to the same HERA data. III. DESCRIPTION OF THE HERA REDUCED CROSS SECTION DATA The H1 and ZEUS experiments from HERA have published two combined data sets for the reduced cross section: the first one in Ref. [11] with the charm contribution in Ref. [12] where the HERA-I data are combined, and the latest final combined result for the inclusive reduced cross section including all HERA (HERA IþII) data in Ref. [13]. Recently the latest charm reduced IN DEPTH ANALYSIS OF THE COMBINED HERA DATA …PHYS. REV. D 98, 036002 (2018) 036002-3 cross section data from HERA I þII have been made available [14]. We will perform fits to both HERA I and HERA I þII data sets, but we will consider the fit to HERA I data to be our main result, as the charm cross section from HERA I þII is not yet published and the additional data in the newer data set mostly affects the high-Q2region not included in the analysis. Moreover, the newer total reduced cross section data set has more than twice as many data points in the region of interest of this work, which renders that data less sensitive to the charm quark if one does not introduce artificial weight factors (the HERA I þII charm reduced cross section data has as many points as in the HERA I results). We include data in the region 1.5<Q 2<50 GeV2.The lower limit, which we require to be larger than μ2 0, guarantees that there is a large scale justifying the perturbative calculation. As the validity of the dipole picture becomes questionable at high Q2, we only include data up to Q2¼50GeV2 inourmainfit,thoughwealsoshow resultsforfitsdoneinthe larger virtuality range with Q2 max ¼500 GeV2. The free parameters in this work are Agand λgthat describe the gluon distribution at μ2¼μ2 0, and Cthat controls the momentum space scale corresponding to the given dipole size jrj. In addition, and as mentioned in the previous section, the light quark and charm quark masses are to be determined by the fit. However, as the bottom quark contribution is small, the fit cannot reliably determine the bquark mass, and we set it to mb¼4.75 GeV. The dependence of the fit quality on the light quark mass mlight is shown in Fig. 1when we fit the HERA I data [11,12]. Throughout this work we show results obtained by using the fit to the HERA I data. Here, the charm mass and all the other parameters are allowed to vary freely with the fixed light quark mass. As one moves to lower values of mlight in the IPsat fits, the χ2reaches a plateau, making it hard to determine a best fit extraction of its value, similarly as in Ref. [8]. Therefore, and in order to have a finite quark mass to act as an infrared regulator, we fix mlight ¼ 0.03 GeV for the IPsat case. The situation is different in the IPnonsat fit, where a relatively large light quark mass ∼0.14 GeV is preferred. This can be interpreted as an effective confinement requirement. The effect of a nonzero quark mass is to suppress dipoles larger than ∼1=mlight, which in the IPnonsat model have an unphysically large (unitarity violating) cross section. The final fit quality in both IPsat and IPnonsat models is similar, suggesting that, when describing the inclusive DIS data, the effective confinement effect in the IPnonsat has a comparable effect as the gluon saturation in the IPsat parametrization. Unlike previous works [7,8] we find that the fit clearly prefers a charm quark mass ∼1.35 GeV, with the χ2 presenting a clear minimum. One main difference of our work with that of Ref. [8] is that the charm data are included in the fit, which allows us to constrain the charm mass simultaneously with the other parameters. The quality of the fit as a function of the charm mass is shown in Fig. 2, where all parameters are allowed to vary (except the light quark mass which is fixed to 0.03 GeV in the IPsat model), while keeping mcconstant. In Table Iwe present the fitted parameters for two different values of Q2 max. As can be seen from the χ2, the description of the precise HERA data is excellent, as already noted in previous works [6–8]. What is new here compared to the previous literature is the IPnonsat model parametrization that we find to describe the combined HERA data equally well.1This is demonstrated in Figs. 3 and 4where the reduced cross section and the charm contribution to it are shown and compared to the IPsat and FIG. 1. Fit quality to HERA inclusive [11] and charm [12] reduced cross section data as a function of the light quark mass. FIG. 2. Fit quality to HERA inclusive [11] and charm [12] reduced cross section data as a function of the charm quark mass. In the IPsat fit the light quark mass is fixed to 0.03 GeV (see text). 1In Ref. [6] the IPnonsat model was fitted to older H1 and ZEUS data that have much larger uncertainties than the combined data set used in this work. H. MÄNTYSAARI and P. ZURITA PHYS. REV. D 98, 036002 (2018) 036002-4 IPnonsat model results, the curves being practically one on top of the other. We are also able to determine the quark masses from the fit. For the charm contribution, we note that there is some tension, the HERA data suggesting slightly slower Q2evolution than what is the outcome of our combined fit. This could be due to shortcomings of the model in describing the heavy quarks, as it happens in the collinear factorized framework where higher order corrections are needed for a proper description of the charm and bottom data [41]. It also might be influenced by the fact that the HERA collider was not particularly tuned to measuring heavy quarks, an issue that will be addressed in future colliders such as the EIC. The fits done to HERA I and HERA I þII combined data sets result in comparable parameters (the largest difference being the scale parameter C, on which the results depend only logarithmically). Also, the newer data set prefers a slightly smaller charm quark mass, but the differences are small. Thus the difference at the level of an observable will be negligible between the two different fits performed to different data sets. We will demonstrate this in Appendix B. We consider the top row of Table Ito be our main fit, as it relies on published data sets and only includes measurements in the kinematical domain where the applied model can be considered to be most reliable. The bottom quark reduced cross section included in the latest combination of HERA heavy quark data [14] is discussed in Appendix A. A. Contribution from large dipoles As discussed above, especially in the IPnonsat model, a nonzero effective light quark mass is needed in order to obtain a satisfactory description of the HERA data. This means that the reduced cross section data are sensitive to the contributions from (possibly nonperturbatively) large dipoles, whose formation should be suppressed by confinement scale effects. In the IPsat parametrization, the TABLE I. All dimensionful parameters in GeV. XþYpoints means Xpoints for σrand Ypoints for σr;charm. Bottom mass is mb¼4.75 GeV, and Bp¼4GeV−2. In the IPsat fit the light quark mass is fixed to prevent numerical instability. The starting scale for the DGLAP evolution is μ2 0¼1.1GeV2. Fit results with HERA I data [11,12] and HERA I þII data [13,14] are shown separately. Type HERA χ2=N NQ2 min Q2 max mlmcCA gλg IPsat I 1.0978 156 þ33 1.5 50 0.03 1.3528 2.2894 2.1953 0.08289 IPsat I þII 1.2781 410 þ33 1.5 50 0.03 1.3210 1.8178 2.0670 0.09575 IPsat I 1.2634 229 þ42 1.5 500 0.03 1.3296 2.6477 2.2097 0.07795 IPsat I þII 1.3014 609 þ42 1.5 500 0.03 1.3113 2.3700 2.1394 0.08388 IPnonsat I 1.122 156 þ33 1.5 50 0.1516 1.3504 4.2974 3.0391 −0.006657 IPnonsat I þII 1.3023 410 þ33 1.5 50 0.1497 1.3180 3.5445 2.8460 0.008336 IPnonsat I 1.2194 229 þ42 1.5 500 0.1332 1.3187 5.6510 3.2820 −0.03460 IPnonsat I þII 1.2944 609 þ42 1.5 500 0.1359 1.3047 4.7328 3.0573 −0.01656 FIG. 3. Inclusive reduced cross section from the fit including data up to 50 GeV2, compared with the HERA data [11]. IN DEPTH ANALYSIS OF THE COMBINED HERA DATA …PHYS. REV. D 98, 036002 (2018) 036002-5 imposed unitarity requirement limits the scattering probability not to exceed unity at large dipoles, but large dipoles can still have a numerically significant contribution. The fractional contribution from large dipoles to the total F2and FLis shown in Figs. 5and 6, respectively. These calculations are done using the IPsat fit (first line of Table I); using the IPnonsat fit would result in very similar rmax dependence. For F2, even at relatively large Q2∼500 GeV2, 10% of the total structure function originates from dipoles larger than 1 fm. On the other hand, the HERA reduced cross section data have relative uncertainties at the percentage level, much smaller than the contribution we obtain from the nonperturbatively large dipoles. The reason for this large contribution is that there is a large so called aligned jet contribution: in the limits z→0 and z→1the large dipole contribution to the transverse photon-proton cross section is only suppressed by the light quark mass as ∼e−mlightr. This can be seen from the virtual photon wave function, Eq. (2). On the other hand, as can be seen from Fig. 6, in the case of FLat moderate Q2the contribution from the region r≳1fm is negligible. This is due to the extra factor z2ð1−zÞ2in the longitudinal photon wave function, Eq. (3), which suppresses the end point FIG. 4. Charm reduced cross section calculated from the fit that includes data up to 50 GeV2, compared with the HERA data [12]. FIG. 5. Contribution to F2at x¼0.01 from the IPsat model from dipoles smaller than rmax at different Q2. Results for IPnonsat are similar. FIG. 6. Contribution to FLat x¼0.01 from dipoles smaller than rmax at different Q2. H. MÄNTYSAARI and P. ZURITA PHYS. REV. D 98, 036002 (2018) 036002-6 contributions. Thus, we would prefer to fit the FLdata instead of the reduced cross section measurements which are dominated by F2. However, the HERA FLdata [42,43] are not precise enough for a detailed comparison with the dipole model calculations. Future DIS facilities EIC and LHeC have plans to measure proton structure functions (including FL)atan unprecedented accuracy [19–21]. Similarly, studying only the charm contribution to the total cross section limits the contribution from large dipoles as demonstrated in Fig. 7. In the case of the F2;charm, even at small Q2contribution from dipoles larger than ∼0.6fm is negligible (but very small dipoles are not sensitive to the saturation effects either). In general we find that F2,FL, and F2;charm are sensitive to different dipole sizes, and future DIS data covering all these structure functions will provide much more precise constraints. In order to further study how much large dipoles affect the fit result, we perform the fits to HERA I inclusive and charm cross section data up to Q2 max ¼50 GeV2with different cutoffs for large dipoles rmax. The resulting fit quality is shown in Fig. 8. We find that in order to obtain a good fit to the HERA data, we have to include dipoles up to ∼2…2.5fm in the IPsat model. In the case of the IPnonsat parametrization, the fit can compensate the effect of the rmax cutoff as the light quark mass is also a fit parameter; thus the fit quality is more stable with respect to the infrared cutoff. The IPnonsat fit drives the light quark mass to zero when rmax becomes ∼1.6fm, and it is not possible to fit the HERA data with a much smaller cutoff. The dependence of the light quark mass on rmax is shown in Fig. 9which further demonstrates that in the IPnonsat model the inclusion of large dipoles requires a larger light quark mass to suppress the contribution from this unphysical region. In the IPsat model, the fits prefer a zero light quark mass at all rmax. IV. DIPOLE SCATTERING AMPLITUDE The resulting dipole amplitude at b¼0is shown in Fig. 10, where we compare it with the previous results from [8] labeled as “RSKV.”Even though our study has incorporated some refinements (variable flavor number scheme in the DGLAP evolution, quark masses determined by the fit, inclusion of the charm reduced cross section data), the base model is essentially the same, and by fitting similar data sets one is expected to obtain compatible dipole amplitudes, despite some numerical differences in the fitted parameters. In the IPnonsat model the evolution at the initial scale in xis very slow (λgbeing close to 0); thus there is basically no evolution at large r, in the region where the IPsat parametrization has already reached unity (and where the IPnonsat model gives unphysical results). FIG. 7. Contribution to charm structure function F2;c at x¼0.01 from dipoles smaller than rmax at different Q2. Note that the scales are different than in Figs. 5and 6. FIG. 8. Fit quality with different cutoffs for large dipoles. In the IPsat model the light quark mass is fixed to mlight ¼0.03 GeV. FIG. 9. Light quark mass obtained as a fit result with the IPnonsat model as a function of the infrared cut for the large dipoles. IN DEPTH ANALYSIS OF THE COMBINED HERA DATA …PHYS. REV. D 98, 036002 (2018) 036002-7 To demonstrate the evolution of the gluon distribution function we plot xgðx; μ2¼μ2 0þC=r2Þin Fig. 11 as a function of rusing both the IPsat and IPnonsat fitted parameters to initialize the DGLAP evolution. At large scales the two parametrizations have small differences, as the effect of the initial condition is washed out in the evolution. Close to the initial scale μ2 0there is basically no evolution if the IPnonsat model parametrization is used (λg≈0, which forces the dipole scattering amplitude not to grow in the region where it is already violating unitarity), unlike in the case of the IPsat model. At large scales and at sufficiently large x≳10−3it is also visible that the gluon density starts to decrease as the scale μ2∼C=r2increases. This is expected, as the momentum conservation in the DGLAP evolution removes the larger-x gluons as they are splitting into the smaller-xones. Similar results were already found in Ref. [6]. This effect is strong close to x∼10−2, where the decreasing gluon density is probed already by dipoles that have large enough sizes to contribute significantly on F2. The point at which the nonlinear effects become relevant is characterized by the saturation scale Q2 s. To determine it we use the definition Nðr2¼2=Q2 s;x;bÞ¼1−e−1=2:ð12Þ The extracted saturation scale as a function of xis shown in Fig. 12. Here, Q2 sis extracted at the central impact parameter b¼0and at the average hbidefined such that Zhbi 0 dbbTpðbÞ¼1 2Z∞ 0 dbbTpðbÞ:ð13Þ This definition gives hbi≈0.46 fm. The difference between the IPsat and IPnonsat parametrizations remains small at all values of Bjorken-x, the IPnonsat model having in general slightly faster evolution. As expected based on previous analyses (e.g. [8,32]), the saturation scale of the proton is at the ΛQCD range in the region x∼10−2, and the region of Q2 sbeing perturbative is reached below x≲ 10−4…10−5[note that the absolute value of Q2 sdepends on the definition (12)]. V. EXCLUSIVE VECTOR MESON PRODUCTION Additional information about the proton structure can be obtained by studying exclusive processes. They are particularly powerful in probing the gluonic structure, as at leading order in collinear factorization the vector meson production cross section is proportional to the squared gluon distribution [44]. In addition to being sensitive to the total gluonic density, in an exclusive process the FIG. 10. The obtained dipole amplitudes at x¼10−6(red lines), x¼10−4(blue lines), and x¼10−2(black lines). RSKV refers to the previous fit [8]. FIG. 11. Gluon density xg as a function of the dipole size rfor x¼10−4(red lines), x¼10−3(blue lines), and x¼10−2(black lines) from top to bottom. FIG. 12. Saturation scale at the center (thick black lines) of the proton and at average impact parameter hbi≈0.46 fm (thin blue lines). H. MÄNTYSAARI and P. ZURITA PHYS. REV. D 98, 036002 (2018) 036002-8 cross section [14]. Due to the relatively large uncertainties and limited number of data points, we did not include this unpublished data set in our fit. Instead, we can use it to check that the bottom quark contribution included in the calculation of the inclusive cross section is compatible with the current measurements. The predicted bottom quark reduced cross section compared with the HERA data is shown in Fig. 21. As noted for the inclusive and charm reduced cross sections, the IPsat and IPnonsat parametrizations give approximately equal results in the HERA kinematics. The description of the data is also good, the χ2=N being 1.81 (1.90) for the IPsat (IPnonsat) model, when comparing with data points at Q2≤500 GeV2. Slightly faster Q2 evolution is obtained compared with the HERA measurements, similarly as in case of the charm reduced cross section (see Fig. 4). APPENDIX B: EFFECT OF THE HERA I+II DATA In this work we have considered the fit performed to HERA I data to be the main result of this work. As was shown in Table I, the best fit parameters change slightly when one fits the final HERA I þII data set which, in addition to having smaller uncertainties, is more dominated by the inclusive reduced cross section. In order to quantify the effect of different data sets on the fit result, we show in Fig. 22. The obtained dipole amplitudes are found to be very similar over a broad range in x. Consequently, the structure function F2obtained with both parametrizations is practically identical. This is demonstrated by showing in Fig. 23 the F2obtained using the fit result to HERA I þII data (second line in Table I) normalized by the result obtained by applying the HERA I fit result (first line in Table I). Similar results are found in the case of FL. Thus, both fits’results can be considered to be equivalent. [1] V. N. Gribov and L. N. Lipatov, Deep inelastic epscattering in perturbation theory, Yad. Fiz. 15, 781 (1972) [Sov. J. Nucl. Phys. 15, 438 (1972)]. [2] V. N. Gribov and L. N. Lipatov, eþe−pair annihilation and deep inelastic ep scattering in perturbation theory, Yad. Fiz. 15, 1218 (1972) [Sov. J. Nucl. Phys. 15, 675 (1972)]. [3] G. Altarelli and G. Parisi, Asymptotic freedom in parton language, Nucl. Phys. B126, 298 (1977). [4] Y. L. Dokshitzer, Calculation of the structure functions for deep inelastic scattering and eþe−annihilation by perturbation theory in quantum chromodynamics., Sov. Phys. JETP 46, 641 (1977). [5] F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, The color glass condensate, Annu. Rev. Nucl. Part. Sci. 60, 463 (2010). [6] H. Kowalski and D. Teaney, An impact parameter dipole saturation model, Phys. Rev. D 68, 114005 (2003). [7] H. Kowalski, L. Motyka, and G. Watt, Exclusive diffractive processes at HERA within the dipole picture, Phys. Rev. D 74, 074016 (2006). [8] A. H. Rezaeian, M. Siddikov, M. Van de Klundert, and R. Venugopalan, Analysis of combined HERA data in the impact-parameter dependent saturation model, Phys. Rev. D 87, 034002 (2013). FIG. 22. Comparison of the IPsat dipole amplitudes obtained by fits to HERA I and HERA I þII data. The used parametrizations are the first and second lines of Table I. Bjorkenxvalues decrease to the left. FIG. 23. Ratio of the F2structure functions obtained by fits to HERA I þII and HERA I data using the IPsat model. The used parametrizations are the first and second lines of Table I. IN DEPTH ANALYSIS OF THE COMBINED HERA DATA …PHYS. REV. D 98, 036002 (2018) 036002-15 [9] A. Luszczak and H. Kowalski, Dipole model analysis of high precision HERA data, Phys. Rev. D 89, 074051 (2014). [10] A. Luszczak and H. Kowalski, Dipole model analysis of highest precision HERA data, including very low Q2’s, Phys. Rev. D 95, 014030 (2017). [11] F. Aaron et al. (H1 and ZEUS Collaborations), Combined measurement and QCD analysis of the inclusive ep scattering cross sections at HERA, J. High Energy Phys. 01 (2010) 109. [12] H. Abramowicz et al. (H1 and ZEUS Collaborations), Combination and QCD analysis of charm production cross section measurements in deep-inelastic ep scattering at HERA, Eur. Phys. J. C 73, 2311 (2013). [13] H. Abramowicz et al. (H1 and ZEUS Collaborations), Combination of measurements of inclusive deep inelastic epscattering cross sections and QCD analysis of HERA data, Eur. Phys. J. C 75, 580 (2015). [14] H1 and ZEUS Collaborations, Combination and QCD analysis of charm and beauty production cross-section measurements in deep inelastic ep scattering at HERA, Eur. Phys. J. C 78, 473 (2018). [15] R. D. Ball, V. Bertone, M. Bonvini, S. Marzani, J. Rojo, and L. Rottoli, Parton distributions with small-xresummation: Evidence for BFKL dynamics in HERA data, Eur. Phys. J. C 78, 321 (2018). [16] K. J. Golec-Biernat and M. Wusthoff, Saturation effects in deep inelastic scattering at low Q2and its implications on diffraction, Phys. Rev. D 59, 014017 (1998). [17] K. J. Golec-Biernat and M. Wusthoff, Saturation in diffractive deep inelastic scattering, Phys. Rev. D 60, 114023 (1999). [18] H. Kowalski, T. Lappi, C. Marquet, and R. Venugopalan, Nuclear enhancement and suppression of diffractive structure functions at high energies, Phys. Rev. C 78, 045201 (2008). [19] A. Accardi et al., Electron Ion Collider: The next QCD frontier—Understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016). [20] E. C. Aschenauer, S. Fazio, J. H. Lee, H. Mäntysaari, B. S. Page, B. Schenke, T. Ullrich, R. Venugopalan, and P. Zurita, The electron-ion collider: Assessing the energy dependence of key measurements, arXiv:1708.01527. [21] J. Abelleira Fernandez et al. (LHeC Study Group Collaboration), A large hadron electron cCollider at CERN: Report on the physics and design concepts for machine and detector, J. Phys. G 39, 075001 (2012). [22] F. Zimmermann, M. Benedikt, D. Schulte, and J. Wenninger, Challenges for highest energy circular colliders, Report No. IPAC-2014-MOXAA01, 2014. [23] Y. V. Kovchegov and E. Levin, Quantum Chromodynamics at High Energy (Cambridge University Press, Cambridge, UK, 2012). [24] A. Aktas et al. (H1 Collaboration), Elastic J=Ψproduction at HERA, Eur. Phys. J. C 46, 585 (2006). [25] S. Chekanov et al. (ZEUS Collaboration), Exclusive photoproduction of J=Ψmesons at HERA, Eur. Phys. J. C 24, 345 (2002). [26] J. Berger and A. M. Stasto, Small xnonlinear evolution with impact parameter and the structure function data, Phys. Rev. D84, 094022 (2011). [27] J. Berger and A. M. Stasto, Exclusive vector meson production and small-xevolution, J. High Energy Phys. 01 (2013) 001. [28] S. Schlichting and B. Schenke, The shape of the proton at high energies, Phys. Lett. B 739, 313 (2014). [29] Y. V. Kovchegov, Small-xF 2structure function of a nucleus including multiple pomeron exchanges, Phys. Rev. D 60, 034008 (1999). [30] I. Balitsky, Operator expansion for high-energy scattering, Nucl. Phys. B463, 99 (1996). [31] J. L. Albacete, N. Armesto, J. G. Milhano, P. Quiroga-Arias, and C. A. Salgado, AAMQS: A non-linear QCD analysis of new HERA data at small-xincluding heavy quarks, Eur. Phys. J. C 71, 1705 (2011). [32] T. Lappi and H. Mäntysaari, Single inclusive particle production at high energy from HERA data to protonnucleus collisions, Phys. Rev. D 88, 114020 (2013). [33] I. Balitsky and G. A. Chirilli, Next-to-leading order evolution of color dipoles, Phys. Rev. D 77, 014019 (2008). [34] T. Lappi and H. Mäntysaari, Direct numerical solution of the coordinate space Balitsky-Kovchegov equation at next to leading order, Phys. Rev. D 91, 074016 (2015). [35] T. Lappi and H. Mäntysaari, Next-to-leading order BalitskyKovchegov equation with resummation, Phys. Rev. D 93, 094004 (2016). [36] T. Lappi and R. Paatelainen, The one loop gluon emission light cone wave function, Ann. Phys. (N.Y.) 379, 34 (2017). [37] R. Boussarie, A. V. Grabovsky, D. Yu. Ivanov, L. Szymanowski, and S. Wallon, Next-to-Leading Order Computation of Exclusive Diffractive Light Vector Meson Production in a Saturation Framework, Phys. Rev. Lett. 119, 072002 (2017). [38] G. Beuf, Dipole factorization for DIS at NLO: Combining the q¯ qand q¯ qg contributions, Phys. Rev. D 96, 074033 (2017). [39] B. Duclou´e, H. Hänninen, T. Lappi, and Y. Zhu, Deep inelastic scattering in the dipole picture at next-to-leading order, Phys. Rev. D 96, 094017 (2017). [40] H. Hänninen, T. Lappi, and R. Paatelainen, One-loop corrections to light cone wave functions: The dipole picture DIS cross section, Ann. Phys. (N.Y.) 393, 358 (2018). [41] R. D. Ball et al. (NNPDF Collaboration), Parton distributions from high-precision collider data, Eur. Phys. J. C 77, 663 (2017). [42] F. Aaron et al. (H1 Collaboration), Measurement of the proton structure function FLðx; Q2Þat low x,Phys. Lett. B 665, 139 (2008). [43] S. Chekanov et al. (ZEUS Collaboration), Measurement of the longitudinal proton structure function at HERA, Phys. Lett. B 682, 8 (2009). [44] M. Ryskin, Diffractive J=Ψelectroproduction in LLA QCD, Z. Phys. C 57, 89 (1993). [45] J. Bartels, K. J. Golec-Biernat, and K. Peters, On the dipole picture in the nonforward direction, Acta Phys. Pol. B 34, 3051 (2003). H. MÄNTYSAARI and P. ZURITA PHYS. REV. D 98, 036002 (2018) 036002-16 [46] A. G. Shuvaev, K. J. Golec-Biernat, A. D. Martin, and M. G. Ryskin, Off-diagonal distributions fixed by diagonal partons at small xand ξ,Phys. Rev. D 60, 014015 (1999). [47] T. Lappi and H. Mäntysaari, Incoherent diffractive J=Ψproduction in high energy nuclear DIS, Phys. Rev. C 83, 065202 (2011). [48] C. A. Bertulani, S. R. Klein, and J. Nystrand, Physics of ultra-peripheral nuclear collisions, Annu. Rev. Nucl. Part. Sci. 55, 271 (2005). [49] T. Lappi and H. Mäntysaari, J=Ψproduction in ultraperipheral Pb þPb and pþPb collisions at LHC energies, Phys. Rev. C 87, 032201 (2013). [50] H. Mäntysaari and B. Schenke, Evidence of Strong Proton Shape Fluctuations from Incoherent Diffraction, Phys. Rev. Lett. 117, 052301 (2016). [51] H. Mäntysaari and B. Schenke, Revealing proton shape fluctuations with incoherent diffraction at high energy, Phys. Rev. D 94, 034042 (2016). [52] C. Alexa et al. (H1 Collaboration), Elastic and protondissociative photoproduction of J=Ψmesons at HERA, Eur. Phys. J. C 73, 2466 (2013). [53] N. Armesto and A. H. Rezaeian, Exclusive vector meson production at high energies and gluon saturation, Phys. Rev. D90, 054003 (2014). [54] B. B. Abelev et al. (ALICE Collaboration), Exclusive J=ψ Photoproduction off Protons in Ultra-Peripheral p-Pb Collisions at ffiffiffiffiffiffiffiffi sNN p¼5.02TeV, Phys. Rev. Lett. 113, 232504 (2014). [55] F. Aaron et al. (H1 Collaboration), Diffractive electroproduction of ρand ϕmesons at HERA, J. High Energy Phys. 05 (2010) 032. [56] P. Amaudruz et al. (New Muon Collaboration), A reevaluation of the nuclear structure function ratios for D, He, Li-6, C and Ca, Nucl. Phys. B441, 3 (1995). [57] M. R. Adams et al. (E665 Collaboration), Shadowing in inelastic scattering of muons on carbon, calcium and lead at low xBj,Z. Phys. C 67, 403 (1995). [58] M. Arneodo et al. (New Muon Collaboration), The structure function ratios F2ðLiÞ=F2ðDÞand F2ðCÞ=F2ðDÞat small x, Nucl. Phys. B441, 12 (1995). [59] H. Kowalski, T. Lappi, and R. Venugopalan, Nuclear Enhancement of Universal Dynamics of High Parton Densities, Phys. Rev. Lett. 100, 022303 (2008). [60] H. Mäntysaari and R. Venugopalan, Systematics of strong nuclear amplification of gluon saturation from exclusive vector meson production in high energy electron-nucleus collisions, Phys. Lett. B 781, 664 (2018). IN DEPTH ANALYSIS OF THE COMBINED HERA DATA …PHYS. REV. D 98, 036002 (2018) 036002-17