High energy evolution for Gribov-Zwanziger confinement: solution to the equation

In this paper we solved the new evolution equation for high energy scattering amplitudethat stems from the Gribov-Zwanziger approach to the confinement of quarks and gluons. We found that (1) the energy dependence of the scattering amplitude turns out to be the same as for QCD BFKL evolution; (2) the spectrum of the new equation does not depend on the details of the Gribov-Zwanzinger approach and (3) all eigenfunctions coincide with the eigenfunctions of the QCD BFKL equation at large transverse momenta $\kappa\,\geq\,1$. The numerical calculations show that there exist no new eigenvalues with the eigenfunctions which decrease faster than solutions of the QCD BFKL equation at large transverse momenta. The structure of the gluon propagator in Gribov-Zwanziger approach, that stems from the lattice QCD and from the theoretical evaluation, results in the exponential suppression of the eigenfunctions at long distances and in the resolution of the difficulties, which the Colour Glass Condensate (CGC) and some other approaches, based on perturbative QCD, face at large impact parameters. We can conclude that the confinement of quark and gluons, at least in the form of Gribov-Zwanziger approach, does not influence on the scattering amplitude except solving the long standing theoretical problem of its behaviour at large impact parameters.


INTRODUCTION
It is well known that perturbative QCD suffers a fundamental problem: the scattering amplitude decreases at large impact parameters (b) as a power of b. Such behaviour contradicts the Froissart theorem [1] and, hence, perturbative QCD cannot lead to an effective theory at high energy.
In particular, the CGC/saturation approach (see Ref. [2] for a review) which is based on perturbative QCD, is confronted by this problem [3,4]. At large b the scattering amplitude is small and, therefore, only the linear BFKL (Balitsky, Fadin, Kuraev and Lipatov) equation [5] describes the scattering amplitude in perturbative QCD. It is known that the eigenfunction of this equation (the scattering amplitude of two dipoles with sizes r and R) has the following form [6] φ γ (r, R, b) = One can see that φ γ (r, R, b) at large impact parameter b decreases as a power of b. In particular, such a decrease leads to the growth of the radius of interaction as a power of the energy [3,4], resulting in the violation of Froissart theorem. Since it was proven in Ref. [6] that the eigenfunction of any kernel with conformal symmetry has the form of Eq. (1), we can only change the large b behaviour by introducing a new dimensional scale in the kernel of the equation. A variety of ideas to overcome this problem have been suggested in Refs. [4,[7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. In our previous paper [27] we used the Gribov-Zwanziger approach [28][29][30][31][32][33][34][35][36][37][38] for the confinement of quarks and gluons to fix this non-perturbative scale. We derived the generalized BFKL evolution equation, which incorporates this new dimensional scale, and demonstrated that this equation leads to the exponential decrease of the scattering amplitude at large b. We will discuss both the equation and large b behaviour of the solution in the next section which has a review character.
The goal of this paper is to find the solution to this new equation. In section III we consider the general properties of the spectrum and eigenfunctions, which follow from an analytical approach. In particular, we prove that the eigenfunctions of Eq. (1) describe the eigenfunctions of the new equation at short distances r R. In section IV we concentrate our efforts on the numerical solution of the equation. We show that all eigenvalues of the new equation, that generate the power energy increase of the scattering amplitude, coincide with the massless BFKL eigenvalues. However the eigenfunctions have quite a different behaviour in comparison with the eigenfunction of the massless BFKL equation and they crucially depend on the input from Gribov-Zwanziger confinement approach. Finally, in Section V we discuss our results and future prospects.

II. BFKL EVOLUTION EQUATION FOR GRIBOV-ZWANZIGER CONFINEMENT -A RECAP
A. Gribov -Zwanziger confinement: gluon propagator As we have alluded that in Ref. [6] it is proven, that eigenfunctions of Eq. (1) have the same form for all kernels with conformal symmetries. Hence we have to modify the kernel of the BFKL equation introducing a new dimensional scale of the non-perturbative origin. In other words, we need an approach which models the confinement of quarks and gluons. Among numerous approaches to confinement, the one proposed by Gribov, [28][29][30][31][32][33][34][35][36][37][38] has special advantages, which makes it most suitable for discussion of the BFKL equation in the framework of this hypotheses. First, it is based on the existence of Gribov copies [28] -multiple solutions of the gauge-fixing conditions, which are the principle properties of non-perturbative QCD. Second, the main ingredient is the modified gluon propagator, which can be easily included in the BFKL-type of equations. Third, in Ref. [24] (see also ref. [64]) it is demonstrated that the Gribov gluon propagator originates naturally from the topological structure of non-perturbative QCD in the form: where χ top = µ 4 is the topological susceptibility of QCD, which is related to the η mass by the Witten-Veneziano relation [65,66]. This allows us to obtain the principal non-perturbative dimensional scale, directly from the experimental data. However, it is shown in Ref. [27] that propagator of Eq. (2), which vanishes at q = 0, does not lead to the exponential suppression of the scattering amplitude at large impact parameters (b). Fortunately, the lattice calculation of the gluon propagator generates the gluon propagator with G (q = 0) = 0 (see Refs. [39,43,48] and references therein), in explicit contradiction with Eq. (2).
As we have mentioned, at high energies q is a two dimensional vector, which corresponds to transverse momentum carried by the gluon. Introducing we can re-write Eq. (3) in the form: where we use notations: B. The BFKL equation in momentum representation.
The BFKL equation for Gribov-Zwanziger gluon propagator has been derived in our previous paper [27], using the procedure that has been described in Ref. [22].
It has two parts: the gluon reggeization and the emission of gluons. The first one has a general form [5]: where G(q) is given by Eq. (3). The analytical expression for Eq. (7) we will discuss below (see also appendix A of Ref. [27]). The emission kernel has been calculated in Ref. [27] using the decomposition of Eq. (2). Indeed, using this decomposition we can treat the production of the gluon as the sum of two sets of the diagrams (see Fig. 1) withM 2 = iµ 2 and withM 2 = − iµ 2 .
We sum the first diagrams of the gluon emission shown in Fig. 1 to find the vertex Γ µ (q, q ) for the kernel of the BFKL equation. It is easy to see that the sum shown in Fig. 1, leads to the Lipatov vertex that has the following form [22]: where p 1,µ and p 2,µ are the momenta of incoming particles (see Fig. 1 for all notations). Using Eq. (7) and Eq. (8), the BFKL equation for Gribov-Zwanziger confinement takes the form (for Q T = 0 † ): This equation looks similar to the BFKL equation for a massive gluon [22] in the non-abelian Yang-Mills theories with a Higgs particle, which is responsible for mass generation. However, we do not have a contact term in Eq. (9). As we have discussed in [27] the absence of a contact term in our equation is a direct indication that Gribov-Zwanziger confinement does not lead to a massive gluon.
Assuming that φ(q) depends only on |q|, we can integrate the emission kernel over the angle and in terms of the variable of Eq. (6), Eq. (9) takes the form: In Eq. (10a) -Eq. (11e) we use the variables which are given in Eq. (6). † Q T is the momentum transferred by the BFKL Pomeron, a conjugate variable to the impact parameter.

III. THE BASICS OF THE SPECTRUM FOR THE MASTER EQUATION
A.
The equation for the eigenfunctions of the massless BFKL equation As has been mentioned the eigenfunctions of the massless BFKL equation form the complete and orthogonal set of functions. Hence, we can expect that the solution to the master equation can be written as the sum over these functions. By this reason we find instructive to consider how the emission kernel of our master equation (see Eq. (10a)) acts on the eigenfunctions of Eq. (12): where the kernel of massless BFKL χ (γ) has the form [5]: where ψ(z) is the Euler ψ-function (formula 8.36 of Ref. [57]). In Eq. (13b) the region of integration over κ is divided in two: κ ≤ κ and κ ≥ κ. In the first region the new variable is introduced t = κ /κ, while in the second the new variable is t = κ/κ . In this way we have both t's in the region (0, 1). In addition, we subtracted in Eq. (13b) (terms with 1 in the numerators of the equation) the contribution from the Regge trajectory (see Eq. (7)). Using formula 3.211 of Ref. [57] we can express these integral over t through the Appel F 1 function (see Ref. [57] formulae 9.180-9.184): κ +m m m(4κ +m) + 3κ +m From Eq. (13b) and Eq. (15) (see also Fig. 2-a) we can see that P (κ, γ) is rather small and decreases at large positive l = ln κ.
Since in Eq. (13b) we subtract the reggeization term we have re-defined the kinetic term in Eq. (10a), subtracting from T (κ) of Eq. (11a) function L (κ) which is equal We denote T (l) is plotted in Fig. 2 Using function P (l, γ) andT (l) we see that our equation for the function e (γ−1)l has the form In Fig. 3 we plot fixing γ = 1 2 + i ν. Fig. 3-a gives Eq. (19) for m = m 0 = 0, which corresponds to the Gribov gluon propagator while in Fig. 3-b the energy is plotted for the lattice QCD gluon propagator with m = 1.27 and m 0 = 3.76 [39].
One can see that the wave functions of the massless BFKL show up as the eigenfunctions of the master equation in the kinematic region of large l = ln κ 1. Generally speaking, it means that the eigenvalues of Eq. (10a) could be (1) the same as the massless one or (2) could be selected out due to behaviour at small κ, leading to the set of the eigenvalues, which is more restricted than the massless BFKL one. In addition, of course, could be some discrete states, whose wave functions decrease steeper than κ γ−1 at large values of κ. From the numerical solution (see below) we see that there is no selection and all energies of the massless BFKL equations occur as the eigenvalues of the master equation. We can understand this, since the BFKL eigenvalues of Eq. (12) are twice degenerate. One can see this since the eigenvalues of the massless BFKL equation do not depend on the sign of ν. At l < 0 we expect the the eigenfunction of the massive BFKL equation should be constant. Replacing this behaviour by the boundary condition φ(l = 0; ν) = Const we see that we can satisfy this boundary condition choosing φ(l) = C 1 φ BFKL (l, ν) + C 2 φ BFKL (l, ν) = sin (l ν + ϕ) where ϕ is the phase. One can see that we do not bring any selection with this procedure.
However, the eigenfunctions of the massless BFKL equation appear the eigenfunctions of Eq. (18) also for κ 1 (l = ln κ −1) (see Fig. 3-b) but only for the case of m = 0 and m 0 = 0 with the eigenvalue E 0 = T (κ = 0) = 0.866 for any value of ν(brown line in Fig. 3-b). The independence of ν means that the eigenvalue ω 0 is infinitely degenerate. Fig. 2 Fig. 2-c), while T (l) approaches a constant (see Fig. 2 In principle, such solutions could be rejected for the master equation if the behaviour at small κ cannot be matched with the behaviour at large κ. However, it looks very unlikely. Indeed, any function of the following type: φ(κ) = P n (2 κ−1) Θ(1−κ) where P n (z) is the Legendre polynomial (see Ref. [57] formulae 8.91), is orthogonal to φ(κ) = Const at l < 0 (for n > 1) and satisfies Eq. (18). The numerical calculations, which we will discuss below, confirm that E = E 0 appears as the eigenvalue of the generalization of the BFKL equation (see Eq. (10a)).  [39]. The brown line in Fig. 3

B. General features of the spectrum
Following the general pattern of Ref. [22] we can re-write Eq. (10a) in the coordinate space, introducing The equation takes the following form with whereκ = −∇ 2 r is the momentum operator and G(r) is equal to For large r, G(r) exponentially decreases as one can see from Eq. (62). Hence, at large r Eq. (21) takes the following form: with the eigenfunctions that have the following form Denoting the large asymptotic behaviour of the eigenfunction as Ψ(r) , we see that the energy is equal to On the other hand, in the region of small r Eq. (22) reduces to the massless QCD BFKL equation (see the previous section and [5,6]: where [6] H 0 = ln p 2 + ln |r| 2 − 2ψ(1) The eigenfunctions of Eq. (27) are Ψ(r) = r 2(1−γ) , and the eigenvalues of Eq. (27) can be parametrized as a function of γ (see Eq. (14)). Therefore, for r → 0 we have the eigenvalue which is equal to From Eq. (26) and Eq. (29) we can conclude, that the value of a and γ are correlated, since Based on Eq. (29) (see also the previous section) we expect that the minimum eigenvalue is equal to χ( 1 2 ) = − 4 ln 2.  On the other hand the estimates of Eq. (26) contradict the result of Ref. [27] that the eigenfunction of the master equation with the Gribov's gluon propagator exhibit the power-like decrease at long distances. We believe that a resolution of this inconsistency is intimately related to the definition of Ψ(r). In particular, instead of Eq. (61) we suggest to introduce the following transform to the coordinate space. First we introduce a newφ(κ) = G −1 (κ) φ(κ). For this function Eq. (9) takes the form: The eigenfunction in the coordinate space has the form: and the master equation has the form of Eq. (21) with the potential energy, which has the different form: with One can see that for m = m 0 = 0 the potential energy U (r, r ) ∝ ln (r) and Eq. (24) turns out to be incorrect. For m = 0 and m 0 = 0 the potential energy decreases exponentially at long distances. Hence, in this case Eq. (24) holds.
C. Eigenfunctions in the vicinity of E0 = T (κ = 0) As we have discussed above, there is a possibility, that the master equation has the eigenvalues in addition to the eigenvalues of the QCD BFKL equation. These states should have the wave functions that decrease much steeper that the eigenfunction of Eq. (12). From Fig. 3 one can expect that the vicinity of E = E 0 = T (κ = 0) can provide such states. Indeed, in the vicinity E → E 0 T (κ) takes the form Eq. (10a) takes the following form in vicinity of κ → 0 one can see from Eq. (37) that φ (κ) has a singularity: or, in other words the wave function has the form: where φ bg is the function which has no singularities. Since E = E 0 is multiple degenerate eigenvalue, the sum of functions of Eq. (39) is also an eigenfunction. It is instructive to note that the eigenfunctions of Eq. (39) does not appear for the QCD BFKL equation. As we have seen, the origin of such eigenfunctions is in the fact that typical κ in Eq. (37) are about the values of mass and not equal to zero.

D. Resume
Concluding this section we wish to emphasize two results that we have proved. First, the eigenvalues of the massless BFKL equation, generally speaking, are expected to be the eigenvalues of the master equation. In principle, it is possible that the behaviour of the wave functions at small values of κ could select out some of the eigenvalues of the BFKL equation in QCD. However, due to double degeneracy of each of the massless BFKL eigenvalues (see, that Eq. (14) has symmetry ν → −ν) the boundary conditions at κ → 0 does not lead to a loss of the eigenvalues of the master equation in comparison with the massless BFKL equation.
Second, it is possible that the eigenvalues of the master equation have a reacher structure than the eigenvalues of the BFKL equation in QCD. Indeed, could be states with the wave functions that are suppressed at large κ: φ(κ) κ − 1 2 +i ν . An example of such function could be Eq. (38). As we see from The separate problem is the state with the wave function that decreases steeper than the eigenfunction of the massless BFKL equation but with the eigenvalue which is smaller than E min = −4 ln 2.
At the moment we cannot answer this question without finding the numerical solution to the master equation.

A. General approach
Generally speaking we need to solve the equation which has the following structure: where K is defined in Eq. (10b). The advantage of using Eq. (10b) in comparison with Eq. (10a), have been discussed in Appendix B of Ref. [27]. For numerical solution we discretize the continuous variables κ and κ using the logarithmic grid {κ n } with N + 1 nodes where the values of κ min and κ max are fixed. In the most details we consider the case with κ min = 10 −10 , κ max = 10 65 and N = 2000, but we investigated the dependence of the solution on the values of κ min , κ max and N .
In the discrete variables Eq. (40) can be approximated in the form Introducing the notations: φ (κ n ) ≡ φ n and κ m ∆ κ K (κ n , κ m ) ≡ K nm we can re-write Eq. (42) in the matrix form where vector φ has N + 1 components φ n and K is (N + 1) × (N + 1) matrix. We need to find the roots of the characteristic polynomial p (E) of the matrix K − E I where I is the identity matrix. Hence, we need to solve the secular equation We believe that we need to compare our numerical procedure, with the equation which has the analytical solution, both to check the accuracy of our numerical estimates and to evaluate our transition to the continuous limit. Recall, that the numerical solution gives the discrete spectrum of the eigenvalues instead of the continuous one. In addition we solve the equation, which was derived in Ref. [22] for non-abelian gauge theory with the Higgs mechanism for mass generation. This theory is not QCD, since it has no confinement of quarks and gluons. However, it has the same colour structure as QCD and introduces the dimensional scale: the mass of Higgs boson. Solving the BFKL equation for this theory we could find out what is more essential: the new dimensional scale or specifics related to the confinement. We will call below this approach "the model" and the main BFKL equation for this model takes the form: This equation has been investigated in detail. In Ref. [22] it was proven that the solution of this equation coincide with the solution to the massless BFKL equation, which is known analytically and can be used as a control of the accuracy of our numerical calculations.

B. Eigenvalues: the general characteristics
The eigenvalues of these four equations are shown in Fig. 5. One can see that (1) numerical estimates do not show the discrete eigenvalues with energy E < E min = −4 ln 2, where E min is the minimal energy of the massless BFKL equation; and (2) none of the eigenvalues of the massless BFKL equation has been selected out in accord with our expectations in section III-A. From Fig. 5 we see that all eigenvalues can be divided in three regions: the eigenvalues E n ≤ E 0 = T (κ = 0), the multiply degenerate eigenvalue E 0 and E n ≥ E 0 .
For E min ≤ E n ≤ E 0 there are no other eigenvalues except the massless BFKL ones. Indeed, we can describe these eigenvalues using the following formulae: β(n) = a β (n + 1), a β = c β / ln κ max /m 2 β For the QCD BFKL equation c β = 3.015 and m 2 β = κ min , while for all other three equations we can put c β = 3.140 and m 2 β = 0.0042 (see Fig. 9-a). Hence, Fig. 5-b and Table I demonstrate a new phenomenon: the eigenvalues of all three equations, which introduce a dimensional scale in the BFKL approach, turns out to be the same. Actually, they are the eigenvalues of the QCD BFKL equation, as it shows Eq. (47) and Fig. 5-b. In Fig. 5 the eigenvalues from Eq. (47) are shown by the solid lines and one can see that all these values for E n ≤ E 0 can be perfectly described by this equation. Eq. (47) can be interpreted as an indication that the transition to the continuous limit reduces to replacement 1 2 + iβ(n) → 1 2 + i ν ≡ γ. In these new variables the eigenvalues of the QCD BFKL looks familiar [2]:  (47)). This indicates that our method of numerical solving provides a good accuracy. As we can see, in both cases the lowest eigenvalue becomes quite close to E min = −4 ln 2 and difference is negligibly small (of the order of 5 × 10 −3 for κ min = 10 −10 and κ max = 10 65 ). Fig. 6 shows the dependence of the first 7 roots versus the value of κ max . One can see that when κ max grows (κ max → ∞) the distance between neighboring roots decreases rapidly, inferring the smooth transition to the continuous limit.
As we can see from    (45)), for the model, developed in Ref. [22] (Higgs, see Eq. it is so large that we can expect something like Bose-Einstein condensation at this energy. The general structure of the eigenfunction at this value of energy we have discussed in section III-C and will consider below. For the scattering amplitude all these eigenfunctions correspond to the cross sections that decrease as a power of energy and, because of this, do not show up in the high energy scattering processes.  For the QCD BFKL equation(see Eq. (45)) the eigenfunctions are given by Eq. (12) and for the numerical solutions they take the form: The eigenfunction for Eq. (46) have been discussed in Ref. [22] and can be described as follows: where  Several examples of the eigenfunctions for Eq. (10b) are shown in Fig. 7. One can see that the number of zeros follows the usual pattern of a quantum mechanical approach: the minimum energy state has no zeros. The next has one and so on. At large κ φ n (κ) ∝ sin (α β n ln κ) or in other words φ n (κ) = C 1 φ BFKL κ; 1 2 + iβ n + C 2 φ BFKL κ; 1 2 − iβ n , where φ BFKL (κ; γ) are given by Eq. (12).
For κ ≥ 1 all eigenfunctions can be parameterized in the following way: where Parameter β n has the simple form defined by Eq. (47) β n = a β (n + 1), a β = 3.140 ln(κ max /m 2 β ) It was found that for the Gribov's propagator (m = 0, m 0 = 0) and for the propagator with m = 0, m 0 = 0 the same m 2 β = 0.0042 can be used (see Fig. 9-a). While ϕ n needs a bit more complicated parametrization ϕ n = a ϕ,0 + a ϕ,1 n + a ϕ,2 (n − n ϕ ) 3 For κ max = 10 65 we obtain the following values for parameters a ϕ,i : a ϕ,0 = 0.486 (1.520); a ϕ,1 = 0.0350 (−0.0223); a ϕ,2 = 0.425 × 10 −4 (1.211 × 10 −4 ); n ϕ = 21.71 (21.76) (56) In Eq. (56) the values of a ϕ,i for the case m = 0, m 0 = 0 are given in parentheses. Fig. 8 shows the n dependence of β n and φ n for n ≤ 40. One can see that the linear dependence of Eq. (54) holds for β n , but ϕ n shows a more complicated pattern: ϕ n ∝ n with variations described by Eq. (55).   As far as dependence on κ max of β n , ϕ n and a n is concerned, we found that they have quite a simple scaling property, which allows to relate their values obtained with different κ max : values of some parameter P n = P (n; κ max1 ) corresponding to κ max1 fit well to results P (n; κ max2 ) after simple change of the scale of n P (n; κ max1 ) ≈ P (nS; κ max2 ), S = ln(κ max2 ) ln(κ max1 ) Such scaling behaviour for β n follows from Eq. (54), while Fig. 10 illustrates such κ max scaling for parameter ϕ n . Figure shows parameters obtained at different κ max in the range from 10 20 to 10 65 . Sets of ϕ were scaled on n with proper coefficient (see Eq. (57)) to κ max = 10 65 . One can see that "scaled" points are in good agreement with original values for κ max = 10 65 . For the Gribov's propagator the eigenfunction φ n (κ) ∝ κ in the region of small κ. In other words, the eigenfunctions in the coordinate space exhibit the power-like decrease at long distances.  (Fig. 11-a) demonstrates scaling property of this parameter for the case with Gribos's propagator. Right plot (Fig. 11-b) shows an for both variants (i.e. Gribov's and QCD propagators) corresponding to eigenfunctions with κmin = 10 −10 and κmax = 10 65 .
It turns out that for small κ we can use Eq. (52), which introduce the dependence of parameter a n on n. Fig. 11-a shows the scaling behaviour of Eq. (57) for a n in the case of Gribov's propagator. It should be stressed that the same scaling behaviour holds for the case of m = 0. Fig. 11-b shows the dependence of a n on n. From this figure we see that a n ∝ n 2 for the Gribov's propagator with a n = 0.3 at n = 0. In other words, the typical κ turns out to be in the region of κ = 0.27 ÷ 0.75. It should be stressed that Eq. (52) describes quite well the behaviour of the eigenfunctions both at large κ ≥ 1 and at small κ ≤ a n , but at κ ∼ a n Eq. (52) does not lead to a good fit of the eigenfunctions. For the case of the lattice QCD propagator (Eq. (11e) with m = 1.27 and m 0 =3.76) Fig. 11-b shows that a n decreases with n approaching a n ≈ 0.4 at n → 40.
One can see that for small n the typical values of κ ∼ 1. and the range of typical κ is 0.4 ÷ 1. We would like to stress that the value of a n cannot be viewed as a typical scale for the κ dependence of the eigenfunctions. Indeed, one can see directly from Fig. 7 that the typical κ is about κ ∼ 1 for both cases. D. Eigenfunctions for En = E0 = T (κ = 0) As we have discussed above, the solution leads to the multiple degenerate state at E n = E 0 = T (κ = 0). Using Eq. (47) we can estimate the value of β 0 : viz. E 0 = −2 ψ(1) + ψ 1 2 + i β 0 + ψ 1 2 − i β 0 . Corresponding eigenfunction index n, where the first degenerate eigenvalue appears, can be estimated using Eq. (54): degenerate eigenfunction sequence starts when β n reaches value β 0 . For the Gribov's propagator β 0 ≈ 0.85 (see Fig. 8-a) leads to the value n = 41 in the case with κ max = 10 65 , while for RGZ gluon propagator β 0 ≈ 0.92 gives n = 45 (see Fig. 8-a).
At such values of n the κ behaviour of the eigenfunctions shows a discontinuity in κ: of course, numerical values of eigenfunctions at ln (κ/κ min ) = n∆ (see Eq. (41)) are still finite, but the values in the neighbouring nodes have a different sign and derivative, indicating that eigenfunctions have pole in κ located somewhere between nodes.
The structure of the eigenfunctions with this eigenvalue is rather simple for the lattice (RGZ) QCD gluon propagator (see Eq. (11e) with m = 1.27 and m 0 = 3.76) and it is close to one, that has been discussed in Ref. [22] for Eq. (46). The eigenfunctions for E n = E 0 have poles in κ = κ p,n as it has been shown in section III-C. Actually the minimal value of κ p,n is equal to κ min (strictly speaking the first pole is located somewhere between the first node κ 0 = κ min and the next node κ 1 = κ min exp(∆ κ ) ≈ κ min (1 + ∆ κ ). With each increase of n pole moved exactly to the next interval on κ (i.e. second pole is located between κ 1 and κ 2 and so on). This sequence terminates when the pole reaches the maximal value of κ = κ 0 ≈ 3, where κ 0 is the location of the first zero of the eigenfunction of Eq. (52). All eigenfunctions with E n = E 0 have the same number of zeros. All these features can be seen from Fig. 12.  Fig. 12-a shows φn (κ) versus κ with n = 50, 150, 250 for which En = E0 and φ400 (κ) with E400 > E0. One can see that all eigenfunctions with En = E0 have the same number of zeros. The same eigenfunctions are shown in Fig. 12-b but in the region of small κ ≤ 100. It is clearly seen that φn (κ) have the pole whose position moves from κmin to κ ≈ 1. Function φ400 (κ) has a different number of zeroes, which corresponds to increase of the value of β in Eq. (58).
Generally for m = 0 and m 0 = 0, the eigenfunction can be approximated by the following expression: Eq. (58), having β 0 (see Fig. 9-b), which does not depend on n, reflects the fact that for all these states the behaviour of the eigenfunctions at large κ ≥ 1 can be described by one function with the same number of zeros. Fig. 13 shows the n-dependence of other parameters of φ (approx) n (κ) (see Eq. (58)). One can see that both the position of the pole κ p,n and its residue a p,n are proportional to n (ln (κ p,n ) ∝ n, ln (a p,n ) ∝ n), while the parameters of φ n (κ; Eq. (52)): a n and ϕ n do not depend on n in the range n = 50 − 250 which corresponds to E n = E 0 . The value of a n = 1 gives us the typical transverse momentum q = µ (see Eq. (6)). For Gribov's gluon propagator the structure of the eigenfunction with E n = E 0 is much more complex. First, one can notice from Fig. 14 that the number of zeros are not the same for these eigenfunction but β 0 in Eq. (59) changes with n rather slowly (see Fig. 8-a for n in the region n = 41 ÷ 70).
Second, we see that φ n (κ) with E n = E 0 have two poles. First pair of poles κ p,1 < κ p,2 appear near κ ≈ 1. With increase of n the smaller one (κ p,1 decreases, while κ p,2 increases. On each increment of n only one of κ p,i moves to the neighbouring κ interval between nodes. So the distance between these poles (in terms of index of κ-nodes) each time increases exactly on 1. The contribution of each of these poles vanishes at κ → 0 and residues of these poles can have the same or opposite sign.
Third, the position of the poles are in the region κ = 0.1÷10 and they exist also in the eigenfunctions with E n > E 0 . Fourth, two poles in the eigenfunctions have close positions. Eq. (59) reflects the main features that we have discussed but the actual structure of the eigenfunction turns out to be much more complex.
However, for the Gribov's propagator (m = 0, m 0 = 0) the eigenfunction vanishes at κ = 0 and can be approximated as follows: The appearance of two poles in Eq. (59) looks natural (see section III-C) due to multiple degeneracy of this eigenvalue. Indeed, due to this the sum of two functions with one pole in each, is also the eigenfunction.  demonstrates that the region of n for the degenerate states with E n = E 0 is very narrow but the structure of the eigenfunction with two poles lasts for E n > E 0 . For E n > E 0 the eigenfunctions take the following form for κ > 1: with β n ∝ n for both Gribov's and lattice QCD gluon propagators as it is seen from Fig. 9-b. One can see from Fig. 5-a, that forκ 1 E n tends to the same asymptotic values for all four equations that we have studied in this paper. For κ < 1 the eigenfunctions have the similar structure as for E n = E 0 , i.e. they have one pole (see Fig. 13-a) for the lattice gluon propagator and two poles for the Gribov's one (see Fig. 15-b). Since these eigenvalues corresponds to the Pomeron intercept ω I P = −ᾱ S E n < ω 0 = −ᾱ S E 0 < 0, they do not contribute to the high energy behaviour of the scattering amplitude.

V. THE SCATTERING AMPLITUDE
A. Green's function of the BFKL Pomeron for the Gribov-Zwanziger confinement The Green's function of the BFKL equation on Y = ln(1/x) representation takes the general form At high energies the main contribution stems from the minimal energy. However, we cannot restrict ourselves by calculating only one term in Eq. (61). To demonstrate this, we use Eq. (52) for the approximate eigenfunctions, which can be written in the following form: where Φ n (κ) = α n (κ + m) The eigenvalues of Eq. (47) we calculate in diffusion approximation in which: where ω BFKL = 4 ln 2ᾱ S ; D = 14 ζ(3)ᾱ S . Therefore, in this approximation the Green function takes the form Taking the integral over β in Eq. (65) we obtain the following Green's function at large values of Y : One can see that at large Y , G (Y, κ fin |0, κ in ) ∝ (D Y ) −3/2 e ω BFKL Y , which should be compared with the massless BFKL case for which G (Y, κ fin |0, κ in ) ∝ (D Y ) −1/2 e ω BFKL Y . These estimates show that we need to sum the contributions of the eigenvalues in the vicinity of n = 0. The source of such contribution can be seen from the first two components of the sum over n in Eq. (61), which can be re-written as follows: where For κ max → ∞ ∆ω 1 → 0, however, the product ∆ω 1 Y at large κ max and Y is undefined. Hence, we have to perform numerical estimates for the sum of Eq. (61) to determine the answer. We will discuss such kind of estimate below. At the moment we wish to emphasize that since in the vicinity of n = 0 the spectrum of the master equation coincide with the QCD BFKL equation, one can see that the influence on the asymptotic behaviour of the scattering amplitude due to Gribov-Zwanziger confinement is rather small. Indeed, as we have seen from the above estimates, we obtain extra suppression of the scattering of the order of 1/(DY ) for our case. Using Eq. (61) we can find the scattering amplitude which will be equal to where N (Y = 0, κ in ) is the initial condition for the scattering amplitude at Y = 0. In Fig. 16 N (Y = 0, κ in ) is taken to be equal to 1/(κ + 1) 2 . We plot in Fig. 16 the contours on which function κ N (Y, κ) (see Eq. (68) is constant. In QCD we have transverse momentum distribution, which depends on | ln κ |,. The QCD evolution results in the increase of | ln κ | with Y . Such an increase leads to two possible branches (depending on different sign of ln κ) with increasing and decreasing average transverse momentum. Such a behaviour is seen in Fig. 16.
For our master equation one can see that the confinement cuts the small transverse momenta and the average κ T are larger than the values of κ T in initial conditions , which we consider κ in = 1, and they grow with Y . Therefore, introducing the Gribov-Zwanziger confinement in the framework of the BFKL equation we obtain the transverse momentum distribution, which is determined by the behaviour of the scattering amplitude at large transverse momenta (at short distances), where we can trust the perturbative QCD approach.

VI. CONCLUSIONS
In this paper we solved the new evolution equation for high energy scattering amplitude that stems from the Gribov-Zwanziger approach to the confinement of quarks and gluons (see Eq. (10a)). The results of this solution we find quite surprising and instructive for future development of high energy physics.
First, the energy dependence of the scattering amplitude turns out to be the same as for QCD BFKL evolution. In particular, the eigenvalues of the new equation, which exceed ω 0 = −ᾱ S E 0 = −ᾱ S T (κ = 0), coincide with the QCD BFKL equation. Second, the spectrum of the new equation does not depend on the details of the Gribov-Zwanzinger approach and coincides with the set of the eigenvalues of the model: non-abelian gauge theories with the Higgs mechanism for mass generation, developed in Ref. [22]. This model has no relation to a QCD approach except having the same colour structure. These features support the ideas, that come out from the analytical analysis of the equation: the main influence of the confinement is in taking off the double degeneration of the QCD BFKL equation, which shows up in independence of the spectrum of the QCD BFKL equation on the sign of ν (see Eq. (14)). Third, all eigenfunctions coincide with the eigenfunctions of the QCD BFKL equation at large transverse momenta κ ≥ 1.
The numerical estimates show that there exist no new eigenvalues with the eigenfunctions that decreases faster than the eigenfunction of the QCD BFKL equation at large transverse momenta.
The eigenfunctions of the master equation with the Gribov's gluon propagator tends to zero at small transverse momenta. In the coordinate representation it means that the eigenfunctions exhibit the power-like decrease at long distances, leading to the power-like decrease in the impact parameters and, therefore, to the severe problem with Froissart theorem and s-channel unitarity (see Refs. [3,4,27]). In other words, the gluon propagator which tends to zero as the Gribov's propagator does, cannot solve the problem with large b dependence of the scattering amplitude in the CGC approach. However, the structure of the gluon propagator in Gribov-Zwanziger approach that stems from the lattice QCD estimates and from the theoretical evaluation (see Refs. [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56]), leads to the gluon propagator which tends to a finite value at zero transverse momentum (G (q → 0) = 0). This results in the exponential suppression of the eigenfunction at long distances and in the resolution of the difficulties, that the CGC approach as well as other approaches, based on perturbative QCD, faces at large impact parameters.
For the intercept ω = −ᾱ S T (κ = 0) we have the multiple degeneration of this eigenvalue, which is strongly correlated with the new dimensional parameter that we introduced to the theory from the confinement. This degeneration looks as Bose-Einstein condensation but it does not contribute to the scattering amplitude at high energy.
We calculate the momentum distributions of the scattering amplitude and found that the typical transverse momentum increases with energy and become independent of the typical confinement scales that we have introduced in our equation.
Therefore, to our surprise, we have to conclude that the confinement of quark and gluons, at least in the form of Gribov-Zwanziger approach, does not influence on the scattering amplitude except solving the long standing theoretical problem of the large impact parameter behaviour of the scattering amplitude. This is a very optimistic message for the CGC approach, but before coming to the strong conclusions we have to check the solution to the non-linear equation with the new kernel. This will be our next problem.