Twist decomposition of non-linear effects in Balitsky–Kovchegov
evolution of proton structure functions

Leszek Motyka1 and Mariusz Sadzikowski2

Jagiellonian University, Institute of Theoretical Physics,
Łojasiewicza 11, 30-348 Kraków, Poland

1<EMAIL_ADDRESS>
2<EMAIL_ADDRESS>

October 31, 2023

Effects of non-linear small-$x$ evolution of the gluon distribution given by
the Balitsky–Kovchegov equation are analyzed within the collinear
approximation framework. We perform a twist decomposition of the proton
structure functions $F_{2}$ and $F_{\mathrm{L}}$ obtained from the
Balitsky–Kovchegov equation using the Mellin representation of the scattering
cross-sections at high energies. In both the structure functions we find
strong corrections coming from the non-linear effects in the gluon evolution
at twist 2, and strongly suppressed higher twist effects. This implies that
unitarization effects of high energy scattering amplitudes are mostly the
leading twist effect. Furthermore we consider the double logarithmic limit of
the Balitsky–Kovchegov equation for the collinear gluon distribution, and
compare the result to the Gribov–Levin–Ryskin equation. We find that these two
equations differ by two powers of the hard scale logarithm for the large
scales.

##  1 Introduction and conclusions

The Electron–Ion Collider will probe deep inelastic scattering (DIS) of
electrons on large nuclei at high collision energies [1, 2] and this will
allow to probe partons that carry a small fraction $x$ of the nucleus
momentum. The parton distribution functions grow steeply with decreasing
values of $x$, and this growth is driven mostly by gluons. Specifically, the
gluon distribution function $g(x,Q^{2})$ at the scale $Q^{2}$ enters the DIS
cross sections as $xg(x,Q^{2})$, that grows at small $x$ as $x^{-\lambda}$
with $\lambda>0$. For sufficiently small values of $x$ and sufficiently low
scales $Q^{2}$ this growth has to be slowed down and finally tamed by
unitarity corrections, see e.g. [3, 4, 5, 6, 7]. These corrections enter into
QCD evolution equations as non-linear terms, that may be interpreted as
effects of parton recombination at high parton density regime. At very high
densities the gluon production by parton splittings and the gluon
recombination are expected to balance, leading a phenomenon called gluon
saturation.

Non-linear corrections to linear evolution equation were studied within the
two main QCD frameworks. Within the collinear
Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) approach a non-linear
evolution equation was proposed by Gribov, Levin and Ryskin (GLR) [8], and by
Mueller and Qiu [9]. The unitarity effects, however, are more pronounced when
viewed from the perspective of small $x$ evolution as given by the
Balitsky–Fadin–Kuraev–Lipatov (BFKL) framework [10, 11, 12, 13]. It is so
because the BFKL framework is more sensitive to lower momentum scales, while
in the DGLAP framework the unitarity effects can be greatly reduced by
choosing sufficiently large factorization scale. The BFKL equation is linear.
The non-linearity enters to this equation in the form of triple BFKL ladder
interaction, called the triple pomeron vertex [3, 4]. The evolution equation
equivalent to BFKL with the triple pomeron interaction was derived by Balitsky
[14] and Kovchegov [5, 6] — the Balitsky–Kovchegov (BK) equation. In results
of further studies of gluon color correlations in the high density regime at
small $x$, the framework of Color Glass Condensate (CGC) was developed [15,
16, 17], and JIMWLK equation was derived for the target wave function [18, 19,
20, 21, 22].

Theoretical understanding of non-linear evolution equations in QCD is pretty
good by now, and a large number of phenomenological studies were perfomed to
constrain these effects. Probably the most successful analysis showing strong
arguments for gluon saturation was performed by Golec–Biernat and Wüsthoff
(GBW) [23, 24]. It describes HERA data from the inclusive DIS down to
photoproduction limit, and the diffractive DIS data in a simple unified
framework, called the GBW saturation model. This model allowed for a simple
and efficient extensions to include scale evolution effects [25] and to
describe also other observables, in particular for exclusive diffractive
processes [26, 27]. In a more formal approach predictions of non-linear
evolution equations were multiply tested against the data, mostly by fitting
the solutions of BK equations to data on proton structure functions, see e.g.
Refs. [28, 29, 30, 31, 32]. Thus, the small $x$ non-linear evolution equations
are quite successful in phenomenological applications. Besides this, however,
theoretical questions arise on their relation to the powerful collinear
framework. In particular, it is interesting to understand how the multiple
scattering and gluon recombination effects in the small $x$ formulation appear
in the Operator Product Expansion (OPE) framework — how they map onto the
twist expansion of hadron structure functions, and how they modify the DGLAP
evolution equation. The aim of the present paper is to address and study these
questions.

Multiple parton exchange in QCD induces higher twist terms in the OPE. The
canonical scaling of twist $\tau$ contribution to hadron structure function is
$(\Lambda/Q)^{\tau-2}$, where $\Lambda$ is a low hadronic scale, and $Q$ is
the large DIS scale. Hence the higher twist terms are power suppressed and for
sufficiently large $Q^{2}$ they may be safely neglected. At small $x$,
however, the evolution of higher twist terms is more rapid than of the leading
twist term. For instance, for the dominant small $x$ gluon exchange, the twist
4 contribution to the structure functions at small $x$ may be estimated as
$\sim[xg(x,Q^{2})]^{2}\Lambda^{2}/Q^{2}$ see e.g. [8, 9, 33, 34, 35] to be
compared with $\sim xg(x,Q^{2})$ twist 2 behavior, where $g(x,Q^{2})$ is the
collinear gluon distribution. It follows that the strong growth of
$xg(x,Q^{2})$ at small $x$ may partially compensate the $1/Q^{2}$ suppression
of the twist 4 term, and the higher twist corrections may become significant,
which would affect the quality of twist 2 DGLAP fits of the structure
functions. In fact, DGLAP fits to proton structure functions measured at HERA
down to $x\sim 2\cdot 10^{-5}$ for $Q^{2}>1$ GeV2 deteriorate for $Q^{2}<5$
GeV2 both for diffractive and inclusive structure functions [36, 37, 38, 39,
40], and inclusion of higher twist corrections at small $x$ was shown to
improve the quality of data description [37, 38, 39, 40].

Theoretical analysis of higher twist corrections to the structure functions is
not easy. First, the number of relevant operators grows quickly with
increasing twist. So does the complexity of evolution equations. Furthermore,
in the standard DGLAP approach the initial conditions are fitted to data, and
with more operators at higher twist, more information would be needed to
constrain the initial conditions for their evolution. Certainly, the currently
available data do not provide such information. Hence, it is justified to use
simplified models of higher twist contributions, that incorporate crucial
features coming from QCD. A very useful guidelines for such models comes from
saturation models or from small $x$ non-linear evolution equations. In the
present analysis we use the approach developed in Refs. [41, 42, 37], at first
for a twist decomposition of the proton structure functions in the saturation
model, and later extended to the twist decomposition of the BFKL cross
sections [43]. The method relies on Mellin transforms of the scattering cross
sections and relating singularities in the complex Mellin moment plane to
contributions with definite twists. Here we apply this approach to the
Balitsky–Kovchegov equation and we perform the twist decomposition of the non-
linear correction.

In the analysis we use a solution of the BK equation in the form of a series
in non-linearity, proposed by Kovchegov [6]. We aim to estimate the
corrections from non-linearity to the linear approximation in the region where
they are moderate, so we truncate the expansion to the first term, that is
from a single contribution of non-linearity. This is sufficient to obtain the
key conclusions: (1) large corrections to proton structure functions from the
non-linearity, that enter at twist 2, and (2) small higher twist corrections
induced by the non-linear term. These results lead to, perhaps surprising,
overall conclusion that non-linear BK corrections in the proton structure
functions are basically the leading twist effect. To be specific, we find that
the twist 2 correction from BK non-linearity to the BFKL result reaches
$-50\%$ in proton structure functions for $Q^{2}=5$ GeV2 already at
$x=10^{-3}$. On the other hand, the higher twist corrections from both BFKL
and BK are found to be at 1% level for $F_{2}$ and below 10% in
$F_{\mathrm{L}}$. It should be kept in mind, however, that the multiple
scattering effects resumed by the BK equation may be probed in a different way
in other processes, leading to a possibly different picture of saturation
effects.

On the top of this, by taking the double logarithmic limit of the BK equation
we obtain a non-linear evolution equation for the collinear gluon distribution
$xg(x,Q^{2})$. This equation resembles the GLR equation, but we find that it
takes a different form of the non-linear part. For a meaningful comparison it
is necessary to consider the double logarithmic limit. The two classes od
logarithms that are resummed correspond to powers of $t=\log(Q^{2})$ and
powers of $y=\log(x_{0}/x)$. For the hierarchy
$t\gg\bar{\alpha}_{\mathrm{s}}y\gg 1$ we find that the non-linear term we
obtain from the BK equation is weaker by two powers of $\log(Q^{2})$ than the
corresponding contribution in the GLR equation. This feature may be traced
back to vanishing of the triple BFKL ladder vertex when the collinear ordering
is imposed, as found in Ref. [44]. This is in full consistency with the
collinear evolution equation for quasi-partonic higher twist operators [33],
in which the gluonic ladder merging vertex is absent. Hence, we conclude that
at high $Q^{2}$ effects of the BK non-linearity are smaller than expected,
both in the gluon evolution and in the higher twist contributions. In the
region where $t\sim\bar{\alpha}_{\mathrm{s}}y\gg 1$, however, the powers of
logarithms agree in the GLR equation and in the double logarithmic limit of
the BK equation.

Furthermore, we investigate the origin of the strong effects of non-linearity
at twist 2. We get the most clear answer by performing an analysis of the
unitegrated gluon distribution $f(x,k^{2})\simeq k^{2}(\partial/\partial
k^{2})xg(x,k^{2})$, that depends on the gluon transverse momentum squared
$k^{2}$. In $f(x,k^{2})$ emerging from the BK equation unitarity effects are
very strong for $k^{2}<Q_{\mathrm{s}}^{2}(x)$, where the characteristic
$x$-dependent momentum scale, that is naturally interpreted as the saturation
scale. In the region $k^{2}<Q_{\mathrm{s}}^{2}(x)$, the distribution
$f(x,k^{2})$ is strongly suppressed. Hence, $Q_{\mathrm{s}}^{2}(x)$ plays the
rôle of a lower cut-off in the integral
$xg(x,Q^{2})\simeq\int^{Q^{2}}\frac{dk^{2}}{k^{2}}f(x,k^{2})$, connecting the
collinear gluon distribution with its unintegrated form. This is a sizable
correction that enters at twist 2. This implies, that the non-linear effects
are concentrated at low momentum scales and it is possible to factor them out
and absorb in the input value of the gluon distribution
$xg(x,\mu_{\mathrm{F}}^{2})$, provided that the factorization scale
$\mu_{\mathrm{F}}^{2}$ is reasonably larger than the saturation scale
$Q_{\mathrm{s}}^{2}(x)$. Hence, in order to clearly see non-linear BK effects
one should probe the region of $\mu_{\mathrm{F}}^{2}<Q_{\mathrm{s}}^{2}$. This
region is not accessible for DIS on the proton, but may be probed in DIS on
large nuclei, e.g. at the Electron–Ion Collider.

##  2 Proton structure functions from Balitsky–Kovchegov equation

The total inclusive cross-section $\sigma^{\gamma^{\ast}A}$ of the virtual
photon $\gamma^{\ast}$ scattering on large nucleus $A$, in the high energy
limit, can be described in terms of the structure functions
$F_{\mathrm{L}},F_{2}$

$F_{\mathrm{T,L}}(x,Q^{2})=\frac{Q^{2}}{4\pi^{2}\alpha_{\mathrm{em}}}\sigma^{\gamma^{\ast}A}_{\mathrm{T,L}},\;\;\;F_{2}(x,Q^{2})=F_{\mathrm{T}}(x,Q^{2})+F_{\mathrm{L}}(x,Q^{2})$
(1)

where

$\sigma^{\gamma^{\ast}A}_{\mathrm{T,L}}(x,Q^{2})=\int\frac{d^{2}\mathbf{r}}{4\pi}\int_{0}^{1}dz|\psi_{\mathrm{T,L}}(z,r,Q^{2})|^{2}\sigma_{q\bar{q}}(x,\mathbf{r}).$
(2)

The transverse ($T$) and longitudinal ($L$) photon wavefunctions
$\psi_{\mathrm{T,L}}(z,r,Q^{2})$ describe probability amplitudes of the
fluctuation of the virtual photon into a quark-antiquark dipole of the
transverse size $r=|\mathbf{r}|$ and a fraction $z$ of the longitudianal
light-cone momentum carried by the quark [45]. The color dipole scatter over a
nucleus with the cross section $\sigma_{q\bar{q}}$ given by an imaginary part
of the forward dipole - nucleon scattering amplitude
$N(x,\mathbf{r},\mathbf{b})$

$\sigma_{q\bar{q}}(x,\mathbf{r})=2\int
d^{2}\mathbf{b}\;N(x,\mathbf{r},\mathbf{b})\equiv\sigma_{0}N(y,\mathbf{r}),$
(3)

where $y=\log(x_{\mathrm{in}}/x)$ is a rapidity variable developed with
respect to some initial value $x_{\mathrm{in}}$ [5]. The cross section
parameter $\sigma_{0}$ is related to the effective nucleus radius
$\sigma_{0}=2\pi R_{A}^{2}$ via integration over the impact parameter
$\mathbf{b}$. In the above equations it is already assumed that the main
contribution to scattering comes from perturbative dipoles located far from
the edges of the nucleus, namely that a size of the dipole $r$ in the
transverse space is much smaller then the nucleus radius $R_{A}$. In this way
one neglects a non-trivial dependence of the amplitude on the impact parameter
$\vec{b}$, limiting only to a simple cylindrical geometry of the high energy
nucleus.

The rapidity evolution of the amplitude $N(y,\mathbf{r})$ is described by the
Balitsky–Kovchegov equation [14, 5], which in momentum space reads

$\frac{\partial\phi(y,k^{2}_{\perp})}{\partial
y}=\bar{\alpha}_{s}\int_{0}^{\infty}\frac{dq^{2}_{\perp}}{q^{2}_{\perp}}\left\\{\frac{q^{2}_{\perp}\phi(y,q^{2}_{\perp})-k^{2}_{\perp}\phi(y,k^{2}_{\perp})}{|q^{2}_{\perp}-k^{2}_{\perp}|}+\frac{k^{2}_{\perp}\phi(y,k^{2}_{\perp})}{\sqrt{4q^{4}_{\perp}+k^{2}_{\perp}}}\right\\}-\bar{\alpha}_{s}\phi^{2}(y,k^{2}_{\perp})$
(4)

where $\bar{\alpha}_{s}=\alpha_{\mathrm{s}}N_{\mathrm{c}}/\pi$ and
$\phi(y,k^{2}_{\perp})$ is the Fourier transform of the amplitude

$\phi(y,k^{2}_{\perp})=\int\frac{d^{2}r}{2\pi}e^{-i\mathbf{k}_{\perp}\cdot\mathbf{r}}\frac{N(y,\mathbf{r})}{r^{2}}.$
(5)

The first term in the BK equation is given by the BFKL kernel that can be
solved exactly using its eigenfunctions [12]

$\phi_{0}(y,k^{2}_{\perp})=\sum_{n=-\infty}^{\infty}\int_{c-i\infty}^{c+i\infty}\frac{d\gamma}{2\pi
i}\left(\frac{k^{2}_{\perp}}{Q_{0}^{2}}\right)^{-\gamma}C_{n}(\gamma)\exp\\{\bar{\alpha}_{s}\chi(n,\gamma)y+in\varphi\\}$
(6)

where $\chi(n,\gamma)$ are eigenvalues of the BK equation and coefficients
$C_{n}(\gamma)$ are determined by the initial condition. Note that we use the
convention in which the standard Mellin moment corresponds to $-\gamma$, and
the reference scale in the Mellin transform is $Q_{0}^{2}$. The fundamental
Mellin strip is located in the interval $c\in(0,1)$. At large rapidity $y$,
which is of main interest in this article, the dominant contribution is given
by $n=0$ eigenvalue

$\chi(\gamma)\equiv\chi(0,\gamma)=2\psi(1)-\psi(1-\gamma)-\psi(\gamma)$ (7)

and we adopt this approximation throughout the whole paper. The non-linear BK
equation can be solved using iterative procedure [6]. In the Mellin space one
can write

$\displaystyle\tilde{\phi}(y,\gamma)$ $\displaystyle=$
$\displaystyle\frac{1}{Q_{0}^{2}}\int_{0}^{\infty}dk^{2}_{\perp}\left(\frac{k^{2}_{\perp}}{Q_{0}^{2}}\right)^{\gamma-1}\phi(y,k^{2}_{\perp})$
$\displaystyle\tilde{\phi}(y,\gamma)$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{\infty}\tilde{\phi}_{i}(y,\gamma),\;\;$ (8)

and the functions $\tilde{\phi}_{i}$ satisfy the infinite set of the coupled
equations:

$\displaystyle\frac{\partial\tilde{\phi}_{0}(y,\gamma)}{\partial y}$
$\displaystyle=$
$\displaystyle\bar{\alpha}_{s}\chi(\gamma)\tilde{\phi}_{0}(y,\gamma),$ (9)
$\displaystyle\frac{\partial\tilde{\phi}_{1}(y,\gamma)}{\partial y}$
$\displaystyle=$
$\displaystyle\bar{\alpha}_{s}\chi(\gamma)\tilde{\phi}_{1}(y,\gamma)-2\pi
i\,\bar{\alpha}_{s}\int_{c_{1}-i\infty}^{c_{1}+i\infty}\frac{d\gamma_{1}}{2\pi
i}\int_{c_{2}-i\infty}^{c_{2}+i\infty}\frac{d\gamma_{2}}{2\pi
i}\delta(\gamma-\gamma_{1}-\gamma_{2})\tilde{\phi}_{0}(y,\gamma_{1})\tilde{\phi}_{0}(y,\gamma_{2}),$
$\displaystyle\frac{\partial\tilde{\phi}_{2}(y,\gamma)}{\partial y}$
$\displaystyle=$
$\displaystyle\bar{\alpha}_{s}\chi(\gamma)\tilde{\phi}_{2}(y,\gamma)$
$\displaystyle-$ $\displaystyle 2\pi
i\,\bar{\alpha}_{s}\int_{c_{1}-i\infty}^{c_{1}+i\infty}\frac{d\gamma_{1}}{2\pi
i}\int_{c_{2}-i\infty}^{c_{2}+i\infty}\frac{d\gamma_{2}}{2\pi
i}\delta(\gamma-\gamma_{1}-\gamma_{2})(\tilde{\phi}_{0}(y,\gamma_{1})\tilde{\phi}_{1}(y,\gamma_{2})+\tilde{\phi}_{1}(y,\gamma_{1})\tilde{\phi}_{0}(y,\gamma_{2})),$
$\displaystyle.....$

and $0<\mathrm{Re}(c_{1},c_{2})<1$. From the structure of the equations one
can infer that $\tilde{\phi}_{i}\sim\tilde{\phi}_{0}^{i+1}$, where the term
$\tilde{\phi}_{0}\sim\exp(\alpha\chi y)$. Therefore the amplitude
$\tilde{\phi}_{i}$ describes the process with $i+1$ pomeron exchanges included
in the DIS diagram. The series (2) is convergent if
$\tilde{\phi}_{0}(y,k^{2}_{\perp})\ll 1$ which locates the quark transverse
momenta within the perturbative domain. The thorough analysis in [6] shows
that the series is convergent for values of quark momenta above the saturation
scale with possible extension using analytical continuation. Such requirement
is consistent with the twist decomposition which assumes that the virtuality
of the quark–antiquark pair is large with respect to the typical hadronic
scale. The procedure is not applicable below the saturation scale. However,
our main goal in this paper is to estimate the influence of the non-linear
corrections from the BK equation on the lowest twists of the structure
functions. Therefore, we limit our calculation to the first order correction.
It is a simple exercise to solve the second equation of (9) with the result

$\displaystyle\tilde{\phi}_{1}(y,\gamma)$ $\displaystyle=$
$\displaystyle\int\frac{d\gamma_{1}}{2\pi i}\frac{d\gamma_{2}}{2\pi i}\,2\pi
i\delta(\gamma-\gamma_{1}-\gamma_{2})C_{0}(\gamma_{1})C_{0}(\gamma_{2})$ (10)
$\displaystyle\frac{\exp{(\bar{\alpha}_{s}y\chi(\gamma))}-\exp{(\bar{\alpha}_{s}y\chi(\gamma_{1})+\bar{\alpha}_{s}y\chi(\gamma_{2}))}}{\chi(\gamma_{1})+\chi(\gamma_{2})-\chi(\gamma)},$

whereas the first equation gives the BFKL solution for $n=0$ eigenvalue. For
further analysis we adopt the exponential form of the initial conditions

$N(y=0,\mathbf{r})=1-\exp\left(-r^{2}Q_{0}^{2}\right)$ (11)

which gives

$C_{0}(\gamma)=-\frac{2^{2\gamma-1}\Gamma(\gamma)}{\Gamma(1-\gamma)}\Gamma(-\gamma)$
(12)

and the BFKL solution takes the form

$\tilde{\phi}_{0}(y,\gamma)=-\frac{2^{2\gamma-1}\Gamma(\gamma)}{\Gamma(1-\gamma)}\Gamma(-\gamma)e^{\bar{\alpha}_{s}\chi(\gamma)y},$
(13)

whereas the Mellin strip is limited to $0<\mathrm{Re}(c)<3/4$. Note that we
changed convention for the $Q_{0}$ parameter of the input fuction
$N(y=0,\mathbf{r})$ w.r.t. Ref. [43]: the present $Q_{0}$ equals $1/2$ of
$Q_{0}$ used in [43].

##  3 Twist decomposition of the DIS cross section

Figure 1: Location of singularities in a complex Mellin space of $\gamma$
variable. Blue points (on the horizontal axis) correspond to singularities
located at integer values, whereas the red ones (above the axis) to integer
$\gamma-\gamma_{1}$. The integration over the large contour is equal to the
infinite sum of integrals over small contours encircling singular points with
a minus sign. One such contour is depicted around point $\gamma=2$.

The twist structure of the photon – nucleus scattering can be obtained from
the Mellin transform of the cross section with respect to the virtuality
scale. The twist decomposition is performed by isolating contributions of
singularities in the complex Mellin moment plane.

The starting point of our analysis is the iterative solution of the BK
equation with respect to the non-linear interaction term, as described in the
previous section. In this framework the cross section can be described as a
series

$\sigma^{\gamma^{\ast}A}_{\mathrm{T,L}}=\sum_{i=0}^{\infty}\sigma^{(i)\,\gamma^{\ast}A}_{\mathrm{T,L}},$

where

$\sigma^{(i)\,\gamma^{\ast}A}_{\mathrm{T,L}}(x,Q^{2})=\sigma_{0}\int_{c-i\infty}^{c+i\infty}\frac{d\gamma}{2\pi
i}\left(\frac{4Q_{0}^{2}}{Q^{2}}\right)^{-\gamma}\tilde{H}_{\mathrm{T,L}}(-\gamma)\frac{\Gamma(1+\gamma)}{2^{-2\gamma-1}\Gamma(-\gamma)}\tilde{\phi}_{i}(y,-\gamma).$
(14)

The Mellin fundamental strip is located in $-3/4<\mathrm{Re}(c)<0$ and
$\tilde{\phi}_{i}$ are solutions of the equations (9). The functions
$\tilde{H}_{\mathrm{T,L}}$ are Mellin transforms of the photon wave functions
that can be found in [41, 42]. The leading BFKL contribution is given by the
formula

$\sigma^{(0)\,\gamma^{\ast}A}_{\mathrm{T,L}}(x,Q^{2})=-\sigma_{0}\int_{-1/2-i\infty}^{-1/2+i\infty}\frac{d\gamma}{2\pi
i}\left(\frac{4Q_{0}^{2}}{Q^{2}}\right)^{-\gamma}\tilde{H}_{\mathrm{T,L}}(-\gamma)\Gamma(\gamma)e^{\bar{\alpha}_{s}y\chi(-\gamma)},$
(15)

and was described in [43]. The lowest order correction
$\sigma^{(1)\,\gamma^{\ast}A}_{\mathrm{T,L}}$ to the BFKL cross section
follows from the solution (10) and decomposition (14)

$\displaystyle\sigma^{(1)\,\gamma^{\ast}A}_{\mathrm{T,L}}$ $\displaystyle=$
$\displaystyle\sigma_{0}\int_{c-i\infty}^{c+i\infty}\frac{d\gamma}{2\pi
i}\left(\frac{4Q_{0}^{2}}{Q^{2}}\right)^{\gamma}\tilde{H}_{\mathrm{T,L}}(\gamma)\frac{\Gamma(1-\gamma)}{\Gamma(\gamma)}\int_{c_{1}-\infty}^{c_{1}+\infty}\frac{d\gamma_{1}}{2\pi
i}\frac{\Gamma(\gamma_{1})\Gamma(\gamma-\gamma_{1})}{2\gamma_{1}(\gamma-\gamma_{1})}\times$
(16) $\displaystyle\times$
$\displaystyle\frac{e^{(\bar{\alpha}_{s}y\chi(\gamma))}-e^{(\bar{\alpha}_{s}y\chi(\gamma_{1})+\bar{\alpha}_{s}y\chi(\gamma-\gamma_{1}))}}{\chi(\gamma_{1})+\chi(\gamma-\gamma_{1})-\chi(\gamma)}$

where $0<\mathrm{Re}(c),\mathrm{Re}(c_{1})<3/4$. It is important to note that
in the above expression both exponent factors are important to maintain
correct analytical structure. Indeed, the multiples zeros of the denominator
in the complex plain are exactly canceled by the contribution from the
exponents in the numerator. The essential singularities are located at integer
values of $\gamma,\gamma_{1}$ and $\gamma-\gamma_{1}$ (see Fig. 1), therefore
the integration over $\gamma$ variable can be decompose into two sums

$\displaystyle\sigma^{(1)\,\gamma^{\ast}A}_{\mathrm{T,L}}\approx\sum_{n=1}^{N}\int_{c-i\infty}^{c+i\infty}\frac{d\gamma_{1}}{2\pi
i}\int_{C_{n}}\frac{d\gamma}{2\pi
i}\tilde{I}(\gamma,\gamma_{1})+\sum_{n=1}^{\infty}\int_{c-i\infty}^{c+i\infty}\frac{d\gamma_{1}}{2\pi
i}\int_{C_{n+\gamma_{1}}}\frac{d\gamma}{2\pi i}\tilde{I}(\gamma,\gamma_{1}),$
(17)
$\displaystyle{I(\gamma,\gamma_{1})=\sigma_{0}\left(\frac{4Q_{0}^{2}}{Q^{2}}\right)^{\gamma}\tilde{H}_{\mathrm{T,L}}(\gamma)\frac{\Gamma(1-\gamma)}{\Gamma(\gamma)}\frac{\Gamma(\gamma_{1})\Gamma(\gamma-\gamma_{1})}{2\gamma_{1}(\gamma-\gamma_{1})}\frac{e^{(\bar{\alpha}_{s}y\chi(\gamma))}-e^{(\bar{\alpha}_{s}y\chi(\gamma_{1})+\bar{\alpha}_{s}y\chi(\gamma-\gamma_{1}))}}{\chi(\gamma_{1})+\chi(\gamma-\gamma_{1})-\chi(\gamma)},}$

where $C_{w}$ is a small clockwise contour located around point $w$ in the
complex $\gamma$ plane. The Cauchy theorem that bridges equations (16, 17) is
satisfied due to the exponential suppression of the large imaginary values
$\mathrm{Im}\gamma$ brought by the photon wavefunctions $H_{\mathrm{T,L}}$.
The asymptotics of $\chi(\gamma)$ does not spoil that in any way, as it
consists of two pieces: the $\,-\log(\gamma)$ term that improves the
convergence and the $-\pi\cot(\pi\gamma)$ that contributes to the (only) poles
of $\chi(\gamma)$ at the integer values of $\gamma$, but is subleading in the
other regions. The series (17) is asymptotic [41, 42] and the optimal number
$N$ depends on the values of $(x,Q)$. In particular, one can check numerically
that for both the transverse and longitudinal cross sections the first two
poles provide the agreement with the full result (16) at the level better than
one per mile at $Q^{2}=5$ GeV2 and $x=10^{-3}$. Additionally, inclusion of the
second pole on the top of the first improve the result by more than the order
of magnitude, which shows the convergence of the expansion for the first two
twists.

The integral $d\gamma$ over $C_{n}$ from the first term in (17) gives a direct
contribution to twist $\tau=2n$ after integration over $d\gamma_{1}$ along the
line parallel to imaginary axis within the Mellin strip. A similar integration
over $C_{n+\gamma_{1}}$ from the second term gives contributions to all twists
of order $\tau\geq 2n$ starting with the lowest value $\tau=4$. This fact
follows from the integration over $d\gamma_{1}$ variable. However, the twists
higher then $\tau=4$ are strongly suppressed, therefore one can assign the
contribution of the second integral to twist $\tau=4$ only, with a small error
of order 1 per cent or less. Summarising, in numerical calculations, the BK
correction to the leading twist was calculated using the first term of (17)
with $n=1$ only, whereas the correction to twist $\tau=4$ by the first term
with $n=2$ and the second term with $n=1$.

##  4 Results

###  4.1 Numerical results

Figure 2: Effects of non-linear corrections in proton structure functions
$F_{2}(x,Q^{2})$ (left) and $F_{\mathrm{L}}(x,Q^{2})$ (right) at
$Q^{2}=5~{}\mbox{GeV}^{2}$ in twist 2 approximation. We show the twist 2
components obtained from the BFKL evolution, the BK correction and their sum.

With the framework described above we calculate higher twist corrections to
the unpolorized proton structure functions $F_{2}$ and $F_{\mathrm{L}}$ at
small $x$. The primary goal is to estimate deviations from linear evolution
(BFKL) regime due to non-linearity introduced by the triple Pomeron
interaction. We compute the first correction in non-linearity in the iterative
solution of Balitsky–Kovchegov equation — which we shall call the BK
correction. Moreover, we perform an explicit twist decomposition of the
structure functions for BFKL result with the BK correction and compare it to
the known results computed earlier in the BFKL approach.

As the reference we take the structure functions obtained from a solution of
the leading logarithmic BFKL equation with parameters obtained in Ref. [43].
Let us remind that the input for the BFKL evolution in the dipole
representation at $x_{\mathrm{in}}=0.1$ is assumed to take the GBW form:
$\sigma(x_{\mathrm{in}},r)=\sigma_{0}[1-\exp(-r^{2}Q_{0}^{2})],$ with
$\sigma_{0}=17.04$ mb and $Q_{0}=0.255$ GeV. The value of the strong coupling
constant $\bar{\alpha}_{\mathrm{s}}$ in the BFKL kernel is set to $0.087$.
This should be understood as an effective value of $\bar{\alpha}_{\mathrm{s}}$
that partially absorbs the higher order corrections to the BFKL kernel, known
to reduce strongly the BFKL Pomeron intercept. This is consistent with
application of the Brodsky–Lepage–Mackenzie scale fixing procedure [46] to the
NLL BFKL kernel [47]. We stress, however, that the primary goal of the present
study is to understand importance of non-linear corrections to BFKL evolution
and its twist decomposition and fine details of the model should not affect
the key, general features of the results.

We start the numerical analysis from evaluating the BK correction to the BFKL
evolution. We compute the leading twist 2 contributions to structure functions
$F_{2}(x,Q^{2})$ and $F_{\mathrm{L}}(x,Q^{2})$. We choose the DIS reference
scale $Q^{2}=5~{}\mbox{GeV}^{2}$, below which the DGLAP fit deteriorates of
the final HERA data on structure functions. As it will be clear from the next
part of the analysis, the higher twist corrections are small and do not change
the conclusions of this part. In Fig. 2 we display the twist 2 contributions
to the structure functions from: the BFKL equation, the BK correction and from
the sum of BFKL and BK parts. Both the BFKL part and the BK correction are
obtained with the same input. The BK corrections to both structure function
are large and negative. The magnitude of the corrections, both absolute and
relative, grows with decreasing $x$. Clearly, when the relative correction is
not small, the higher order corrections of the non-linearity would be
necessary to achieve a good approximation of the complete solution of the
Balitsky–Kovchegov equation. For the present study it is sufficient to
conclude that the non-linear corrections to BFKL results at twist 2 are large
already at $x=0.001$.

Figure 3: Relative effects of twist 4 corrections in proton structure
functions $F_{2}(x,Q^{2})$ (left) and $F_{\mathrm{L}}(x,Q^{2})$ (right) at
$Q^{2}=5~{}\mbox{GeV}^{2}$. We show the ratios of the twist 4 corrections
obtained from the BFKL evolution, the BK correction and their sum to the twist
2 BFKL + BK result.

Next we turn to the analysis of higher twist corrections to the proton
structure functions. The most important measure of higher twist effects is its
relative magnitude to the twist 2 approximation. We choose as a reference the
twist 2 estimates of the structure functions obtained from BFKL equation with
the BK correction. We restrict the analysis to twist $\tau=4$ effects. As it
is clearly seen in Fig. 3, the relative twist 4 effects are small or moderate,
so it is not necessary to consider the higher, $\tau>4$, corrections. The
magnitudes of relative corrections are different, but the general pattern is
similar for both the structure functions. The twist 4 corrections coming from
BFKL are negative, and the BK twist 4 contributions are positive, but with
smaller absolute value than the leading BFKL contribution. The overall higher
twist corrections are negative in both structure functions. The relative
higher twist corrections are found to be larger for $F_{\mathrm{L}}$, where
they reach up to about (negative) 10%. For $F_{2}$ we find corrections up to
about (negative) $1.5\%$. This is expected, as the coefficient function for
$F_{\mathrm{L}}$ generates a scale logarithm for twist 4, and does not for
twist 2, while for $F_{\mathrm{T}}$, that is dominant in $F_{2}$, the twist 4
contribution carries one less power of the scale logarithm than the leading
twist 2 term [41, 42]. One should keep in mind, however, that the twist
content of the cross sections depends not only on the evolution equation, but
also on the form of the input, and this dependence is stronger when $x$ is not
very small.

The presented results are obtained taking the first order corrections in non-
linearity. Most likely this leads to overestimating the non-linear effects. It
is expected to happen because the expansion in non-linearity produces the
alternating series for the dipole scattering amplitude, as demonstrated in
[6]. Effects of the higher orders in non-linearity should turn in when the
relative first order correction stops being much smaller than one, and this is
certainly the case when the first order correction reaches 50%. We expect that
at higher orders the correction should be somewhat smaller, but this not
changes the overall pattern.

We performed a similar analysis also for $Q^{2}=2$ GeV2 and for $Q^{2}=10$
GeV2 and found that the general pattern is very similar to the case $Q^{2}=5$
GeV2. Therefore, we do not depict them in separate figures. Instead, the
numerical values of non-linear corrections at twist 2 are given in Table 1 for
these three values of $Q^{2}$ at $x=0.01$ and $x=0.001$. To be specific, we
provide the numerical values of BK corrections $\Delta F_{\mathrm{2,BK}}$,
$\Delta F_{\mathrm{L,BK}}$ compared to the BFKL results in the leading twist
for $F_{2}$ and $F_{\mathrm{L}}$ structure functions.

$x=0.01$ |  $Q^{2}=2$ GeV2 |  $Q^{2}=5$ GeV2 |  $Q^{2}=10$ GeV2
---|---|---|---
$\Delta F_{\mathrm{2,BK}}/F_{\mathrm{2,BFKL}}$ | -26$\%$ | -24$\%$ | -22$\%$
$\Delta F_{\mathrm{L,BK}}/F_{\mathrm{L,BFKL}}$ | -21$\%$ | -19$\%$ | -18$\%$
$x=0.001$ |  |  |
$\Delta F_{\mathrm{2,BK}}/F_{\mathrm{2,BFKL}}$ | -73$\%$ | -64$\%$ | -59$\%$
$\Delta F_{\mathrm{L,BK}}/F_{\mathrm{L,BFKL}}$ | -60$\%$ | -53$\%$ | -48$\%$
Table 1: Contribution of the BK corrections to structure functions
$F_{2},F_{\mathrm{L}}$ at the leading twist.

In Table 2 one can find the ratio of twist 4 and twist 2 for $F_{2}$
($R_{F_{2}}$) and $F_{\mathrm{L}}$ ($R_{F_{\mathrm{L}}}$) structure functions:
$R_{F_{2}}=(F_{\mathrm{2,BFKL}}^{(4)}+\Delta
F_{\mathrm{2,BK}}^{(4)})/(F_{\mathrm{2,BFKL}}^{(2)}+\Delta
F_{\mathrm{2,BK}}^{(2)})$,
$R_{F_{\mathrm{L}}}=(F_{\mathrm{L,BFKL}}^{(4)}+\Delta
F_{\mathrm{L,BK}}^{(4)})/(F_{\mathrm{L,BFKL}}^{(2)}+\Delta
F_{\mathrm{L,BK}}^{(2)})$, as well as, the ratio of BK correction twist 4 to
twist 2 for both structure functions: $R_{\Delta F_{2}}=\Delta
F_{\mathrm{2,BK}}^{(4)}/(F_{\mathrm{2,BFKL}}^{(2)}+\Delta
F_{\mathrm{2,BK}}^{(2)})$, $R_{\Delta F_{\mathrm{L}}}=\Delta
F_{\mathrm{L,BK}}^{(4)}/(F_{\mathrm{L,BFKL}}^{(2)}+\Delta
F_{\mathrm{L,BK}}^{(2)})$.

$x=0.01$ |  $Q^{2}=2$ GeV2 |  $Q^{2}=5$ GeV2 |  $Q^{2}=10$ GeV2
---|---|---|---
$R_{F_{2}}$ | -2.2$\%$ | -1.2$\%$ | -0.7$\%$
$R_{\Delta F_{2}}$ | 0.4$\%$ | 0.15$\%$ | 0.08$\%$
$R_{F_{\mathrm{L}}}$ | -17.4$\%$ | -8.9$\%$ | -5.1$\%$
$R_{\Delta F_{\mathrm{L}}}$ | 2.1$\%$ | 1.0$\%$ | 0.5$\%$
$x=0.001$ |  |  |
$R_{F_{2}}$ | -0.2$\%$ | -0.5$\%$ | -0.3$\%$
$R_{\Delta F_{2}}$ | 1.6$\%$ | 0.6$\%$ | 0.3$\%$
$R_{F_{\mathrm{L}}}$ | -8.6$\%$ | -4.3$\%$ | -2.4$\%$
$R_{\Delta F_{\mathrm{L}}}$ | 10.3$\%$ | 4.0$\%$ | 2.0$\%$
Table 2: Ratio of twist 4 to twist 2 for structure functions
$F_{2},F_{\mathrm{L}}$ and non-linear corrections. Definitions can be found in
the main text.

The overall pattern seen in both the tables agrees with expectations. The
effects of non-linearity are strongest at the leading twist, they grow with
decreasing $x$ and slightly decrease with $Q^{2}$. The higher twist
corrections due to non-linearity behave in a similar way, but they are much
smaller and their decrease with $Q^{2}$ is much faster. Let us mention that
the results for $Q^{2}=2$ GeV2 are presented here rather for completeness than
for phenomenological purposes. At such low values of $Q^{2}$ and low $x$ the
non-linear corrections are of the same order as the linear part. Therefore,
accurate numerical predictions for sure require higher terms of the expansion
in non-linearity.

###  4.2 Discussion

The obtained results show a stronger dependence on $x$ of the twist 4 BK
correction than the twist $4$ BFKL contribution. This is an expected results.
In an earlier study [43] we found that in the double logaritmic saddle point
approximation the rapidity dependence of the twist 4 BFKL term is governed by
$\exp\left(2\sqrt{\bar{\alpha}_{\mathrm{s}}y\log(Q^{2}/Q_{0}^{2})}-2\alpha_{\mathrm{s}}y\right)$
(up to power factors of $y$), to be compared with twist 2 BFKL in the same
approximation
$\sim\exp\left(2\sqrt{\bar{\alpha}_{\mathrm{s}}y\log(Q^{2}/Q_{0}^{2}}\right)$.
This means that in the leading logarithmic BFKL evolution terms corresponding
to a double gluon ladder exchange in the total cross section, expected to grow
as
$\sim\left[\exp\left(2\sqrt{\bar{\alpha}_{\mathrm{s}}y\log(Q^{2}/Q_{0}^{2})}\right)\right]^{2}$,
are absent. In other words the multiple elementary $t$-channel gluons present
in the reggeized gluons of the BFKL formalism have zero projection on the
leading twist 4 exchange in the double logarithmic approximation. The BK
correction is different because it is generated by the triple gluon ladder
vertex, that couples the genuine two gluon ladder contribution to the BFKL
evolved $\gamma^{*}$ scattering. Therefore one expects the strong
$\sim\left[\exp\left(2\sqrt{\bar{\alpha}_{\mathrm{s}}y\log(Q^{2}/Q_{0}^{2})}\right)\right]^{2}$
growth of the BK twist 4 correction at asymptotically large $y$ and its
dominance in total twist 4 at very small $x$. The presented results indicate
however, that this asymptotic regime is not reached in HERA kinematics.

The overall picture emerging from the numerical analysis is quite clear. It
turns out that the non-linear evolution, as given by the BK equation, have
strong effects in the leading twist 2 component of the structure functions,
and the higher twist components coming from the BK equation are strongly
suppressed w.r.t. the leading twist terms. It should be stressed that the
large non-linear corrections found in BK at twist 2 affect the results
obtained in the BFKL framework, and the corresponding effect in the DGLAP
framework depends on the factorization scale, as discussed below. Combining
the large BK effects in twist 2 and weak at higher twists, we conclude that
the non-linear corrections are concentrated at low scales. Hence we expect
that in DGLAP framework the BK corrections can be mostly absorbed into the
input for the twist 2 gluon evolution, provided the initial scale of DGLAP
evolution $\mu_{\mathrm{F,0}}$ is big enough. Clearly, the scale for the BK
effects is the saturation scale $Q_{\mathrm{s}}$, so we conclude that with
$\mu_{\mathrm{F,0}}\gg Q_{\mathrm{s}}$ the DGLAP description should not be
significantly affected by non-linear evolution effects. This conclusion may
change when $\mu_{\mathrm{F,0}}<Q_{\mathrm{s}}(x)$ for some range of $x$
probed by the data. Then one would expect a significant modification the twist
2 DGLAP evolution due to non-linearity. Fits to the proton structure functions
at HERA assuming saturation effects indicate that in HERA kinematics
$Q_{\mathrm{s}}^{2}<1$ GeV2. Hence, with DGLAP fits assuming typically
$\mu_{\mathrm{F,0}}^{2}=2$ GeV2 or a higher value, they should not be affected
by non-linearity. The situation may change, however, for DIS on a large
nucleus with the mass number $A$, for which the saturation scale
$Q_{\mathrm{s}}^{2}$ is enhanced by $A^{1/3}$ w.r.t. the saturation scale in
proton. In order to account for the possibility of $Q_{\mathrm{s}}>Q_{0}$, in
the next section we shall consider the effects of non-linearity on the $Q^{2}$
evolution in this regime.

##  5 Non-linear evolution in the collinear approximation

###  5.1 The double logarithmic regime

In the numerical analysis we found strong effects of the BK correction in
twist 2 components of the structure functions. In order to better understand
the origin of this effect let us consider the double logarithmic limit of the
BK equation. It is convenient to start from the BK evolution for the
unintegrated gluon distribution, $f(x,k^{2})$, related to the collinear gluon
distribution by the LL formula,
$xg(x,Q^{2})=\int^{Q^{2}}dk^{2}f(x,k^{2})/k^{2}$. In what follows we shall
also use notation $f(x,k^{2})\to f(y,k^{2})$ with $y=\log(x_{\mathrm{in}}/x)$.
We keep the convention for the Mellin transform used in the previous sections:

$\tilde{f}(y,\gamma)=\int{dk^{2}\over
k^{2}}f(y,k^{2})(k^{2}/Q_{0}^{2})^{\gamma},\qquad\qquad
f(y,k^{2})=\int_{c-i\infty}^{c+i\infty}{d\gamma\over 2\pi
i}\tilde{f}(y,\gamma)(k^{2}/Q_{0}^{2})^{-\gamma}.$ (18)

The unintegrated gluon distribution is related to the dipole amplitude through
the formula [48]

$f(x,k^{2})=\frac{N_{c}S_{\perp}}{4\alpha_{s}\pi^{2}}(k^{2})^{2}\nabla^{2}_{k}\phi(y,k^{2})|_{y=\log(1/x)},$
(19)

where the function $\phi(y,k^{2})$ is the Fourier transform of
$N(y,\mathbf{r})/r^{2}$, see Eq. (5). The BK equation for the unintegrated
gluon distribution reads [28, 48]:

${\partial f(y,k^{2})\over\partial
y}=\bar{\alpha}_{\mathrm{s}}k^{2}\int{da^{2}\over
a^{2}}\left[{f(y,a^{2})-f(y,k^{2})\over|a^{2}-k^{2}|}+{f(y,k^{2})\over\sqrt{4a^{4}+k^{4}}}\right]$
$-{2\pi\alpha_{\mathrm{s}}^{2}\over
S_{\perp}}\left[k^{2}\,\int_{k^{2}}{da^{2}\over
a^{4}}f(y,a^{2})\int_{k^{2}}{db^{2}\over
b^{4}}f(y,b^{2})\,+\,f(y,k^{2})\int_{k^{2}}{da^{2}\over
a^{4}}\log(a^{2}/k^{2})\,f(y,a^{2})\right],$ (20)

where $S_{\perp}$ is the transverse target area, for a uniform target with
radius $R_{A}$, $S_{\perp}=\pi R_{A}^{2}$. The first line represents the
linear part (the BFKL equation), and the second line is a non-linear
correction corresponding to the triple pomeron interaction in the BK equation.
Note that the non-linear term corresponds to integrals with a strict anti-
collinear ordering, $a^{2},b^{2}>k^{2}$, so it vanishes in the collinear
limit. In the Mellin representation this equation reads:

${\partial\tilde{f}(y,\gamma)\over\partial
y}=\bar{\alpha}_{\mathrm{s}}\chi(-\gamma)\tilde{f}(y,\gamma)$
$-{2\pi\alpha_{\mathrm{s}}^{2}\over S_{\perp}Q_{0}^{2}}\int{d\gamma_{1}\over
2\pi i}\int{d\gamma_{2}\over 2\pi i}\,2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,\left[{1\over(\gamma_{1}+1)(\gamma_{2}+1)}+{1\over(\gamma_{1}+1)^{2}}\right]\,\tilde{f}(y,\gamma_{1})\tilde{f}(y,\gamma_{2}),$
(21)

which, using the symmetry between $\gamma_{1}$ and $\gamma_{2}$ can be
rewritten as

${\partial\tilde{f}(y,\gamma)\over\partial
y}=\bar{\alpha}_{\mathrm{s}}\chi(-\gamma)\tilde{f}(y,\gamma)$
$-{\pi\alpha_{\mathrm{s}}^{2}\over S_{\perp}Q_{0}^{2}}\int{d\gamma_{1}\over
2\pi i}\int{d\gamma_{2}\over 2\pi i}\,2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,{(\gamma+1)^{2}\over(\gamma_{1}+1)^{2}(\gamma_{2}+1)^{2}}\,\tilde{f}(y,\gamma_{1})\tilde{f}(y,\gamma_{2}).$
(22)

This implies the evolution equation for collinear gluon distribution:

${\partial\tilde{g}(y,\gamma)\over\partial
y}=\bar{\alpha}_{\mathrm{s}}\chi(-\gamma)\tilde{g}(y,\gamma)$
$+{\pi\alpha_{\mathrm{s}}^{2}\over S_{\perp}Q_{0}^{2}}\int{d\gamma_{1}\over
2\pi i}\int{d\gamma_{2}\over 2\pi i}\,2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,{(\gamma+1)^{2}\over(\gamma_{1}+1)^{2}(\gamma_{2}+1)^{2}}\,{\gamma_{1}\tilde{g}(y,\gamma_{1})\,\gamma_{2}\tilde{g}(y,\gamma_{2})\over\gamma},$
(23)

where $\tilde{g}(y,\gamma)=\left.\int
dQ^{2}/Q^{2}\,xg(x,Q^{2})(Q^{2}/Q_{0}^{2})^{\gamma}\right|_{x=x_{\mathrm{in}}\exp(-y)}$
is the collinear gluon distribution $xg(x,Q^{2})$ in the Mellin
representation. Note that the relation
$f(x,k^{2})=k^{2}\partial_{k^{2}}xg(x,k^{2})$ leads to the
$\tilde{f}(y,\gamma)=(-\gamma)\tilde{g}(y,\gamma)$ relation in the
$(y,\gamma)$ variables. In the double logaritmic limit, which corresponds to
the leading powers of $\gamma$ around $\gamma=0$, we obtain:

$(-\gamma){\partial\tilde{g}(y,\gamma)\over\partial
y}=\bar{\alpha}_{\mathrm{s}}\tilde{g}(y,\gamma)$
$-{\pi\alpha_{\mathrm{s}}^{2}\over S_{\perp}Q_{0}^{2}}\int{d\gamma_{1}\over
2\pi i}\int{d\gamma_{2}\over 2\pi i}\,2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,{(\gamma+1)^{2}\over(\gamma_{1}+1)^{2}(\gamma_{2}+1)^{2}}\,{\gamma_{1}\tilde{g}(y,\gamma_{1})\,\gamma_{2}\tilde{g}(y,\gamma_{2})},$
(24)

where we approximated $\chi(-\gamma)\simeq-1/\gamma+{\cal O}(1)$ around
$\gamma\to 0$. Before further analysis of this equation let us compare it to
the GLR equation [8, 49]:

$\frac{\partial^{2}xg(x,Q^{2})}{\partial
y\partial\log(Q^{2}/Q_{0}^{2})}=\bar{\alpha}_{s}xg(x,Q^{2})-{N_{\mathrm{c}}\over
2C_{\mathrm{F}}}{\pi\alpha_{\mathrm{s}}^{2}\over
S_{\perp}Q^{2}}(xg(x,Q^{2}))^{2}$ (25)

which in the Mellin representation $(y,\gamma)$ takes the form:

$(-\gamma){\partial\tilde{g}(y,\gamma)\over\partial
y}=\bar{\alpha}_{\mathrm{s}}\tilde{g}(y,\gamma)$ $-{N_{\mathrm{c}}\over
2C_{\mathrm{F}}}{\pi\alpha_{\mathrm{s}}^{2}\over
S_{\perp}Q_{0}^{2}}\int{d\gamma_{1}\over 2\pi i}\int{d\gamma_{2}\over 2\pi
i}\,2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,\tilde{g}(y,\gamma_{1})\tilde{g}(y,\gamma_{2}),$
(26)

with $C_{\mathrm{F}}=(N_{\mathrm{c}}^{2}-1)/2N_{\mathrm{c}}$. Let us notice
that after the Mellin transform the non-linear term becomes the convolution as
expected, and the $1/Q^{2}$ factor changes into the $1/Q_{0}^{2}$ which
carries the physical dimension. We find differences between the non-linear
equation (24) and the GLR equation. The color prefactor of GLR
$N_{\mathrm{c}}/2C_{\mathrm{F}}=N_{\mathrm{c}}^{2}/(N_{\mathrm{c}}^{2}-1)$
differs from our by ${\cal O}(1/N_{\mathrm{c}}^{2})$ and this is beyond the
leading $N_{\mathrm{c}}$ accuracy of the BK equation. There is, however, a
more important difference in the integral kernel in the Mellin space. GLR
gives

$\tilde{K}_{\mathrm{GLR}}=2\pi i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma),$
(27)

while we obtain:

$\tilde{K}_{\mathrm{BK}}=2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,{\gamma_{1}\gamma_{2}(\gamma+1)^{2}\over(\gamma_{1}+1)^{2}(\gamma_{2}+1)^{2}}.$
(28)

In order to compare the essential properties of these two kernels and
resulting evolution equations let us consider their large $Q^{2}$ behavior. In
this region the evolution is dominated by linear term. Further, we want to
determine the leading powers of the logarithm $\log(Q^{2}/Q_{0}^{2})$ that
emerge from both the kernels.

At first, let us focus on the limit
$t=\log(Q^{2}/Q_{0}^{2})\gg\bar{\alpha}_{\mathrm{s}}y\gg 1$. In this region
the saddle point solution of the linear evolution equation is dominated by the
anomalous dimension $\gamma_{s}=\sqrt{\bar{\alpha}_{\mathrm{s}}y/t}\ll 1$.
Note that in the convention for the Mellin moments applied in this paper, they
equal to minus anomalous dimensions.

As already said, in the large $Q^{2}$ regime, the solution for the gluon
distribution $\tilde{g}(y,\gamma_{i})$ is dominated by $\gamma_{i}\sim 0$, so
the Dirac $\delta$ imposes $\gamma\simeq 1$ in the non-linear correction term
(28). Therefore it may be approximated by

$\tilde{K}_{\mathrm{BK}}\simeq 2\pi
i\,\delta(\gamma_{1}+\gamma_{2}+1-\gamma)\,{4\gamma_{1}\gamma_{2}}.$ (29)

It is important to notice that the factor $4\gamma_{1}\gamma_{2}$ in the
integral kernel does not change the dominance of the $\gamma_{i}\sim 0$ region
in the integrand. It is because in this region the collinear gluon
distributions are strongly enhanced, they behave as:
$\tilde{g}(y,\gamma_{i})\sim\exp(-\bar{\alpha}_{\mathrm{s}}y/\gamma_{i})$, and
the prefactors $\gamma_{i}$ are subleading with respect to the exponentiated
$\gamma_{i}$ pole part. The aforementioned condition $\gamma\simeq 1$ yields
the leading $1/Q^{2}$ dependence of the non-linear correction term, as in the
case of the GLR equation, while the leading part of the gluon distribution is
given by the linear evolution and remains localized in the region of
$\gamma\simeq 0$. The integral operator corresponding to this kernel has the
same structure as the GLR correction term induced by (27), but with the
replacement: $\tilde{g}(y,\gamma_{i})\to\gamma_{i}\tilde{g}(y,\gamma_{i})$. In
the double logarithmic regime we get the following integral representation of
the linear rapidity evolution equation:

$-\gamma\tilde{g}(y,\gamma)\simeq\bar{\alpha}_{\mathrm{s}}\int^{y}dy^{\prime}g(y^{\prime},\gamma).$
(30)

This implies that the factors $\gamma_{i}\tilde{g}(y,\gamma_{i})$ that appear
in the non-linear term of Eq. (24) are suppressed by one order of
$\bar{\alpha}_{\mathrm{s}}$ in the resummations of scale logarithms. Another
way to see this is to use the saddle point solution to the gluon evolution
equation, the solution,

$xg(x,Q^{2})=\left.A(\bar{\alpha}_{\mathrm{s}}y/\log^{3}(Q^{2}))^{1/4}\exp(2\sqrt{\bar{\alpha}_{\mathrm{s}}y\log(Q^{2})})\right|_{y=\log(x_{\mathrm{i}n}/x)}.$

Multiplication by the Mellin moment $\gamma$ corresponds to the
differentiation w.r.t. $\log(Q^{2})$. It leads to lowering the power of
logarithm. Explicitly:

${\partial
xg(x,Q^{2})\over\partial\log(Q^{2})}=\left(\sqrt{\bar{\alpha}_{\mathrm{s}}y\over\log(Q^{2})}-{3\over
4\log(Q^{2})}\right)xg(x,Q^{2}).$ (31)

In the counting of leading logarithms the relative powers of
$\alpha_{\mathrm{s}}$ and $\log(k^{2})$ are important, and both the factors
$\sqrt{\bar{\alpha}_{\mathrm{s}}y/\log(Q^{2})}$ and $-3/(4\log(Q^{2}))$ lower
the power of logarithm by one w.r.t. the power of $\bar{\alpha}_{\mathrm{s}}$.
It follows that one iteration of the non-linear term from the BK equation
expressed in terms of the collinear gluon distribution, comes at the
$\alpha_{\mathrm{s}}^{4}\log(Q^{2}/Q^{2}_{0})/Q^{2}$ order, to be compared
with the GLR non-linear term,
$\sim\alpha_{\mathrm{s}}^{2}\log(Q^{2}/Q^{2}_{0})/Q^{2}$.

The last conclusion holds true for the collinear gluon distribution resulting
from the resummation of term enhanced by powers of $\log(Q^{2}/Q^{2}_{0})$,
and for $t=\log(Q^{2}/Q^{2}_{0})\gg\bar{\alpha}_{\mathrm{s}}y$. For the other
asymptotic regions: $\bar{\alpha}_{\mathrm{s}}y\gg t\gg 1$, and
$\bar{\alpha}_{\mathrm{s}}y\sim t\gg 1$ , the dominant value of the anomalous
dimension $\gamma_{s}=\sqrt{\bar{\alpha}_{\mathrm{s}}y/t}$ is not bounded to
be much smaller than one, and there is no significant value reduction in the
transition from the gluon collinear gluon distribution function
$\tilde{g}(y,\gamma)$ to the unintegrated one,
$\tilde{f}(y,\gamma)\simeq\gamma\tilde{g}(y,\gamma)$, as the loss of one power
of the scale logarithm is compensated by a similar or larger enhancement by
the factor of $\bar{\alpha}_{\mathrm{s}}y$. In these regions we recover the
logaritmic scaling of the GLR equation.

The conclusion about the strong suppression of the non-linear BK term in the
double logarithmic approximation with $t\gg\bar{\alpha}_{\mathrm{s}}y$
hierarchy can be checked by a direct analysis of this term in momentum space.
The transverse momentum integrals in Eq. (20) have lower boundary of $k^{2}$ —
it corresponds to the anticollinear ordering. Hence the non-linear term,
proportional to $\alpha_{\mathrm{s}}^{2}$ cannot produce the logarithm
$\log(k^{2}/Q_{0}^{2})$, to be contrasted with the linear term, that yields
the leading $\alpha_{\mathrm{s}}\log(k^{2}/Q_{0}^{2})$ contribution. Thus we
get the same conclusion as from the analysis in the Mellin moments space: the
non-linear correction in the BK equation enters at lower order of the
logaritmic $\log(Q^{2})$ resummation of the perturbative series, than the
corresponding term in the GLR equation. This difference may be traced back to
the logaritmic integration needed to obtain the collinear gluon distribution
from the unintegrated one. As a direct consequence, we expect that for
$t\gg\bar{\alpha}_{\mathrm{s}}y\gg 1$, the non-linear corrections from BK
equation are significantly weaker than in the GLR equation.

It should be interesting to revisit the original argument about the connection
between the GLR equation and the BK equation in the double logartithmic
approximation (DLA) given in Ref. [5], that has lead to a different conclusion
than ours. That analysis was performed in the transverse position
representation. Two key approximations were applied: (i) the dipole kernel was
approximated by the leading behavior for large daughter dipoles, $r\ll
r^{\prime},|\mathbf{r}^{\prime}-\mathbf{r}|$,

${\mathbf{r}^{2}\over{\mathbf{r}^{\prime}}^{2}(\mathbf{r}^{\prime}-\mathbf{r})^{2}}\;\longrightarrow\;\theta(r^{\prime}-r){\mathbf{r}^{2}\over{\mathbf{r}^{\prime}}^{4}},$
(32)

where the parent and daughter size dipole vectors are given by $\mathbf{r}$
and $\mathbf{r}^{\prime},\mathbf{r}^{\prime}-\mathbf{r}$ respectively,
$\theta$ is the Heaviside function. In addition, consistently,
$N(x,\mathbf{r}^{\prime}-\mathbf{r})$ was approximated by
$N(x,\mathbf{r}^{\prime})$. Furthermore, (ii): the DLA relation between the
dipole cross section and the collinear gluon distribution,

${N(x,\mathbf{r})\over r^{2}}={\alpha_{\mathrm{s}}\pi^{2}\over
2N_{\mathrm{c}}S_{\perp}}xg(x,1/r^{2}),$ (33)

was employed. This lead to the following approximate form of the BK equation:

$x{\partial\over\partial x}xg(x,1/r^{2})={\bar{\alpha}_{\mathrm{s}}\over
2}\int_{r^{2}}^{1/\Lambda^{2}}{d{r^{\prime}}^{2}\over{r^{\prime}}^{2}}\left[2xg(x,1/{r^{\prime}}^{2})-{\pi^{2}\alpha_{s}\over
2N_{c}S_{\perp}}{r^{\prime}}^{2}(xg(x,1/{r^{\prime}}^{2}))^{2}\right],$ (34)

where $\Lambda$ is a nonperturbative energy scale of QCD or alternatively the
saturation scale. From this equation the GLR equation was obtained by taking
the logartithmic derivative, $r^{2}\partial/\partial r^{2}$. Unfortunately,
step (i) is not accurate enough to provide the correct dependence on $r^{2}$
of the non-linear term. In Eq. (34) the $r^{2}$ dependence of the r.h.s. is
coming entirely from the lower end-point region of the ${r^{\prime 2}}$
integration. This is justified for the linear term in the DLA due to its
$1/{r^{\prime}}^{2}$ leading behavior in the integrand, but is not correct for
the non-linear term. It is clear when one considers the non-linear term in the
BK equation rewritten in terms of $xg(x,1/r^{2})$, using (33). It reads:

$-{\alpha_{\mathrm{s}}^{2}\over
4S_{\perp}}\int^{r^{\prime}<1/\Lambda}d^{2}\mathbf{r}^{\prime}\,xg(x,1/{r^{\prime}}^{2})\,xg(x,1/(\mathbf{r}^{\prime}-\mathbf{r})^{2}).$
(35)

When this more accurate expression is differentiated with respect to
$\log(r^{2})$, the result:

$-{\alpha_{\mathrm{s}}^{2}\over
4S_{\perp}}\int^{r^{\prime}<1/\Lambda}d^{2}\mathbf{r}^{\prime}\,xg(x,1/{r^{\prime}}^{2})\,r^{2}{\partial\over\partial
r^{2}}xg(x,1/(\mathbf{r}^{\prime}-\mathbf{r})^{2})$ (36)

receives contributions from the whole integration region of
$\mathbf{r}^{\prime}$ with no enhancement for $r^{\prime}\to r$. Hence, in
this integral all scales of the collinear gluon distributions are probed down
to $\Lambda^{2}$, with no enhancement of a particular region from the
integration kernel. So, the contributions of gluon distributions at different
scales $1/{r^{\prime}}^{2}$ are weighted only by the integration volume
effects, and they enhance the large $r^{\prime}$ region, corresponding to
small $1/{r^{\prime}}^{2}\sim\Lambda^{2}$ scales111In fact, the evolution of
gluon distributions introduces some enhancement of the large scales
$1/{r^{\prime}}^{2}$ and $1/(\mathbf{r}^{\prime}-\mathbf{r})^{2}$ scales (i.e.
the small $r^{\prime 2}$ and $(\mathbf{r}^{\prime}-\mathbf{r})^{2}$), but
these effects are enhanced only moderately by powers of logarithms and cannot
fully compensate the strong power enhancement of large $r^{\prime}$ values due
to the integration volume effects.. This is different than in the case of
analogous differentiation of approximate integral (34), where the scales in
both gluon distributions go to large $1/r^{2}$. Of course, such modification
of the scales in the gluon distributions changes the evolution lenght in the
scale and consequently, also the powers of logarithms of the hard scale
included in the final equation of the GLR type with respect to the exact BK
equation. In particular, when the scale approaches $\Lambda^{2}$, the effects
of evolution are not present and there are no contributions of the large scale
logarithms. This point, however rather subtle, is essential. The scheme that
we implemented is free from inaccuracies introduced by approximation (32).

The BK equation describes the evolution of color dipole scattering amplitude.
From this amplitude we obtain the underlying unintegrated gluon distribution
that belongs to a wider class of Transverse Momentum Distributions (TMDs), see
e.g. [50, 51]. More precisely, the distribution $f(x,k^{2})$ that we use is
called the dipole TMD (in the TMD notation: $xG^{(2)}(x,k^{2})$). In general
TMDs are defined by expectation values of gauge link contours, that correspond
to parton configurations, and therefore depend on the scattering process. An
important TMD is called the Weizsäcker–Williams distribution
($xG^{(1)}(x,k^{2})$ in the TMD notation), relevant for many physical
processes, e.g. for dijet production in DIS. The Weizsäcker–Williams and
dipole gluon TMDs have the same behavior for $k^{2}\gg Q_{\mathrm{s}}^{2}$,
but differ dramatically for $k^{2}\sim Q_{\mathrm{s}}^{2}$ or
$k^{2}<Q_{\mathrm{s}}^{2}$. They also obey different evolution equations:
while the dipole gluon TMD is governed by the BK equation, in the evolution of
the Weizsäcker–Williams gluon TMD the quadrupoles play an important role and
the evolution is more complicated [52]. Hence, in general, one may expect that
the non-linear evolution equation of the integrated versions of different TMDs
are also different. The present analysis is focused only on the dipole gluon
TMD and the conclusions may be not applicable to other gluon TMDs and their
integrated counterparts. Let us add that in the McLerran–Venugopalan model
[15, 16] the the Weizsäcker–Williams TMD is directly related to
$\phi(x,k^{2})$ (see e.g. [53]), and in the double logarithmic approximation
the evolution equation of $\phi(x,k^{2})$ takes the GLR form. The equivalence
of $xG^{(1)}(x,k^{2})$ and $\phi(x,k^{2})$, however, does not hold when the
QCD evolution effects are taken into account.

###  5.2 Effects of the high gluon density regime

Now we turn to effects of the non-linearity when the gluon density is large,
in particular to the gluon saturation regime. We shall perform a heuristic
analysis of the impact of saturation effects. We choose a simple model of
$f(x,k^{2})$, that provides clear analytic insight. The key feature of the
model is presence of $x$-dependent saturation scale, that plays the role of
lower cutoff on $k$, separating the linear evolution region from the region
with strong suppression effects due to unitarity corrections, which impose
fundamental constraints on the dipole scattering amplitude. In addition, the
geometric scaling property of the dipole cross section is assumed to hold. The
geometric scaling was initially discovered in the HERA data for the total
$\gamma^{*}p$ cross section [54], and it holds with a good accuracy for color
dipole cross section obtained from the BK equation at larger rapidities [55,
56, 57]. The scaling is equivalent to universality of shape of the dipole
cross section at different rapidities, and is also known as the “travelling
wave” phenomenon [58, 59]. The physics picture and generic features of the
obtained results do not depend on details of the shape of $f(x,k^{2})$. The
dipole scattering amplitude obtained from the simple model we assume has the
geometric scaling built in, and it is consistent with the unitatity
constraints. The way the amplitude approaches the unitarity limit is slightly
different from the Levin–Tuchin law [55] that follows from the BK equation,
but this difference does not affect the conclusions of the analysis. A similar
reasoning lead to an proposal to impose unitarity effects on QCD evolution at
small $x$ as an absorptive boundary [60, 61], which is not sensitive to a
particular way to approach the unitarity limit. The geometric scaling implies
the following form of the unintegrated gluon density

$f(x,k^{2})=S_{\perp}Q_{\mathrm{s}}^{2}(x)\,h(k^{2}/Q_{\mathrm{s}}^{2}(x)),$
(37)

where $Q_{\mathrm{s}}(x)$ is the $x$-dependent saturation scale, and the
function $h$ is the universal profile of the BK solution. BK equation leads to
an approximate power dependence of the saturation scale,
$Q_{\mathrm{s}}^{2}(x)\simeq Q_{0}^{2}(x_{\mathrm{in}}/x)^{\lambda}$, with
$\lambda\simeq 0.3$. Unitarity of the color dipole cross section scattering of
a very dense target implies the asymptotic behavior

$h(\xi)\simeq\xi^{2}\quad\mbox{for}\;\;\xi\to 0,$ (38)

that corresponds to $f(y,k^{2})\sim k^{4}$ for $k^{2}\ll
Q_{\mathrm{s}}^{2}(y)$ [57]. For $\xi\gg 1$ the behavior of $h(\xi)$ is driven
mostly by the linear evolution. For the purpose of this analysis it is
sufficient to approximate the large $\xi$ behavior of $h(\xi)$ by
$\xi^{\gamma_{c}}$, where $\gamma_{c}$ is a positive number, much smaller than
1, related to the anomalous dimension of the gluon distribution function. The
simplest model of $h(\xi)$ that incorporates both the features is

$h(\xi)=A[\xi^{2}\theta(1-\xi)+\xi^{\gamma_{c}}\theta(\xi-1)],$ (39)

where $\theta$ is the Heaviside function and $A$ is a numerical constant.

We apply this model to estimate the effect of the non-linear term in the BK
equation on the collinear gluon distribution for various hierarchy of scales.
The regime of $Q^{2}\gg Q^{2}_{s}(y)$ was studied above in the double
logarithmic limit. Using the model of the BK solution for $f(y,k^{2})$ we get

$xg(x,Q^{2})\simeq\int^{Q^{2}}{dk^{2}\over
k^{2}}f(x,k^{2})=A(x_{\mathrm{in}}/x)^{\lambda}\left[\int^{Q_{\mathrm{s}}^{2}(x)}{k^{2}dk^{2}\over
Q_{\mathrm{s}}^{4}(x)}+\int_{Q_{\mathrm{s}}^{2}(x)}^{Q^{2}}{dk^{2}\over
k^{2}}(k^{2}/Q_{\mathrm{s}}^{2}(x))^{\gamma_{c}}\right]$
$=A(x_{\mathrm{in}}/x)^{\lambda}\,\left[1/2+{(Q^{2}/Q_{\mathrm{s}}^{2}(x))^{\gamma_{c}}-1\over\gamma_{c}}\right].$
(40)

Applying the expansion in $\gamma_{c}$ around zero up to the first order, and
keeping only the leading logarithmic term we get

$xg(x,Q^{2})\simeq(x_{\mathrm{in}}/x)^{\lambda}\log(Q^{2}/Q_{\mathrm{s}}^{2}(x)).$
(41)

In the absence of non-linear correction the saturation scale in
$\log(Q^{2}/Q_{\mathrm{s}}^{2}(x))$ should be replaced by a much smaller scale
$\mu_{0}\ll Q_{\mathrm{s}}(x)$, giving
$xg(x,Q^{2})|_{\mathrm{linear}}\simeq(x_{\mathrm{in}}/x)^{\lambda}\log(Q^{2}/\mu_{0}^{2})$.
Hence the relative correction due to non-linearity reads

${xg(x,Q^{2})-xg(x,Q^{2})|_{\mathrm{linear}}\over
xg(x,Q^{2})|_{\mathrm{linear}}}\simeq{\log(Q^{2}/Q_{\mathrm{s}}^{2}(x))-\log(Q^{2}/\mu_{0}^{2})\over\log(Q^{2}/\mu_{0}^{2})}=-{\log(Q_{\mathrm{s}}^{2}(x)/\mu_{0}^{2})\over\log(Q^{2}/\mu_{0}^{2})}\,.$
(42)

This correction enters without a suppressing power factor of $1/Q^{2}$, hence
at the leading twist, and due to logarithmic dependencies on the scales it is
not small. This is consistent with our findings of strong non-linear
correction in the proton structure functions at twist 2. The presented
estimate is rather crude, but it clearly shows how the non-linear corrections
contribute to twist observables. It happens because the gluon recombination /
unitarity leads to a strong suppression of unintegrated gluon distribution
$f(x,k^{2})$ in the region $k^{2}<Q_{\mathrm{s}}^{2}(x)$, and this imposes an
effective lower cut-off on the logaritmic integration in
$xg(x,Q^{2})=\int_{Q_{\mathrm{s}}^{2}(x)}dk^{2}/k^{2}f(x,k^{2})$.

In these consideration we assumed that the saturation scale
$Q_{\mathrm{s}}\gg\mu_{0}$, where $\mu_{0}$ should be interpreted as an
intrinsic hadronic scale of the proton, for instance it could be related to
the inverse proton size. This is not in the perturbative domain, but this does
not endanger the conclusions as $\mu_{0}$ enters only as a lower cut-off of
logarithmic integrations. It should be clearly distinguished from an initial
scale of the DGLAP evolution $\mu_{\mathrm{F,0}}$, which is typically set to
be greater than 1 GeV, and greater than the saturation scale in the proton.
For the case $\mu_{\mathrm{F,0}}\gg Q_{\mathrm{s}}(x)$ the leading twist non-
linear correction enter mostly as the input of the DGLAP evolution, with weak
correction terms, as described in Sec. 5.1.

###  5.3 The intermediate region

Above we discussed the non-linear effects for $Q^{2}\gg Q_{\mathrm{s}}^{2}(x)$
and in the high gluon density regime. The intermediate region of $Q^{2}\sim
Q_{\mathrm{s}}^{2}(x)$ is hardest to analyze analytically, as in this region
$Q_{\mathrm{s}}^{2}(x)/Q^{2}\simeq 1$ and the linear and non-linear effects
have similar size. Also, the relevant scale logarithms are of order 1, and the
double logarithmic approximation is expected to have a very limited accuracy
there. Hence, we believe that for a reliable numerical predictions in this
region should be obtained within a complete non-linear evolution framework as
given e.g. by the complete BK equation or the JIMWLK equation. Measurements of
proton structure functions from HERA and their phenomenological analysis
suggest, that the bulk of HERA data, say for $Q^{2}>2$ GeV2 are driven by the
linear evolution regime, as the saturation scale at HERA is below 1 GeV. The
situation is expected to be different for deep inelastic scattering on heavy
nuclei, that is going to be studied at the Electron–Ion Collider and in the
possible future Large Electron–Hadron Collider (LHeC) [62, 63]. For a nucleus
with mass number $A$ one expects the saturation scale $Q_{\mathrm{s}}^{2}$ to
be enhanced by a factor of $A^{1/3}$ with respect to the proton case. Moreover
at the LHeC one expects to probe the range of $x$ that may extend down to
$2\cdot 10^{-6}$ at $Q^{2}=1$ GeV2, see e.g. [64]. Hence one may reach
$Q_{\mathrm{s}}^{2}$ of few GeV2 and there are good prospects to probe the
non-linear effects well in the perturbative regime — at the lowest $x$ the
LHeC data on the proton and nuclear structure functions should have
$Q^{2}<Q_{\mathrm{s}}^{2}$. In this region one expects that the dominance of
twist 2 does not hold and the DGLAP framework may be not sufficient to
describe the data. Perhaps an even more remarkable scenario may be realized —
given the fact that the twist expansion in our approach leads to an asymptotic
series with the expansion parameter $\sim Q_{\mathrm{s}}^{2}/Q^{2}$, one may
expect that the twist expansion ceases to make sense for
$Q^{2}<Q_{\mathrm{s}}^{2}$. Therefore it should be extremely interesting to
probe this kinematic region experimentally and investigate the theoretical
side within the framework of non-linear evolution equations. The problem of
potential impact of non-linear effects on the proton and nucleus structure
functions at LHeC was addressed in several studies [64, 65, 66, 67]. In
particular, in a recent analysis [67] careful matching of the DGLAP and BK
description was perfomed for $(x,Q^{2})$ range where both approaches are
valid, and the differences between the predictions at small $x$ and moderate
$Q^{2}$ were studied. On a general level, the numerical results found in Ref.
[67] have a similar pattern to the one following from our study, which,
however, is more focused on identifying the structural features, and not yet
on precision fits to the data.

## Acknowledgements

This research was supported by the Polish National Science Centre (NCN) grant
no. 2017/27/B/ST2/02755.

## References

  * [1] A. Accardi et al., _Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all_ , _Eur. Phys. J. A_ 52 (2016) 268, [1212.1701].
  * [2] 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].
  * [3] J. Bartels, _Unitarity corrections to the Lipatov pomeron and the small x region in deep inelastic scattering in QCD_ , _Phys. Lett. B_ 298 (1993) 204–210.
  * [4] J. Bartels, _Unitarity corrections to the Lipatov pomeron and the four gluon operator in deep inelastic scattering in QCD_ , _Z. Phys. C_ 60 (1993) 471–488.
  * [5] Y. V. Kovchegov, _Small x F(2) structure function of a nucleus including multiple pomeron exchanges_ , _Phys. Rev. D_ 60 (1999) 034008, [hep-ph/9901281].
  * [6] Y. V. Kovchegov, _Unitarization of the BFKL pomeron on a nucleus_ , _Phys. Rev. D_ 61 (2000) 074018, [hep-ph/9905214].
  * [7] J. Bartels and C. Ewerz, _Unitarity corrections in high-energy QCD_ , _JHEP_ 09 (1999) 026, [hep-ph/9908454].
  * [8] L. V. Gribov, E. M. Levin and M. G. Ryskin, _Semihard Processes in QCD_ , _Phys. Rept._ 100 (1983) 1–150.
  * [9] A. H. Mueller and J.-w. Qiu, _Gluon Recombination and Shadowing at Small Values of x_ , _Nucl. Phys. B_ 268 (1986) 427–452.
  * [10] E. A. Kuraev, L. N. Lipatov and V. S. Fadin, _The Pomeranchuk Singularity in Nonabelian Gauge Theories_ , _Sov. Phys. JETP_ 45 (1977) 199–204.
  * [11] I. I. Balitsky and L. N. Lipatov, _The Pomeranchuk Singularity in Quantum Chromodynamics_ , _Sov. J. Nucl. Phys._ 28 (1978) 822–829.
  * [12] L. N. Lipatov, _Small x physics in perturbative QCD_ , _Phys. Rept._ 286 (1997) 131–198, [hep-ph/9610276].
  * [13] V. S. Fadin and L. N. Lipatov, _BFKL pomeron in the next-to-leading approximation_ , _Phys. Lett. B_ 429 (1998) 127–134, [hep-ph/9802290].
  * [14] I. Balitsky, _Operator expansion for high-energy scattering_ , _Nucl. Phys. B_ 463 (1996) 99–160, [hep-ph/9509348].
  * [15] L. D. McLerran and R. Venugopalan, _Computing quark and gluon distribution functions for very large nuclei_ , _Phys. Rev. D_ 49 (1994) 2233–2241, [hep-ph/9309289].
  * [16] L. D. McLerran and R. Venugopalan, _Gluon distribution functions for very large nuclei at small transverse momentum_ , _Phys. Rev. D_ 49 (1994) 3352–3355, [hep-ph/9311205].
  * [17] 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].
  * [18] J. Jalilian-Marian, A. Kovner and H. Weigert, _The Wilson renormalization group for low x physics: Gluon evolution at finite parton density_ , _Phys. Rev. D_ 59 (1998) 014015, [hep-ph/9709432].
  * [19] E. Iancu, A. Leonidov and L. D. McLerran, _Nonlinear gluon evolution in the color glass condensate. 1._ , _Nucl. Phys. A_ 692 (2001) 583–645, [hep-ph/0011241].
  * [20] E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran, _Nonlinear gluon evolution in the color glass condensate. 2._ , _Nucl. Phys. A_ 703 (2002) 489–538, [hep-ph/0109115].
  * [21] E. Iancu, A. Leonidov and L. D. McLerran, _The Renormalization group equation for the color glass condensate_ , _Phys. Lett. B_ 510 (2001) 133–144, [hep-ph/0102009].
  * [22] A. Kovner, M. Lublinsky and Y. Mulian, _Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov, Kovner evolution at next to leading order_ , _Phys. Rev. D_ 89 (2014) 061704, [1310.0378].
  * [23] K. J. Golec-Biernat and M. Wusthoff, _Saturation effects in deep inelastic scattering at low Q**2 and its implications on diffraction_ , _Phys. Rev. D_ 59 (1998) 014017, [hep-ph/9807513].
  * [24] K. J. Golec-Biernat and M. Wusthoff, _Saturation in diffractive deep inelastic scattering_ , _Phys. Rev. D_ 60 (1999) 114023, [hep-ph/9903358].
  * [25] J. Bartels, K. J. Golec-Biernat and H. Kowalski, _A modification of the saturation model: DGLAP evolution_ , _Phys. Rev. D_ 66 (2002) 014001, [hep-ph/0203258].
  * [26] H. Kowalski and D. Teaney, _An Impact parameter dipole saturation model_ , _Phys. Rev. D_ 68 (2003) 114005, [hep-ph/0304189].
  * [27] H. Kowalski, L. Motyka and G. Watt, _Exclusive diffractive processes at HERA within the dipole picture_ , _Phys. Rev. D_ 74 (2006) 074016, [hep-ph/0606272].
  * [28] K. Kutak and J. Kwieciński, _Screening effects in the ultrahigh-energy neutrino interactions_ , _Eur. Phys. J. C_ 29 (2003) 521, [hep-ph/0303209].
  * [29] K. Kutak and A. M. Staśto, _Unintegrated gluon distribution from modified BK equation_ , _Eur. Phys. J. C_ 41 (2005) 343–351, [hep-ph/0408117].
  * [30] J. L. Albacete, N. Armesto, J. G. Milhano and C. A. Salgado, _Non-linear QCD meets data: A Global analysis of lepton-proton scattering with running coupling BK evolution_ , _Phys. Rev. D_ 80 (2009) 034031, [0902.1112].
  * [31] E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez and D. N. Triantafyllopoulos, _Collinearly-improved BK evolution meets the HERA data_ , _Phys. Lett. B_ 750 (2015) 643–652, [1507.03651].
  * [32] G. Beuf, H. Hänninen, T. Lappi and H. Mäntysaari, _Color Glass Condensate at next-to-leading order meets HERA data_ , _Phys. Rev. D_ 102 (2020) 074028, [2007.01645].
  * [33] A. P. Bukhvostov, G. V. Frolov, L. N. Lipatov and E. A. Kuraev, _Evolution equations for quasi-partonic operators_ , _Nucl. Phys. B_ 258 (1985) 601–646.
  * [34] J. Bartels and M. G. Ryskin, _The Analytic structure of the anomalous dimension of the four gluon operator in deep inelastic scattering_ , _Z. Phys. C_ 62 (1994) 425–430.
  * [35] J. Bartels and C. Bontus, _An Estimate of twist four contributions at small x (B) and low Q**2_ , _Phys. Rev. D_ 61 (2000) 034009, [hep-ph/9906308].
  * [36] 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) 580, [1506.06042].
  * [37] L. Motyka, M. Sadzikowski and W. Słomiński, _Evidence of strong higher twist effects in diffractive DIS at HERA at moderate $Q^{2}$_, _Phys. Rev. D_ 86 (2012) 111501, [1203.5461].
  * [38] L. A. Harland-Lang, A. D. Martin, P. Motylinski and R. S. Thorne, _The impact of the final HERA combined data on PDFs obtained from a global fit_ , _Eur. Phys. J. C_ 76 (2016) 186, [1601.03413].
  * [39] I. Abt, A. M. Cooper-Sarkar, B. Foster, V. Myronenko, K. Wichmann and M. Wing, _Study of HERA ep data at low Q 2 and low $x_{Bj}$ and the need for higher-twist corrections to standard perturbative QCD fits_, _Phys. Rev. D_ 94 (2016) 034032, [1604.02299].
  * [40] L. Motyka, M. Sadzikowski, W. Słomiński and K. Wichmann, _Evidence of quasi-partonic higher-twist effects in deep inelastic scattering at HERA at moderate $Q^{2}$_, _Eur. Phys. J. C_ 78 (2018) 80, [1707.05992].
  * [41] J. Bartels, K. J. Golec-Biernat and K. Peters, _An Estimate of higher twist at small x(B) and low Q**2 based upon a saturation model_ , _Eur. Phys. J. C_ 17 (2000) 121–128, [hep-ph/0003042].
  * [42] J. Bartels, K. Golec-Biernat and L. Motyka, _Twist expansion of the nucleon structure functions, F(2) and F(L), in the DGLAP improved saturation model_ , _Phys. Rev. D_ 81 (2010) 054017, [0911.1935].
  * [43] L. Motyka and M. Sadzikowski, _Twist decomposition of proton structure from BFKL and BK amplitudes_ , _Acta Phys. Polon. B_ 45 (2014) 2079, [1411.7774].
  * [44] J. Bartels and K. Kutak, _A Momentum Space Analysis of the Triple Pomeron Vertex in pQCD_ , _Eur. Phys. J. C_ 53 (2008) 533–548, [0710.3060].
  * [45] N. N. Nikolaev and B. G. Zakharov, _Color transparency and scaling properties of nuclear shadowing in deep inelastic scattering_ , _Z. Phys. C_ 49 (1991) 607–618.
  * [46] S. J. Brodsky, G. P. Lepage and P. B. Mackenzie, _On the Elimination of Scale Ambiguities in Perturbative Quantum Chromodynamics_ , _Phys. Rev. D_ 28 (1983) 228.
  * [47] S. J. Brodsky, V. S. Fadin, V. T. Kim, L. N. Lipatov and G. B. Pivovarov, _The QCD pomeron with optimal renormalization_ , _JETP Lett._ 70 (1999) 155–160, [hep-ph/9901229].
  * [48] S. Bondarenko and L. Motyka, _Solving effective field theory of interacting QCD pomerons in the semi-classical approximation_ , _Phys. Rev. D_ 75 (2007) 114015, [hep-ph/0605185].
  * [49] Y. V. Kovchegov and E. Levin, _Quantum Chromodynamics at High Energy_ , vol. 33.  Oxford University Press, 2013, 10.1017/9781009291446.
  * [50] R. Angeles-Martinez et al., _Transverse Momentum Dependent (TMD) parton distribution functions: status and prospects_ , _Acta Phys. Polon. B_ 46 (2015) 2501–2534, [1507.05267].
  * [51] A. van Hameren, P. Kotko, K. Kutak, C. Marquet, E. Petreska and S. Sapeta, _Forward di-jet production in p+Pb collisions in the small-x improved TMD factorization framework_ , _JHEP_ 12 (2016) 034, [1607.03121].
  * [52] F. Dominguez, A. H. Mueller, S. Munier and B.-W. Xiao, _On the small- $x$ evolution of the color quadrupole and the Weizsäcker–Williams gluon distribution_, _Phys. Lett. B_ 705 (2011) 106–111, [1108.1752].
  * [53] F. Dominguez, C. Marquet, B.-W. Xiao and F. Yuan, _Universality of Unintegrated Gluon Distributions at small x_ , _Phys. Rev. D_ 83 (2011) 105005, [1101.0715].
  * [54] A. M. Staśto, K. J. Golec-Biernat and J. Kwieciński, _Geometric scaling for the total gamma* p cross-section in the low x region_ , _Phys. Rev. Lett._ 86 (2001) 596–599, [hep-ph/0007192].
  * [55] E. Levin and K. Tuchin, _Solution to the evolution equation for high parton density QCD_ , _Nucl. Phys. B_ 573 (2000) 833–852, [hep-ph/9908317].
  * [56] E. Levin and K. Tuchin, _New scaling at high-energy DIS_ , _Nucl. Phys. A_ 691 (2001) 779–790, [hep-ph/0012167].
  * [57] K. J. Golec-Biernat, L. Motyka and A. M. Staśto, _Diffusion into infrared and unitarization of the BFKL pomeron_ , _Phys. Rev. D_ 65 (2002) 074037, [hep-ph/0110325].
  * [58] S. Munier and R. B. Peschanski, _Geometric scaling as traveling waves_ , _Phys. Rev. Lett._ 91 (2003) 232001, [hep-ph/0309177].
  * [59] S. Munier and R. B. Peschanski, _Traveling wave fronts and the transition to saturation_ , _Phys. Rev. D_ 69 (2004) 034008, [hep-ph/0310357].
  * [60] A. H. Mueller and D. N. Triantafyllopoulos, _The Energy dependence of the saturation momentum_ , _Nucl. Phys. B_ 640 (2002) 331–350, [hep-ph/0205167].
  * [61] D. N. Triantafyllopoulos, _The Energy dependence of the saturation momentum from RG improved BFKL evolution_ , _Nucl. Phys. B_ 648 (2003) 293–316, [hep-ph/0209121].
  * [62] 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].
  * [63] LHeC, FCC-he Study Group collaboration, P. Agostini et al., _The Large Hadron–Electron Collider at the HL-LHC_ , _J. Phys. G_ 48 (2021) 110501, [2007.14491].
  * [64] J. Rojo and F. Caola, _Parton distributions and small-x QCD at the Large Hadron Electron Collider_ , in _17th International Workshop on Deep-Inelastic Scattering and Related Subjects_ , (Berlin, Germany), p. 222, Science Wise Publ., 2009.  0906.2079.
  * [65] A. Accardi, V. Guzey and J. Rojo, _Nuclear parton distributions and deviations from DGLAP at an Electron Ion Collider_ , 1106.3839.
  * [66] C. Marquet, M. R. Moldes and P. Zurita, _Unveiling saturation effects from nuclear structure function measurements at the EIC_ , _Phys. Lett. B_ 772 (2017) 607–614, [1702.00839].
  * [67] N. Armesto, T. Lappi, H. Mäntysaari, H. Paukkunen and M. Tevio, _Signatures of gluon saturation from structure-function measurements_ , _Phys. Rev. D_ 105 (2022) 114017, [2203.05846].