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

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 F2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_2$$\end{document} and FL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{\textrm{L}}$$\end{document} 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.


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 −λ with λ > 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 τ contribution to hadron structure function is (Λ/Q) τ −2 , where Λ 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 ∼ [xg(x, Q 2 )] 2 Λ 2 /Q 2 see e.g.[8,9,33,34,35] to be compared with ∼ 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 ∼ 2 • 10 −5 for Q 2 > 1 GeV 2 deteriorate for Q 2 < 5 GeV 2 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 GeV 2 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 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 ≫ ᾱs y ≫ 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 ∼ ᾱs y ≫ 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 ) ≃ k 2 (∂/∂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 2 s (x), where the characteristic x-dependent momentum scale, that is naturally interpreted as the saturation scale.In the region k 2 < Q 2 s (x), the distri-bution f (x, k 2 ) is strongly suppressed.Hence, Q 2 s (x) plays the rôle of a lower cut-off in the integral , 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, µ 2 F ), provided that the factorization scale µ 2 F is reasonably larger than the saturation scale Q 2 s (x).Hence, in order to clearly see non-linear BK effects one should probe the region of µ 2 F < Q 2 s .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.

Proton structure functions from Balitsky-Kovchegov equation
The total inclusive cross-section σ γ * A of the virtual photon γ * scattering on large nucleus A, in the high energy limit, can be described in terms of the structure functions where The transverse (T ) and longitudinal (L) photon wavefunctions ψ 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 = |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 σ q q given by an imaginary part of the forward dipole -nucleon scattering amplitude N (x, r, b) where y = log(x in /x) is a rapidity variable developed with respect to some initial value x in [5].The cross section parameter σ 0 is related to the effective nucleus radius σ 0 = 2πR 2 A via integration over the impact parameter 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 ⃗ b, limiting only to a simple cylindrical geometry of the high energy nucleus.
The rapidity evolution of the amplitude N (y, r) is described by the Balitsky-Kovchegov equation [14,5], which in momentum space reads where ᾱs = α s N c /π and ϕ(y, k 2 ⊥ ) is the Fourier transform of the amplitude The first term in the BK equation is given by the BFKL kernel that can be solved exactly using its eigenfunctions [12] ϕ 0 (y, where χ(n, γ) are eigenvalues of the BK equation and coefficients C n (γ) are determined by the initial condition.Note that we use the convention in which the standard Mellin moment corresponds to −γ, and the reference scale in the Mellin transform is Q 2 0 .The fundamental Mellin strip is located in the interval c ∈ (0, 1).At large rapidity y, which is of main interest in this article, the dominant contribution is given by n = 0 eigenvalue 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 and the functions φi satisfy the infinite set of the coupled equations: From the structure of the equations one can infer that φi ∼ φi+1 0 , where the term φ0 ∼ exp(αχy).Therefore the amplitude φi describes the process with i + 1 pomeron exchanges included in the DIS diagram.The series ( 8) is convergent if φ0 (y, k 2 ⊥ ) ≪ 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 whereas the first equation gives the BFKL solution for n = 0 eigenvalue.For further analysis we adopt the exponential form of the initial conditions Figure 1: Location of singularities in a complex Mellin space of γ variable.Blue points (on the horizontal axis) correspond to singularities located at integer values, whereas the red ones (above the axis) to integer γ − γ 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 γ = 2.
which gives and the BFKL solution takes the form whereas the Mellin strip is limited to 0 < Re(c) < 3/4.Note that we changed convention for the Q 0 parameter of the input fuction N (y = 0, r) w.r.t.Ref. [43]: the present Q 0 equals 1/2 of Q 0 used in [43].

Twist decomposition of the DIS cross section
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 nonlinear interaction term, as described in the previous section.In this framework the cross section can be described as a series , where The Mellin fundamental strip is located in −3/4 < Re(c) < 0 and φi are solutions of the equations ( 9).The functions HT,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 and was described in [43].The lowest order correction σ to the BFKL cross section follows from the solution (10) and decomposition ( 14) where 0 < Re(c), 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 γ, γ 1 and γ −γ 1 (see Fig. 1), therefore the integration over γ variable can be decompose into two sums where C w is a small clockwise contour located around point w in the complex γ plane.The Cauchy theorem that bridges equations (16,17) is satisfied due to the exponential suppression of the large imaginary values Imγ brought by the photon wavefunctions H T,L .The asymptotics of χ(γ) does not spoil that in any way, as it consists of two pieces: the − log(γ) term that improves the convergence and the −π cot(πγ) that contributes to the (only) poles of χ(γ) at the integer values of γ, 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 GeV 2 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γ over C n from the first term in (17) gives a direct contribution to twist τ = 2n after integration over dγ 1 along the line parallel to imaginary axis within the Mellin strip.A similar integration over C n+γ 1 from the second term gives contributions to all twists of order τ ≥ 2n starting with the lowest value τ = 4.This fact follows from the integration over dγ 1 variable.However, the twists higher then τ = 4 are strongly suppressed, therefore one can assign the contribution of the second integral to twist τ = 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 τ = 4 by the first term with n = 2 and the second term with n = 1.
GeV 2 in twist 2 approximation.We show the twist 2 components obtained from the BFKL evolution, the BK correction and their sum.

Numerical results
With the framework described above we calculate higher twist corrections to the unpolorized proton structure functions F 2 and F 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 in = 0.1 is assumed to take the GBW form: with σ 0 = 17.04 mb and Q 0 = 0.255 GeV.The value of the strong coupling constant ᾱs in the BFKL kernel is set to 0.087.This should be understood as an effective value of ᾱ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 L (x, Q 2 ).We choose the DIS reference scale Q 2 = 5 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 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.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.
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 τ = 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, τ > 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 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 L generates a scale logarithm for twist 4, and does not for twist 2, while for F 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 nonlinearity 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 GeV 2 and for Q 2 = 10 GeV 2 and found that the general pattern is very similar to the case Q 2 = 5 GeV 2 .Therefore, we do not depict them in separate x = 0.01 Table 1: Contribution of the BK corrections to structure functions F 2 , F L at the leading twist.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 ∆F 2,BK , ∆F L,BK compared to the BFKL results in the leading twist for F 2 and F L structure functions.
In Table 2 one can find the ratio of twist 4 and twist 2 for F 2 (R F 2 ) and F L (R F L ) structure functions: R F 2 = (F L,BK ).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 GeV 2 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.

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 2 ᾱs y log(Q 2 /Q 2 0 ) − 2α s y (up to power factors of y), to be compared with twist 2 BFKL in the same approximation ∼ exp 2 ᾱs y log(Q 2 /Q 2 0 .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 ∼ exp 2 ᾱs y log(Q 2 /Q 2 0 ) 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 γ * scattering.Therefore one expects the strong ∼ exp 2 ᾱs y log(Q 2 /Q 2 0 ) 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 µ F,0 is big enough.Clearly, the scale for the BK effects is the saturation scale Q s , so we conclude that with µ F,0 ≫ Q s the DGLAP description should not be significantly affected by non-linear evolution effects.This conclusion may change when µ F,0 < Q 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 2 s < 1 GeV 2 .Hence, with DGLAP fits assuming typically µ 2 F,0 = 2 GeV 2 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 2 s is enhanced by A 1/3 w.r.t. the saturation scale in proton.In order to account for the possibility of Q 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

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, In what follows we shall also use notation f (x, k 2 ) → f (y, k 2 ) with y = log(x in /x).We keep the convention for the Mellin transform used in the previous sections: The unintegrated gluon distribution is related to the dipole amplitude through the formula [48] f where the function ϕ(y, k 2 ) is the Fourier transform of N (y, r)/r 2 , see Eq. ( 5).The BK equation for the unintegrated gluon distribution reads [28,48]: where S ⊥ is the transverse target area, for a uniform target with radius R A , S ⊥ = πR 2 A .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: (21) which, using the symmetry between γ 1 and γ 2 can be rewritten as This implies the evolution equation for collinear gluon distribution: ∂g(y, γ) ∂y = ᾱs χ(−γ)g(y, γ) where g(y, is the collinear gluon distribution xg(x, Q 2 ) in the Mellin representation.Note that the relation f (x, k 2 ) = k 2 ∂ k 2 xg(x, k 2 ) leads to the f (y, γ) = (−γ)g(y, γ) relation in the (y, γ) variables.In the double logaritmic limit, which corresponds to the leading powers of γ around γ = 0, we obtain: where we approximated χ(−γ) ≃ −1/γ + O(1) around γ → 0. Before further analysis of this equation let us compare it to the GLR equation [8,49]: which in the Mellin representation (y, γ) takes the form: with C F = (N 2 c − 1)/2N 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 2 0 which carries the physical dimension.We find differences between the non-linear equation (24) and the GLR equation.The color prefactor of GLR N c /2C F = N 2 c /(N 2 c − 1) differs from our by O(1/N 2 c ) and this is beyond the leading N c accuracy of the BK equation.There is, however, a more important difference in the integral kernel in the Mellin space.GLR gives while we obtain: 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 2 0 ) that emerge from both the kernels.At first, let us focus on the limit t = log(Q 2 /Q 2 0 ) ≫ ᾱs y ≫ 1.In this region the saddle point solution of the linear evolution equation is dominated by the anomalous dimension γ s = ᾱs y/t ≪ 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 g(y, γ i ) is dominated by γ i ∼ 0, so the Dirac δ imposes γ ≃ 1 in the non-linear correction term (28).Therefore it may be approximated by It is important to notice that the factor 4γ 1 γ 2 in the integral kernel does not change the dominance of the γ i ∼ 0 region in the integrand.It is because in this region the collinear gluon distributions are strongly enhanced, they behave as: g(y, γ i ) ∼ exp(−ᾱ s y/γ i ), and the prefactors γ i are subleading with respect to the exponentiated γ i pole part.The aforementioned condition γ ≃ 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 γ ≃ 0. The integral operator corresponding to this kernel has the same structure as the GLR correction term induced by ( 27), but with the replacement: g(y, γ i ) → γ i g(y, γ i ).In the double logarithmic regime we get the following integral representation of the linear rapidity evolution equation: This implies that the factors γ i g(y, γ i ) that appear in the non-linear term of Eq. ( 24) are suppressed by one order of ᾱ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, .
Multiplication by the Mellin moment γ corresponds to the differentiation w.r.t.log(Q 2 ).It leads to lowering the power of logarithm.Explicitly: In the counting of leading logarithms the relative powers of α s and log(k 2 ) are important, and both the factors ᾱs y/ log(Q 2 ) and −3/(4 log(Q 2 )) lower the power of logarithm by one w.r.t. the power of ᾱ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 α 4 s log(Q 2 /Q 2 0 )/Q 2 order, to be compared with the GLR non-linear term, 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 ) ≫ ᾱs y.For the other asymptotic regions: ᾱs y ≫ t ≫ 1, and ᾱs y ∼ t ≫ 1 , the dominant value of the anomalous dimension γ s = ᾱ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 g(y, γ) to the unintegrated one, f (y, γ) ≃ γg(y, γ), as the loss of one power of the scale logarithm is compensated by a similar or larger enhancement by the factor of ᾱ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 ≫ ᾱ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 α 2 s cannot produce the logarithm log(k 2 /Q 2 0 ), to be contrasted with the linear term, that yields the leading α s log(k 2 /Q 2 0 ) 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 ≫ ᾱs y ≫ 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 ≪ r ′ , |r ′ − r|, where the parent and daughter size dipole vectors are given by r and r ′ , r ′ − r respectively, θ is the Heaviside function.In addition, consistently, N (x, r ′ − r) was approximated by N (x, r ′ ).Furthermore, (ii): the DLA relation between the dipole cross section and the collinear gluon distribution, was employed.This lead to the following approximate form of the BK equation: where Λ 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 ∂/∂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 ′2 integration.This is justified for the linear term in the DLA due to its 1/r ′ 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: When this more accurate expression is differentiated with respect to log(r2 ), the result: receives contributions from the whole integration region of r ′ with no enhancement for r ′ → r.Hence, in this integral all scales of the collinear gluon distributions are probed down to Λ 2 , with no enhancement of a particular region from the integration kernel.So, the contributions of gluon distributions at different scales 1/r ′ 2 are weighted only by the integration volume effects, and they enhance the large r ′ region, corresponding to small 1/r ′ 2 ∼ Λ 2 scales1 .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 Λ 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 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 ϕ(x, k 2 ) (see e.g.[53]), and in the double logarithmic approximation the evolution equation of ϕ(x, k 2 ) takes the GLR form.The equivalence of xG (1) (x, k 2 ) and ϕ(x, k 2 ), however, does not hold when the QCD evolution effects are taken into account.

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 γ * 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 where Q 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 2 s (x) ≃ Q 2 0 (x in /x) λ , with λ ≃ 0.3.Unitarity of the color dipole cross section scattering of a very dense target implies the asymptotic behavior h(ξ) ≃ ξ 2 for ξ → 0, that corresponds to f (y, k 2 ) ∼ k 4 for k 2 ≪ Q 2 s (y) [57].For ξ ≫ 1 the behavior of h(ξ) is driven mostly by the linear evolution.For the purpose of this analysis it is sufficient to approximate the large ξ behavior of h(ξ) by ξ γc , where γ c is a positive number, much smaller than 1, related to the anomalous dimension of the gluon distribution function.The simplest model of h(ξ) that incorporates both the features is where θ 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 ≫ 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 Applying the expansion in γ c around zero up to the first order, and keeping only the leading logarithmic term we get In the absence of non-linear correction the saturation scale in log(Q 2 /Q 2 s (x)) should be replaced by a much smaller scale µ 0 ≪ Q s (x), giving xg(x, Q 2 )| linear ≃ (x in /x) λ log(Q 2 /µ 2 0 ).Hence the relative correction due to non-linearity reads 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 2 s (x), and this imposes an effective lower cut-off on the logaritmic integration in xg(x, Q 2 ) = Q 2 s (x) dk 2 /k 2 f (x, k 2 ).In these consideration we assumed that the saturation scale Q s ≫ µ 0 , where µ 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 µ 0 enters only as a lower cut-off of logarithmic integrations.It should be clearly distinguished from an initial scale of the DGLAP evolution µ F,0 , which is typically set to be greater than 1 GeV, and greater than the saturation scale in the proton.For the case µ F,0 ≫ Q 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.

The intermediate region
Above we discussed the non-linear effects for Q 2 ≫ Q 2 s (x) and in the high gluon density regime.The intermediate region of Q 2 ∼ Q 2 s (x) is hardest to analyze analytically, as in this region Q 2 s (x)/Q 2 ≃ 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 GeV 2 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 2 s 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 • 10 −6 at Q 2 = 1 GeV 2 , see e.g.[64].Hence one may reach Q 2 s of few GeV 2 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 2 s .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 ∼ Q 2 s /Q 2 , one may expect that the twist expansion ceases to make sense for Q 2 < Q 2 s .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.

Figure 2 :
Figure 2: Effects of non-linear corrections in proton structure functionsF 2 (x, Q 2 ) (left) and F L (x, Q 2 ) (right) at Q 2 = 5GeV 2 in twist 2 approximation.We show the twist 2 components obtained from the BFKL evolution, the BK correction and their sum.

Figure 3 :
Figure 3: Relative effects of twist 4 corrections in proton structure functionsF 2 (x, Q 2 ) (left) and F L (x, Q 2 ) (right) at Q 2 = 5 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.
), as well as, the ratio of BK correction twist 4 to twist 2 for both structure functions: R ∆F 2 = ∆F

Table 2 :
Ratio of twist 4 to twist 2 for structure functions F 2 , F L and non-linear corrections.Definitions can be found in the main text.