Proton PDFs with non-linear corrections from gluon recombination

We present numerical studies of the leading non-linear corrections to the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi evolution equations of parton distribution functions (PDFs) resulting from gluon recombination. The effect of these corrections is to reduce the pace of evolution at small momentum fractions $x$, while slightly increasing it at intermediate $x$. By implementing the non-linear evolution in the xFitter framework, we have carried out fits of proton PDFs using data on lepton-proton deep inelastic scattering from HERA, BCDMS and NMC. While we find no evidence for the presence of non-linearities, they cannot be entirely excluded either and we are able to set limits for their strength. In terms of the recombination scale $Q_r$, which plays a similar role as the saturation scale in the dipole picture of the proton, we find that $Q_r \lesssim 2.5\,{\rm GeV}$. We also quantify the impact of the non-linear terms for the longitudinal structure function at the Electron-Ion Collider and the Large Hadron Electron Collider and find that these measurements could provide stronger constraints on the projected effects.


I. INTRODUCTION
In the framework of Quantum Chromodynamics (QCD) and collinear factorization [1], the cross sections for hard processes involving intial-state hadrons are calculated as convolutions of perturbative, process-dependent coefficient functions and non-perturbative, process-independent parton distribution functions (PDFs) f i (x, Q 2 ).The PDFs describe the density of partons (gluons and quarks) of flavor i carrying a fraction x of the hadron's momentum at a resolution scale Q.While the x dependence cannot be calculated by the perturbative means of QCD, its Q 2 dependence is governed by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [2][3][4][5].These evolution equations originate from the renormalization of divergences induced by collinear radiation from the initial-state partons.As a result of the radiation, the density of partons that carry a momentum smaller than that from which the radiation tree started, grows, leading to a rise of PDFs at low momentum fractions x as the scale Q 2 increases.This effect is especially pronounced for the gluon distribution.However, if the parton densities become sufficiently high, also the inverse processes in which several partons recombine, will begin to play a part, moderating the growth of PDFs at small x as Q 2 increases.These contributions are formally of higher-twist, i.e., they are suppressed by inverse powers of Q 2 .Since the Q 2 -dependence of PDFs is only logarithmic, the effects of recombination should disappear as Q 2 grows.
In the language of PDFs, the first recombination terms in the evolution can be derived from diagrams, where two initial-state gluons merge into collinear gluons/quarks that eventually participate in the hard scattering.The leading contributions at small x were originally calculated in Refs.[6][7][8][9] giving rise to the Gribov-Levin-Ryskin-Mueller-Qiu (GLR-MQ) corrections to the DGLAP evolution equations.Phenomenological applications of this approach have been studied in a number of publications, including its impact on observables in deep inelastic scattering (DIS) [10], comparisons with HERA data [11,12], its relation to the Balitsky-Kovchegov (BK) evolution equation [13,14], and the evolution of nuclear PDFs [15].However, the GLR-MQ equations violate the momentum conservation.This issue was fixed in the calculation of Zhu and Ruan [16,17], which is based on the same class of diagrams, but keeps also the non-leading contributions, see also Ref. [18].The resulting evolution equation is the one whose implications we will here study the framework of global fits of PDFs.The effect of the non-linear corrections is to moderate the Q 2 growth of PDFs at small x, while somewhat increasing the speed of evolution at moderate values of x such that the total momentum is conserved.This warrants to speak about dynamically generated shadowing and antishadowingterms which are perhaps more familiar from the context of nuclear PDFs [19].
For dimensional reasons, the higher-twist effects imply a presence of a dimensionful parameter R, which controls the strength of the non-linear effects and which can be associated with the size of the area in the transverse plane, where the partonic overlap leading to the non-linear corrections becomes important.The inverse of this parameter 1/R = Q r gives correspondingly the momentum scale -the recombination scale as we will call it -below which the effects of recombination will be considerable.Once Q 2 is so low that the non-linear terms in the evolution are as important as the linear ones, also terms suppressed by higher inverse powers of Q 2 become presumably important and the resummation of all these terms can be thought to give rise to the phenomenon of partonic saturation [9].This regime is traditionally discussed in the framework of the Color Glass Condensate (CGC) effective field theory [20,21] and is implemented in the dipole picture of DIS.In this case, the saturation effects are characterized by the saturation momentum scale Q s , which can be related to the critical dipole transverse size.It would be tempting to link the two emergent scales Q r and Q s .However, in the dipole language Q s also depends on x roughly as Q s ∼ x −b , whereas in the picture discussed here it is a constant, with the x dependence being entirely dynamical.Thus the two cannot be compared directly and the comparison is at best only indicative.
In the presented work, we have implemented the first momentum-conserving non-linear corrections into the HOPPET evolution code [22] and interfaced it with the xFitter framework [23,24].This allows us to perform global fits of proton PDFs with a range of different values for R. Scanning through various values of R and inspecting the resulting goodnesses of the fits and comparisons with the data, we are able to place lower limits for R, or equivalently upper limits for the recombination scale Q r .We have used BCDMS [25], HERA [26] and NMC data [27] on lepton-proton DIS in our analysis.Apart from searching for signatures of and placing limits on the effects of recombination from these DIS data, our motivation is to qualitatively compare the systematics of non-linearities with those resulting from the Balitsky-Fadin-Kuraev-Lipatov (BFKL) resummation [28][29][30][31] of log x-type terms relevant at small x.The BFKL resummation was incorporated into the DGLAP equations in Ref. [32] and evidence in favor of these dynamics has thereafter been seen in fits to the HERA data [33,34].It could be expected that the general features of the partonic recombination and resumming BFKL logarithms are, however, qualitatively similar at small-x as both tend to tame the rapid Q 2 growth of the gluon PDF.
While in the present paper we discuss the recombination in the case of free protons, we note that the non-linear corrections and partonic saturation are expected to be enhanced in the case of heavier nuclear targets [6,8,15,35,36], the recombination scale increasing naively as A 1/6 as a function of the nuclear mass number A. As of today, however, no clear signal has been seen.To discover this enhancement is also one of the cornerstones of the envisioned physics programs of the Electron-Ion Collider (EIC) [37,38] and Large Hadron Electron Collider (LHeC) [39,40] and a strong motivation to run the colliders with both light and heavy nuclei.The construction of a Future Circular Collider (FCC) [41], which would reach kinematics even further beyond those of the LHeC in electron-ion collisions, would naturally be able to put even stronger constraints on the nature and strength of non-linear corrections to the evolution of PDFs.
The paper is organized as follows: First, Sec.II provides an overview of the theoretical background and the implementation of the Zhu and Ruan non-linear corrections to the evolution of PDFs, including numerical comparisons between the linear and non-linear evolution equations.In Sec.III we present our new global fits of proton PDFs with non-linear corrections.Section IV discusses the impact on the longitudinal structure function F L (x, Q 2 ) in the HERA, EIC and LHeC kinematics.Finally, we summarize our findings and outline the way forward in Sec.V.

II. NON-LINEAR EVOLUTION EQUATIONS
The DGLAP evolution equations are a set of renormalization group equations arising from the renormalization of collinear divergences in ladder-type Feynman graphs describing parton emission in QCD.They take the following standard form [2][3][4][5], where q i (x, Q 2 ) are the PDFs for quark of flavor i, G(x, Q 2 ) is the gluon PDF and P ij (x/y) are the partonic splitting functions.These equations are clearly linear in PDFs.The first non-linear corrections stem from diagrams like the one in Fig. 1 showing a contribution to the DIS process with two gluons from the proton merging into one, which eventually splits into a quark-antiquark pair coupling to the virtual photon.Here, we will consider the corrections to the evolution of the gluon and quark-singlet distributions, which can be written as [16,17], where "linear terms" refer to the right-hand side of the usual DGLAP evolution, see Eq. (1).The square of the gluon PDF G 2 (x 1 , Q 2 ) arises due to the modelling of the four-gluon correlation function, see the lower part of Fig. 1.The parameter R with the dimension of length is introduced on the basis of dimensional analysis.In a sense, the factor 1/(RQ) 2 can be thought to be proportional to the overlap probability of two partons and R regarded as the length scale on which it takes place.A smaller value for the parameter R corresponds to stronger non-linear effects.The recombination functions in Eqs. ( 3) and (4) are given by i i with C A = N = 3 and T f = 1/2.Unlike the GLR-MQ equations, these evolution equations conserve the PDF momentum sum rule, because the non-linear terms satisfy the following conditions separately for the gluon and singlet quark distributions, Note that the linear parts of the full evolution equations also conserve the total momentum, but only in the sum of the gluon and quark singlet PDFs.The method employed in Refs.[16,17] can also be used to derive the quark-antiquark and quark-gluon recombination, but these effects are omitted in the present analysis.We also note that the gluon and quark singlet PDFs form a closed set of evolution equations, and that the quark flavor dependence enters through different non-singlet combinations of quarks, where the gluon recombination terms cancel.It is therefore consistent to consider the non-linear terms only in the evolution equations for gluon and quark singlet PDFs, Eqs. ( 3) and (4), while still solving for the full flavor-dependent evolution of PDFs.
We have implemented the non-linear corrections as an extension to the evolution toolkit HOPPET [22].This code is optimized to rapidly solve the linear DGLAP equations by converting the numerical convolution of the PDFs with the splitting functions to multiplication of a matrix and a vector, which can be computed extremely fast on modern hardware.To do this, HOPPET first changes the variables from x and Q 2 to y = ln(1/x) and t = ln Q 2 , which helps with numerical accuracy.The PDF of each flavor can be represented numerically by the sum of a set of PDF values at fixed y values y α , weighted with the interpolation weights w α (y): Using Eq. ( 10), the convolution (P ⊗ q) can then be written in an analogous way, where HOPPET precomputes the values of the matrix P αβ (t) once and reuses it for all future calculations, To apply this formalism to the non-linear terms, we treat G 2 (x, Q 2 ) as an independent flavor, which then enters the evolution as a linear term and is updated to be equal to [G(x, Q 2 )] 2 at each iteration.The integration limits are handled by substituting x with x/2 in the first non-linear term and by using a modified integration method, which integrates only up to y = ln(2) for the second term.This method allows us to compute the evolution with non-linear corrections, while increasing the computing time by no more than ≈ 20%.We have also numerically checked our implementation by verifying the momentum sum rule.In Fig. 2 we illustrate the effects of non-linearities in comparison to the regular DGLAP evolution and show the ratios between the PDFs evolved using the non-linear and linear evolution equations as a function of x at different values of Q.Here, we have initialized the evolution with the boundary condition given by the CJ15 proton PDFs [42] at Q 0 = 1.3 GeV, and performed the evolution in Q 2 according to the linear and non-linear equations up to different target values of Q.In both cases the linear terms in the evolution are included up to next-to-next-to-leading order (NNLO) accuracy in the strong coupling constant α s .The parameter R is set to 0.5, 1 and 2 GeV −1 in the top, center and bottom panels, respectively -we will see in the next section that these are reasonable choices.The Q values are chosen as Q 0 + ∆Q with Q 0 = 1.3 GeV and ∆Q increasing by factors of 100 from 10 −4 to 10 4 GeV.The dynamically generated suppression at small x ("shadowing") and an enhancement at intermediate x ("antishadowing") are clearly visible for both quark singlet (right panels) and gluon (left panels) PDFs.The effects are generally more significant for the gluon PDF than for the quark singlet, which, in part, results from the fact that the gluon distribution at the initial scale is smaller than the quark singlet one -very close to zero at small x -and therefore even small changes are more significant in the ratios.As a result, even a tiny step to higher Q causes visible differences in the case of gluons, whereas for the quark singlet much larger step has to be taken in order to see a difference.While the relative effects of non-linearities first increase as Q increases, there is eventually a turnaround and the differences begin to decrease.This was to be expected as so that the non-linear terms die away at asymptotically large Q 2 .For gluons this turnaround happens already around ∆Q ≈ 1 GeV, but for quarks it takes place only at the electroweak scale ∆Q ≈ 100 GeV.The dependence on the choice of R works as expected -a smaller R leads to more pronounced effects although qualitatively there are no significant differences between different choices of R.

III. PDF FITS WITH RECOMBINATION EFFECTS
As was observed in Fig. 2, the effects of recombination can persist up to rather high values of Q 2 , where the PDFs are constrained by experimental data.As a result, in order to be consistent with these data constraints, the initial condition at Q 0 has to be iterated.To get a more complete and consistent picture of the impact of non-linearities it is thus necessary to perform global PDF fits including the non-linear terms in the evolution.In addition to the corrections to the evolution equations, higher-twist corrections also exist for the DIS coefficient functions like the ones for the structure function F 2 [16][17][18].However, we neglect these effects in the present study.

A. Fitting Framework
The framework used to perform our global PDF fits with non-linear corrections is based on the xFitter package [23,24], which we have extended by adding our modified version of HOPPET as a new module, see Sec.II.We have verified that this reproduces the results obtained with the default QCDNUM [43,44] evolution in the limit R → ∞.The fits are performed at the NNLO accuracy in both evolution and scattering coefficients.We use the RTOPT general-mass variable-flavor-number scheme [45] for the heavy quark coefficient functions.In order to study the dependence of the non-linear effects on the assumed parameterization, we have prepared PDF fits with two different parameterizations for gluons and also varied the initial scale.We refer to these as parameterization 1 and 2. For parameterization 1 the PDFs are parameterized at the initial scale Q 0 = 1.0 GeV using the same parametric form as employed in the HERAPDF2.0PDF analysis [26]: For the valence distributions u v , d v , and the sea-quark densities ū, d, we use the form and the gluon PDF is parameterized as When using the parameterization in Eq. ( 14), the gluon distribution tends to become negative at small values of x and Q 2 .In the case of parameterization 2 we apply the same parameterization for gluons as we have for quarks in Eq. ( 13), which guarantees that the small-x gluons will remain positive at the initial scale.To partly compensate for this restriction we increase the initial scale to Q 0 = 1.3 GeV for parameterization 2. The strange quark is set to xs = xs = (2/3)x d at Q 0 in both cases.The normalization parameters A uv , A dv and A g are fixed through the quark-number and momentum sum rules.For simplicity and for lack of constraints for the flavor separation, we fix A ū = A d and B ū = B d.With parameterization 1, we fix C ′ g = 25, where the exact value is unimportant as long as C ′ g ≫ C g such that the negative term does not contribute beyond low x.The full set of 14 open PDF parameters in the case of parameterization 1 is, and in the case of parameterization 2, Any parameters not explicitly mentioned are set to zero.The parameter R is kept fixed in each fit and the procedure is repeated for R values from 0.2 GeV −1 to 3.0 GeV −1 in steps of 0.01 GeV −1 .For R > ∼ 3.0 GeV −1 , the non-linear modifications were found to be negligible.
We find the optimal set of parameters by minimizing the χ 2 function defined as The theoretical predictions for data point i and a given set of PDF parameters are represented by T i .The measured data values are given by D i and δ i,stat , δ i,uncorr and γ i α are the statistical, uncorrelated and correlated uncertainties associated with each data point.Minimizing this χ 2 function with respect to the 14 fit parameters and the nuisance parameters b α determines the central set of PDFs and systematic shifts, − α γ i α D i b α , which give the best description of the data.The uncertainties of the obtained PDFs are then determined using the Hessian method [46], which is based on a quadratic expansion of the χ 2 function around its minimum, χ 2 = χ 2 min +∆χ 2 .The choice of the maximum displacement ∆χ 2 , the tolerance, depends on the definition of uncertainties and in PDF fits it is usually conservatively chosen to be ∆χ 2 ≫ 1.Indeed, assuming that the data points are Gaussianly distributed around the "truth", the expected χ 2 /N dof should be less than ∼1.05 in 90% of the cases with ∼1600 data points, which is the amount of points we are now considering.However, the minimum values we find are χ 2 /N dof ∼1.18, i.e., the ideal choice ∆χ 2 = 1 is, perhaps, an underestimate.Assuming a Gaussian probability density for the fit parameters results in ∆χ 2  max ≈ 20 at the 90% confidence level [47].In practice, we determine the PDFs with the tolerance set to ∆χ 2 = 1 and then rescale the resulting uncertainties by √ 20 ≈ 4.5 to match the actual tolerance since this leads to a better convergence of the xFitter algorithms.

B. Data selection
Our analysis, as the ones presented in Refs.[48] and [49], includes the lepton-proton DIS data from the BCDMS, HERA and NMC experiments, which are summarized in Table I.The data from BCDMS and NMC are given in terms of the structure function F p 2 (x, Q 2 ) measured in neutral-current (NC) deep inelastic muon scattering.The HERA data sets provide the reduced cross section of DIS in electron-proton and positron-proton collisions, both for charged-current (CC) and NC processes.We do not include the HERA data on the longitudinal structure function F L [50] since it is derived from the other HERA data sets, and including it would to some extent double-count some constraints.We will, however, discuss F L separately in Sec.IV.The data on heavy-quark production from HERA are omitted as well [51][52][53].
In the first round of fits, we apply the same Q 2 > 3.5 GeV 2 cut as, e.g., in the HERAPDF2.0analysis, but with parameterization 1 we also perform a second round of fits with the cut lowered to Q 2 > 1.0 GeV 2 to see whether the inclusion of the 1/Q 2 suppressed terms in the evolution equation can provide an improved description of the data in the low-Q 2 region.Indeed, it is known that in NNLO fits, χ 2 /N dof tends to get systematically worse as the minimum Q 2 is lowered from 15 GeV 2 to 2 GeV 2 [33,34] -a tension that the BFKL resummation appears to ease.The two rightmost columns of Table I indicate the number of data points per experiment, which pass these cuts.The total number of data points is 1568 for the higher-Q 2 cut and 1636 for the lower-Q 2 one, respectively.

C. Fit results
The resulting values of χ 2 /N dof minima for all the fits with different R values are shown in Fig. 3 and the contributions of different data sets to the total χ 2 /N dof are presented in Appendix A, Tables II, III and IV.In Fig. 3, the orange (parameterization 1) and green (parameterization 2) crosses correspond to the Q 2 > 3.5 GeV 2 TABLE I: Summary of DIS data sets used in out fits of proton PDFs with non-linear corrections.The last two columns show the number of data points passing the Q 2 > 3.5 GeV 2 and Q 2 > 1 GeV 2 cuts, respectively.The energies listed for the HERA measurements refer to the proton beam energy, while those listed for BCDMS refer to the muon beam energy.NMC uses the same energies as BCDMS but combines them into one dataset.

Experiment
Data set Year Ref. cut and the blue ones (parameterization 1) are for the fits with Q 2 > 1.0 GeV 2 .We have not performed a fit with parameterization 2 with the lower cut, as it lies below the initial scale of the PDF evolution in this case.
In the case of parameterization 1, the resulting values of χ 2 /N dof generally decrease with increasing R, i.e., with a decreasing strength of the non-linear corrections, but remain constant for R > 1.5 GeV −1 .Lowering the value of R down to ≈ 0.7 GeV −1 causes only a small increase, but below that χ 2 /N dof rapidly starts diverging.Comparing the curves between the two different cuts, there is no qualitative difference in the R dependence, but the lower cut shifts χ 2 /N dof upwards by approximately 0.04 units across all values of R. While this increase in χ 2 /N dof from the lowered cut might not seem too significant at first sight, one has to consider the small number of added data points that causes it.Namely, the 68 added data points increase the total χ 2 by 145, i.e., more than 2 units per point.Therefore, the non-linear corrections do not help to alleviate the tensions with lower-Q 2 data.In the case of parametrization 2 (green crosses), i.e., with a positive-defined gluon distribution at the initial scale, the best descriptions are obtained in the interval 0.5 < R < 0.7 GeV −1 although the difference with respect to R = 3 GeV −1 is small.
Tables II, III and IV in Appendix A break down the obtained χ 2 into contributions from different data sets.Tables II and IV, which correspond to the fits with parameterization 1 at Q 2 cut = 3.5 GeV 2 and 1 GeV 2 , respectively, show that the preference towards higher values of R is a feature generally shared by all the data sets.At the same time, some data sets are more sensitive to R than others with the NMC data set showing the largest difference in χ 2 /N dof .Table III shows that with parametrization 2, the χ 2 /N dof profiles as a function of R are flatter than in the case of parametrization 1.
In the following, we will discuss proton PDFs determined from the fits with the Q 2 > 3.5 GeV 2 cut for a representative selection of R values of {0.4,0.6, 0.8, 1, 2, 3} GeV −1 .Figures 4 and 5 show the resulting gluon (left) and quark singlet (right) distributions at Q = {Q 0 , 2 GeV, 10 GeV, 100 GeV} for parameterization 1 and 2, respectively.They are compared to the HERAPDF2.0results given by the black dashed curve.To avoid visual clutter, the uncertainties (90% confidence level) are only shown for HERAPDF2.0and R = 3.0 GeV −1 case, but the remaining uncertainties of the rest closely resemble the latter one.The 90% confidence level for HERAPDF2.0 is obtained by multiplying the uncertainties, which are given at 68% confidence level with a factor of 1.645.
In the case of parameterization 1, the fit with the highest value of R, i.e., the one with smallest non-linear effects, behaves similarly as the HERAPDF2.0reference but, despite the added data from BCDMS and NMC, has significantly larger uncertainties due to the higher Hessian tolerance -20 for our fit vs. 2.71 for HERAPDF2.0at 90% confidence level.The differences in the central values can, at least partly, be attributed to the fact the we included these additional data, which may affect the PDFs even at low-x due to the sum rules.As the value of R is lowered, the shapes of the fitted PDFs change progressively further away from the baseline with R = 3.0 GeV −1 .In particular, one observes a significant suppression of the gluon and quark-singlet PDFs at small x and Q ≥ 2 GeV as R decreases.This suppression becomes weaker as the scale Q is increased due to the weakening of the non-linear effects by a factor of 1/Q 2 in the evolution equations, but it is still visible even at Q = 100 GeV.As expected, the largest differences between the PDFs corresponding to different values of R are visible at the initial scale Q 0 = 1 GeV, though they carry little meaning since there are no data that constrain this region due to the Q 2 cut > 3.5 GeV 2 cut.A comparison of the results presented Figs. 4 and 5 shows that the systematics of the R and Q dependence of the PDFs obtained using parameterizations 1 and 2 are very similar even though there are large differences around the initial scale Q 0 .This gives us confidence that our conclusions are not very sensitive to the details of the parameterization.The parameters of the fits are summarized in Tables V, VI and VII in Appendix B.

IV. IMPACT ON THE LONGITUDINAL STRUCTURE FUNCTION
The longitudinal structure function F L carries a direct sensitivity to the gluon PDFs and should therefore have an increased sensitivity to the non-linear effects.We note that since the HERA data entered our analysis as reduced cross sections σ r = F 2 − (y 2 /Y + )F L , where Y + = 1 + (1 − y) 2 , y = s/xQ 2 , s being the squared center-of-mass energy of the lepton-proton collision, they already carry information of F L .In addition, the Q 2 dependence of F 2 is directly sensitive to F L [54], particularly at small x.In this sense the discussion below is not entirely independent of the fits presented above.

A. The longitudinal structure function at HERA
To gain a deeper understanding of the sensitivity of the HERA data to non-linear corrections, we compare the values for the longitudinal structure function F L (x, Q 2 ) that result from our PDFs with the data taken by H1 [50] and ZEUS [55]. Figure 6 presents the HERA F L measurements at different values of Q 2 and integrated over x bins.The left panel shows the data compared with the predictions using the PDFs fitted at Q 2 cut = 3.5 GeV 2 with parameterization 1 and the right panel with parameterization 2. For Q 2 > ∼ 100 GeV 2 the predictions converge and become indistinguishable.Towards lower Q 2 values the predictions spread out with higher R values corresponding to slightly larger values of F L .At low Q 2 , the behaviour depends on the parameterization.In the case of parameterization 2, the predictions converge again towards F L → 0 at the PDFs parameterization scale Q 2 0 = 1.69 GeV 2 .The predictions for parameterization 1 do not converge towards F L → 0 since the parameterization allows the PDFs to be The same as Fig. 4, but showing PDFs fitted with parameterization 2. FIG.6: Structure function F L calculated with PDFs resulting from different choices of R compared with the measurements by H1 [50] and ZEUS [55].We only show predictions for negative and therefore yield negative GeV 2 , the NNLO predictions generally lie below the values for F L measured by H1, as has been observed in other studies [33,34].On the other hand, the ZEUS measurements generally lie below the predictions, but come with larger uncertainties.Interestingly, the systematics of non-linearities we find goes in the oppostite direction in comparison to the effects of small-x BFKL resummation: While the non-linear effects tend to decrease F L , the small-x resummation has a tendency to increase it.This can be traced to the small-x behaviour of gluon and quark-singlet PDFs which are, at perturbative values of Q 2 , increased in the BFKL approach in comparison to the fixed-order NNLO result.In the case of non-linearities discussed here, the gluon and quark-singlet PDFs are generally below the standard NNLO result at small x.At very low values of Q 2 -very close to the parametrization scale -they go in the opposite way, but the evolution is slower in the case with non-linear effects and leads to suppressed gluon and quark singlet PDFs very quickly as Q 2 increases.As a result, the systematics of BFKL and non-linear effects are quite different.
Figure 7 shows F L as a function of x at different values of Q 2 .The theoretical calculations are performed with the PDFs corresponding to fits with different values of R using parameterization 2. Predictions for F L made with PDFs from parameterization 1 result in negative values of F L at Q 2 ≤ 10 GeV 2 due to the gluon PDFs not being restrained to positive values at the initial scale.Since the negativity of F L seen in the previous figure indicates that the corresponding theoretical setup is not entirely consistent, e.g., too low Q 2 or too low R, we omit the PDFs fitted with parameterization 1 from the following discussion.For clarity, we show the uncertainty band only for the R = 3 GeV −1 case.One can see from the figure that for x > ∼ 0.3 × 10 −2 , the predictions with different R are almost indistinguishable.However for Q 2 > 5 GeV 2 and towards lower x values, the calculations with a smaller R lead to significantly reduced F L .This was to be expected as both gluon and quark-singlet PDFs at small x were found to be increasingly suppressed as R decreased, see Figs. 4 and 5.At lower values of Q 2 , where the relative differences between predictions corresponding to different R values are expected to be the largest, the absolute values are too close to zero to see any meaningful differences.As already observed in Sec.II, the effects of non-linearities die out rather slowly at small x as Q 2 increases.Also here, we see in Fig. 7 that the differences between calculations with different R actually still increase at x ≈ 10 −5 until the highest considered Q 2 .We note that at Q 2 > 8.5 GeV 2 the F L data are well described by most of the predictions.However, with R ≤ 0.4 GeV −1 the description begins to visibly deteriorate even at higher values of Q 2 .

B. Prospects for the longitudinal structure function at EIC and LHeC
The observations made above naturally raise a question of how much better future DIS experiments such as EIC [37] and LHeC [39,40] could constrain the presence of non-linearities.To this end, we compare the spread of F L predictions by taking the ratio of the prediction resulting from each set of PDFs determined at a different value of R over that at R = 3 GeV −1 : These ratios are compared to the projected relative uncertainties for F L taken from the EIC development documents [56] and the LHeC White Paper [39] in Figs. 8 and 9, respectively, which are shown as black error bars.Again, we only show the results for parameterization 2 due to the negativity of F L when calculated with the PDFs from parameterization 1.This allows us to estimate whether F L measurements at these future experiments can help to constrain the size of non-linear effects under the assumption that F L can be accurately reproduced by the theory, e.g., by adding small-x resummation as in Refs.[33,34].The projected EIC data only reach the x region, where the different R values are clearly distinguishable at Q < 2.5 GeV 2 , but even at Q = 2.47 GeV 2 the uncertainties of the pseudodata are barely smaller than the the spread between predictions with R > 0.6 GeV −1 .In the case of LHeC, however, the reach towards lower x values should lead to more stringent constraints for the non-linear effects, since the different R values lead to vastly different predictions in this kinematic region.x x ) The same as Fig. 8, but for LHeC [39].

V. CONCLUSIONS AND OUTLOOK
We have presented numerical studies of the non-linear gluon recombination corrections to the DGLAP evolution as derived by Zhu and Ruan, which improve upon the commonly used GLR-MQ equations and restore the PDF momentum sum rule.We have extended the HOPPET and xFitter toolkits to account for these corrections.With these extensions we performed several fits of proton PDFs to BCDMS, HERA and NMC data on DIS at NNLO accuracy varying the dimensionful parameter R, which controls the strength of the recombination effects.We examined two different parametrizations for the gluon PDF -one corresponding to HERAPDF2.0,which allows the gluon PDF to go negative at the parametrization scale, and another which is positive definite at the parametrization scale.We found that both parameterizations result in a similarly good description of the data, χ 2 /N dof ∼ 1.19, for R > 1 GeV −1 , but that the description begins to eventually deteriorate towards small R. For the parametrization allowing for a negative gluon PDF this deterioration sets in for R < ∼ 0.7 GeV −1 , but for the positive-definite parametrization at somewhat lower values, R < ∼ 0.5 GeV −1 .In both cases R < 0.4 GeV −1 seems to be excluded, which translates to an upper limit for the recombination scale Q r = 1/R < ∼ 2.5 GeV.We also studied the sensitivity of our fits to the choice of Q 2 cut : lowering Q 2 cut from 3.50 GeV 2 to 1 GeV 2 leads to an increase of χ 2 /N dof , which the non-linear corrections are unable to tame.In other words, we do not find evidence for the presence of non-linear effects in the used DIS data, but can set an upper limit for the parameter controlling their strength, R > ∼ 0.4 GeV −1 or Q r = 1/R < ∼ 2.5 GeV.Our conclusions about HERA data not supporting non-linear effects is in line with a recent study made in Ref. [57].We also compared our PDF fits with the longitudinal structure function F L available from HERA.While the ZEUS data lies below the predictions from our fits, the large uncertainties still make the two compatible.However, the H1 data, which reaches lower Q 2 values and has smaller uncertainties, lies well above the predictions at low Q 2 .Non-linear effects in the evolution do not help to resolve this tension between NNLO calculations and the data, indicating that other theoretical ingredients, such as small-x resummation are necessary to describe these data properly.Assuming that this tension can be resolved, we show that the improved statistics of the EIC, and particularly the reach towards lower x values at the LHeC could provide significantly stronger constraints on non-linear effects.
In general, after a short interval of Q 2 evolution, the non-linear effects tend to reduce both gluon and quark PDFs at small values of x in comparison to the case with no recombination.This suppression becomes more pronounced as the parameter R decreases and the effects can persist up to Q > ∼ 100 GeV.This behavior was not sensitive to the form of the initial parameterization and thereby appears to be a general feature of the non-linear effects.This is opposite to the BFKL resummation which tends to increase the small-x PDFs at large Q 2 .Nevertheless, the fact that the Q 2 evolution does not wash away all the possible remnants of non-linearities even at Q = O(100 GeV) means that some effects could be seen in proton-proton collisions at the LHC as well.For example, the direct photon and heavy-flavor production at forward direction are sensitive to PDFs at small x and some visible differences could be induced from the fits with different R. To chart such effects is one possible way to make use of the results of the present analysis.One can also extend the analysis itself by, e.g., including the contributions of two-gluon initiated processes in the coefficient functions as well.For example, the recombination contributions enter the coefficient functions of F L at O(α 2 s ), i.e., they are of the same perturbative order as the leading-twist NLO terms.As these leading-twist NLO terms turn out rather significant for F L , the inclusion of higher-twist terms could potentially change the picture in the case of F L .For F 2 the higher-twist terms are of the same perturbative order as the leading-twist NNLO terms.In addition, one could also include a wider range of data sets in the fit, e.g., from the LHC, and/or generalize the framework to the case of heavier nuclei with the expectation of observing enhanced effects of recombination.
The tools developed in this work and LHAPDF files [58] for the fitted PDFs are available here: https://research.hip.fi/qcdtheory/nuclear-pdfs/.

FIG. 1 :
FIG.1:A typical cut diagram representing the recombination of two gluon ladders into one.Note that the gluon ladders in black and red do not interact.The momentum fractions x and x 1 appearing in Eqs.(3) and (4) are indicated.

x 1 FIG. 2 :
FIG.2: Ratios between the PDFs evolved using the non-linear and linear evolution equations with the CJ15 initial condition at Q 0 = 1.3 GeV as a function of x and at different values of Q = Q 0 + ∆Q.The left column corresponds to the gluon PDF, while the right column shows the quark singlet distribution.The top, middle and bottom rows correspond to R = 1/2, 1 and 2 GeV −1 , respectively.

FIG. 4 :
FIG.4: PDFs resulting from fits with parameterization 1 at various R values.The left column shows the gluon PDF, while the right side shows the quark singlet PDF.The rows correspond to different scales Q.

2 )Q 2 =60.0 GeV 2 FIG. 7 :
FIG.7: Structure function F L resulting from PDFs fitted with parameterization 2 and different R values compared with the HERA F L measurements[50].We only show predictions for Q 2 > Q 2 0 .

2 FIG. 8 : 2 FIG. 9 :
FIG. 8:The ratios of the longitudinal structure functions corresponding to different values of R to that with R = 3 GeV −1 , see Eq. (18), as a function of x.The calculations use parametrization 2 and are carried out at different values of Q 2 .The vertical bars represent the expected statistical uncertainty of the corresponding measurements at the EIC[56].We only show predictions for Q 2 > Q 2 0 .

TABLE III :
The same as TableII, but with parameterization 2.

TABLE VI :
The same as TableV, but for PDFs fitted with the lower cut Q 2 cut = 1.0 GeV 2 .

TABLE VII :
The same as TableV, but for PDFs fitted with parameterization 2.