# Proton PDFs with non-linear corrections from gluon recombination

P. Duwentäster<EMAIL_ADDRESS>University of Jyvaskyla, Department of
Physics, P.O. Box 35, FI-40014 University of Jyvaskyla, Finland  Helsinki
Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland
V. Guzey<EMAIL_ADDRESS>University of Jyvaskyla, Department of Physics,
P.O. Box 35, FI-40014 University of Jyvaskyla, Finland  Helsinki Institute of
Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland  I. Helenius
<EMAIL_ADDRESS>University of Jyvaskyla, Department of Physics, P.O.
Box 35, FI-40014 University of Jyvaskyla, Finland  Helsinki Institute of
Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland  H. Paukkunen
<EMAIL_ADDRESS>University of Jyvaskyla, Department of Physics, P.O.
Box 35, FI-40014 University of Jyvaskyla, Finland  Helsinki Institute of
Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland

###### Abstract

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
Collins:1989gx , 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 Dokshitzer:1977sg ; Gribov:1972ri ; Gribov:1972rt ; Altarelli:1977zs
. 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. Qiu:1986wh ;
Kovchegov:2012mbw ; Mueller:1985wy ; Gribov:1983ivg  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) Bartels:1990zk , comparisons with HERA data
Prytz:2001tb ; Eskola:2002yc , its relation to the Balitsky-Kovchegov (BK)
evolution equation Boroun:2013mgv ; Boroun:2009zzb , and the evolution of
nuclear PDFs Rausch:2022nkv . However, the GLR-MQ equations violate the
momentum conservation. This issue was fixed in the calculation of Zhu and Ruan
Zhu:1998hg ; Zhu:1999ht , which is based on the same class of diagrams, but
keeps also the non-leading contributions, see also Ref. Blumlein:2001bk . 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 antishadowing – terms which are perhaps more familiar
from the context of nuclear PDFs Klasen:2023uqj .

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 Gribov:1983ivg .
This regime is traditionally discussed in the framework of the Color Glass
Condensate (CGC) effective field theory Gelis:2010nm ; Morreale:2021pnn  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}\sim 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 Salam:2008qg  and interfaced
it with the xFitter framework xFitter:2022zjb ; Alekhin:2014irh . 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 BCDMS:1989qop , HERA H1:2015ubc  and NMC data
NewMuon:1996fwh  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 Fadin:1975cb ; Kuraev:1976ge ; Kuraev:1977fs ;
Balitsky:1978ic  of $\log x$-type terms relevant at small $x$. The BFKL
resummation was incorporated into the DGLAP equations in Ref. Bonvini:2017ogt
and evidence in favor of these dynamics has thereafter been seen in fits to
the HERA data Ball:2017otu ; xFitterDevelopersTeam:2018hym . 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 Mueller:1985wy
; Qiu:1986wh ; Eskola:1993mb ; Eskola:2003gc ; Rausch:2022nkv , 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) Accardi:2012qut ;
AbdulKhalek:2021gbh  and Large Hadron Electron Collider (LHeC) LHeC:2020van ;
LHeCStudyGroup:2012zhm  and a strong motivation to run the colliders with both
light and heavy nuclei. The construction of a Future Circular Collider (FCC)
FCC:2018byv , 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 Dokshitzer:1977sg ; Gribov:1972ri ; Gribov:1972rt ;
Altarelli:1977zs ,

$\displaystyle Q^{2}\frac{d}{dQ^{2}}\begin{pmatrix}q_{i}(x,Q^{2})\\\\[5.0pt]
G(x,Q^{2})\end{pmatrix}=\frac{\alpha_{s}(Q^{2})}{2\pi}\sum_{j}\int_{x}^{1}\frac{dy}{y}\begin{pmatrix}P_{q_{i}q_{j}}\left(\frac{x}{y}\right)&P_{q_{i}g}\left(\frac{x}{y}\right)\\\\[5.0pt]
P_{gq_{j}}\left(\frac{x}{y}\right)&P_{gg}\left(\frac{x}{y}\right)\end{pmatrix}\begin{pmatrix}q_{j}(y,Q^{2})\\\\[5.0pt]
G(y,Q^{2})\end{pmatrix}\,,$ (1)

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,

$p$$(proton1)+(0cm,0.1cm)$$p$$(proton1)+(0cm,-0.1cm)$$p$$(proton1)+(2cm,0cm)$$(pBlob1)+(0cm,0.1cm)$$(pBlob1)+(0cm,-0.1cm)$$(proton1)+(2.5cm,0.5cm)$$(proton1)+(2.5cm,-0.5cm)$$(proton1)+(7.5cm,0.5cm)$$(proton1)+(7.5cm,-0.5cm)$$(proton1)+(8cm,0cm)$$(pBlob2)+(0cm,0.1cm)$$(pBlob2)+(0cm,-0.1cm)$$(proton1)+(10cm,0cm)$$p$$(proton2)+(0cm,0.1cm)$$p$$(proton2)+(0cm,-0.1cm)$$p$$(pBlobCircleLeft1)+(0cm,1cm)$$(pBlobCircleLeft1)+(0cm,2cm)$$(pBlobCircleLeft1)+(0cm,3.5cm)$$(pBlobCircleRight1)+(0cm,1cm)$$(pBlobCircleRight1)+(0cm,2cm)$$(pBlobCircleRight1)+(0cm,3.5cm)$$(ladderLeftBlueV3)+(0cm,1cm)$$(ladderRightBlueV3)+(0cm,1cm)$$(box1TopLeft)+(0cm,1cm)$$(box2BotLeft)+(0cm,1cm)$$(box1TopRight)+(0cm,1cm)$$(box2BotRight)+(0cm,1cm)$$(box2TopLeft)+(-1cm,1cm)$$(box2TopRight)+(1cm,1cm)$$(leptonGammaLeft)+(-1cm,1cm)$$l$$(leptonGammaRight)+(1cm,1cm)$$l$$(leptonIn)+(4.5cm,0.5cm)$$(cutTop)+(0cm,-11cm)$$(pBlobCircleLeft1)+(1cm,0cm)$$(pBlobCircleLeft1)+(1cm,1.5cm)$$(pBlobCircleLeft1)+(1cm,2.5cm)$$(pBlobCircleLeft1)+(1cm,3.5cm)$$(pBlobCircleLeft1)+(4cm,0cm)$$(pBlobCircleLeft1)+(4cm,1.5cm)$$(pBlobCircleLeft1)+(4cm,2.5cm)$$(pBlobCircleLeft1)+(4cm,3.5cm)$$x_{1}$$x_{1}$$x$$x$$\gamma$$\gamma$$x_{1}$$x_{1}$
Figure 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.

$\displaystyle\Sigma(x,Q^{2})=\sum_{i}\left[q_{i}(x,Q^{2})+\bar{q}_{i}(x,Q^{2})\right],\
\ \mathrm{with}\ i\in\\{u,d,s,c,b,t\\}\,,$ (2)

which can be written as Zhu:1999ht ; Zhu:1998hg ,

$\displaystyle x\frac{dG(x,Q^{2})}{d\ln{Q^{2}}}\ =$ linear terms
$\displaystyle+\frac{9}{32\pi^{2}}\left(\frac{1}{RQ}\right)^{2}\int_{x/2}^{1/2}dx_{1}xx_{1}G^{2}(x_{1},Q^{2})\sum_{i}P_{i}^{GG\rightarrow{G}}(x_{1},x)$
(3)
$\displaystyle-\frac{9}{16\pi^{2}}\left(\frac{1}{RQ}\right)^{2}\int_{x}^{1/2}dx_{1}xx_{1}G^{2}(x_{1},Q^{2})\sum_{i}P_{i}^{GG\rightarrow{G}}(x_{1},x),$
$\displaystyle\equiv$ $\displaystyle\quad\text{linear terms}+\Delta
G(x,Q^{2})\,,$ $\displaystyle x\frac{d\Sigma(x,Q^{2})}{d\ln{Q^{2}}}\ =$ linear
terms
$\displaystyle+\frac{9}{32\pi^{2}}\left(\frac{1}{RQ}\right)^{2}\int_{x/2}^{1/2}dx_{1}xx_{1}G^{2}(x_{1},Q^{2})\sum_{i}P_{i}^{GG\rightarrow{q\overline{q}}}(x_{1},x)$
(4)
$\displaystyle-\frac{9}{16\pi^{2}}\left(\frac{1}{RQ}\right)^{2}\int_{x}^{1/2}dx_{1}xx_{1}G^{2}(x_{1},Q^{2})\sum_{i}P_{i}^{GG\rightarrow{q\overline{q}}}(x_{1},x),$
$\displaystyle\equiv$ $\displaystyle\quad\text{linear
terms}+\Delta\Sigma(x,Q^{2})\,,$

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

$\displaystyle\sum_{i}P_{i}^{GG\rightarrow{G}}(x_{1},x)$
$\displaystyle=\frac{3\alpha_{s}^{2}}{8}\frac{C_{A}^{2}}{N^{2}-1}\frac{(2x_{1}-x)(-136xx_{1}^{3}-64x_{1}x^{3}+132x_{1}^{2}x^{2}+99x_{1}^{4}+16x^{4})}{xx_{1}^{5}}\,,$
(5) $\displaystyle\sum_{i}P_{i}^{GG\rightarrow{q\overline{q}}}(x_{1},x)$
$\displaystyle=\alpha_{s}^{2}\frac{T_{f}}{N(N^{2}-1)}(2x_{1}-x)^{2}\left[\frac{(4x^{2}+5x_{1}^{2}-6x_{1}x)}{x_{1}^{5}}+\frac{N^{2}}{(N^{2}-1)}\frac{(4x^{2}+4x_{1}^{2}-6x_{1}x)}{x_{1}^{5}}\right]\,,$
(6)

with $C_{A}=N=3$ and $T_{f}=1/2$. Unlike the GLR-MQ equations, these evolution
equations conserve the PDF momentum sum rule,

$\frac{d}{dQ^{2}}\int_{0}^{1}dxx\left[G(x,Q^{2})+\Sigma(x,Q^{2})\right]=0\,,$
(7)

because the non-linear terms satisfy the following conditions separately for
the gluon and singlet quark distributions,

$\displaystyle\frac{d}{d\ln{Q^{2}}}\int_{0}^{1}dx\Delta G(x,Q^{2})=0\,,$ (8)
$\displaystyle\frac{d}{d\ln{Q^{2}}}\int_{0}^{1}dx\Delta\Sigma(x,Q^{2})=0\,.$
(9)

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. Zhu:1999ht ; Zhu:1998hg  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 Salam:2008qg . 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_{\alpha}$, weighted with the interpolation weights $w_{\alpha}(y)$:

$\displaystyle q(y,t)\ \ $ $\displaystyle=\ \
\sum_{\alpha}w_{\alpha}(y)q(y_{\alpha},t)\ \ \equiv\ \
\sum_{\alpha}w_{\alpha}(y)q_{\alpha}(t)\,.$ (10) Using Eq. (10), the
convolution $(P\otimes q)$ can then be written in an analogous way,
$\displaystyle(P\otimes q)(y,t)\ \ $ $\displaystyle=\ \
\sum_{\alpha}w_{\alpha}(y)(P\otimes q)_{\alpha}(t)=\ \
\sum_{\beta}\sum_{\alpha}w_{\alpha}(y)P_{\alpha\beta}(t)q_{\beta}(t)\,,$ (11)
where HOPPET precomputes the values of the matrix $P_{\alpha\beta}(t)$ once
and reuses it for all future calculations, $\displaystyle P_{\alpha\beta}(t)\
\ $ $\displaystyle=\ \
\int_{e^{-y_{\alpha}}}^{1}dzP(z,t)w_{\beta}(y_{\alpha}+\ln z)\,.$ (12)

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 $\approx 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 Accardi:2016qay  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 $\alpha_{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}+\Delta Q$ with $Q_{0}=1.3$ GeV and $\Delta 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
$G^{2}(x,Q^{2})/Q^{2}\rightarrow 0$ as $Q^{2}\rightarrow\infty$ so that the
non-linear terms die away at asymptotically large $Q^{2}$. For gluons this
turnaround happens already around $\Delta Q\approx 1\,{\rm GeV}$, but for
quarks it takes place only at the electroweak scale $\Delta Q\approx 100\,{\rm
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$.

Figure 2:  Ratios between the PDFs evolved using the non-linear and linear
evolution equations with the CJ15 initial condition at $Q_{0}=1.3\,{\rm GeV}$
as a function of $x$ and at different values of $Q=Q_{0}+\Delta 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.

##  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}$ Zhu:1999ht ; Zhu:1998hg ; Blumlein:2001bk .
However, we neglect these effects in the present study.

###  III.1 Fitting Framework

The framework used to perform our global PDF fits with non-linear corrections
is based on the xFitter package xFitter:2022zjb ; Alekhin:2014irh , 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 Botje:2010ay ; Botje:2016wbq  evolution in the limit
$R\rightarrow\infty$. The fits are performed at the NNLO accuracy in both
evolution and scattering coefficients. We use the RTOPT general-mass variable-
flavor-number scheme Thorne:2006qt  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.0 PDF analysis H1:2015ubc : For the valence
distributions $u_{v},d_{v}$, and the sea-quark densities $\bar{u},\bar{d}$, we
use the form

$\displaystyle
xf(x)=A_{f}x^{B_{f}}(1-x)^{C_{f}}\left(1+D_{f}x+E_{f}x^{2}\right)\,,$ (13)

and the gluon PDF is parameterized as

$\displaystyle
xg(x)=A_{g}x^{B_{g}}(1-x)^{C_{g}}-A^{{}^{\prime}}_{g}x^{B^{{}^{\prime}}_{g}}(1-x)^{C^{{}^{\prime}}_{g}}\,.$
(14)

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~{}\text{GeV}$ for parameterization 2.
The strange quark is set to $xs=x\bar{s}=(2/3)x\bar{d}$ at $Q_{0}$ in both
cases.

The normalization parameters $A_{u_{v}},A_{d_{v}}$ 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_{\bar{u}}=A_{\bar{d}}$ and
$B_{\bar{u}}=B_{\bar{d}}$. With parameterization 1, we fix
$C^{\prime}_{g}=25$, where the exact value is unimportant as long as
$C^{\prime}_{g}\gg 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,

$\displaystyle\\{B_{g},C_{g},A^{\prime}_{g},B^{\prime}_{g},B_{u_{v}},C_{u_{v}},E_{u_{v}},B_{d_{v}},C_{d_{v}},C_{\bar{u}},D_{\bar{u}},A_{\bar{d}},B_{\bar{d}},C_{\bar{d}}\\}\,,$
(15)

and in the case of parameterization 2,

$\displaystyle\\{B_{g},C_{g},D_{g},E_{g},B_{u_{v}},C_{u_{v}},E_{u_{v}},B_{d_{v}},C_{d_{v}},C_{\bar{u}},D_{\bar{u}},A_{\bar{d}},B_{\bar{d}},C_{\bar{d}}\\}\,.$
(16)

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\gtrsim 3.0$ GeV-1, the non-
linear modifications were found to be negligible.

We find the optimal set of parameters by minimizing the $\chi^{2}$ function
defined as

$\displaystyle\chi^{2}=\sum_{i}\frac{[T_{i}-\sum_{\alpha}\gamma_{\alpha}^{i}D_{i}b_{\alpha}-D_{i}]^{2}}{\left(\delta_{i,\mathrm{stat}}\sqrt{D_{i}T_{i}}\right)^{2}+\left(\delta_{i,\mathrm{uncorr}}T_{i}\right)^{2}}+\sum_{\alpha}b_{\alpha}^{2}.$
(17)

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 $\delta_{i,\mathrm{stat}},\delta_{i,\mathrm{uncorr}}$ and
$\gamma_{\alpha}^{i}$ are the statistical, uncorrelated and correlated
uncertainties associated with each data point. Minimizing this $\chi^{2}$
function with respect to the 14 fit parameters and the nuisance parameters
$b_{\alpha}$ determines the central set of PDFs and systematic shifts,
$-\sum_{\alpha}\gamma_{\alpha}^{i}D_{i}b_{\alpha}$, which give the best
description of the data. The uncertainties of the obtained PDFs are then
determined using the Hessian method Pumplin:2001ct , which is based on a
quadratic expansion of the $\chi^{2}$ function around its minimum,
$\chi^{2}=\chi^{2}_{\mathrm{min}}+\Delta\chi^{2}$. The choice of the maximum
displacement $\Delta\chi^{2}$, the tolerance, depends on the definition of
uncertainties and in PDF fits it is usually conservatively chosen to be
$\Delta\chi^{2}\gg 1$. Indeed, assuming that the data points are Gaussianly
distributed around the ”truth”, the expected $\chi^{2}/N_{dof}$ should be less
than $\sim$$1.05$ in 90% of the cases with $\sim$$1600$ data points, which is
the amount of points we are now considering. However, the minimum values we
find are $\chi^{2}/N_{dof}\sim$$1.18$, i.e., the ideal choice
$\Delta\chi^{2}=1$ is, perhaps, an underestimate. Assuming a Gaussian
probability density for the fit parameters results in
$\Delta\chi^{2}_{\mathrm{max}}\approx 20$ at the 90% confidence level
Eskola:2009uj . In practice, we determine the PDFs with the tolerance set to
$\Delta\chi^{2}=1$ and then rescale the resulting uncertainties by
$\sqrt{20}\approx 4.5$ to match the actual tolerance since this leads to a
better convergence of the xFitter algorithms.

###  III.2 Data selection

Our analysis, as the ones presented in Refs. Walt:2019slu  and
Helenius:2021tof , includes the lepton-proton DIS data from the BCDMS, HERA
and NMC experiments, which are summarized in Table 1. The data from BCDMS and
NMC are given in terms of the structure function $F_{2}^{p}(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}$
H1:2013ktq  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 H1:2009uwa ; H1:2012xnw ; ZEUS:2014wft .

In the first round of fits, we apply the same $Q^{2}>3.5$ GeV2 cut as, e.g.,
in the HERAPDF2.0 analysis, but with parameterization 1 we also perform a
second round of fits with the cut lowered to $Q^{2}>1.0$ GeV2 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, $\chi^{2}/N_{dof}$ tends to get systematically
worse as the minimum $Q^{2}$ is lowered from $15\,{\rm GeV}^{2}$ to $2\,{\rm
GeV}^{2}$ Ball:2017otu ; xFitterDevelopersTeam:2018hym  – a tension that the
BFKL resummation appears to ease. The two rightmost columns of Table 1
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.

Table 1: 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$ GeV2 and $Q^{2}>1$ GeV2 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. | $N_{\text{data}}^{Q^{2}>3.5\,\mathrm{GeV^{2}}}$ | $N_{\text{data}}^{Q^{2}>1.0\,\mathrm{GeV^{2}}}$
---|---|---|---|---|---
BCDMS | NC $\mu$ $F_{2}^{p}$ 100 GeV | 1996 | BCDMS:1989qop  | 83 | 83
NC $\mu$ $F_{2}^{p}$ 120 GeV | 91 | 91
NC $\mu$ $F_{2}^{p}$ 200 GeV | 79 | 79
NC $\mu$ $F_{2}^{p}$ 280 GeV | 75 | 75
HERA | NC e+ 920 GeV | 2015 | H1:2015ubc  | 377 | 425
NC e+ 820 GeV | 70 | 78
NC e+ 575 GeV | 254 | 260
NC e+ 460 GeV | 204 | 210
NC e- 920 GeV | 159 | 159
CC e+ 920 GeV | 39 | 39
CC e- 920 GeV | 42 | 42
NMC | NC $\mu$ $F_{2}^{p}$ | 1997 | NewMuon:1996fwh  | 95 | 95
Total | 1568 | 1636

###  III.3 Fit results

Figure 3: Minimum $\chi^{2}/N_{dof}$ as a function of $R$ for parameterization
1 with $Q^{2}>3.5$ GeV2 (orange) and $Q^{2}>1$ GeV2 (blue), and
parameterization 2 with $Q^{2}>3.5$ GeV2 (green).

The resulting values of $\chi^{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 $\chi^{2}/N_{dof}$ are presented in Appendix A, Tables
2, 3 and 4. In Fig. 3, the orange (parameterization 1) and green
(parameterization 2) crosses correspond to the $Q^{2}>3.5$ GeV2 cut and the
blue ones (parameterization 1) are for the fits with $Q^{2}>1.0$ GeV2. 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 $\chi^{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 $\approx 0.7$ GeV-1 causes only a small increase, but
below that $\chi^{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 $\chi^{2}/N_{dof}$ upwards by
approximately 0.04 units across all values of $R$. While this increase in
$\chi^{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 $\chi^{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 2, 3 and 4 in Appendix A break down the obtained $\chi^{2}$ into
contributions from different data sets. Tables 2 and 4, which correspond to
the fits with parameterization 1 at $Q^{2}_{\mathrm{cut}}=3.5$ GeV2 and 1
GeV2, 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 $\chi^{2}/N_{dof}$. Table 3 shows that with
parametrization 2, the $\chi^{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$ GeV2 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\,{\rm
GeV},10\,{\rm GeV},100\,{\rm GeV}\\}$ for parameterization 1 and 2,
respectively. They are compared to the HERAPDF2.0 results given by the black
dashed curve. To avoid visual clutter, the uncertainties (90% confidence
level) are only shown for HERAPDF2.0 and $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.0 reference 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.0 at 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\geq 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}_{\mathrm{cut}}>3.5$ GeV2
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 5, 6 and 7 in Appendix B.

Figure 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$.  Figure 5: The
same as Fig. 4, but showing PDFs fitted with parameterization 2.

##  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 $\sigma_{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}$
Lappi:2023lmi , particularly at small $x$. In this sense the discussion below
is not entirely independent of the fits presented above.

###  IV.1 The longitudinal structure function at HERA

Figure 6: Structure function $F_{L}$ calculated with PDFs resulting from
different choices of $R$ compared with the measurements by H1 H1:2013ktq  and
ZEUS ZEUS:2014thn . We only show predictions for $Q^{2}>Q_{0}^{2}$.  Figure 7:
Structure function $F_{L}$ resulting from PDFs fitted with parameterization 2
and different $R$ values compared with the HERA $F_{L}$ measurements
H1:2013ktq . We only show predictions for $Q^{2}>Q_{0}^{2}$.

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
H1:2013ktq  and ZEUS ZEUS:2014thn . 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}_{\mathrm{cut}}=3.5$ GeV2 with parameterization 1 and the right panel
with parameterization 2. For $Q^{2}\gtrsim 100$ GeV2 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}\rightarrow 0$ at the PDFs parameterization scale $Q_{0}^{2}=1.69$ GeV2.
The predictions for parameterization 1 do not converge towards
$F_{L}\rightarrow 0$ since the parameterization allows the PDFs to be negative
and therefore yield negative $F_{L}$ values for $Q^{2}\lesssim 5$ GeV2. At
$Q^{2}\lesssim 10$ GeV2, the NNLO predictions generally lie below the values
for $F_{L}$ measured by H1, as has been observed in other studies Ball:2017otu
; xFitterDevelopersTeam:2018hym . 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}\leq 10$ GeV2 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\gtrsim 0.3\times 10^{-2}$, the predictions with different
$R$ are almost indistinguishable. However for $Q^{2}>5$ GeV2 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\approx 10^{-5}$ until the highest considered $Q^{2}$. We
note that at $Q^{2}>8.5$ GeV2 the $F_{L}$ data are well described by most of
the predictions. However, with $R\leq 0.4$ GeV-1 the description begins to
visibly deteriorate even at higher values of $Q^{2}$.

###  IV.2 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 Accardi:2012qut  and LHeC LHeC:2020van ;
LHeCStudyGroup:2012zhm  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:

$\displaystyle
R_{F_{L}}(x,Q^{2})\equiv\frac{F_{L}(x,Q^{2})}{F_{L}^{R=3\mathrm{GeV}^{-1}}(x,Q^{2})}.$
(18)

These ratios are compared to the projected relative uncertainties for $F_{L}$
taken from the EIC development documents Aschenauer:9999abc  and the LHeC
White Paper LHeC:2020van  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. Ball:2017otu ;
xFitterDevelopersTeam:2018hym . The projected EIC data only reach the $x$
region, where the different $R$ values are clearly distinguishable at $Q<2.5$
GeV2, but even at $Q=2.47$ GeV2 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.

Figure 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
Aschenauer:9999abc . We only show predictions for $Q^{2}>Q_{0}^{2}$. Figure 9:
The same as Fig. 8, but for LHeC LHeC:2020van .

##  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, $\chi^{2}/N_{dof}\sim 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\lesssim 0.7$ GeV-1, but for the positive-definite
parametrization at somewhat lower values, $R\lesssim 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\lesssim 2.5\,{\rm GeV}$. We also studied the
sensitivity of our fits to the choice of $Q^{2}_{\mathrm{cut}}$: lowering
$Q^{2}_{\mathrm{cut}}$ from $3.50\,{\rm GeV}^{2}$ to $1\,{\rm GeV}^{2}$ leads
to an increase of $\chi^{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\gtrsim 0.4$ GeV-1 or
$Q_{r}=1/R\lesssim 2.5\,{\rm GeV}$. Our conclusions about HERA data not
supporting non-linear effects is in line with a recent study made in Ref.
Mantysaari:2018nng . 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\gtrsim 100\,\mathrm{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={\cal O}(100\,\rm 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
$\mathcal{O}(\alpha_{s}^{2})$, 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 Buckley:2014ana  for the
fitted PDFs are available here: https://research.hip.fi/qcdtheory/nuclear-
pdfs/.

## Acknowledgements

Our work has been supported by the Academy of Finland, projects 331545 and
330448, and was funded as a part of the Center of Excellence in Quark Matter
of the Academy of Finland, project 346326. This research is part of the
European Research Council project ERC-2018-ADG-835105 YoctoLHC. We acknowledge
grants of computer capacity from the Finnish Grid and Cloud Infrastructure
(persistent identifier urn:nbn:fi:research-infras-2016072533).

## References

  * (1) J. C. Collins, D. E. Soper, and G. F. Sterman, “Factorization of Hard Processes in QCD,” Adv. Ser. Direct. High Energy Phys. 5 (1989) 1–91, hep-ph/0409313.
  * (2) Y. L. Dokshitzer, “Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics.,” Sov. Phys. JETP 46 (1977) 641–653.
  * (3) V. N. Gribov and L. N. Lipatov, “Deep inelastic e p scattering in perturbation theory,” Sov. J. Nucl. Phys. 15 (1972) 438–450.
  * (4) V. N. Gribov and L. N. Lipatov, “e+ e- pair annihilation and deep inelastic e p scattering in perturbation theory,” Sov. J. Nucl. Phys. 15 (1972) 675–684.
  * (5) G. Altarelli and G. Parisi, “Asymptotic Freedom in Parton Language,” Nucl. Phys. B 126 (1977) 298–318.
  * (6) J.-w. Qiu, “Nuclear Shadowing at Small Values of x,” Nucl. Phys. B 291 (1987) 746–764.
  * (7) Y. V. Kovchegov and E. Levin, Quantum Chromodynamics at High Energy, vol. 33.  Oxford University Press, 2013.
  * (8) A. H. Mueller and J.-w. Qiu, “Gluon Recombination and Shadowing at Small Values of x,” Nucl. Phys. B 268 (1986) 427–452.
  * (9) L. V. Gribov, E. M. Levin, and M. G. Ryskin, “Semihard Processes in QCD,” Phys. Rept. 100 (1983) 1–150.
  * (10) J. Bartels, G. A. Schuler, and J. Blumlein, “A Numerical study of the small x behavior of deep inelastic structure functions in QCD,” Z. Phys. C 50 (1991) 91–102.
  * (11) K. Prytz, “Signals of gluon recombination in deep inelastic scattering,” Eur. Phys. J. C 22 (2001) 317–321, hep-ph/0112095.
  * (12) K. J. Eskola, H. Honkanen, V. J. Kolhinen, J.-w. Qiu, and C. A. Salgado, “Nonlinear corrections to the DGLAP equations in view of the hera data,” Nucl. Phys. B 660 (2003) 211–224, hep-ph/0211239.
  * (13) G. R. Boroun and S. Zarrin, “An approximate approach to the nonlinear DGLAP evaluation equation,” Eur. Phys. J. Plus 128 (2013) 119, 1402.0678.
  * (14) G. R. Boroun, “Shadowing correction to the gluon distribution behavior at small x,” Eur. Phys. J. A 42 (2009) 251–255, 1402.0864.
  * (15) J. Rausch, V. Guzey, and M. Klasen, “Numerical evaluation of the nonlinear Gribov-Levin-Ryskin-Mueller-Qiu evolution equations for nuclear parton distribution functions,” Phys. Rev. D 107 (2023), no. 5, 054003, 2211.14094.
  * (16) W. Zhu, “A New approach to parton recombination in a QCD evolution equation,” Nucl. Phys. B 551 (1999) 245–274, hep-ph/9809391.
  * (17) W. Zhu and J.-h. Ruan, “A New modified Altarelli-Parisi evolution equation with parton recombination in proton,” Nucl. Phys. B 559 (1999) 378–392, hep-ph/9907330.
  * (18) J. Blumlein, V. Ravindran, J.-h. Ruan, and W. Zhu, “Twist four gluon recombination corrections for deep inelastic structure functions,” Phys. Lett. B 504 (2001) 235–240, hep-ph/0102025.
  * (19) M. Klasen and H. Paukkunen, “Nuclear PDFs After the First Decade of LHC Data,” 2311.00450.
  * (20) F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, “The Color Glass Condensate,” Ann. Rev. Nucl. Part. Sci. 60 (2010) 463–489, 1002.0333.
  * (21) A. Morreale and F. Salazar, “Mining for Gluon Saturation at Colliders,” Universe 7 (2021), no. 8, 312, 2108.08254.
  * (22) G. P. Salam and J. Rojo, “A Higher Order Perturbative Parton Evolution Toolkit (HOPPET),” Comput. Phys. Commun. 180 (2009) 120–156, 0804.3755.
  * (23) xFitter Collaboration, H. Abdolmaleki et al., “xFitter: An Open Source QCD Analysis Framework. A resource and reference document for the Snowmass study,”  6, 2022.  2206.12465.
  * (24) S. Alekhin et al., “HERAFitter,” Eur. Phys. J. C 75 (2015), no. 7, 304, 1410.4412.
  * (25) BCDMS Collaboration, A. C. Benvenuti et al., “A High Statistics Measurement of the Proton Structure Functions F(2) (x, Q**2) and R from Deep Inelastic Muon Scattering at High Q**2,” Phys. Lett. B 223 (1989) 485–489.
  * (26) H1, ZEUS Collaboration, H. Abramowicz et al., “Combination of measurements of inclusive deep inelastic ${e^{\pm}p}$ scattering cross sections and QCD analysis of HERA data,” Eur. Phys. J. C 75 (2015), no. 12, 580, 1506.06042.
  * (27) New Muon Collaboration, M. Arneodo et al., “Measurement of the proton and deuteron structure functions, F2(p) and F2(d), and of the ratio sigma-L / sigma-T,” Nucl. Phys. B 483 (1997) 3–43, hep-ph/9610231.
  * (28) V. S. Fadin, E. A. Kuraev, and L. N. Lipatov, “On the Pomeranchuk Singularity in Asymptotically Free Theories,” Phys. Lett. B 60 (1975) 50–52.
  * (29) E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, “Multi - Reggeon Processes in the Yang-Mills Theory,” Sov. Phys. JETP 44 (1976) 443–450.
  * (30) E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, “The Pomeranchuk Singularity in Nonabelian Gauge Theories,” Sov. Phys. JETP 45 (1977) 199–204.
  * (31) I. I. Balitsky and L. N. Lipatov, “The Pomeranchuk Singularity in Quantum Chromodynamics,” Sov. J. Nucl. Phys. 28 (1978) 822–829.
  * (32) M. Bonvini, S. Marzani, and C. Muselli, “Towards parton distribution functions with small-$x$ resummation: HELL 2.0,” JHEP 12 (2017) 117, 1708.07510.
  * (33) R. D. Ball, V. Bertone, M. Bonvini, S. Marzani, J. Rojo, and L. Rottoli, “Parton distributions with small-x resummation: evidence for BFKL dynamics in HERA data,” Eur. Phys. J. C 78 (2018), no. 4, 321, 1710.05935.
  * (34) xFitter Developers’ Team Collaboration, H. Abdolmaleki et al., “Impact of low-$x$ resummation on QCD analysis of HERA data,” Eur. Phys. J. C 78 (2018), no. 8, 621, 1802.00064.
  * (35) K. J. Eskola, J.-w. Qiu, and X.-N. Wang, “Perturbative gluon shadowing in heavy nuclei,” Phys. Rev. Lett. 72 (1994) 36–39, nucl-th/9307025.
  * (36) K. J. Eskola, H. Honkanen, V. J. Kolhinen, J.-w. Qiu, and C. A. Salgado, “Nonlinear corrections to the DGLAP equations: Looking for the saturation limits,” hep-ph/0302185.
  * (37) A. Accardi et al., “Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all,” Eur. Phys. J. A 52 (2016), no. 9, 268, 1212.1701.
  * (38) R. Abdul Khalek et al., “Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report,” Nucl. Phys. A 1026 (2022) 122447, 2103.05419.
  * (39) LHeC, FCC-he Study Group Collaboration, P. Agostini et al., “The Large Hadron–Electron Collider at the HL-LHC,” J. Phys. G 48 (2021), no. 11, 110501, 2007.14491.
  * (40) LHeC Study Group Collaboration, J. L. Abelleira Fernandez et al., “A Large Hadron Electron Collider at CERN: Report on the Physics and Design Concepts for Machine and Detector,” J. Phys. G 39 (2012) 075001, 1206.2913.
  * (41) FCC Collaboration, A. Abada et al., “FCC Physics Opportunities: Future Circular Collider Conceptual Design Report Volume 1,” Eur. Phys. J. C 79 (2019), no. 6, 474.
  * (42) A. Accardi, L. T. Brady, W. Melnitchouk, J. F. Owens, and N. Sato, “Constraints on large-$x$ parton distributions from new weak boson production and deep-inelastic scattering data,” Phys. Rev. D 93 (2016), no. 11, 114017, 1602.03154.
  * (43) M. Botje, “QCDNUM: Fast QCD Evolution and Convolution,” Comput. Phys. Commun. 182 (2011) 490–532, 1005.1481.
  * (44) M. Botje, “Erratum for the time-like evolution in QCDNUM,” 1602.08383.
  * (45) R. S. Thorne, “A Variable-flavor number scheme for NNLO,” Phys. Rev. D 73 (2006) 054019, hep-ph/0601245.
  * (46) J. Pumplin, D. Stump, R. Brock, D. Casey, J. Huston, J. Kalk, H. L. Lai, and W. K. Tung, “Uncertainties of predictions from parton distribution functions. 2. The Hessian method,” Phys. Rev. D 65 (2001) 014013, hep-ph/0101032.
  * (47) K. J. Eskola, H. Paukkunen, and C. A. Salgado, “EPS09: A New Generation of NLO and LO Nuclear Parton Distribution Functions,” JHEP 04 (2009) 065, 0902.4154.
  * (48) M. Walt, I. Helenius, and W. Vogelsang, “Open-source QCD analysis of nuclear parton distribution functions at NLO and NNLO,” Phys. Rev. D 100 (2019), no. 9, 096015, 1908.03355.
  * (49) I. Helenius, M. Walt, and W. Vogelsang, “NNLO nuclear parton distribution functions with electroweak-boson production data from the LHC,” Phys. Rev. D 105 (2022), no. 9, 094031, 2112.11904.
  * (50) H1 Collaboration, V. Andreev et al., “Measurement of inclusive $ep$ cross sections at high $Q^{2}$ at $\sqrt{s}=$ 225 and 252 GeV and of the longitudinal proton structure function $F_{L}$ at HERA,” Eur. Phys. J. C 74 (2014), no. 4, 2814, 1312.4821.
  * (51) H1 Collaboration, F. D. Aaron et al., “Measurement of the Charm and Beauty Structure Functions using the H1 Vertex Detector at HERA,” Eur. Phys. J. C 65 (2010) 89–109, 0907.2643.
  * (52) H1, ZEUS Collaboration, H. Abramowicz et al., “Combination and QCD Analysis of Charm Production Cross Section Measurements in Deep-Inelastic ep Scattering at HERA,” Eur. Phys. J. C 73 (2013), no. 2, 2311, 1211.1182.
  * (53) ZEUS Collaboration, H. Abramowicz et al., “Measurement of beauty and charm production in deep inelastic scattering at HERA and measurement of the beauty-quark mass,” JHEP 09 (2014) 127, 1405.6915.
  * (54) T. Lappi, H. Mäntysaari, H. Paukkunen, and M. Tevio, “Evolution of structure functions in momentum space,” 2304.06998.
  * (55) ZEUS Collaboration, H. Abramowicz et al., “Deep inelastic cross-section measurements at large y with the ZEUS detector at HERA,” Phys. Rev. D 90 (2014), no. 7, 072002, 1404.6376.
  * (56) E. Aschenauer et al., “The longitudinal structure function $f_{l}$ at an eic.” [https://www.phenix.bnl.gov/WWW/publish/elke/EIC/EIC-R&D-Tracking/Meetings/fl.pdf; accessed 03-May-2023].
  * (57) H. Mäntysaari and P. Zurita, “In depth analysis of the combined HERA data in the dipole models with and without saturation,” Phys. Rev. D 98 (2018) 036002, 1804.05311.
  * (58) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, “LHAPDF6: parton density access in the LHC precision era,” Eur. Phys. J. C 75 (2015) 132, 1412.7420.

##  Appendix A Obtained $\chi^{2}$ values for each data set

The following three Tables 2-4 show the values of $\chi^{2}/N_{dof}$ obtained
in the different sets of PDF fits for a selection of $R$ values
$R=\\{0.4,0.6,0.8,1.0,2.0,3.0\\}$ GeV-1. Table 2 shows the values obtained in
the fits with parameterization 1 and $Q^{2}_{\mathrm{cut}}=3.5$ GeV2. Table 4
shows the same for $Q^{2}_{\mathrm{cut}}=1.0$ GeV2 and Table 3 for
parameterization 2 with $Q^{2}_{\mathrm{cut}}=3.5$ GeV2.

Table 2: Values of $\chi^{2}/N_{dof}$ obtained in the fits with parameterization 1 and $Q^{2}_{\mathrm{cut}}=3.5$ GeV2 for various values for $R$ for each individual data set. The final row lists the total $\chi^{2}/N_{dof}$ for each $R$. $R$ [GeV-1] | 0.4 | 0.6 | 0.8 | 1.0 | 2.0 | 3.0
---|---|---|---|---|---|---
BCDMS | NC $\mu$ $F_{2}^{p}$ 100 GeV | 1.257 | 1.273 | 1.236 | 1.229 | 1.202 | 1.197
NC $\mu$ $F_{2}^{p}$ 120 GeV | 0.770 | 0.766 | 0.753 | 0.754 | 0.748 | 0.746
NC $\mu$ $F_{2}^{p}$ 200 GeV | 1.118 | 1.106 | 1.091 | 1.093 | 1.087 | 1.086
NC $\mu$ $F_{2}^{p}$ 280 GeV | 0.949 | 0.947 | 0.946 | 0.951 | 0.948 | 0.948
HERA | NC e+ 920 GeV | 1.324 | 1.248 | 1.218 | 1.210 | 1.200 | 1.200
NC e+ 820 GeV | 1.088 | 1.040 | 1.029 | 1.008 | 0.991 | 0.988
NC e+ 575 GeV | 0.924 | 0.895 | 0.883 | 0.879 | 0.875 | 0.874
NC e+ 460 GeV | 1.137 | 1.109 | 1.094 | 1.095 | 1.095 | 1.094
NC e- 920 GeV | 1.495 | 1.469 | 1.457 | 1.458 | 1.454 | 1.455
CC e+ 920 GeV | 1.330 | 1.421 | 1.349 | 1.412 | 1.373 | 1.365
CC e- 920 GeV | 1.505 | 1.434 | 1.501 | 1.505 | 1.494 | 1.491
NMC | NC $\mu$ $F_{2}^{p}$ | 1.179 | 1.038 | 0.975 | 0.972 | 0.972 | 0.972
$\chi^{2}/N_{dof}$ | 1.244 | 1.205 | 1.194 | 1.188 | 1.182 | 1.182
Table 3: The same as Table 2, but with parameterization 2. $R$ [GeV-1] | 0.4 | 0.6 | 0.8 | 1.0 | 2.0 | 3.0
---|---|---|---|---|---|---
BCDMS | NC $\mu$ $F_{2}^{p}$ 100 GeV | 1.045 | 1.006 | 1.012 | 1.013 | 1.016 | 1.017
NC $\mu$ $F_{2}^{p}$ 120 GeV | 0.755 | 0.749 | 0.749 | 0.750 | 0.751 | 0.752
NC $\mu$ $F_{2}^{p}$ 200 GeV | 1.090 | 1.083 | 1.082 | 1.082 | 1.082 | 1.082
NC $\mu$ $F_{2}^{p}$ 280 GeV | 0.934 | 0.934 | 0.935 | 0.935 | 0.937 | 0.937
HERA | NC e+ 920 GeV | 1.211 | 1.188 | 1.188 | 1.191 | 1.197 | 1.198
NC e+ 820 GeV | 1.066 | 1.041 | 1.029 | 1.024 | 1.021 | 1.021
NC e+ 575 GeV | 0.883 | 0.878 | 0.877 | 0.877 | 0.878 | 0.878
NC e+ 460 GeV | 1.074 | 1.066 | 1.068 | 1.070 | 1.072 | 1.073
NC e- 920 GeV | 1.455 | 1.455 | 1.456 | 1.456 | 1.457 | 1.457
CC e+ 920 GeV | 1.006 | 0.952 | 0.946 | 0.944 | 0.942 | 0.941
CC e- 920 GeV | 1.218 | 1.263 | 1.262 | 1.261 | 1.262 | 1.262
NMC | NC $\mu$ $F_{2}^{p}$ | 1.092 | 1.082 | 1.081 | 1.083 | 1.086 | 1.087
$\chi^{2}/N_{dof}$ | 1.206 | 1.187 | 1.187 | 1.188 | 1.191 | 1.192
Table 4: The same as Table 2, but with the lower cut $Q^{2}_{\mathrm{cut}}=1.0$ GeV2.  $R$ [GeV-1] | 0.4 | 0.6 | 0.8 | 1.0 | 2.0 | 3.0
---|---|---|---|---|---|---
BCDMS | NC $\mu$ $F_{2}^{p}$ 100 GeV | 1.272 | 1.262 | 1.227 | 1.206 | 1.175 | 1.167
NC $\mu$ $F_{2}^{p}$ 120 GeV | 0.797 | 0.762 | 0.755 | 0.751 | 0.747 | 0.746
NC $\mu$ $F_{2}^{p}$ 200 GeV | 1.133 | 1.099 | 1.092 | 1.089 | 1.087 | 1.087
NC $\mu$ $F_{2}^{p}$ 280 GeV | 0.955 | 0.949 | 0.952 | 0.952 | 0.952 | 0.952
HERA | NC e+ 920 GeV | 1.778 | 1.331 | 1.256 | 1.231 | 1.213 | 1.212
NC e+ 820 GeV | 1.816 | 1.218 | 1.213 | 1.196 | 1.179 | 1.177
NC e+ 575 GeV | 1.017 | 0.919 | 0.907 | 0.901 | 0.898 | 0.898
NC e+ 460 GeV | 1.229 | 1.111 | 1.078 | 1.081 | 1.076 | 1.075
NC e- 920 GeV | 1.560 | 1.495 | 1.472 | 1.466 | 1.463 | 1.463
CC e+ 920 GeV | 1.420 | 1.332 | 1.362 | 1.361 | 1.363 | 1.361
CC e- 920 GeV | 2.039 | 1.598 | 1.504 | 1.494 | 1.484 | 1.480
NMC | NC $\mu$ $F_{2}^{p}$ | 1.648 | 1.068 | 1.008 | 1.008 | 1.023 | 1.028
$\chi^{2}/N_{dof}$ | 1.507 | 1.259 | 1.226 | 1.218 | 1.213 | 1.212

##  Appendix B Obtained PDF parameters

The following three Tables 5-7 show the parameters obtained in the different
sets of PDF fits for a selection of $R$ values
$R=\\{0.4,0.6,0.8,1.0,2.0,3.0\\}$ GeV-1. Table 5 shows the parameters obtained
in the fits with parameterization 1 and $Q^{2}_{\mathrm{cut}}=3.5$ GeV2. Table
6 shows the same for $Q^{2}_{\mathrm{cut}}=1.0$ GeV2 and Table 7 for
parameterization 2 with $Q^{2}_{\mathrm{cut}}=3.5$ GeV2.

Table 5: Parameters obtained in the fits with parameterization 1 and $Q^{2}_{\mathrm{cut}}=3.5$ GeV2 at various values for $R$. Parameters in parentheses are fixed at the given value, while parameters not shown are set to zero. $R$ [GeV-1] | 0.4 | 0.6 | 0.8 | 1.0 | 2.0 | 3.0
---|---|---|---|---|---|---
$A_{\bar{d}}$ | 0.404469 | 0.161368 | 0.104539 | 0.061467 | 0.056025 | 0.054686
$A^{\prime}_{g}$ | 2.218085 | 0.435239 | 0.410244 | 0.201921 | 0.494136 | 0.428424
$B_{\bar{d}}$ | 0.010765 | -0.085343 | -0.108657 | -0.128949 | -0.144103 | -0.145589
$B_{d_{v}}$ | 1.422865 | 1.454721 | 1.363376 | 1.323526 | 1.314441 | 1.313151
$B_{g}$ | -0.211063 | -0.342516 | -0.174697 | -0.042353 | -0.123857 | -0.069253
$B^{\prime}_{g}$ | 0.256574 | -0.449499 | -0.309159 | -0.302806 | -0.260254 | -0.235
$B_{u_{v}}$ | 0.914925 | 0.921869 | 0.934284 | 0.921569 | 0.923402 | 0.923452
$C_{\bar{d}}$ | 18.358827 | 13.091428 | 6.798707 | 2.759383 | 2.15534 | 2.020401
$C_{d_{v}}$ | 5.876824 | 5.987926 | 5.560861 | 5.694443 | 5.733379 | 5.735318
$C_{g}$ | 2.208999 | 2.120541 | 3.683178 | 4.354789 | 4.851641 | 5.230455
$C^{\prime}_{g}$ | (25.0) | (25.0) | (25.0) | (25.0) | (25.0) | (25.0)
$C_{\bar{u}}$ | 15.873966 | 17.198562 | 16.234897 | 16.742077 | 15.787451 | 15.521866
$C_{u_{v}}$ | 2.48766 | 2.496799 | 2.529641 | 2.489838 | 2.489735 | 2.487229
$D_{\bar{u}}$ | 16.866158 | 28.489481 | 22.125291 | 23.503082 | 22.000836 | 21.205419
$E_{u_{v}}$ | -0.997742 | -0.988268 | -0.971517 | -1.003805 | -1.020271 | -1.027135
Table 6: The same as Table 5, but for PDFs fitted with the lower cut $Q^{2}_{\mathrm{cut}}=1.0$ GeV2. $R$ [GeV-1] | 0.4 | 0.6 | 0.8 | 1.0 | 2.0 | 3.0
---|---|---|---|---|---|---
$A_{\bar{d}}$ | 0.150472 | 0.105445 | 0.094215 | 0.082058 | 0.06272 | 0.061345
$A^{\prime}_{g}$ | -0.481786 | 0.104096 | 0.142766 | 0.234982 | 0.24655 | 0.268521
$B_{\bar{d}}$ | -0.089444 | -0.091097 | -0.100384 | -0.104039 | -0.109834 | -0.110594
$B_{d_{v}}$ | 1.026321 | 1.222286 | 1.314041 | 1.297563 | 1.266633 | 1.265011
$B_{g}$ | 0.112897 | -0.090972 | 0.012384 | 0.019837 | 0.103054 | 0.1072
$B^{\prime}_{g}$ | -0.148734 | -0.378519 | -0.278353 | -0.215802 | -0.153492 | -0.140128
$B_{u_{v}}$ | 0.872954 | 0.913072 | 0.927043 | 0.924927 | 0.919602 | 0.91878
$C_{\bar{d}}$ | 4.207006 | 5.524392 | 6.20659 | 5.084081 | 2.925174 | 2.766838
$C_{d_{v}}$ | 4.871045 | 5.154293 | 5.408388 | 5.359225 | 5.395968 | 5.409781
$C_{g}$ | 2.001991 | 2.939454 | 4.006964 | 4.490556 | 5.487986 | 5.642055
$C^{\prime}_{g}$ | (25.0) | (25.0) | (25.0) | (25.0) | (25.0) | (25.0)
$C_{\bar{u}}$ | 22.296374 | 16.795233 | 15.487592 | 14.595522 | 13.624455 | 13.376877
$C_{u_{v}}$ | 2.365049 | 2.476635 | 2.499119 | 2.481032 | 2.461876 | 2.458407
$D_{\bar{u}}$ | 30.234495 | 17.061131 | 15.205939 | 12.889824 | 10.462935 | 9.919356
$E_{u_{v}}$ | -1.056546 | -1.003855 | -0.993606 | -1.007731 | -1.02491 | -1.028614
Table 7: The same as Table 5, but for PDFs fitted with parameterization 2. $R$ [GeV-1] | 0.4 | 0.6 | 0.8 | 1.0 | 2.0 | 3.0
---|---|---|---|---|---|---
$A_{\bar{d}}$ | 0.248371 | 0.238772 | 0.229188 | 0.226517 | 0.224377 | 0.224177
$B_{\bar{d}}$ | -0.095279 | -0.108632 | -0.115726 | -0.118081 | -0.120341 | -0.120637
$B_{d_{v}}$ | 0.989654 | 1.05047 | 1.050583 | 1.053705 | 1.057556 | 1.058181
$B_{g}$ | -0.063721 | 0.176504 | 0.300771 | 0.368229 | 0.46931 | 0.489468
$B_{u_{v}}$ | 0.330037 | 0.10342 | 0.10933 | 0.108711 | 0.108755 | 0.108843
$C_{\bar{d}}$ | 9.978846 | 10.592777 | 10.523469 | 10.595957 | 10.745633 | 10.781037
$C_{d_{v}}$ | 4.363475 | 4.575691 | 4.5757 | 4.591326 | 4.61234 | 4.61603
$C_{g}$ | 19.131659 | 21.235982 | 21.791519 | 21.993852 | 22.228964 | 22.267116
$C_{\bar{u}}$ | 8.828193 | 9.34704 | 9.525786 | 9.554065 | 9.563123 | 9.560328
$C_{u_{v}}$ | 3.24987 | 2.501855 | 3.573471 | 3.572397 | 3.568641 | 3.567745
$D_{g}$ | -9.671812 | -8.386368 | -7.830814 | -7.681217 | -7.689173 | -7.720356
$D_{\bar{u}}$ | 9.759016 | 10.223122 | 11.255959 | 11.466205 | 11.578213 | 11.574223
$D_{u_{v}}$ | 31.344142 | 232.987381 | 215.930099 | 217.666728 | 217.330601 | 216.978654
$E_{g}$ | 465.912581 | 399.008994 | 335.978167 | 302.820564 | 257.896509 | 249.69638
$E_{u_{v}}$ | -20.809607 | -238.165353 | 11.060703 | 11.049585 | 9.585398 | 9.200095


# Segmentation of Maya hieroglyphs through fine-tuned foundation models

FNU Shivam1, Megan Leight2, Mary Kate Kelly3, Claire Davis2, Kelsey
Clodfelter2, Jacob Thrasher1,
Yenumula Reddy1, Prashnna Gyawali1
1 Lane Department of Computer Science and Electrical Engineering, West
Virginia University, USA
2 Department of Art History, School of Art and Design, West Virginia
University, USA
3 Mount Royal University, Canada
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

The study of Maya hieroglyphic writing unlocks the rich history of cultural
and societal knowledge embedded within this ancient civilization’s visual
narrative. Artificial Intelligence (AI) offers a novel lens through which we
can translate these inscriptions, with the potential to allow non-specialists
access to reading these texts and to aid in the decipherment of those
hieroglyphs which continue to elude comprehensive interpretation. Toward this,
we leverage a foundational model to segment Maya hieroglyphs from an open-
source digital library dedicated to Maya artifacts. Despite the initial
promise of publicly available foundational segmentation models, their
effectiveness in accurately segmenting Maya hieroglyphs was initially limited.
Addressing this challenge, our study involved the meticulous curation of image
and label pairs with the assistance of experts in Maya art and history,
enabling the fine-tuning of these foundational models. This process
significantly enhanced model performance, illustrating the potential of fine-
tuning approaches and the value of our expanding dataset. We plan to open-
source this dataset for encouraging future research, and eventually to help
make the hieroglyphic texts legible to a broader community, particularly for
Maya heritage community members.

##  1 Introduction

Figure 1: Examples of Maya vessels from Kerr’s Maya Vase Database (K505 and
K3433). A blue box highlights the shape of the glyph block, although K505 is
pseudoglyphic. The segmentation of these glyph blocks, to their individual
hieroglyphic characters or symbols, is the primary focus of this study.

The ancient Maya civilization, one of the most sophisticated and vibrant pre-
Columbian societies in the Americas, has captivated scholars and enthusiasts
alike with its rich tapestry of cultural achievements, particularly in art and
writing [4]. Renowned for their advanced knowledge in mathematics, astronomy,
and architecture, the Maya archaeological record also includes a vast array of
intricate art with hieroglyphic inscriptions that hold additional information
to understanding their complex societal structures, history, rituals, and
cosmology. Texts found on temple walls, stone monuments, pottery, and
manuscripts all serve as a visual narrative of the ancient Maya people,
offering invaluable insights into their daily lives, religious beliefs, and
historical events. However, despite their significance, the intricate nature
of the writing system is not easily parsed and presents a formidable challenge
to traditional methods of study; it also serves as a barrier to interested
non-specialists, and the modern Maya heritage community who are stakeholders
in this historical knowledge.

Though the artifact record is robust, it is not nearly to the scale necessary
for comprehensive decipherment of all known hieroglyphs. Many of them appear
too infrequently in the corpus to have enough contexts to securely identify
them. This is where artificial intelligence offers a groundbreaking
opportunity in understanding Maya art, as demonstrated by previous research
[2, 5]. By leveraging AI and machine learning (ML) technologies, researchers
can develop frameworks for the segmentation and interpretation of Maya
hieroglyphs, thereby allowing for more comprehensive, thorough analysis of the
corpus. Importantly, none of the existing AI frameworks are able to recognize
nor segment the Maya glyph character blocks at the center of our study. Thus,
this technological approach offers a promising solution to bridge the gap
between the expert knowledge restricted to a small number of specialists and
the interested public who could gain access to reading texts from this ancient
culture.

Current state-of-the-art AI frameworks like deep learning predominantly rely
on extensive datasets, often in the magnitude of millions of well-annotated
examples [13]. Such annotation is typically achieved through crowd-sourcing
platforms, such as Mechanical Turk, where numerous individuals contribute to
the labor-intensive process of tagging and categorizing data [14]. This
approach facilitates the development of AI systems capable of recognizing
patterns and making informed predictions or analyses; however, the application
of this methodology to the domain of ancient writing encounters significant
obstacles.

The specialized nature of Maya hieroglyphs in ancient art means that
accurately annotating these images requires a depth of expertise in the
subject matter, a criterion that cannot be easily met through general crowd-
sourcing methods. Experts in Maya epigraphy are relatively scarce, and the
nuanced understanding required to identify and categorize the myriad of
logograms (word signs) and syllabograms (sound signs) is not something that
can be distilled into a simple task for the average person. Moreover, unlike
the fields where AI has seen extensive application, the dataset of Maya
artifacts is not vast but rather limited. Many of these artifacts are unique,
with only a finite number available for study, further complicating the task
of assembling a comprehensive and representative dataset for training AI
models. This scarcity of data, coupled with the need for expert-level
annotation, presents a substantial barrier to leveraging AI for the analysis
and interpretation of Maya inscriptions, underscoring the need for innovative
solutions tailored to the challenges of working with ancient and specialized
datasets.

Foundational models are recent advancements within deep learning, and they are
trained on vast datasets spanning millions of examples across various domains
[3, 8]. These models have demonstrated remarkable capabilities, including
zero-shot learning, which enables them to perform tasks on data they have
never seen before, such as segmenting a wide array of images with precision,
including complex medical images [10]. The success of foundational models in
these areas, especially their ability to generalize and adapt to new
challenges without direct training, has inspired us to explore their potential
in the field of ancient writing.

Toward this end, we first explored the application of foundational models,
specifically the Segment Anything (SAM) model [8], for segmenting glyph blocks
on Maya ceramic vessels. Despite its success in various domains, the
performance of such foundational segmentation models was initially sub-optimal
in this context. Secondly, recognizing the critical importance of expertly
annotated data in addressing the challenges posed by the limited and
specialized nature of Maya hieroglyphs, we enlisted the help of several
individuals with deep expertise in Maya art and epigraphy. These experts
dedicated their efforts to meticulously annotating images sourced from Justin
Kerr’s Maya Vase Database [6], which serves as a photographic repository for
ancient Maya artifacts. In focusing on the vase database, each vessel was
carefully annotated around each glyph block to ensure the highest level of
accuracy. Finally, we fine-tuned the SAM model using our specifically curated
dataset, which led to a demonstrable improvement in performance. Overall, our
contributions can be summarized as follows:

  * •

For the first time, we explored the capabilities of publicly available
foundational models, specifically Segment Anything (SAM), in segmenting
hieroglyphs from Maya vessels.

  * •

Our study involved the manual annotation of hieroglyphs within publicly
available images of Maya art, which we plan to open source to empower other
researchers with a valuable resource for exploring and understanding Maya
hieroglyphs.

  * •

We fine-tuned SAM, demonstrating improved segmentation performance.

##  2 Dataset preparation and Pre-processing

###  2.1 Data Sources

The dataset utilized in our project is derived from the MayaVase database, an
extensive collection of images curated by Justin Kerr from 1966-2012, recently
acquired by Dumbarton Oaks Research Library and Collection [1]. The database
comprises images of vessels and objects housed in museums and private
collections, the majority of which have been gathered from across the United
States, Europe, Mexico, and Central America. The MayaVase database is
particularly notable for its publication of Kerr’s pioneering rollout
photographs of Maya vessels, which allowed cylindrical objects to be viewed as
a single narrative image [7]. In addition to cylindrical vessels, the database
includes images of plates, bowls, archaeological sites, various art objects,
and items from private collections, providing a comprehensive visual catalog
of Maya material culture in the public domain. We have provided two examples
of vessels from the Maya Vase Database in Fig. 1. Note the presence of
pseudoglyphic text on K505, the vessel with the dog peering into a vase. The
text is placed correctly and the blocks are glyph-shaped; however, they do not
bear legible meaning, which confirms the importance of having Maya epigraphers
on these types of AI data research projects.

###  2.2 Data Preparation

Annotation of the images represented a pivotal contribution to our project,
encompassing extensive efforts in data preparation by a dedicated team
comprising faculty and students from art historical and anthropological
disciplines. These team members, who have devoted years to the study of these
objects, played a crucial role in the meticulous process of data curation and
annotation. The process began with a comprehensive web-scraping of the entire
Maya Vase Database to collect our sample images. Kerr, in collaboration with
other specialists, analyzed each vessel and assigned similar iconographic
themes (e.g. bees, canoe, dog) across all the images in the database. We
focused on vessels from approximately 27 of the 88 different iconographic
themes.

Subsequent to the image collection, the team embarked on the manual annotation
of each image, carefully selecting those containing text. Hieroglyphs were
present in each image, resulting in a curated collection of over 117 images,
varying in shape and size, specifically chosen for their presence of glyph
blocks or running hieroglyphic text. Then, more detailed annotation took place
by carefully identifying and labeling each block, a task that required a deep
understanding of Maya iconography and script. To facilitate this process, we
utilized labelme.com, an open-source annotation tool, which allowed for
efficient and accurate annotation work. For each image, polygons were drawn
around each hieroglyph block and given different file names, and a
corresponding JSON file was created, containing all the annotations for a
single image. These JSON files were used in generating ground truth binary
masks for each image. An example of input image and its corresponding binary
mask is presented in Fig. 2.

(a) Input image

(b) Binary mask

Figure 2: Binary mask (right) generated through manual annotation of glyphs
from the original input image (left).

##  3 Methodology

We present a novel approach for segmenting Maya characters using a
foundational segmentation model. Specifically, we consider the Segment
Anything (SAM) [8] model and fine-tuned this model with our dataset.

###  3.1 Preliminary: Segment Anything (SAM)

The Segment Anything (SAM) model is a foundational segmentation network,
comprising an image encoder, a prompt encoder, and a mask decoder [8]. SAM
works by first encoding the image into a high dimensional vector
representation using the image encoder, which is a Masked Autoencoder pre-
trained Vision Transformer. The prompt is encoded into a separate vector
representation by the prompt encoder. These two vector representations are
then combined and passed to a mask decoder, which gives the output specified
by the prompt.

Even though SAM has shown great promise in segmenting objects that it has
never seen before, there are several issues when it is applied to Maya
hieroglyph segmentation. Maya images have too much noise and unwanted
characters, which makes it difficult for SAM to segment the Maya hieroglyphs
from the other elements of the visual scenes. Despite the fact that SAM is
trained on using a large and diverse number of images, they are not specific
to the ancient Maya image datasets, which makes segmentation of Maya
characters more inaccurate than what we expected. So, in this regard, fine-
tuning the SAM on these ancient Maya images is important to move forward.

###  3.2 Finetuned-SAM

We propose fine-tuning SAM to enhance segmentation performance for Maya
hieroglyph recognition. By refining SAM, we aim to adapt it specifically to
the nuances of Maya inscriptions. For fine-tuning, we utilize the pre-trained
weights of the image encoder and prompt encoder, ensuring the model benefits
from prior learning. However, we freeze the parameters of the image and prompt
encoder, calculating gradients only for the mask decoder.

Towards learning the Finetuned-SAM, we use the combination of Dice loss and
Cross-Entropy loss as the objective function. The overall loss can be
represented as:

$L=\alpha\cdot L_{\text{CE}}+\beta\cdot L_{\text{Dice}}$

where $\alpha$ and $\beta$ are the weights that determine the contribution of
each loss component to the total loss. In our experiments, both components are
given equal importance.

Figure 3: Qualitative analysis of the results from SAM and Finetuned-SAM for
selected samples in our test dataset. For both SAM and Finetuned-SAM, we
display the result when two random points are chosen as prompts (denoted by
▲).

##  4 Experiments and Results

###  4.1 Setup

Data: Our experiment utilized a dataset of 117 total images with corresponding
annotations for the Maya hieroglyphs in each image, generated in the form of a
binary mask. Each image and its corresponding ground truth mask was resized to
256 by 256 pixels in order to standardize the input size. This dataset is
partitioned randomly into three subsets: training (75 images), validation (19
images), and test (23 images). All the reported results, including evaluation
metrics, are based on the test split.

Model: For modeling, we utilized the SAM ViT Base version of the SAM model.
For fine-tuning the SAM model (Finetuned-SAM), we trained with the following
hyperparameters: 300 epochs, a learning rate of 1e-3, a batch size of 6, and
used the Adam optimizer. A validation split was employed to identify the best-
performing model during training. We tracked the model’s performance on the
validation set and selected the model demonstrating the highest performance
for testing. The compute time on a NVIDIA RTX 3090 24 GB for fine-tuning was
approximately 2 hours. We assessed our model’s performance using two standard
metrics: Intersection-over-Union (IoU) and DICE score.

###  4.2 Results

We provide the main results of our study in Table 1. As shown, we also
compared our models against standard segmentation approaches like UNet [12]
and autoencoder-based segmentation [11, 9]. Compared to these, the foundation-
based approach demonstrates improved performance. Within the foundation-based
approach, our proposed Finetuned-SAM model outperforms the off-the-shelf SAM
model. Given that these foundation models rely on prompts, we explored various
prompt options in our inference analysis. Additionally, to ensure unbiased
results, we conducted our inference analysis five times with randomly selected
prompts.

Model | IoU | DICE
---|---|---
UNet | 0.179 | 0.289
Autoencoder | 0.174 | 0.277
SAM [1 random point] | 0.118 | 0.199
SAM [2 random point] | 0.115 | 0.196
SAM [0.5 Bbox] | 0.128 | 0.219
SAM [0.75 Bbox] | 0.097 | 0.173
Finetuned-SAM [1 random point] | 0.395 | 0.532
Finetuned-SAM [2 random point] | 0.397 | 0.538
Finetuned-SAM [0.5 Bbox] | 0.327 | 0.471
Finetuned-SAM [0.75 Bbox] | 0.321 | 0.453
Table 1: The comparison of the presented Finetuned-SAM with standard segmentation methods and regular SAM evaluated by IoU and Dice score. For SAM and Finetuned-SAM, we report the average of five random inference runs.  Model | IoU | DICE
---|---|---
SAM [1 random point/block] | 0.316 | 0.438
SAM [2 random points/block] | 0.387 | 0.522
SAM [3 random points/block] | 0.382 | 0.517
Finetuned-SAM [1 random point/block] | 0.417 | 0.555
Finetuned-SAM [2 random points/block] | 0.439 | 0.580
Finetuned-SAM [3 random points/block] | 0.410 | 0.543
Table 2:  Quantitative analysis of Maya glyph segmentation with individualized
input prompts demonstrating improved accuracy.

Further, we analyzed the effect of providing prompts for individual glyph
blocks in our database. We present the results in Table 2, where compared to
Table 1, both SAM and Finetuned-SAM’s performances have improved. This
improvement is anticipated, as providing block-specific prompts offers more
spatial context. Nonetheless, even in this scenario, Finetuned-SAM outperforms
SAM, demonstrating the effectiveness of fine-tuning.

In Figure 3, we provide qualitative results comparing the performance of SAM
and Finetuned-SAM. As observed, Finetuned-SAM outperforms SAM in segmenting
hieroglyphs from Maya images across all examples. Notably, in some cases
(third and fourth row), the SAM model fails significantly, whereas Finetuned-
SAM, though not perfect, still focuses on the correct region.

###  4.3 Conclusion

We present Maya hieroglyph segmentation using a foundational segmentation
model, which could lay the groundwork for using AI in the recognition of these
hieroglyphs. We employed the foundational segmentation model, Segment Anything
(SAM), as a baseline to segment Maya texts, and introduced Finetuned-SAM, a
version of the SAM model fine-tuned with a dataset meticulously annotated with
the assistance of experts in Maya art and writing. The results from the
Finetuned-SAM model underscore its superior ability to segment inscriptions
accurately, even when training is conducted with a relatively limited dataset.
Our immediate future work includes expanding the dataset, as it is evident
that a more annotated dataset would potentially improve the results further.
The public release of these data could potentially ignite interest among
researchers in using recent advancements in AI to study ancient Maya writing,
as well as provide access to the inscriptions and their translations to
interested publics - both for the Maya heritage community to gain increased
access to their ancestral textual records, as well as for interested non-
specialists.

###  4.4 Acknowledgement

We are grateful to the Humanities Center at West Virginia University for their
generous support of this research through a humanities collaboration grant,
which will enhance our project and its contributions to the academic
community.

## References

  * Archive [2012] Justin Kerr Maya Archive.  Dumbarton oaks research library and collection.  _https://www.doaks.org/research/library-archives/icfa/collections/justin-kerr-maya-archive-1966-2012-1_ , 2012.
  * Bogacz et al. [2018] Bartosz Bogacz, Felix Feldmann, Christian Prager, and Hubert Mara.  Visualizing networks of maya glyphs by clustering subglyphs.  In _Eurographics Workshop on Graphics and Cultural Heritage_ , pages 105–111. The Eurographics Association, 2018.
  * Bommasani et al. [2021] Rishi Bommasani, Drew A Hudson, Ehsan Adeli, Russ Altman, Simran Arora, Sydney von Arx, Michael S Bernstein, Jeannette Bohg, Antoine Bosselut, Emma Brunskill, et al.  On the opportunities and risks of foundation models.  _arXiv preprint arXiv:2108.07258_ , 2021.
  * Fash and Fash [Spring-Autumn 1996] William Fash and Barbara Fash.  Building a world-view: Visual communication in classic maya architecture.  _RES: Anthropology and Aesthetics_ , (29/30):127–147, Spring-Autumn 1996.
  * Hu et al. [2017] Rui Hu, Jean-Marc Odobez, and Daniel Gatica-Perez.  Extracting maya glyphs from degraded ancient documents via image segmentation.  _Journal on Computing and Cultural Heritage (JOCCH)_ , 10(2):1–23, 2017.
  * Kerr [1966–2012] Justin Kerr.  Maya vase database: An archive of rollout photographs.  https://www.mayavase.com/, 1966–2012.
  * Kerr [2007] Justin Kerr.  A short history of rollout photography.  _www.famsi.org/research/kerr/rollout.html_ , 2007.
  * Kirillov et al. [2023] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C Berg, Wan-Yen Lo, et al.  Segment anything.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 4015–4026, 2023.
  * Li [2019] Lianfa Li.  Deep residual autoencoder with multiscaling for semantic segmentation of land-use images.  _Remote Sensing_ , 11(18):2142, 2019.
  * Ma et al. [2024] Jun Ma, Yuting He, Feifei Li, Lin Han, Chenyu You, and Bo Wang.  Segment anything in medical images.  _Nature Communications_ , 15(1):654, 2024.
  * Ohgushi et al. [2021] Toshiaki Ohgushi, Kenji Horiguchi, and Masao Yamanaka.  Road obstacle detection method based on an autoencoder with semantic segmentation.  In _Proceedings of the Asian conference on computer vision (ACCV), Kyoto, Japan, November 30-December 4, 2020_ , 2021.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical image computing and computer-assisted intervention–MICCAI 2015: 18th international conference, Munich, Germany, October 5-9, 2015, proceedings, part III 18_ , pages 234–241. Springer, 2015.
  * Sun et al. [2017] Chen Sun, Abhinav Shrivastava, Saurabh Singh, and Abhinav Gupta.  Revisiting unreasonable effectiveness of data in deep learning era.  In _Proceedings of the IEEE international conference on computer vision_ , pages 843–852, 2017.
  * Vaughan [2018] Jennifer Wortman Vaughan.  Making better use of the crowd: How crowdsourcing can advance machine learning research.  _Journal of Machine Learning Research_ , 18(193):1–46, 2018.