Non-linear equation in the re-summed next-to-leading order of perturbative QCD: the leading twist approximation

In this paper, we use the re-summation procedure, suggested in Refs.\cite{DIMST,SALAM,SALAM1,SALAM2}, to fix the BFKL kernel in the NLO. However, we suggest a different way to introduce th non-linear corrections in the saturation region, which is based on the leading twist non-linear equation. In the kinematic region:$\tau\,\equiv\,r^2 Q^2_s(Y)\,\leq\,1$ , where $r$ denotes the size of the dipole, $Y$ its rapidity and $Q_s$ the saturation scale, we found that the re-summation contributes mostly to the leading twist of the BFKL equation. Assuming that the scattering amplitude is small, we suggest using the linear evolution equation in this region. For $\tau \,>\,1$ we are dealing with the re-summation of $\Lb \bas \,\ln \tau\Rb^n$ and other corrections in NLO approximation for the leading twist.We find the BFKL kernel in this kinematic region and write the non-linear equation, which we solve analytically. We believe the new equation could be a basis for a consistent phenomenology based on the CGC approach.


INTRODUCTION
The Colour Glass Condensate(CGC) approach is the only candidate for an effective theory at high energies, which is based on our microscopic theory: QCD (see Ref. [1] for a review). However, it has been known for a long time, that to describe the scattering amplitude in the framework of CGC [2][3][4][5][6][7][8] we need to include at least the next-to-leading order corrections to the non-linear equations. Indeed, the two essential parameters, that determine the high energy scattering, are the BFKL Pomeron [9] intercept, which is equal to 2.8ᾱ S , which leads to the energy behaviour of the scattering amplitude N ∝ exp 2.8ᾱ S ln( 1 x ) , and 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.
The non-linear equations in the NLO has been written in Refs. [11][12][13][14][15][16][17][18], but their use in high energy phenomenology is marred by the instabilities due to the presence of large and negative NLO corrections enhanced by double collinear logarithms (see Ref. [18] for discussions and references). These problems have been found in Refs. [19,20], identified and solved in Refs. [21][22][23] for the linear BFKL equation. It turns out that instabilities are closely related to the wrong choice of the energy(rapidity) scale, and we need to introduce the re-summed NLO corrections to cure this problem.
In this paper, we use the re-summation procedure, suggested in Refs. [21][22][23], to fix the BFKL kernel in the NLO, which coincides with the kernel used in Ref. [18]. However, we suggest a different way to introduce the non-linear corrections, than in Ref. [18], which is based on the approach, developed in Ref. [24] for the leading twist non-linear equation in the LO. Our first observation is that the re-summation contributes mostly to the leading twist of the BFKL equation in the kinematic region τ ≡ r 2 Q 2 s (Y ) < 1 , where r denotes the size of the dipole, Y its rapidity and Q s the saturation scale. We assume that at τ = 1 the scattering amplitude is small and we can neglect the non-linear corrections. Therefore, for τ < 1 we can restrict ourselves by the linear BFKL equation. For τ > 1 we are dealing with the re-summation of (ᾱ S ln τ ) n and other corrections in NLO approximation for the leading twist. In this paper we find the BFKL kernel in this kinematic region, and write the non-linear equation. We found the analytical solution of this equation, and believe the new equation could be a basis for a consistent phenomenology based on the CGC approach.

II. BK NON -LINEAR EQUATION
The BK evolution equation for the dipole-target scattering amplitude N (x 10 , b, Y ; R) has in the leading order (LO) of perturbative QCD [1,7,[25][26][27] 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 which in the leading order has the following form: K LO (x 02 , x 12 ; x 10 ) = x 2 10 x 2 02 x 2 12 (2) In Eq. (1) R denotes the size of the target dipole andᾱ S = N c α S /π where N c is the number of colours. For the kernel of the LO BFKL equation (see Eq. (2)) the eigenvalues take the form [9,10]: where ψ(z) denotes the Euler psi-function ψ (z) = d ln Γ(z)/dz. In the next-to-leading order (NLO), the non-linear equation has a more complicated form [11][12][13][14]: and N f and N c are the number of fermions and colours, respectively. S ij denotes the S-matrix for scattering of a dipole of size x ij with the target, which can be written using the scattering amplitude N ij , as follows S ij = 1 − N ij . Eq. (4) gives the explicit form of the BFKL kernel in the NLO, but as it has been alluded to we need to re-sum the NLO corrections to avoid instabilities. We re-sum in the approximation, that was suggested in Ref. [18], which we will discuss below. It turns out, that in the framework of this re-summation we can neglect the contribution, which is proportional to S 13 S 32 − S 13 S 32 in Eq. (4); and reduce Eq. (4) to Eq. (1) with the kernel, which has to be found in the re-summed NLO. An additional argument for such simplification stems from the fact that deep in the saturation region, where S ij → 0, all terms except the first one, which is proportional to S 12 , are small and can be neglected. In other words, deep in the saturation region Eq. (4) reduces to Eq. (1) without addressing the specific form of re-summation.
In the next-to-leading order the kernel is derived in Refs. [19,20] and has the following form: The explicit form of χ N LO (γ) is given in Ref. [19]. However, χ N LO (γ) turns out to be singular at γ → 1, Such singularities indicate, that we have to calculate higher order corrections to obtain a reliable result. The procedure to re-sum high order corrections is suggested in Ref. [21][22][23]28]. The resulting spectrum of the BFKL equation in the NLO, can be found from the solution of the following equation [21][22][23] where and Functions χ N LO (γ) and A T (ω) as well as the constants (F and b), are defined in Refs. [21][22][23]. In Ref. [28] Khoze, Martin, Ryskin and Stirling (KMRS) sugested an economic form of χ 1 (ω, γ), which coincides with Eq. (8) to within 7%, and, therefore, gives reasonable estimates of all constants and functions in Eq. (8). The equation for ω takes the form One can see that γ(ω) → 0 when ω → 1 as follows from energy conservation.
The general solution to Eq. (10) can be written as follows where φ γ (r, R, b) is the eigenfunction of the BFKL equation and φ in (γ, R) can be found from the initial condition at Y = 0. In Eq. (11) r ≡ x 10 denotes the size of the scattering dipole, while R the size of the target. In Ref. [10] it was proved that the eigenfunction of the BFKL equation has the following form with ξ = ln for any kernel, which satisfies the conformal symmetry. Using Eq. (12) we can re-write the general solution in the form: where ξ = ln r 2 Q 2 s (Y = 0; b, R) . Comparing Eq. (13) with this definition of ξ, one can see that for r R.

B. Eigenvalues of the BFKL equation
As has been mentioned, the re-summation, that has been suggested in Ref. [18] is determined by the anomalous dimension γ in the vicinity of the eigenvalues at γ → 1. The singular part of the general kernel in the NLO (see Eq. (6)) has the following form: with the solution: It is instructive to note that Eq. (9) gives All other terms in Eq. (9) vanish at γ = 1. Solving Eq. (18) we obtain Eq. (19) gives a good description of the eigenvalues of Eq. (9) ω KMRS for γ > 1 − γ cr =γ (see Fig. 1). This is the reason that we denote this eigenvalue as ω > .

C. Saturation and geometric scaling behaviour of the scattering amplitude in the NLO
It is well known [1,24,25,[30][31][32], that one does not need to know the precise structure of the non-linear corrections for finding the saturation momentum, as well as for discussing the behaviour of the scattering amplitude in the vicinity of the saturation scale. We only need to find the solution of the linear BFKL equation, which is a wave package that satisfies the condition, that phase and group velocities are equal [25]. This solution determines the new dimensional scale of the problem: the saturation momentum Q s . The equation takes the form For the eigenvalues of the BFKL equation in the NLO given by Eq. (17), the solution to Eq. (20) takes the form: The equation determines the saturation momentum ξ s = − ln r 2 1 Actually Q 2 s (Y = Y 0 ) depends on other variables as well (see Eq. (15)), but we will omit these variable in the further presentation. For Eq. (19) it turns out that From Fig. 2 one can see that Eq. (23) leads to the value of λ ≈ 0.25, which describes the experimental data at α S ≈ 0.2.
In the vicinity of the saturation scale where r 2 Q 2 s (Y ) = τ → 1, the scattering amplitude has the following form [31,32]: This equation gives the initial conditions for the scattering amplitude in the saturation domain. These conditions have the forms: D. Double log approximation (DLA) in the re-summed NLO approximation

DLA in the LO approximation
In the leading order BFKL equation, the DLA stems from the eigenfunction In the coordinate representation (see the linear part of Eq. (1)) the DLA comes from the distances x 12 ∼ x 02 x 01 . For such distances the BFKL equation has the form In the derivation of Eq. (27) we used Eq. (2) for the kernel and assumed that b x 12 . Note, that the gluon reggeization term (N (x 10 , b, Y )) in Eq. (27) does not contribute in the DLA.
IntroducingÑ (x ij ) = N (x ij )/x 2 ij and changing ξ → ξ = −ξ we obtain Eq. (27) in the traditional form: Eq. (28) follows naturally from the general solution to the BFKL equation 2 : using the eigenvalue of Eq. (26). Indeed, plugging Eq. (26) into Eq. (29) and taking the integral over γ, closing the contour of integration we obtaiñ Taking the integral over ω we obtain the solution which is a series, with a general term which is proportional (ᾱ S Y ξ ) n . The particular form of the solution is determined by the functionñ in (ω) which we can find from the initial condition at Y = 0.

DLA in the NLO for BFKL equation
We wish to find the DLA in the NLO , using Eq. (29) in which we substitute the solution 3 to Eq. (16) with respect to γ: Plugging Eq. (31) into Eq. (29) we obtaiñ where we have introduced new variables: η = Y − ξ and ξ = −ξ. From Eq. (32) one can see that the solution is a series with the general term ∝ (ᾱ S η ξ ) n , whose exact form is determined by the initial condition at Y = 0. The solution of Eq. (32) can be re-written in more economical form: Therefore, we see that the only difference of the DLA in the NLO, stems from a new variable η, which replaces Y = ln(1/x). The physical meaning of this replacement has been examined in detail in Ref. [18]. The point of it is the fact, that the partons (gluons) in the wave function of the fast hadrons are ordered in accord of the lifetimes of the partons, which are proportional t i ∼ x i P/k 2 i,T , where x i denotes the fraction of the longitudinal momentum of the fast hadron, with the moment P , and k i,T is the transverse momentum of the i-th parton. The ordering that leads to the logs in η, has the form It is instructive to note, that the form of the kernel and the differential Eq. (33) are quite different, from the ones that has been discussed in Refs. [18,34]. Indeed, in Ref. [18] it is suggested that the linear equation has the form: where ρ ≡ L x02,x01 L x12,x01 and L xi2,x01 ≡ ln(x 2 i2 /x 2 01 ). The solution to Eq. (34) has the following form for η = Y − ξ > 0 [18,34] Assuming that solution to Eq. (33) depends on one variable ζ = 2 √ᾱ S ξ η we obtain that this equation takes the form: which has the solution N (ζ) = C 1 I 0 (ζ) (see formula 8.494(1) in Ref. [35]). Summing double logs and comparing these two solutions, one can see that both are function ofᾱ S η ξ . Both solutions atᾱ S η ξ 1 have the asymptotic behaviour The function H depends on Y and ξ logarithmically. The expression for the saturation momentum stems from the equation which gives the condition when N ∼ H (Y, ξ) ∼ Const, since H is slowly changing function.
It is easy to see that Eq. (38) leads to the saturation momentum of Eq. (22). However, to understand better the difference between these two equations, we would like to re-write Eq. (34) in the form of Eq. (33). The general solution to Eq. (34) has the following form: where we use Eq. (17) for ω(γ), and replace 1 − γ by γ . The function φ in (γ , R) is determined by the initial condition at Y = 0 and R denotes the size of the target dipole. Using Eq. (39) we see that N (ξ , Y ) satisfies the following equation: , can be reduced to Eq. (33). Indeed, one can see that the solution of Eq. (32) satisfies both these equations.

Non-linear evolution in LO DLA
Gluon reggeization: The linear evolution equation in the leading order DLA is given by Eq. (28). To include non-linear corrections, we need to consider the general BK equation (see Eq. (1)), with the kernel in the DLA approximation. The first question which arises, is how to treat the reggeization term in Eq. (27), which is neglected in the DLA. Indeed, this term does not contribute in the DLA and including it, is a particular way to make estimates beyond those of the DLA. We believe that it is necessary to do this, to provide the correct behaviour of the scattering amplitude which should approach 1 (N → 1) for large Y . Bearing this in mind we need to preserve the gluon reggeization term in Eq. (27), which has the form we can re-write Eq. (41) in the form: This equation corresponds to the eigenvalue ω(γ) for the amplitudeÑ (ξ , Y ): For the dipole amplitude, N (ξ , Y ), Eq. (41) takes the following form: which leads to the eigenvalue ω(γ) : Applying the procedure, that has been discussed in subsection III-C we obtain We would like to mention that the value of λ turns out to be 30% less that the LO DLA value, λ = 4ᾱ S . Non-linear equation: Finally, we need to re-write the general BK equation (see Eq. (1)) and account for the non-linear term. Using the BFKL kernel in the DLA we obtain For the dipole amplitude N = e −ξ Ñ , Eq. (47) has the following form 4 In definition of ξ we assumed that b R. For b R, ξ = ln b 2 /x 2 01 .
Eq. (48) has the traveling wave solution, which is a function of one variable z = α Y + β ξ . Indeed Eq. (48) takes the form From Eq. (43) the natural choice for the parameters α and β is α = 1 2ᾱ S 3 + 2 √ 2 and β = −1. We postpone the discussion of the nonlinear equation to the next section. For small N 1, we can neglect the non-linear term, and Eq. (49) leads to the scattering amplitude N which is equal to this coincides with the general form of Eq. (24) where λ andγ are determined by Eq. (46).

Equation:
As we have discussed in section III-D-2, the difference between LO DLA and NLO DLA (see Eq. (31)) stems only from the new energy variable: η = Y − ξ . Therefore, the non-linear equations have the form: The last term in Eq. (51) comes from the gluon reggeization (the third term in Eq. (1)). It does not contribute to the DLA approximation, but has to be taken into account together with the non-linear corrections , since it provides the correct asymptotic behaviour of the scattering amplitude N → 1 at Y 1. The linear equation for the amplitudeÑ (ξ , Y ) , whose origin is Eq. (51) neglecting the non-linear term, can be written, using Eq. (42), as the following equation for the eigenvalue ω(γ ): Eq. (52) differs from Eq. (31), by the last term in the r.h.s. of this equation, which takes into account the gluon reggeization contribution.
The solution which has the following form: where ω(γ ) is taken from Eq. (52). In Fig. 3 we plot ω m , which is the solution to the equation: This equation is the generalization of Eq. (52) in where we took into account, energy conservation. Note, that this equation introduces conservation of energy into Eq. (52). The solution to Eq. (54) has the form: . ω< is the same as in Fig. 1, which we discuss below in section IV-A . ω KMRS is given by Eq. (9).
For the dipole amplitude N = e −ξ Ñ Eq. (51) has the following form For Eq. (56), as well as for Eq. (48), we can find a solution, which is the function of one variable z = α η + β ξ . Indeed, the equation Eq. (56) takes the form As we have seen in the discussion of Eq. (48), the natural choice of parameters α and β is α = 1 2ᾱ S 3 + 2 √ 2 ≡ aᾱ S and β = −1.
Solution: generalities: We have not found the explicit solution to Eq. (57). This is the reason for discussing the general features of the solution in this subsection, based on the phase portrait (see Ref. [36]). The equation can be re-written in the matrix form as Eq. (58) has two critical points: (0, 0) and (1, 0). Near these critical points, Eq. (58) can be re-written in the matrix form: where ∆N denotes a small deviation of N in the vicinity of the critical point N crit . Matrices DF have the following forms: The eigenvalues of these matrices are equal to The path for this solution in the plane (N (z), N z (z)) is shown in Fig. 4 by a red line. Using Eq. (61), we can find solutions in the vicinity of the critical points. Indeed, at z → −∞ the eigenfunctions of Eq. (59) in the vicinity of the critical point (0,0) have the form: Therefore, the general solution at z → −∞ has the form Note that λ 1 = λ 2 , similarly, at z → ∞ we have in the vicinity of the critical point (1,0): and Since in vicinity of the critical point (1,0) λ 1 is positive, we need to put C 1 = 0 and therefore, we have Analytical solutions in the bounded regions: The previous analysis can be improved in two different kinematic regions: at small N < 1, when non-linear corrections are small; and at N → 1, where we can develop the semiclassical approach.
Eq. (57) can be solved, by looking for the solution d d z N (z) = F (N ). The equation takes the form: For N 1 the solution of Eq. (68) takes the form: As we have seen (see Eq. (48) and Eq. (56)), the natural choice of parameters α and β is, α = 1 2ᾱ S 3 + 2 √ 2 ≡ aᾱ S and, β = −1, which leads to C 1 = √ 2 − 1 . The dependence of N on z can be recovered from the equation: which has the solution: Note that Eq. (71) gives N (z), as a function of z = α η + ξ (see Eq. (57)) in the vicinity of the saturation scale. As we have discussed, the energy variable η leads to the Eq. (51), and Eq. (20) determines Q 2 s (η) = Q 2 s ( 0) exp λ η η. Therefore, to obtain N = N 0 r 2 Q 2 s (Y ) γ we need to re-calculate the variable z in terms of Y and ξ: z = λ η (Y + ξ) + ξ. Hence, in Eq. (71) Eq. (74) leads too the final expression in Eq. (71) . It should be noted, that Eq. (71) is the same as the solution of Eq. (61) and Eq. (63), which we obtained in a different way.
In Fig. 5 we plot the dependence λ on the values ofᾱ S . One can see that Eq. (72) leads to λ ≈ 0.3 at rather large values ofᾱ S ≈ 0.2, which is needed for the description of the HERA experimental data In the kinematic region, where N → 1 it is more convenient to solve Eq. (57) directly, by looking for the solution in the form N = 1 − exp (−Ω(z)), assuming that Ω (z) is a smooth function of z (Ω zz (z) Ω 2 z (z)). For such a function Ω (z), Eq. (57) can be re-written in the form: Solving Eq. (75) we obtain Therefore, Ω (z) can be found from the equation We choose the sign + in Eq. (77) since we are looking for Ω(z) which is positive at large z, leading to N ≤ 1.
Using Eq. (77) we can evaluate Ω z (z) and Ω zz (z): The integral over Ω can be taken then, and Eq. (77)has the form: where Ω 0 is the initial value of Ω at z = 0. At large z, Ω (z) → √ 2 − 1 z, while at small z, Ω (z) = exp 1 2 √ 2 − 1 z . Note, that we obtain similar behaviour of the amplitude at z → 0 as in Eq. (71), but with a different power of z.
In Fig. 6 we plot the solution of Eq. (54). Note, that Fig. 6-c shows that Ω zz Ω 2 z only for z > 12 ÷ 15. Therefore, for z < 12 ÷ 15 we have to solve the following equation: Practical way to obtain estimates for the scattering amplitude: The difficulties in solving the general equations (see Eq. (57) and Eq. (80)), are that the critical point (1,0) is not a stable point, but only a saddle point. This means that there exists only one path which gives the solution with N → 1 at large z (see Fig. 4). It is difficult to find this trajectory analytically or/and numerically. These difficulties are illustrate by Fig. 7 where N z /N and N for the semiclassical solution of Eq. (75) are shown in green, while the values of N z /N for small N , are plotted in red. One can see that it is not possible to find a function which provides a smooth matching to these two solutions. The only possibility is to find the exact solution, that covers a large region of z .
The one way out, comes from the hope that the solution at small N , is sufficient to cover the region of γ fromγ to γ = 1 2 , where we have to use Eq. (57). We believe that for smaller γ we need to use the different expression for the eigenvalues ω (γ), which is given by Eq. (91).

Energy conservation in DLA
We need to consider Eq. (18), which follows from Eq. (9) in the region γ → 1, for taking energy conservation into account. Resolving this equation with respect to γ we obtain for the amplitudeÑ (ξ , Y ): which leads to the following equation instead of Eq. (33). Repeating the same procedure as in section III-C-3, the following non-linear equation can be written: Neglecting the N 2 term in Eq. (83), we have a linear equation which determines the value of the saturation scale Q 2 s (η) = Q 2 s ( 0) exp λ η η andγ η : Tedious but simple calculations following the procedure of section III-C lead to The values of λ(ᾱ S ) andγ(ᾱ S ) are plotted in Fig. 5-b. Introducing a new variable z = α η + β ξ , this equation has the form: A natural choice of α and β is α = λ η (ᾱ S ) and β = −1. Considering N z (z) = F (N ) we can re-write Eq. (86) in the form For N 1, the solution of Eq. (87) can be obtained in the same way as the solution of Eq. (69): The solution to Eq. (88) can be written as follows: with λ andγ given in Eq. (85).
In the previous section we specified how we changed the kernel in the perturbative QCD region, taking into account the NLO corrections. We now desire to find how we need to change the kernel 1/γ for τ = rQ s > 1 (in the saturation region). For this purpose we re-write Eq. (9) in the vicinity of γ → 0, where it has the form: This eigenvalue describes the behaviour of ω KMRS for γ ≤ 1 2 (see Fig. 3). We start from the simplified expression for Eq. (91): Resolving Eq. (92) with respect of γ we obtain γ =ᾱ 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.

B. 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 [24]. 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. (93), Eq. (95) takes the form:

C. The solution
For solving this equation we introduce function Ω (Y ; ξ, b) [24] N (Y, ξ) = 1 − exp (−Ω (Y, ξ)) Substituting Eq. (97) into Eq. (96) we reduce it to the form where λ is given by Eq. (72). The variable ξ s is defined as The use of this variable indicates the main idea of our approach: we wish to match the solution of the non-linear Eq. (98b) with the solution of the non-linear Eq. (87). However, we assume that forγ ≥ 1 2 we need the solution for Eq. (87) for N < 1 , where it has the form N = N 0 exp (γ z) with z z = ξ s + ξ with t ± = ξ s ± ξ, the solution takes the form: where all constants have to be determined from the initial and boundary conditions of Eq. (25). First we see that C 2 = 0 and κ = 0. From the condition Ω z /Ω =γ at t + = 0 we can find C 1 . Indeed, differentiating Eq. (102) with respect to t + one can see that at t + = 0 we have: From Eq. (103) one can see that choosing we satisfy the initial condition d ln(Ω) dz | t+=0 =γ of Eq. (72). Finally, the solution of Eq. (102) can be re-written in the following form for Ω 0 1: For Ω → Ω 0 and if Ω 0 1, Eq. (105) can be solved explicitly giving Eq. (106) gives the solution which depends only on one variable z = ξ s + ξ, and satisfies the initial conditions of Eq. (25). At large z we obtain the solution [24]: or in terms of the amplitude We wish to stress that Eq. (108) reproduces the asymptotic solution to Eq. (4), which has been derived in Refs. [38,39], for fixedᾱ S . It should be noted that both solutions of Eq. (106) and Eq. (107) can be derived directly from Eq. (101) assuming 1 − exp (−Ω) → Ω and 1 − exp (−Ω) → 1 for small z and large z, respectively.

D. Non-linear equation: NLO + energy conservation
In this section we discuss the general form of the eigenvalues which is given by Eq. (91). We obtain the following solution for ω from this equation.
where α S =ᾱ S /(1 +ᾱ S ) and ρ = ξ − ξ = ln x 2 01 Q 2 s (Y ) − ln x 2 02 Q 2 s (Y ) (see Eq. (34)). Therefore, the general solution for linear equation takes the form: One can see that for e α S ξ N (Y, ξ) =N (Y, ξ) we have: We assume that the main contributions stem from the region x 02 x 01 and/or x 12 x 01 Bearing this in mind we can use Eq. (94) for the kernel , which can be re-written in the form: where ξ ik = ln x 2 ik Q 2 s (Y ) . Using Eq. (112) we can re-write the general BK equation (see Eq. (1)) in the form: Using Eq. (111) we can re-write Eq. (113) in the form: Introducing N (Y, ξ) = 1 − exp (−Ω (Y, ξ)) we reduce Eq. (114) to the following expression: Differentiating Eq. (115) with respect to ξ we obtain: Plugging the last term in Eq. (116) from Eq. (115) we have: Looking for the solution which has the geometric scaling behaviour, we re-write Eq. (117) in new variable z = ξ s + ξ and it takes the form One can see that deep in the saturation region, where Ω 1, Ω (z) = 1 λ z. However, for z < 1/λ we reproduce the solution of Eq. (107). Therefore, we can infer that the energy conservation crucially changed the behaviour of the scattering amplitude deep in the saturation region. In Fig. 8 we plot the solution to Eq. (98b) and to Eq. (118). One can see that the difference is large for z > 6. At z ≤ 5, which corresponds the kinematic region of HERA, this difference is not so large and can be neglected.
However, accounting of energy conservation is beyond the NLOBK, which is the subject of this paper. We believe, that it is necessary to find the corrections of next-to-next-to leading order to treat the energy conservation on a theoretical footing.

V. CONCLUSIONS
Concluding, we wish to formulate clearly the three stages of our approach: 1. For Q 2 s r 2 < 1 (perturbative QCD region) we suggest to use the linear evolution equation(see Eq. (56): We have discussed this equation in the section III-D-4. 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. (119). The value of the saturation momentum is given by Eq. (85). 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 .
3. For Q 2 s r 2 1 (saturation region), we propose to use the solution to the non-linear equation which has the following form: This equation does not take into account the corrections relating to energy conservation, but it is simple, and for large region of z the solution to the non-linear equation of section IV-D-1, is described quite well by Eq. (121).
In general, we developed the approach based on the DLA approximation of perturbative QCD, and on the nonlinear evolution for the leading twist approximation. We considered the non-linear equations, which arise in different kinematic regions , and discussed solutions to them.
We believe that our suggested approach which takes into account the main features of the NLO corrections to the BFKL kernel, and which is likely to describe the experimental data for DIS processes. Indeed, Fig. 5-b shows that the energy dependence of the saturation scale, as well as the value ofγ, are very close to the values that stem from saturation models, that describe the data quite well (see Ref. [40] for example).