Inferring the initial condition for the Balitsky-Kovchegov equation

We apply Bayesian inference to determine the posterior likelihood distribution for the parameters describing the initial condition of the small-x Balitsky-Kovchegov evolution equation at leading logarithmic accuracy. The HERA structure function data is found to constrain most of the model parameters well. In particular, we find that the HERA data prefers an anomalous dimension γ ≈ 1 unlike in previous fits where γ > 1 which led to e.g. the unintegrated gluon distribution and the quark-target cross sections not being positive definite. The determined posterior distribution can be used to propagate the uncertainties in the non-perturbative initial condition when calculating any other observable in the Color Glass Condensate framework. We demonstrate this explicitly for the inclusive quark production cross section in proton-proton collisions and by calculating predictions for the nuclear modification factor for the F 2 structure function in the EIC and LHeC/FCC-he kinematics.


I. INTRODUCTION
Understanding the high-energy structure of protons and nuclei is one of the central goals of the nextgeneration Electron-Ion Collider (EIC) [1].Especially the fact that one can perform nuclear-DIS experiments for the first time in collider kinematics is intriguing: at high center-of-mass energies it is possible to access the high-density part of the nuclear wave function where nonlinear QCD dynamics with emergent saturation phenomena is expected to play a major role.
Precise theory calculations taking into account saturation effects are necessary to probe in detail the non-linear QCD dynamics in hadronic collisions at RHIC and at the LHC, as well as in accurate nuclear-DIS experiments [2] at the EIC or at the proposed future LHeC/FCC-he collider [3].The Color Glass Condensate (CGC) effective theory of QCD [4] provides a natural framework to describe gluon saturation at high energies.In recent years there has been a rapid progress towards next-to-leading order accuracy in CGC calculations, see e.g.Refs.[5][6][7][8][9][10][11][12][13][14][15].In addition to higher order corrections, precise theoretical predictions require a well constrained nonperturbative input describing the proton or nuclear structure at moderately small-x, which is an initial condition to the perturbative evolution equations such as the Balitsky-Kovchegov (BK) [16,17] equation resumming contributions enhanced by large logarithms of center-ofmass energy.
The initial condition for the BK evolution has been fitted to the HERA structure function data [18,19] at leading [20][21][22][23] and recently also at next-to-leading [24,25] order accuracy (see also Ref. [26] for a complementary approach to determine the initial condition from the large-x structure of the proton).This initial condition is a necessary input when calculating any other scattering process at small-x, such as particle production in proton-nucleus collision at the LHC, see e.g.Refs.[8,21,27,28].
All current fits to the DIS data extracting the initial condition for the BK evolution lack an uncertainty analysis, i.e. do not provide any method to propagate the uncertainties in the initial condition parametrization to the calculated observables.When developing the CGC calculations to the precision level, a statistically rigorous treatment of these uncertainties can be seen to be of equal importance as the higher order corrections when comparing predictions to current and future experimental data.
In this work we extract, for the first time with an uncertainty estimate, the non-perturbative initial condition for the BK evolution from the HERA structure function data [19].This is achieved by employing a Bayesian inference setup extensively used in the field, recently when extracting e.g.properties of the Quark Gluon Plasma [29][30][31][32][33][34] and the event-by-event fluctuating proton geometry [35].The determined posterior distribution for the initial condition parameters allows for a rigorous propagation of initial condition parametrization uncertainties to calculations of all other observables in the CGC framework.In this work we present results based on a leading order analysis but note that the developed computationally efficient setup can be extended to the numerically heavy next-to-leading order case, where it has been recently shown that combined analyses including both total cross section and charm production data are feasible [24].
This paper is structured as follows.In Sec.II we review the calculation of the proton structure functions at high energy in the CGC approach.The Bayesian inference setup is described in Sec.III.The determined initial condition for the BK evolution with an uncertainty estimate is presented in Sec.IV.To illustrate the propagation of uncertainties to calculations of various observables, we calculate in Sec.V inclusive quark production in proton-proton collisions and the nuclear modification factor for the structure function F 2 which will be measured at the EIC but for which no predictions that include the BK evolution consistently with the HERA data exist in the literature (see however Ref. [36]).Conclusions are presented in Sec.VI, and the importance of including correlations among the experimental uncertainties is analyzed in Appendix A.

II. DEEP INELASTIC SCATTERING AT HIGH ENERGY
The total cross section in lepton-nucleus scattering is typically expressed in terms of the reduced cross section defined as Here y = Q 2 /(sx) is the inelasticity, x Bj is the Bjorkenx, Q 2 the photon virtuality and √ s the electron-proton center-of-mass energy.The structure functions F 2 and F L are related to the total virtual photon-nucleus cross section as Here the subscripts T and L refer to the transverse and longitudinal virtual photon polarization states, respectively.
In the dipole picture the total cross section for the virtual photon-nucleus scattering at high energy factorizes to a product of the virtual photon wave function ψ γ * →q q describing γ * → q q fluctuation, and the dipole-target scattering amplitude N as [37] Here z is the fraction of the photon longitudinal momentum carried by the quark, r is the transverse size of the q q dipole (higher Fock states would enter beyond leading order) and b is the center-of-mass of the dipole.We sum over light quarks f = u, d, s, as it has been shown in Ref. [22] that in a leading order calculation with the energy dependence obtained from the BK evolution it is not possible to simultaneously describe the total and charm production cross section data.Following Refs.[20][21][22], when calculating the virtual photon-proton cross section we neglect the impact parameter dependence in the dipole-proton amplitude N and replace Here the proton transverse area σ 0 /2 is taken to be a free parameter.
The squared photon wave functions summed over quark helicities read [37] |ψ and with Here m q = 0.14 GeV is the light quark mass and e q is the quark charge.One could also in principle take the light quark mass to be a fit parameter as e.g. in Refs.[22,38].However these previous analyses have found only a small sensitivity on m q and as such we choose a fixed value that is compatible with previous fits.
The Bjorken-x dependence of the dipole amplitude is given by the Balitsky-Kovchegov equation: Here x ij = x i − x j and Y = ln x 0 /x (see also Ref. [39] for a detailed discussion of the evolution variable).In this work we initiate the BK evolution at x 0 = 0.01.When the running coupling corrections are included following the Balitsky prescription [40] the kernel reads In this work we present a leading-order analysis.We note that in addition to next-to-leading order corrections [6,7,13,14], it is also possible to resum corrections enhanced by large transverse logarithms to the BK equation following Refs.[41][42][43] (that in practice approximate the full NLO BK equation accurately [44,45]).Furthermore, it would be worthwhile to use the recently developed resummation procedure from Ref. [39] where the evolution is formulated in terms of the target rapidity instead of the projectile rapidity.These improvements we leave for future work.
For the initial condition of the BK evolution we use a McLerran-Venugopalan model [46] inspired parametriza-tion as e.g. in Ref. [21] N Here the free parameters are Q 2 s,0 related to the initial saturation scale, the anomalous dimension γ controlling the shape of the dipole amplitude at small |x ij | and the infrared regulator e c .The Bayesian inference setup discussed in Sec.III can efficiently handle multidimensional parameter space, which allows us to use this more flexible parametrization for the initial condition compared to functional forms used in previous fits (e.g. in Ref. [21] either γ or e c was fixed).
The strong coupling constant in Eq. ( 9) depends on the transverse distance scale.Following again Refs.[20,21] we write the coordinate space strong coupling as Here we take n f = 3 as we only consider light quarks in this work, and Λ QCD = 0.241 GeV.In the infrared region the coupling is frozen to α s = 0.7.The parameter C 2 controls the strong coupling scale in the coordinate space, and based on Refs.[40,47] its expected value is C 2 = e −2γ E ≈ 0.32.However we consider C 2 to be a free parameter to control the uncertainty of the strong coupling scale and to absorb some missing higher order corrections, for example the fact that the next-to-leading order evolution is known to result in slower evolution speed compared to the leading order case [45].Previous leading order fits have also been found to prefer a much larger C 2 = 7 . . .15 [21].

III. BAYESIAN ANALYSIS SETUP
The non-perturbative parameters describing the initial condition for the BK evolution, the proton transverse area and the scale of the coordinate space running coupling can be estimated using Bayesian parameter inference.In this work we use a similar inference procedure as e.g. in Refs.[29,33,35].The output is a multidimensional likelihood distribution that comprehensively describes the parameter uncertainties and correlations based on prior beliefs and experimental data.
Bayesian inference is based on Bayes' theorem: Here the model output vector with given parametriza- The posterior function, represented by P (y(θ)|y exp ), is the probability of y(θ) being true given the experimental data, y exp , as evidence.
The probability represented by P (y exp |y(θ)) is a measure of the agreement between the observation and the model calculation at the kinematical point x given the model parameters θ.We use a multivariate normal distribution in which case the logarithm of the likelihood function reads log P (y exp |y(θ)) Here Σ = Σ model (θ) + Σ exp is a matrix sum of the model, or emulator as discussed shortly, and the experimental covariance matrices and the difference ∆y(θ) = y(θ) − y exp is a vector with length N .The experimental data considered in this work is the reduced cross section data measured by the H1 and ZEUS collaborations [19], where the published dataset also includes the full covariance matrix describing correlations among the 162 sources of systematic uncertainty.Furthermore we include as correlated uncertainties the 7 different procedural uncertainties originating from the combination of the H1 and ZEUS datasets.These correlations between different systematic uncertainties have not been included in previous dipole model fits, and the effect of these correlations on the final posterior distribution is illustrated in Appendix A. The uncorrelated systematic uncertainties and statistical uncertainties are added in quadrature.The prior, represented by P (θ) in Eq. ( 12), encodes initial knowledge on the range and shape of the probability distribution of the model parameters.Not much information is known about the true distribution of the model parameters other than estimates of previous fits from Refs.[20][21][22].In this study, the prior is chosen to be a flat multivariate distribution with bounds corresponding to guesses of possible parameter ranges, outside of which the posterior is set to zero.Initially, the width of the posterior distribution is unknown, so the study conservatively begins with a wide allowable extent and is later narrowed down.The prior ranges used in the final analysis are shown in Table I.
In order to effectively sample the multidimensional parameter space, we employ Gaussian Process Emulators (GPEs) [48] as computationally efficient surrogates for the actual theory calculation, and the Markov Chain Monte Carlo (MCMC) sampling method.The MCMC importance sampling is implemented via the emcee package [49] that uses an affine-invariant ensemble sampler.This sampler uses an algorithm similar to the Metropolis-Hastings (MH) where each step of a walker from A → B is accepted with a probability of q ∼ P (θ B ) P (θ A ) .Unlike the MH algorithm, emcee employs an ensemble of k walkers where the advancement of the chain of one walker is also based on the current positions of the other walkers; therefore, it converges much faster than MH.Each step is then appended to the MCMC chain to, eventually, form the target posterior distribution.
The GPEs are trained at different points in the parameter space sampled from a latin hypercube which ensures that the design points are optimally distributed in the parameter space.This significantly reduces the number of training points needed to get accurate estimates from GPEs, and in this analysis we use 500 training parametrizations.We include data in the region 2.0 GeV 2 < Q 2 < 50.0 GeV 2 where there are 403 datapoints.The low-Q 2 cut ensures the presence of a perturbative scale, and the upper limit removes contribution from the region dominated by DGLAP dynamics not included in the present analysis.We represent the HERA data in this kinematical domain by 6 principal components that capture 99.99999% variance of the model output within the prior.The GP output at a certain (x, θ) is a mean estimate and a covariance matrix that describes the emulator uncertainties and is used to form the matrix Σ model .
The GPEs are validated by comparing their predictions with the actual model calculations of the reduced cross section.The left panel of Fig. 1 show the comparison between the model and the emulator output at all kinematical (x Bj , Q 2 , y) points in the HERA data computed using separate validation sets of model parameters not included in the training of the GPEs.The emulator accuracy is found to very good, with the average relative difference being 0.047% in our standard setup where Q 2 s,0 , e c , C 2 , σ 0 /2 and γ are free parameters.The righthand side panel of Fig. 1 shows the corresponding z-score defined as where σ GPE is the uncertainty estimate of the GPE.The width of the z distribution is less than unity which means that the emulator uncertainties are typically slightly overestimated.However, it is the experimental uncertainty that dominates in the likelihood function.

IV. INFERRED INITIAL CONDITION FOR THE BK EVOLUTION
We determine the posterior distribution of model parameters in two separate setups.
First, in our main setup we consider all parameters discussed above (Q 2 s,0 , e c , C 2 , σ 0 /2 and γ) to be free and determine the corresponding posterior distribution.We will refer to this setup as the 5-parameter case.For comparison we also consider the case where we set γ = 1 as in the standard MV model, and refer to this as the 4-parameter setup.This parametrization corresponds to the "MV e " one used in Ref. [21].
The obtained posterior distributions are shown in Fig. 2 for both the 4-and 5-parameter cases.The plots along the diagonal present the 1D projections of the posterior showing the likelihood distributions for individual parameters.The off-diagonal plots are 2D histograms that illustrate the correlations between the pairs of model parameters.
The HERA data is found to constrain most of the model parameters well, except the infrared regulator e c for which the likelihood distribution extends to large values of e c in both setups.The broad distribution is not unexpected as the dipole amplitude depends on e c weakly and only in the region where r ≳ 1/Λ QCD , see Eq. (10).Also, in the 5-parameter case, large values for Q 2 s,0 are allowed and compensated by small e c , and in both setups we obtain a clear negative correlation between Q 2 s,0 and e c .
In the 5-parameter case where γ is a free parameter, we find that the HERA data prefers γ = 0.95 . . .1.05.This is in contrast to previous fits performed in Refs.[20][21][22] where the Q 2 -dependence of the HERA data was shown to prefer γ ≈ 1.1 . . .1.2 (in fits with fixed e c = 1).As we will discuss in more detail in Sec.V, having γ ≲ 1 is advantageous e.g. when calculating inclusive particle production cross section in proton-nucleus collisions.The most likely values for the individual parameters in the 4-parameter case are similar to what is reported in the previous leading order fit [21].
In addition to the negative Q 2 s,0 -e c correlation discussed above, we also find a clear negative correlation between Q 2 s,0 and σ 0 /2 especially in the 5-parameter case.This is because, in the region where the virtual photonproton cross section is dominated by small dipoles, one has where we use the fact that r 2 ∼ 1/Q 2 (except in the aligned jet limit where z ≈ 0 or z ≈ 1).This dependence also explains the negative correlation seen between γ and σ 0 /2.As stated earlier, Q s,0 parametrizes the saturation scale at initial x Bj .The saturation scale, Q s , determines the scale at which non-linear dynamics becomes impor-0 .0 5 0 0 .0 7 5 0 . 1 0 0 .
tant.Independent of the functional form of the BK initial condition, we can define (following e.g.[21]) Q 2 s as the solution to N (r 2 = 2/Q 2 s , x Bj ) = 1 − e −1/2 .The determined saturation scales at initial x Bj in both 4-and 5-parameter setups are shown in Table .I, and are found to be identical in both cases.Furthermore, Fig. 4 shows the x Bj dependence of Q 2 s , which is identical in both setups despite the differences in the posterior distributions.
At large Q 2 s,0 , the typical running coupling α s (r) ∼ α s (1/Q s ) is smaller.Consequently in that case there is  no need for as large C 2 to get a slow enough evolution speed compatible with the x-dependence of the HERA data.Hence the negative correlation between C 2 and Q 2 s,0 is obtained, and consequently a positive correlation between C 2 and σ 0 /2.These correlations appear very clearly in both setups.Consequently, it is essential to take into account the correlations between the model parameters when estimating uncertainties when calculating different cross sections.
The median values for the parameters are presented in Table I along with the maximum a posteriori (MAP) estimates corresponding to the maximum of the posterior distribution.
The compatibility with the HERA data is quantified by χ 2 /d.o.f defined as FIG. 4: The saturation scale as a function of x Bj , defined as the solution to N (r) = 1−e −1/2 when r 2 = 2/Q 2 s .This is shown for both the 4-and 5-parameter cases with 2σ uncertainty bands.
where N = 403 is the number of experimental datapoints in the considered kinematical domain and p = 4 or p = 5 is the number of free model parameters.The χ 2 /d.o.f ≈ 1.01 shows that an excellent description of the precise HERA data is obtained, similarly as in previous leading order fits.By calculating the average cross section over many samples from the posterior and not using a single parametrization (e.g.median or MAP), we get χ 2 /d.o.f = 1.02 in both the 4-and 5-parameter cases, showing that both the MAP and median parametrizations are good estimates when calculating the proton structure functions.
The good agreement with the HERA data is also illustrated in Fig. 3 where we show a comparison to the reduced cross section data in a few selected virtuality bins.The central lines are the average values obtained by computing the cross section using many posterior samples, and the uncertainty band corresponds to two stan- x B j = 1 0 FIG. 5: The dipole amplitude at the initial condition and evolved down to x Bj = 10 −5 calculated in the 4-and 5parameter cases.The band shows a 2σ uncertainty.
dard deviation variation.To allow for asymmetric uncertainty estimates, in this work we always calculate separately the 2 standard deviation (2σ) uncertainty band above and below the mean value.We note that the reduced cross section is generically underestimated at low Q 2 .However, we still get a good χ 2 when the correlated systematic uncertainties in the HERA data are taken into account.If the statistical and systematical uncertainties in the HERA data are simply added in quadrature, we get a somewhat larger χ 2 /d.o.f = 2.2 in both the 4-and 5-parameter case.See also discussion in Appendix A for the effect of the correlated uncertainties on the posterior distribution.
The initial condition for the dipole amplitude N (r, x Bj = 0.01) and the evolved amplitude at a much smaller x Bj = 10 −5 are shown in Fig. 5.Note that the initial dipole amplitude only depends on Q s,0 , e c and γ, and furthermore the evolution speed is also sensitive to C 2 .The central value of the dipole amplitude is obtained as an average of the dipole amplitudes computed using different samples from the posterior distribution, and the uncertainty estimate is a 2 standard deviation band.We show the dipole amplitudes both from the 4-and 5-parameter setups, and, as expected, the uncertainty band in the 5-parameter case is slightly larger due to a broader posterior distribution.However the actual (2 standard deviation) uncertainty in the dipole amplitude is small, typically a few percent and maximally ≈ 10%.

V. APPLICATIONS
As the dipole-target scattering amplitude is a convenient degree of freedom at small-x (see e.g.Ref. [50]), the determined posterior distribution is a necessary input to calculations of all observables.The fact that the full posterior distribution, and not just the best fit values, are available, also enables us to rigorously take into account the uncertainties in the non-perturbative parameters describing the initial condition for the BK evolution.In this Section, we demonstrate how the uncertainties in the extracted dipole amplitude propagate to the inclusive quark-target cross section and to the nuclear modification factor for the F 2 structure function that will be measured at the EIC.

A. Inclusive quark production
Let us first consider quark-target scattering which is the relevant subprocess in inclusive forward hadron production in pp or pA collisions.The inclusive forward quark production cross section in proton-proton collisions is directly proportional to the two-dimensional Fourier transform of the initial dipole-amplitude [21,27,51,52]: where the produced quark transverse momentum and rapidity are k and y, respectively.The dipole amplitude is evaluated at x = |k|/ √ se −y where √ s is the center-ofmass energy.Furthermore, xq(x, µ 2 ) is the quark parton distribution function evaluated at scale µ 2 , and The two-dimensional Fourier transform of the dipole amplitude scaled by σ 0 /2 at the initial x Bj = 0.01 is calculated in both the 4-and 5-parameter case, and the results are shown in Fig. 6.Although the median γ obtained is close to 1, the posterior covers parametrizations where anomalous dimensions slightly larger than unity (up to γ ≈ 1.05) are encountered, resulting in negative values of Sp (k) in the large k region.This is because the Fourier transform is not positive definite with γ > 1 [53].At smaller x after the BK evolution the Fourier transform would be positive definite as the BK evolution drives the anomalous dimension towards an asymptotic value γ = 0.6 . . .0.8.
Additionally we observe that although the uncertainties in the determined dipole amplitude are typically around 5% (see Fig. 5), uncertainties in the quark-target cross section (in the Sp ×σ 0 /2) can be significantly larger.Even in the 4-parameter case where the Fourier transform is positive definite, e.g. at k = 2 GeV the upper limit of the 2σ uncertainty band is around 150% larger than the mean value.This highlights the importance of properly taking into account the parametrization uncertainties when applying the dipole amplitude determined from the HERA data when describing particle production processes at the LHC as e.g. in Ref. [21].At next-to-leading order [9,10] the possibility to emit a gluon changes the kinematics from 1 → 1 at LO to 1 → 2 at NLO, which can have a significant effect on the p T distribution (e.g. one can obtain a power-law tail at high p T even with a Gaussian dipole).The importance of properly take into account uncertainties in the initial dipole-target scattering amplitude when calculating inclusive spectra at NLO has been recently emphasize in Ref. [28].

B. Nuclear effects on F2
Much more pronounced saturation effects can be expected to seen in scattering processes where the target is a heavy nuclei instead of the proton, as the saturation scale in the nucleus scales roughly as ∼ A 1/3 .The dipoleproton scattering amplitude determined in this work can be generalized to the dipole-nucleus case e.g. by employing the optical Glauber model.Following Ref. [21], we write the initial condition for the BK evolution of the dipole-nucleus scattering amplitude at fixed impact parameter b as where T A (b) is the transverse thickness function of the nucleus of mass number A. The dipole-nucleus amplitudes at fixed b are then evolved to small-x by solving the BK equation without the impact parameter dependence.This thickness function is obtained by integrating the Woods-Saxon distribution 2 N (r, x), i.e. use a dipole-proton scattering amplitude scaled such that all nuclear effects vanish.
We compute the nuclear modification factor for the structure function F 2 defined as where F 2,A is the F 2 structure function for the nucleus with mass number A. In this work we consider gold, i.e.A = 197, for which this modification factor will be measured at the EIC.The obtained R eA as a function of x Bj is shown in Fig. 7 and as a function of Q 2 in Fig. 8.We predict a significant nuclear suppression in the F 2 structure function already in the EIC kinematics.By construction R eA → 1 at the initial x = 0.01 in the limit Q 2 → ∞, and a stronger suppression is seen both towards small-x Bj and low-Q 2 .Unlike above when studying the inclusive quark production cross section, in the case of R eA the parametrization uncertainties effectively cancel in the structure function ratio and the 2σ uncertainty is at 1-2% level.The uncertainty grows slightly towards small-x Bj where also the C 2 parameter controlling the evolution speed become important.

VI. CONCLUSIONS AND OUTLOOK
We have determined the posterior distribution for the non-perturbative parameters describing the dipole- proton scattering amplitude at the initial condition of the small-x BK evolution.This is a necessary input to calculations of all scattering processes at high energy in the Color Glass Condensate framework.We have also determined, for the first time, uncertainty estimates for these non-perturbative parameters.The obtained posterior distribution makes it possible to rigorously take into account the uncertainties in the initial condition parametrization when looking for signals of gluon saturation from various different observables.The posterior distributions determined in this work can be found from the supplementary material.
We have shown that the HERA structure function data constrains most of the model parameters well and the remaining uncertainties in the obtained dipole amplitude are at a few percent level.The flexible Bayesian inference setup allows us to include more freedom in the initial dipole-proton amplitude than in the previous fits, and in particular we find that the HERA data does not require a large anomalous dimension γ ≈ 1.1 . . .1.2 unlike in the analyses presented in Refs.[20][21][22].Instead we get γ ≈ 1 which is more natural considering that γ ≤ 1 is required to obtain a positive definite quark production cross section and unintegrated gluon distribution [53].
We demonstrate the uncertainty propagation by calculating the inclusive quark production cross section in proton-proton collisions and the nuclear modification factor R eA for the structure function F 2 to be measured at the EIC and at the LHeC/FCC-he.Significant uncertainties up to ∼ 100% were found in the case of the inclusive quark production in proton-proton collisions.On the other hand, the uncertainties were found to mostly cancel in the nuclear modification factor R eA for the F 2 structure function.As such it can be crucial to properly propagate the uncertainties in the non-perturbative input when comparing the CGC predictions describing the gluon saturation effects to e.g.LHC inclusive spectra.
In the future, we will apply the the developed computationally efficient setup to determine the initial condition for the BK evolution at next-to-leading order accuracy.This is especially intriguing as it has been recently shown that at NLO global analyses including both total and heavy quark production cross section data simultaneously become feasible [24,25].We expect the proper treatment of initial condition uncertainties to be of equal importance as higher-order corrections when the Color Glass Condensate calculations are promoted to the precision level.

FIG. 1 :
FIG. 1: Validation for the 5-parameter Bayesian fit (left) comparing the calculated reduced cross section to the emulator prediction.The z-score plot (right) is constructed with 500 validation points over 403 kinematics.The dashed line shows a Gaussian fit to the z-score histogram, for which the mean and the standard deviation (SD) are given in the figure.

FIG. 2 : 1 -
FIG.2:1-and 2-dimensional projections of the posterior probability distributions for both the 4-parameter (red) and 5-parameter case (blue).The dotted green line shows the best fit values from the MV e fit of Ref.[21]

FIG. 7 :
FIG. 7: Nuclear modification factor for the structure function F 2 as a function of x Bj at fixed Q 2 = 10 GeV 2 .

FIG. 9 :
FIG. 9: Posterior distribution in the 5-parameter case obtained with (blue) and without (red) including the correlations among the experimental uncertainties.

TABLE I :
MAP and median values for the 4-and 5-parameter setups.The uncertainty estimates in the 95% credible intervals are also shown.
FIG.8: Nuclear modification factor for the structure function F 2 as a function of photon virtuality Q 2 at fixed x Bj = 0.001.