Signatures of gluon saturation from structure-function measurements

We study experimentally observable signals for nonlinear QCD dynamics in deep inelastic scattering (DIS) at small Bjorken variable $x$ and moderate virtuality $Q^2$, by quantifying differences between the linear Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution and nonlinear evolution with the Balitsky-Kovchegov (BK) equation. To remove the effect of the parametrization freedom in the initial conditions of both equations, we first match the predictions for the DIS structure functions $F_2$ and $F_{\rm L}$ from both frameworks in a region in $x,Q^2$ where both frameworks should provide an accurate description of the relevant physics. The differences in the dynamics are then quantified by the deviations when one moves away from this matching region. For free protons we find that the differences in $F_2$ remain at a few-percent level, while in $F_{\rm L}$ the deviations are larger, up to $10\,\%$ at the EIC and $40\,\%$ at the LHeC kinematics. With a heavy nucleus the differences are up to $10\,\%$ in $F_2$, and can reach $20\,\%$ and $60\,\%$ in $F_{\rm L}$ for the EIC and the LHeC, respectively.


I. INTRODUCTION
The high-energy regime of Quantum Chromodynamics (QCD) is characterised by high partonic densities where non-linear phenomena play a dominant role [1]. In this regime, the predictions of usual linear evolution equations for parton densities like the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation [2][3][4][5] which provides an accurate description of the QCD dynamics at moderate to large values of the momentum fractions x of the probed parton and virtualities Q 2 Λ 2 QCD , or the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [6,7] which accurately describes the partonic dynamics at small x, are expected to become less reliable. Highdensity effects should become relevant and unitarity corrections or, alternatively, gluon recombination processes, tame the growth of parton densities towards small x and make them saturate -thus non-linear dynamics is also referred to as saturation. Non-linear modifications to DGLAP evolution equations were first proposed in Refs. [8,9]. Later on, a formulation suitable for QCD at such high parton densities and energies (or small x) was developed as a weak coupling but non-perturbative effective field theory, the Color Glass Condensate (CGC), see Refs. [1,10]. Such approach resums all higher-twist corrections that are neglected in linear DGLAP evolution. In this framework, the BFKL equation for evolution in x is generalized to the Balitsky-Kovchegov one (BK) [11][12][13]. The BK equation includes a nonlinear term needed to preserve the unitarity of the (dipole) scattering amplitude.
Finding such a novel regime of QCD would be of uttermost importance for understanding QCD, for developing more realistic initial conditions for simulations of heavy ion collisions at high energy and for computing the main backgrounds for signals of physics beyond the Standard Model in hadronic collisions, especially at the high energies of the Future Circular Collider [14,15]. Linear, fixedorder DGLAP evolution equations are known to provide a very good description of the structure of the proton at the available values of x and Q 2 , see e.g. Ref. [16]. However, there have also been claims that the HERA data [17,18] at small x and moderate Q 2 necessitates the inclusion of small-x resummation [19,20] or non-linear evolution [21]. 1 Note also that while DGLAP-based fits provide a good description of available data constraining parton densities inside nuclei (nPDFs), see e.g. [23][24][25][26][27][28], non-linear effects are density effects that should be enhanced by the nuclear size. Therefore, searches for saturation with heavier nuclei are crucial not only because the expected effects are stronger, and thereby more easily detectable, but also to eventually check the theoretical explanation of such non-linear dynamics, if found. For the moment, however, even the D-meson production at the LHC that can probe the nuclei down to x ∼ 10 −6 at low Q 2 shows no visible deviation from the DGLAP approach [28,29].
Searches of saturation are ongoing in hadronic colliders [30] -the Relativistic Heavy Ion Collider (RHIC) at BNL and the Large Hadron Collider (LHC) at CERN, see e.g. [31] for future prospects at the LHC. However, with a cleaner environment and direct experimental access to the kinematical variables x, Q 2 , deep inelastic scattering (DIS) in lepton-proton/nucleus collisions offers the ideal environment to study the partonic structure of protons and heavier nuclei. The currently projected and proposed accelerators, the Electron Ion Collider [32,33] (EIC) at BNL, and the Large Hadron-electron Collider [34,35] (LHeC) and the Future Circular Collider in electronhadron mode [14,15] (FCC-he) at CERN, will provide lepton-hadron collisions in the centre-of-mass range from a few tens of GeV to a few TeV per nucleon. Together with the use of new detector techniques, they promise to revolutionise our understanding of hadron structure and partonic dynamics.
In DIS, several observables have been proposed as being sensitive to non-linear dynamics: particle correlations [34,36], diffraction [32,37], and deviations from linear evolution in inclusive observables, in particular the total cross section. Concerning the latter, it is expected that the EIC and the LHeC will be sensitive to saturation, i.e. to deviations from DGLAP due to highertwist effects. However, conclusively distinguishing these two pictures using experimental data can, in practice, be difficult. This is due to the fact that both non-linear and DGLAP-based calculations require non-perturbative inputs that are obtained by fitting experimental data. When the initial conditions of the evolution are fit to the same data, also the predictions will not deviate dramatically. Thus different theoretical frameworks might be equally capable of describing experimental data a posteriori. Additionally, the DGLAP and BK equations predict the evolution of cross sections in different directions: increasing Q 2 and decreasing x, respectively. Thus one cannot solve both starting from the same initial condition: the perturbative QCD prediction of one equation is the nonperturbative initial condition of the other. Therefore, it is important to be careful to distinguish between genuine effects of the different evolution dynamics on one hand, and the ability of the initial condition parametrizations to adjust to data on the other hand.
Several works have investigated the possibility to disentangle DGLAP linear dynamics from saturation in ep collisions at the LHeC [34,35,38] and eA collisions at the EIC [39]. The basic technique employed in these studies is the inclusion of pseudodata from the projected experiments in a DGLAP-based global fit. Then, some measure of the fit quality is computed. Generically, a worsening is observed when the introduced pseudodata are generated using a model that contains non-linear dynamics in comparison to the case when the corresponding pseudodata are produced using linear DGLAP dynamics. Focusing on the difference between the generated pseudodata and existing DGLAP-evolved parton distributions, as in Ref. [39], addresses potential tensions between predictions from nonlinear dynamics and existing PDF parametrizations. The latter are determined not only by the QCD dynamics, but by existing data and to some extent by parametrization choices. Thus the question addressed by such a comparison has many more facets than the simple one we pose here, namely: how different is the actual evolution dynamics once one allows both fits to adjust themselves to a common set of data?
The purpose of this paper is to evaluate the difference between DGLAP and non-linear evolution in the kinematic regions corresponding to the EIC and the LHeC, both for ep and eA collisions, with a slightly different emphasis than previous studies. Our procedure is the following. We first produce results for the F 2 and F L structure functions using non-linear evolution equations. Then, employing reweighting techniques, a DGLAP evolved set of parton densities is produced that results in structure functions coinciding with the non-linear evolution ones in a line of the x, Q 2 plane, as precisely as possible. Technically this is achieved by a Bayesian reweighting process. Here one takes an existing collection of "replica" PDF sets that all satisfy the DGLAP equation, and calculates for each of them a weight defined by how close it is to the given F 2 and F L structure functions. The reweighted structure functions and PDFs are then defined as a weighted sum over all replicas. We will explain this procedure in more detail below in Sec. II C. Now the difference between the non-linear results and the DGLAP ones for F 2 and F L when moving away from this line provides a quantification of the signatures of saturationof the difference between linear and non-linear dynamics. Although we are interested in the experimentally accessible kinematical regimes, we do not use estimated experimental errors in this procedure.
Our approach can be compared to the one used in the previous studies [34,35,38,39], where a DGLAP fit to pseudodata from BK-evolution in the whole kinematical region was performed. There one fits to the pseudodata in a wide region in the kinematical plane, and the weight given to specific kinematical regions is determined by the projected errors. Our goal, in contrast, is to force the two approaches to agree quite precisely on just a line in the kinematical plane, preferably in a region where both theoretical frameworks should be valid. The deviations as one moves away from this line provide a clearer picture of the genuine differences in the evolution dynamics, independently of experimental errors, kinematical ranges or PDF parametrizations.
The structure of the manuscript is as follows. In Section II A we review the method to compute the structure functions in collinear factorisation that will be employed to produce the DGLAP results. Section II B addresses the same point in the framework of non-linear, saturation dynamics. The PDF reweighting technique to produce a DGLAP fit coincident with the saturation results in a line of the x, Q 2 plane is presented in Section II C. Results for proton and nucleus are shown and discussed in Section III. Finally, we present our conclusions in Section IV.

A. Structure functions in collinear factorization
To compute the structure functions in collinear factorization we use APFEL [40] which accesses the PDFs through the LHAPDF library [41]. For the proton we use the next-to-leading order (NLO) NNPDF3.1 set NNPDF31 nlo as 0118 1000 [42] which has N rep = 1000 replicas describing the PDF uncertainty, and for gold (our representative of a heavier nucleus) we use the NLO nNNPDF2.0 nNNPDF20 nlo as 0118 Au197 [27] set which provides 1000 replicas as well 2 . The central value for any PDF-dependent quantity O is obtained as an average over the individual predictions O[f k ] for each individual replica, and the uncertainty δO is computed as the variance, We calculate the neutral-current DIS structure functions at NLO accuracy in the FONLL-B general-mass scheme [44]. This scheme matches with the one used in the actual global analyses. With the NNPDF3.1 proton PDF set we ensure the consistent treatment of fitted charm quarks by setting the appropriate options in APFEL [45], in particular enabling intrinsic charm by APFEL::EnableIntrinsicCharm(true). All the parameters (masses of heavy quarks, values of strong coupling α s , etc...) are set to be the same as in the corresponding PDF fits.
We note that the NNPDF collaboration has also prepared PDF sets including small-x resummation effects at leading and next-to-leading logarithmic accuracy [19]. In particular the (linear) small-x resummation was found to be necessary to describe the HERA structure function data in the small-x, moderate Q 2 region at next-to-nextto-leading order fits. At NLO accuracy, on the other hand, this resummation does not have a numerically significant effect on structure functions at HERA kinematics and is not considered in this work.

B. Structure functions with gluon saturation
In order to calculate proton and nuclear structure functions taking into account non-linear saturation effect we use the Color Glass Condensate (CGC) framework [1,10,46]. Here it is convenient to describe DIS processes in the dipole picture in a frame where the target is at rest and the photon plus light-cone momentum is large, and write the leading order total photon-nucleus cross section as (3) Here the photon light front wave function Ψ γ * →qq (r, z, Q 2 ) describes the photon splitting to a quark (with flavor f ) dipole with transverse size r and the quark carrying a fraction z of the photon plus momentum [1]. All information about the target is encoded in the dipole-target scattering amplitude N (r, b, x), which describes the eikonal propagation of the quark-antiquark pair through the target color field. The subscripts T and L refer to the photon polarization states. The structure functions can then be written as The energy (or Bjorken-x) evolution of the dipoletarget scattering amplitude N (and thus the x dependence of the structure functions) can be obtained by solving the perturbative Balitsky-Kovchegov (BK) evolution equation [11,12] which resums contributions ∼ α s ln 1/x to all orders. In contrast to the collinear factorization based approach described above, the Q 2 dependence is not calculated perturbatively -at least at leading order. Instead, it is effectively determined, along with other parameters, by the non-perturbative initial condition for the BK evolution. These initial conditions are obtained by fitting the HERA structure function data [17,18]. In this work we use the "MV e " leading order fit with running coupling corrections [47] from Ref. [48] (see also a similar fit in Ref. [49]).
One now needs to generalize the dipole-proton scattering amplitude to the dipole-nucleus case in order to calculate nuclear structure functions. We again do this using the procedure of Ref. [48], where the squared nuclear saturation scale is taken to be proportional to the transverse density . This scaling is used to generate the initial condition for the BK evolution separately at every impact parameter b, and different impact parameters are then evolved independently.
In principle it would be possible to include impact parameter dependence also in the BK evolution, but that would require us to include an additional phenomenological description of confinement effects, see e.g. Refs. [50][51][52][53]. As the structure functions are only sensitive to the integral over the impact parameter, we expect additional effects from the geometry evolution to only have a small effect on our results. In the region where the saturation scale of the nucleus would fall below that of the proton we use, again following Ref. [48], a dipole-proton scattering amplitude scaled such that the total dipole-nucleus cross section would be simply A times the dipole-proton cross section. By construction the nuclear modifications then vanish in the dilute region at the edge of the nucleus. This and similar setups have been successfully used to describe various LHC and RHIC measurements [54][55][56][57][58][59][60].
Equation (3) corresponds to the leading order contribution (with α s ln 1/x contributions resummed to all orders). Currently there is a rapid progress in the field towards next-to-leading order (NLO) accuracy, and necessary ingredients for NLO phenomenology are becoming available. These include, for example, the virtual photon wave function [61][62][63][64][65][66] and the BK evolution equation [67][68][69][70] at NLO. In the first phenomenological structure function calculations at NLO accuracy [62], very small differences were found between running coupling LO and NLO calculations in the HERA kinematics when the nonperturbative initial condition for the BK evolution is fitted to the data. This means that in leading order setups such as the one used in this work the fit parameters can effectively capture most of the higher order corrections. Consequently we expect our leading order framework to be sufficient for the purposes of this work, where the aim is to study differences between the DGLAP and BK dynamics.

C. Matching
As a technical method to match the structure functions in collinear factorization to those calculated with the BK evolution as discussed in Sec. II B, we use the Bayesian reweighting method [71,72]. The matching is achieved by reweighting the PDFs to reproduce structure functions calculated from the BK equation in a line in the x, Q 2 -plane. The preferred region for this matching is the region where both theoretical frameworks should be valid. Thus one wants to be sufficiently above the saturation scale Q 2 s to be in a linear regime described by DGLAP, but not so high that the large logarithms of Q 2 , not resummed by BK evolution, dominate. An additional technical requirement is that the matching points should be above the minimum scale Q 2 0 for which the PDF sets that we use are available. Based on these considerations we have chosen to match the BK and collinear factorization calculations in the region 10Q 2 is the saturation scale of the proton as defined in Ref. [48]. This leads to the matching interval x < 5.6 × 10 −3 for the proton and x < 1 × 10 −2 for the gold nucleus. We perform the matching above a lower limit x > 1.03 × 10 −5 that is chosen to be in the LHeC kinematical regime. In our matching region the saturation effects are expected to be weak, but Q 2 is still not too large. Consequently the DGLAP-and BKevolution based frameworks are expected to deviate from each other less dramatically than if the matching was performed in some more extreme region of phase space. Our results naturally depend on the choice of the matching scale; the results obtained with a matching at a higher value of Q 2 are discussed in Appendix A.
Matching the DGLAP and BK predictions using reweighting amounts to forming an appropriate linear combination of the available PDF replicas such that the predictions given by this linear combination reproduce the structure functions from the BK framework. In practice, for each PDF replica f k we define a figure-of-merit χ 2 k by where y i denote the N data matching values from BK and y i [f k ] are the corresponding values calculated by using a given PDF replica. The order of magnitude of χ 2 can be controlled by the entries of the covariance matrix σ ij which is taken to be of the diagonal form where δ BK is a constant. In a sense, we are thus assigning a constant relative uncertainty on the predictions ("pseudodata") generated from the BK setup.
For each replica we assign the so-called Giele-Keller weight [73], which always sum up to unity, The value for any quantity O using the reweighted PDFs is then obtained as a weighted average, with an uncertainty estimate We note that in the case of reweighting to real experimental data, the appropriate weights to be used in conjunction with the NNPDF replicas include also a term (χ 2 k ) (N data −1)/2 [71,74]. While the Giele-Keller weights favor replicas with χ 2 k /N data ≈ 0, the other weights favor replicas with χ 2 k /N data ≈ 1, as discussed in Ref. [74]. Here we are matching PDFs to smooth values resulting from a theory calculation, without point-by-point statistical fluctuations, not reweighting to experimental data. Thus the preferred values of χ 2 k /N data should be close to 0 and the Giele-Keller weights are better suited to our purpose.  It is useful to define the so-called effective number of replicas N eff , which serves as an proxy to the number of replicas with a significant weight [71,74], In the present analysis we use N eff to choose an appropriate value for the error parameter δ BK : If δ BK is too low, N eff ≈ 1 (for large N rep ) and the procedure picks up only a single replica irrespectively of how well it fits with the matching values. On the other hand, if N eff ∼ N rep all replicas have the equal weight and the reweighting does nothing, i.e., δ BK is too high. We have found that iteratively adjusting δ BK such that N eff ≈ 10 is a good strategy for finding a set of PDFs that matches the given boundary conditions, i.e., structure functions from the BK framework. In order to obtain N eff ≈ 10 we have fixed δ BK = 4.5 for the proton F 2 , δ BK = 11.5 for the proton F L , δ BK = 39.5 for the nuclear F 2 , and δ BK = 46 for the nuclear F L . For the nuclear reweighting we use N data = 138, and for the proton N data = 125.

A. Proton
The structure functions F 2 and F L for the proton before and after the reweighting on the Q 2 = 10Q 2 s (x) line are shown in Figs. 1a and 1b. The reweighting is done separately for F 2 and F L , as also in reality these two quantities will be measured in different kinematical domains and with a different experimental precision. The structure functions obtained after the reweighting can be seen to match very well to the BK results. This was to be expected since the proton PDFs and the initial condition for the BK evolution are fitted to the same precise HERA data at x 10 −4 , and the central NNPDF3.1 results are already very close to the BK values to begin with in this domain. However, a nearly perfect agreement with the BK results is obtained also at x 10 −4 . All in all, the matching procedure is thus found to work extremely well here.
Next we study how the differences in the BK vs. DGLAP dynamics become visible when we move away from the Q 2 ≈ 10Q 2 s (x) line. In Figs. 2a and 2b we show the relative difference as a function of both x and Q 2 , where F Rew 2,L refers to the corresponding structure function calculated using the reweighted PDFs. The points used in the reweighting are also indicated in these figures. One-dimensional projections of the same quantity are plotted at fixed values of x in Fig. 3.
For the F 2 structure function shown in Fig. 2a the differences remain very small, at most at a few-percent level almost everywhere in the studied x, Q 2 range, except in the high-x, high Q 2 and low-x, low Q 2 corners. This is better visible in Fig. 3a where we show the relative differences as a function of virtuality Q 2 at four different x values from x = 5.6 × 10 −3 (largest x for which is the initial scale in the NNPDF3.1 PDF set) to x = 10 −5 . The smallest x values in our plots are beyond reach for the EIC, which will collide electrons with energies 5 − 18 GeV on protons and nuclei with energies 250 and 100 GeV/nucleon respectively, resulting in a kinematic reach (at Q 2 = 10 GeV 2 ) down to x ∼ 10 −3 [33]. Smaller x values could be probed at the LHeC (50 GeV electrons on Z/A × 7 TeV/nucleon protons and nuclei) whose kinematic reach goes down to x ∼ 10 −5 [35] and at the FCC-he [14] (60 GeV electrons on Z/A×50 TeV/nucleon protons and nuclei) whose kinematic coverage extends to even lower x. We see that around x ∼ 10 −4 the Q 2 dependencies are nearly equal in both frameworks. In the higher-x region the BK equation predicts a stronger Q 2 dependence than the DGLAP equation, while in the x 10 −4 region the BK dynamics results with a weaker Q 2 dependence than what the DGLAP equation predicts. As a result, at fixed Q 2 ∼ 10 GeV 2 the relative difference changes sign as a function of x. Since the relative differences remain at a few-percent level, a very precise determination of the proton F 2 is required in order to distinguish between the two physical pictures in a statistically meaningful manner.
The differences between the BK and DGLAP dynamics are more clearly visible in the case of the structure function F L . This can be seen from Fig. 2b and Fig. 3b which show the analogous plots for F L that were above discussed for F 2 . There are now larger differences even within the HERA kinematics as the F L data from HERA are rather scarce. The DGLAP evolved F L shows gener-ically a significantly stronger Q 2 dependence in comparison to the BK evolved F L . The x dependencies at fixed Q 2 are now also clearly different, particularly at Q 2 10 GeV 2 where the BK evolution predicts a stronger x dependence than what the reweighted F L with DGLAP dynamics has. In the EIC kinematics the relative differences are maximally ∼ 10%. In the LHeC kinematics differences can be up to ∼ 40% which would be easier to resolve.
Note that for the proton F 2 in Fig. 3a there is a small b quark threshold effect in the FONLL-B scheme at Q 2 = m 2 b ≈ 24 GeV 2 .

B. Heavy Nucleus
Next we study the structure functions F 2 and F L for the 197 Au nucleus. The procedure here is the same as for protons above: first we match the nuclear PDF predictions with the BK predictions in the region 10Q 2 where Q s is the same proton saturation scale as what was used above. As the saturation scale of a large nucleus is generically just a few times larger than that of the proton (see e.g. Ref. [48]), at Q 2 = 10Q 2 s (x) the non-linear effects should still be relatively small and both frameworks are expected to be applicable.
The nuclear structure functions F 2,L at Q 2 = 10Q 2 s (x) before and after reweighting are shown in Figs. 4a and 4b. As one can see, this time the reweighting is not as successful as it was in the case of proton. The proton and nuclear PDF sets have exactly the same amount of Monte Carlo replicas but, as there are much fewer experimental constraints for the heavy-nucleus structure functions, there are larger variations between the different replicas (as illustrated by the large uncertainty band before reweighting). Consequently, the remaining uncertainty band calculated with N eff = 10 effective replicas can be much larger than with protons, and a more precise matching would require an even larger number of replicas. Thus the larger error band after reweighting is mostly just a technical limitation, not a physics effect in itself. A partial explanation for the difference between the central values before and after reweighting is that the BK prediction is not constrained by experimental data on hard processes involving heavy nuclei, which do affect the nuclear-PDF predictions.
Next we again study the relative differences between the DGLAP-and BK-evolved structure functions, defined by Eq. (13), on the x, Q 2 plane. The results are shown in Fig. 5a for F 2 and in Fig. 5b for F L . In comparison to the proton results shown in Figs. 2a and 2b we find much larger differences with a heavy nucleus. This is expected, as non-linear dynamics not present in the DGLAP evolution is enhanced in heavy nuclei by a factor ∼ A 1/3 . The fact that the results obtained with the matched PDFs do not exactly agree with the BK predictions results leads to the line of agreement between the two approaches (white in the figure) not exactly aligning with the points where the matching is done. However, the differences between the two evolution equations when moving away from the Q 2 = 10Q 2 s line are still very clear. Again, we see more dramatic differences in F L compared to F 2 . For F 2 the cleanest systematical effect is that the BK evolution predicts slower x evolution at all Q 2 , and this difference in the evolution speed is much larger than with proton targets. In the case of F L the DGLAP evolution results in a significantly faster Q 2 evolution except at the highest x ∼ 10 −2 , and the x dependence from BK evolution is again generically slower. The fact that BK evolution generically results in a slower x dependence with heavy nuclei is expected, as the effect of the non-linear contributions in the BK evolution is to slow down the evolution speed.
To quantify more accurately the expected magnitude of the deviations, we show in Figs. 6a and 6b the relative difference from Figs. 5a and 5b for fixed x values in the EIC and LHeC/FCC-he kinematics as a function of Q 2 . The reweighted F 2 values differ from the BK-evolved results less than 10 % in the studied kinematical domain, which implies that also with heavy nuclei the F 2 should be measured at a few-percent precision in order to be sensitive to the differences in the evolution. For F L the differences are again significantly larger and up to 60% in the LHeC/FCC-he kinematics and up to 15% in the EIC kinematics. Similarly to protons, a bottom quark mass threshold effect in the FONLL-B scheme is seen in the nuclear F 2 in Fig. 6a at Q 2 ≈ 24 GeV 2 .

C. Effects on PDFs
The matching process results in a new PDF set that is constrained by the BK evolved structure functions on the Q 2 = 10 . . . 11Q 2 s (x) line. In this section we study how these constraints originating from the non-linear dynamics included in the BK evolution alter the parton distribution functions. The coefficients ω k defined in Eq. (8) determined in the reweighting process for the applied PDF sets are made available in the Supplementary material. Using these coefficients and the publicly available NNPDF3.1 or nNNPDF2.0 sets it is then possible to construct DGLAP evolved parton distribution functions that match the saturation model predictions at moderate Q 2 and are compatible with the neutral-current HERA data.
First we illustrate here the effect of reweighting on the proton parton distribution functions. The results are shown in Figs. 7a and 7b for the gluon and up-quark distributions at two different virtualities. In order to study how complementary the observables F 2 and F L are, we separately show the results with matching to F 2 and matching to F L . The central set values from the NPDF3.1 are also shown for comparison. Comparing Fig. 7a to 7b we find the reweighting to have a slightly stronger effect on the gluon distribution compared to up quarks, but in general the reweighting has only a modest effect on the proton PDFs. This is expected, as the  PDFs are fitted to the same HERA data that is used to constrain the BK boundary conditions. Whether F 2 or F L is used in reweighting has only a small effect on the determined reweighted PDFs. Thus, we do not expect to see strong tensions when measurements from the EIC or LHeC/FCC-he are eventually used to disentangle the BK and DGLAP dynamics.
The reweighted nuclear up-quark and gluon distributions are shown in Figs. 8a and 8b. Comparing to the proton results shown in Figs. 7a and 7b we see that nuclear PDFs are affected much more by the reweighting already in the x 10 −3 region, which is expected, as in nNNPDF2.0 there are only few data constraints in this region. The reweighted nuclear PDFs are suppressed by a large factor compared to the central values from the nNNPDF2.0 set. Again both F 2 and F L pseudodata have similar effects and as such no strong tensions with already existing data included in the nuclear PDF fits are expected in global analyses. In Fig. 8a the nuclear gluon distribution, reweighted with F 2 data, becomes negative at small x 2 · 10 −5 and at Q 2 = 3.1 GeV 2 . However, the gluon distribution is not an observable, and structure functions remain positive.

IV. CONCLUSIONS
Understanding the importance of nonlinear effects is a major motivation for future high-energy scattering experiments. Experimental measurements have, however, a limited precision and kinematical range. It is therefore essential to quantify in a clear way the magnitude of the effects that one is looking for. In this paper we have done this by comparing the DIS structure functions F 2 and F L calculated using the running coupling BK equation, which includes saturation effects, to collinear factorization at NLO. The approaches require a set of non-perturbative initial conditions in some initial line in the (x, Q 2 )-plane and, given these initial conditions, make predictions in other regions. In order to compare the two approaches to evolution it is essential to differentiate between the effects of these initial conditions, and the actual dynamical predictions of the evolution. Our focus here is on the evolution dynamics. Thus we perform this comparison by choosing the nonperturbative initial conditions so that they match in a regime where both frameworks should be valid, at virtualities somewhat above the saturation scale. Technically this is achieved by profiting from large samples of DGLAP-evolved Monte-Carlo  replicas of PDFs available from the NNPDF collaboration, and using reweighting to match them to the structure functions from BK evolution. There is already a significant agreement due to the fact that both calculations have been (for protons) fit to the same HERA data. The deviations of the two calculations, when moving away from the matching line, provide a clear assessment of the quantitative difference in the dynamics between the two.
As a result of this analysis, we see that F 2 has a quite weak sensitivity to the difference between the evolution equations, at least in the region where Q 2 is large enough for the PDF sets to be available. The relative differences in F 2 are at a percent level for protons, and somewhat higher for nuclei. The longitudinal structure function F L is significantly more sensitive, but unfortunately also harder to measure as precisely. Overall the difference, as expected, is strongest at the lowest Q 2 values. Pushing collinear factorization to lower scales would clearly make the difference more manifest. In finer details, the behavior of the difference in the x, Q 2 -plane also has more unexpected aspects.
We have here focused on providing a comparison of the predictions on a purely theoretical level, independently of any assumptions about a specific experiment. To turn this calculation into a more controlled estimate relevant to a specific measurement requires, in addition, a more detailed inclusion of estimated experimental errors, their correlations, the expected kinematical coverage and so on. This we leave to future work. Another interesting avenue for the future would be use this as a starting point for using the two approaches together. It should be possible to use BK evolution to predict the behavior at relatively low Q 2 values and very low x and combine this with DGLAP evolution towards higher Q 2 .
Such a combined framework would significantly extend the predictive power of weak coupling QCD by replacing a part of what are now initial conditions fitted to data by predictions from a first-principles calculation. It would be interesting to perform the matching required in such a calculation analytically. This is, however, beyond the scope of this work. Instead we take a first step in this direction by providing, as a supplementary material to this paper, the DGLAP-evolved PDF sets for protons and gold nuclei that do match the BK-evolved cross section at our chosen matching scale. s (x) for the gold nucleus. These figures should be contrasted with Figs. 4, 5 and 6, where the same quantities are plotted with our default matching scale. One can see that an equally good matching at this higher scale can be achieved. In this case there is much more phase space available below the matching scale, and at the lower Q 2 values the disagreement between the two calculations becomes quite substantial. This indicates that a factor ∼ 2.5 in Q 2 is already sufficient to reveal a substantial difference between what is in some sense a backward DGLAP evolution towards lower scales, and the BK dynamics. This is, however, not relevant for an actual determination based on experimental data. First, experimentally there is more phase space and thus more data available at the lowest scales. Second, and more important, allowing both calculations to adjust nonperturbative initial conditions to the same data will tend to match them at a scale which makes them differ less, which in this case is at smaller values of Q 2 .