Non-linear evolution in the re-summed next-to-leading order of perturbative QCD:\confronting the experimental data

In this paper we compare the experimental HERA data with the next-to-leading order approach (NLO) of Ref.[C.~Contreras, E.~Levin, R.~Meneses and M.~Sanhueza,Eur. Phys. J. C 80 (2020) no.11, 1029). This approach includes the re-summed NLO corrections to the kernel of the evolution equation, the correct asymptotic behaviour in the NLO at $\tau = r^2 Q^2_s \,\gg\,1$; the impact parameter dependence of the saturation scale in accord with the Froissarrt theorem as well as the non-linear corrections. In this paper, we successfully describe the experimental data with the quality, which is not worse, than in the leading order fits with larger number of the phenomenological parameters. It is demonstrated, that the data could be described, taking into account both the diffusion on $\ln(k_T)$, which stems from perturbative QCD, and the Gribov's diffusion in impact parameters. It is shown an ability to describe the data at rather large values of $\alpha_S$.

The goal of this paper is to compare with the experimental (HERA) data the next-to-leading order approach (NLO) of Ref. [1]. In Ref. [1], we develop the approach in which we include the re-summation procedure, suggested in Refs. [2][3][4], to fix the BFKL kernel in the NLO. In particular, we introduce the rapidity variable, which plays the role of the "evolution time", in the same way as in Ref. [5]. However, we suggest a different way to account for the non-linear corrections, than in Ref. [5], which leads to additional change of the NLO kernel of the evolution equation. The advantage of our kernel of the BFKL equation [6,7], is that the scattering amplitude satisfies the high energy limits, which follows from the approach of Ref. [8] (see Refs. [9,10]) to the NLO Balitsky-Kovchegov (BK) [11] evolution [12][13][14][15][16][17][18].
We firmly believe that finding the correct NLO approximation for the non-linear evolution is one of the most important and urgent problem in the theoretical description of the high energy scattering. Indeed, in the Colour Glass Condensate(CGC) approach, which is the only candidate for an effective theory at high energies (see Ref. [19] for a review), the two essential parameters, that determine the high energy scattering, calculated in leading order of perturbative QCD [6,11,[20][21][22][23][24][25][26][27][28] turns out to be in an apparent contradiction with the experimental data. The first one is the BFKL Pomeron [6] intercept, which is equal to 2.8ᾱ S and leads to the energy behaviour of the scattering amplitude N ∝ exp 2.8ᾱ S ln( 1 x ) . The second is the energy behaviour of the new dimensional scale: saturation momentum Q 2 s ∝ exp 4.88ᾱ S ln( 1 x ) . Both show the increase in the leading order CGC approach, which cannot be reconciled with the available experimental data. So, the large NLO corrections appear as the only way out, now as well as two decades ago.
In the next section we will outline the main results of Ref. [1] and will specify our theoretical description of the dipole scattering amplitude. However, the current stage of our theoretical understanding of non-perturbative QCD is such, that we have to build a model. We need to take into account the non-perturbative corrections that will reproduce the correct, exponentially decreasing at large impact parameters (b) scattering amplitude. It has been demonstrated in Refs. [29][30][31][32], that the CGC equations [11,[23][24][25][26][27][28] as well as all other approaches, based on perturbative QCD, lead to the amplitude that increases as a power of energy, resulting in the violation of the Froissart theorem [33] 1 . Unfortunately, without a theoretical control on non-perturbative QCD we have to use a phenomenological approach to model the large b behaviour. In this paper we will exploit two approaches: 1. the non-perturbative behaviour of the saturation scale, which we will parameterize as follows: and the value ofγ we will discuss below. In the vicinity of the saturation scale such b dependance results in the large b-dependence of the scattering amplitude, which is proportional to exp (−m b) at b 1/m, in accordance of the Froissart theorem [33]. In addition, we reproduce the large Q T dependence of this amplitude proportional to Q −4 T which follows from the perturbative QCD calculation [39]. Theoretically the fact that we can absorb the non-perturbative b-dependence in Q s (Y, b) (see Eq. (1) for example), follows from the semi-classical approach to BK equation [40] and has been widely used in all, so called, saturation models [9,[41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60].
The paper is organized as follows. In section II we will give a brief review of the approach that we have developed in Ref. [1]. Next we will discuss the procedure of obtaining the NLO kernel of the BFKL equation based on the anomalous dimensions for re-summed NLO corrections, which is suggested in Refs. [2][3][4]72]. Then we will consider the modification that we need to introduce in the NLO kernel to account for correct behaviour of the scattering amplitude at high energy and the structure of the non-linear equation. In section III we specify our approach, which we use for describing the experimental data. In particular, we introduce the phenomenological parameters, which have to be calculated in the non-perturbative QCD approach, and discuss their physical meaning and the range of possible values. In section IV we collect the results of the fit of the experimental data on DIS. We summarize our results in the conclusion.

II. LEADING TWIST APPROXIMATION FOR NON-LINEAR EVOLUTION IN THE NLO -A RECAP
A.
Re-summed anomalous dimensions in NLO and the kernel of linear evolution The general solution to the linear equation has the following form: where N is the scattering amplitude of the dipole with size r at the impact parameter b. Y is the rapidity of the dipole. φ γ (r, b) is the eigenfunction of the BFKL equation which has the general form: with ξ = ln where R is the size of the target. φ in can be found from the initial condition at Y = 0.

B. Non-linear equation and the feedback to the kernel of the linear evolution
The general structure of the non-linear Balitsky-Kovchegov equation [11] has the following form: where x ik = x i − x k and x 10 ≡ r, x 20 ≡ r and x 12 ≡ r − r . Y is the rapidity of the scattering dipole and b is the impact factor. K (x 02 , x 12 ; x 10 ) is the kernel of the BFKL equation. In our approach we wish to preserve this form, but include the NLO corrections to the kernel. In particular, we would like to include the reggeization term in Eq. (16), which contribute to the linear equation, but has been neglected in DLA (see Eq. (15)). It should be stressed that only keeping this term we can provide the correct asymptotic behaviour at large Y : N → 1. Therefore, the linear equation takes the form: and Eq. (14) can be re-written in the form: For the dipole amplitude, N (ξ , Y ), Eq. (14) takes the following form: In red we putγ = 0.63 which stems from LO estimates and λ = 0.2 which is a typical value that is need to describe DIS data (see for example Ref. [59]).
which leads to the eigenvalue ω(γ): Using Eq. (11) in the vicinity of γ → 1 one can see that Eq. (20) takes the form [1]: Using the general equation to determine the critical anomalous dimension and the energy behaviour of the saturation scale (see review [19]): we obtain:γ Re-writing Eq. (23) for the saturation momentum as a function of Y we obtain C. NLO BFKL kernel in the saturation domain

The kernel in γ-representation
The non-linear corrections are essential in the region of small values of x (large Y ) where τ = r 2 Q 2 s > 1. In Ref. [8] it is shown that in the BFKL equation we have two types of logarithms, which are generated by the following kernels In Ref. [1] we suggest to use Eq. (25) instead of the full expression of Eq. (10).
In the previous section we specified how we changed the kernel in the perturbative QCD region, taking into account the NLO corrections. In the region of large τ we need to find the anomalous dimension at γ → 0, as it follows from Eq. (25). Hence we re-write Eq. (11) in the vicinity of γ → 0, where it has the form: which differs from the behaviour of the LO kernel, only by replacingᾱ S →ᾱ S / (1 +ᾱ S ). Therefore, we can discuss the non-linear equation in the LO, substitutingᾱ S / (1 +ᾱ S ) in place ofᾱ S , at the final stage.

The non-linear equation
In the saturation region where τ > 1, the logarithms originate from the decay of a large size dipole, into one small size dipole and one large size dipole [8]. However, the size of the small dipole is still larger than 1/Q s . This observation can be translated in the following form of the kernel in the LŌ . Inside the saturation region the BK equation of the LO takes the form For the NLO kernel of Eq. (27), Eq. (29) takes the form:

The solution
For solving this equation Substituting Eq. (31) into Eq. (30) we reduce it to the form where λ is given by Eq. (24). The variable ξ s is defined as The use of this variable indicates the main idea of our approach in the region of τ > 1: we wish to match the solution of the non-linear Eq. (32b) with the solution of linear equation (see Eq. (19)) in the kinematic region where it has the form [76]: andγ is determined by Eq. (24). Eq. (34) leads to the initial and boundary conditions: Eq. (32b) has a traveling wave solution (see formula 3.4.1.1 of Ref. [77]). For Eq. (32b) in the canonical form: with t ± = ξ s ± ξ, the solution takes the form: where all constants have to be determined from the initial and boundary conditions of Eq. (35). First we see that C 2 = 0 and κ = 0. From the condition Ω z /Ω =γ at t + = 0 we can find C 1 . Indeed, differentiating Eq. (37) with respect to t + one can see that at t + = 0 we have: From Eq. (38) one can see that choosing we satisfy the initial condition d ln(Ω) dz | t+=0 =γ of Eq. (24). Finally, the solution of Eq. (37) can be re-written in the following form for Ω 0 1: One can see that Eq. (40) gives the geometric scaling solution [8,[78][79][80] which depends only on the variable z. For Ω → Ω 0 and if Ω 0 1, Eq. (40) can be solved explicitly giving Eq. (41) gives the solution which depends only on one variable z = ξ s + ξ, and satisfies the initial conditions of Eq. (35). At large z we obtain the solution [8]: 42) or in terms of the amplitude We wish to stress that Eq. (43) reproduces the asymptotic solution to the BK equation in the NLO, which has been derived in Refs. [9,10], for fixedᾱ S . It should be noted that both solutions of Eq. (41) and Eq. (42) can be derived directly from Eq. (36) assuming 1 − exp (−Ω) → Ω and 1 − exp (−Ω) → 1 for small z and large z, respectively.

A. Generalties
The observables in deep inelastic scattering can be expressed through the following scattering amplitudes (see Fig. 2 and Ref. [19] for the review and references therein) where Y = ln (1/x Bj ) and x Bj is the Bjorken x. z is the fraction of energy carried by quark. Q is the photon virtuality. b denotes the impact parameter of the scattering amplitude.
Eq. (44) shows the main features of the interactions at high energies. They go in two stages. The first is the decay of virtual photon in quark-antiquark pair, described by |Ψ γ * (Q, r, z) | 2 in Eq. (44). For large Q 2 ( Q 2 ≥ Q 2 0 with Q 2 0 ≈ 0.7GeV 2 , see Ref. [81]) the wave function is well known (see Ref. [19] and references therein) where T(L) denotes the polarization of the photon and f is the flavours of the quarks. 2 = m 2 f + Q 2 z(1 − z). However, even for DIS the non-perturbative corrections become essential and we cannot use Eq. (45) and Eq. (46) for Q 2 ≤ Q 2 0 . The second is the amplitude N (r, Y ; b) of the interaction of the dipole with the target. We have discussed this amplitude in the previous section and will specify the way how we will use the finding of this section for the practical application to DIS below.

B. The dipole scattering amplitudes
Based on the approach, briefly discussed in the previous section, we can see three different kinematic regions where we are going to use the scattering amplitude N (Y, ξ ; b) in three different forms.
1. For Q 2 s r 2 < 1 (perturbative QCD region) we suggest to use the linear evolution equation of Eq. (19), which has been discusssed above. Recalling that in the derivation of this equation we use the DLA in which both η = Y − ξ and ξ are considered to be large:ᾱ S η ξ 1. For η > 0 (Y > ξ) we can use the experimental data as the initial condition for Eq. (19). The value of the saturation momentum is given by Eq. (24). Since its value turns out to be much less than Q max , which stems from the condition ξ max = ln Q 2 max /Q 2 0 = η, we can use the DGLAP evolution equation in the next-to-leading order for Q 2 > Q 2 max . 2. For Q 2 s r 2 ∼ 1 (vicinity of the saturation scale) we use the scattering amplitude in the form [76] N (r, withγ from Eq. (24). For τ < 1 we have to solve the linear equation with some initial condition which we need to take from the experimental data. This solution will give us the scattering amplitude in the vicinity of the saturation scale where the scattering amplitude can be estimated using Eq. (48). We use a different way to introduce the parameters from the experimental data: we expand Eq. (48) to the region τ < 1 , replacingγ by following expression:
The mass of the c-quark (about m c = 1.4 GeV) is not small and we took this into account replacing x in the cc scattering by For τ = Q 2 s r 2 1 we need to use the solution of Eq. (40) to the non-linear evolution equation. However, it has been found in Ref. [82] that the following formula: with Ω (z) from Eq. (41) and with a = 0.65 describes the exact solution within accuracy less that 2.5 %. Therefore, we use Eq. (51) in our attempts to describe the HERA data.

C. Phenomenological input
As has been mentioned above, we need to set the initial conditions at Y = 0 for the linear evolution equation for τ ≤ 1. In this section we wish to clarify how we introduce the phenomenological parameters to describe this conditions. The first one we have considered: N 0 which determines the scattering amplitude at τ = 1. Two other parameters describe the saturation momentum at Y = 0: the value of Q s at b = 0 and the behaviour as a function of b. collecting everything that we have discussed about Q s we use two expression for the saturation momentum: where Z and λ are defined in Eq. (2) and Eq. (24), respectively. The values of these parameters have to be found from fitting of the experimental data. The experience with such fitting gives Q Finally, we introduce three phenomenological parameters from the initial conditions whose values have to be found from the fit of the experimental data. The masses of quarks, that determine the wave function of the virtual photon, determine the infra-red behaviour of the wave function. We take two sets of them: the current masses and the masses of light quarks are equal to 140 MeV which is the typical infra-red cutoff in our approach.

IV. RESULTS OF THE FITS
Using the approach, that has been discussed in the previous section, we attempt to describe the most accurate data for the deep inelastic structure function F 2 [83]. The implicit parameters of our approach are the restrictions of the kinematic region of the experimental data that we include in the fit. We chose: 0.85 GeV 2 ≤ Q 2 ≤ 27 GeV 2 and x ≤ 0.01. The lower limit of Q 2 stems from non-perturbative correction to the wave function of the virtual photon, while the upper limit originates from two restrictions: x ≤ 0.01 and contribution of the additional term in γ of Eq. (49) is small. The choice of the largest x is dictated by the needs to have low x for legitimate use of our theoretical formulae and the practical wish to use in the fit as more data as possible.

Dipole amplitude
Wave function    1 and 2), which is the typical infra-red cutoff in our approach. For the sets 2 and 4 the value ofᾱS = 0.2 is fixed and only three parameters were chosen from the fit.
In Table 1 we show the values of parameters from our fits. It should be noted that we did not fit the mass of the quarks but use two sets of them: the masses of light quarks are equal to m = 140 M eV , which we consider as the typical infra-red cutoff in our approach (sets 1 and 2); and the current masses (sets 3 and 4). We used the data for the inclusive F 2 in the kinematic region: 0.85 GeV 2 ≤ Q 2 ≤ 27 GeV 2 and x ≤ 0.01, for finding the parameters while F cc X Bj with the experimental data, to demonstrate our ability to describe the experimental data. The quality of the fit one can see from Fig. 3. In Fig. 4 we show the data on F 2 as function of Q 2 at different values of x. This figure illustrates that we can describe the data for x ≤ 0.013. In addition, we show in Fig. 3 F 2 for 27 ≤ Q 2 ≤ 60 GeV 2 . The quality of the fit is not very good due to the small experimental errors. Indeed, χ 2 /d.o.f. for the range of Q 2 from 0.85 GeV 2 to 60 GeV 2 is equal to 2.16 (set 1), 1.93 (set 3) and 3.19 (set 5). The reason for this stems from the approximate character of Eq. (49). The contributions of the nonlinear corrections are negligible for such large photon virtualities and, therefore, our approach coincides with the one of Ref. [5] (see also Ref. [85]), which describes the data very well.
We believe, that we demonstrated that our approach can describe the experimental data at rather small values of Q 2 where the non-linear corrections give an essential contribution. For large Q 2 our approach coincides with the one of Ref. [5] which takes into account the linear evolution in much better way than we [85]. Fig. 5 demonstrates how our fits describe the experimental data on F cc 2 while Fig. 6 shows that our fit is able to reproduce the data at low values of Q 2 . One can see that the agreement with the experimental data is good even at low Q 2 . It should be stressed that two sets of the light quark mass give the same description illustrating a possibility to use the wave function of the virtual photon in perturbative QCD at rather low values of Q 2 .
In addition to F 2 and F CC 2 we compare our approach with the experimental data on F L [89,90]. In spite of the fact that the data have sufficiently large errors one can see from Fig. 7 that we are able to describe the data quite well.
Sets 2 and 4 need special comments. In these sets we chooseᾱ S = 0.2 to illustrate the ability of our approach to overcome the difficulties of next-to-leading approaches in describing the experimental data: the small value of neededᾱ S . One can see, that we obtain a good χ 2 /d.o.f. but for the specific set of data in the wide kinematic region 0.85 ≤ Q 2 ≤ 60 GeV 2 but choosing a restricted selection of data. Note that these data are marked by the squares in the panel of Fig. 3 with sets 2 and 4.
The general characteristics of our fits are: (1) the possibility to describe the data at rather small values of Q 2 , where the non-linear corrections turn out to be essential, (2)  Unfortunately, in spite of large numbers of the papers in which the experimental data have been compared with the saturation models [9,[41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60] we can confront our approach only with Refs. [9,60] since in all other papers the assumptions were made that contradict the theoretical informations. For example, in the most models Q 2 s is proportional to exp −b 2 /B while this behaviour disagrees with the Froissart theorem [33]; or/and the behaviour for r 2 Q 2 s 1 contradicts the approach of Ref. [8]. Comparing our parameters with Ref. [9] one can see that the typical value ofᾱ S ≈ 0.14 − 0.15 in Ref. [9] is less than in sets 2 and 4. The value of Q 2 0 is about in three times larger than in this paper. The value of m in Ref. [9] is the same as in the electro-magnetic form factor of proton, while in this paper m is larger, than in the electro-magnetic form factor of proton in agreement with the new experimental information on the behaviour of the two gluon form factors [86,87]. Taking into account that χ 2 /d.o.f. are better in this paper than in Ref. [9,60], we believe that the theoretical approach of section II, give reliable description of the current experimental data on DIS.
We wish to draw your attention to the sets 5 and 6, which take into account the different dependence of the saturation momentum on b (see Eq. (52b)). From Fig. 3 and Fig. 5 one can see that we are able to describe the experimental data. Hence we see that the suggested dependence on rapidity of the impact parameter distribution is in accord with the experimental data.

V. CONCLUSIONS
In this paper we compare with the experimental data the NLO approach of Ref. [1]. In this approach we include the re-summation procedure, suggested in Refs. [2][3][4][5], to fix the BFKL kernel in the NL0, but we treat differently the non-linear corrections. The advantage of our treatment is that we reproduce the correct asymptotic behaviour at large τ = r 2 Q 2 s . Fixing our phenomenological parameters from the fit to the experimental data on F 2 DIS structure function in the region: 0.85 ≤ Q 2 ≤ 27 GeV 2 and x ≤ 0.013, we found that can describe these data with very good χ 2 /d.o.f.. It should be stressed, that for the first time it is demonstrated that the impact parameter dependence of the saturation scale, which is given by Eq. (52b), is an agreement with the data. We believe that this is an interesting result since Eq. (52b) follows from the first attempt to take into account both the diffusion on ln(k T ), which stems from perturbative QCD, and the Gribov's diffusion in b [71].
Using the phenomenological parameters, that we found from the fit of F 2 , we demonstrated that our approach is able to describe the experimental data both on F cc 2 and on F 2 at low values of Q 2 = 0.4−0.7 GeV 2 . We also reproduce the experimental data on F L .
It should be pointed out that our description in the NLO is better or of the same quality as the description by LO with more fitting parameters. Two problems, that we faced in describing the data in the NLO approach, have been partly healed: we are able to fit the data with not small value ofᾱ S = 0.2 and the value of the saturation momentum at Y = 0 is about Q 2 0 ≈ 1 GeV 2 instead of 3 GeV 2 of Ref. [9]. However, we are aware that Q 0 = 1 GeV is rather large (see for example Ref. [88]).  Concluding we wish to mentioned that the model that we developed here, is the only one which includes (i) the resummed NLO corrections; (2) the correct theoretical behaviour at large τ ; (3) the large b dependence of the saturation scale in accord with the Froissart theorem as well as the non-linear evolution. We believe that it will be useful for further discussion of the high energy scattering in QCD.

VI. ACKNOWLEDGEMENTS
We thank our colleagues at Tel Aviv university and UTFSM for encouraging discussions. Our special thanks go to E. Gotsman for all his remarks and suggestions on this paper, and to Yuri Ivanov for technical support of the USM