Energy evolution and the Bose-Einstein enhancement for double parton densities

In this paper we found that the Bose-Einstein enhancement generates the strong correlations, which increase with energy in the BFKL evolution. This increase leads to the double parton densities ( $\Phi$), that are much larger than the product of the single parton densities ($\phi$). However, numerically, it turns out that the ratio $\Phi/\phi^2 \propto \Lb 1/x\Rb^{\delta_2}$ with $\delta_2 \sim \bas/\Lb N^2_c - 1\Rb^{2/3}\,\,\ll\,\,1$ and we do not expect a large correction for the accessible range of energies. However, for $N_c=3$ it tuns out that $\delta_2 = 0.07 \Delta_{\rm BFKL}$ where $\Delta_{\rm BFKL}$ is the intercept of the BFKL Pomeron and we can anticipate an substantial increase for the range of rapidities $Y \sim 20$.It is shown that all $1/(N^2_c -1)$ corrections to the double parton densities stem from the Bose-Einstein enhancement.

For a long time the double parton distribution functions (DPDFs) and their DGLAP evolution * have been of interest to the theoretical high energy community, and has been discussed in detail . On the other hand the BFKL evolution † and related to them the double gluon densities (transverse momentum distributions(TMD2s)) have attracted less interest of the theorists, in spite of the fact, that they give the simplest way to estimate the possible correlations in the QCD parton cascade at high energies, where experimental observations of the double parton interactions [32][33][34][35][36][37][38] were made.
In this paper we re-visit the evolution equation in the BFKL kinematic region of small x, where partons are either gluons or colorless dipoles. In the coordinate representation we use colorless dipoles as partons, while in the momentum representation, it is more convenient to discuss the parton cascade in terms of gluons. This evolution equation was written in Ref. [39](see also Refs. [40]) for the double gluon densities Φ (x 1 , p 1,T ; x 2 , p 2,T ) ≡ Φ (Y − y 1 , p 1,T ; Y − y 2 , p 2,T ) with respect to rapidity (Y ) of the initial hadron (projectile). Note that we use notation Φ (x 1 , p 1,T ; x 2 , p 2,T ) for the double gluon density, and φ (x, p T ) for the single gluon density. x i is the fraction of energy of the gluon 'i' while p i,T denotes its transverse momentum. In the reference frame where the initial hadron is fast moving ln(1/x i ) = Y − y i , Y denotes the rapidity of the projectile(hadron) and y i the rapidity of the parton i. This evolution answers the question, what are the multiplicities of two colorless dipoles in one dipole, that moves with rapidity Y . We believe that in the spirit of the BFKL evolution we need to answer a different question, what is the multiplicity of two gluons with rapidities y 1 and y 2 , if we know their multiplicities at y 1 = y 0 1 and y 2 = y 0 2 . Therefore, the first goal of our paper is re-write the evolution equation in a convenient form to answer this question. It turns out that such evolution has been discussed in Refs. [41][42][43] in the framework of the CGC approach (see Ref. [44] for the review of this approach). The evolution equations have been derived in these papers in the desired form, for y 1 = y 2 . It was shown, that it is not necessary to take into account the non-linear corrections for the double (and multi) gluons densities, and that the evolution equations reduce to the BFKL evolution. In Ref. [43] this evolution is written taking into account the corrections of the order of 1/N 2 c , where N c is the number of colors. The linear evolution equations for the double gluon density are closely related to the BKP equations [56] for the scattering amplitude, and 1/N 2 c corrections to these equations have been discussed in Refs. [53][54][55][59][60][61][62].
The second goal is to include the Bose-Einsten enhancement coming from the correlations of identical gluons, in the evolution. Bose-Einstein correlations have drawn considerable attention recently, since they give essential contributions to the azimuthal angle correlations [47][48][49][50][51][52]. It has been shown (see Ref. [50] for example) that the Bose-Einstein enhancement leads to a significant contribution to the measured angle correlations. We believe that this fact calls for a generalization of the evolution equation by taking into account this enhancement. We show that the Bose-Einstein enhancement is responsible for a term in the linear evolution which is suppressed as 1/ N 2 c − 1 , which has been found in Ref. [43]. In other words, we state that all corrections of the order of 1/N 2 c in the evolution equations for the double gluon density, stem from the Bose-Einstein enhancement. In particular, the symmetry between the azimuthal angle ϕ and the angle π − ϕ does not appear. This symmetry, which is not based on the principle features of our CGC approach, reveals itself in the scattering amplitudes, but we do not find any indication of it in the double gluon densities. We wish to stress that, the Bose-Einstein corrections are closely related to the BKP equations [56] and generate the energy behaviour of the twist four operator which increases as s ∆4 , with ∆ 4 > 2∆ 2 , ∆ 2 denotes the intercept of the BFKL Pomeron (see Refs. [53][54][55]).
We have discussed the Bose-Einstein enhancement for the DGLAP evolution [57], and have shown that it changes considerably, the high energy behavior of the DPDFs. In particular it turns out that the widely used assumption: does not hold, even at small x 1 and x 2 , due the Bose-Einstein correlations ‡ .
Before describing the structure of the paper we would like to introduce the observables of interest: the parton density φ (x, p t , q T ) and the double parton density Φ x 1 , p 1,T ; x 2 , p 2,T ; q T . The single gluon density characterizes the multiplicity of gluons with fraction of energy x and transverse momentum p T at q T = 0, and it can be written as follows † The Balitsky, Fadin, Kuraev and Lipatov (BFKL) equation [46] is written for evolution in x (the energy scale of the process) assuming thatᾱ S ln(1/x) ∼ 1 butᾱ S 1 andᾱ S ln(Q 2 ) 1. This evolution has been considered in Ref. [39][40][41][42][43] for the double gluon densities. ‡ In Ref. [58,64] it is shown that Eq. (1) is a good approximation to the solution of the DGLAP evolution equation for the double parton distribution functions at small x and large p T , but this claim is only correct, when neglecting the contributions of the order of 1/(N 2 c −1).
where Ψ denotes the partonic wave function of the fast hadron, |Ω n = n i a + (x i , k i,T , c i )|0 (|0 denotes the vacuum state) and a + (x i , k i,T , c i ) and a(x i , k i,T , c i ) denote the creation and annihilation operators for partons (gluons for small x i ) with fraction of energy x i , transverse momentum k i,T and color c i . The produced gluon has longitudinal momentum x and transverse momentum p T , while b indicates its color.
The double transverse momentum densities describe the number of gluons with (x 1 , p 1,T ) and (x 2 , p 2,T ) in the parton cascade, and it can be written with the aid of the wave function of the produced gluon Ψ ({x i , k i,T }) as follows The paper is organized as follows: In the next section we discuss the BFKL evolution of the double gluon density in the region of low x. For completeness of presentation we review the derivation of the evolution equations for the double gluon densities which was done in Refs. [41][42][43] in the framework of the dipole approach [63]. In spite of the fact that this derivation is only valid for y 1 = y 2 , it shows that these evolution equations are linear BFKL equations which are not affected by non-linear shadowing corrections. In section 3 we re-derive the BFKL equations for the double gluon densities directly in the momentum representation. In this representation, we generalize the equation for the case of y 1 = y 2 and re-write the equations in the form which is suitable for taking into account the Bose-Einstein enhancement (BEE). In this section we find the solutions to the equations without BEE. The interference diagrams that are responsible for the BEE, are discussed in section 4. In section 5 the evolution equations with BEE are proposed, and we find that Φ/φ 2 ∝ (1/x) δ2 with δ 2 ∼ᾱ S / N 2 c − 1 2/3 . Section 5 is devoted to a discussion of the energy behavior of the double gluon densities. In the conclusions we discuss our main results.

II. BFKL EVOLUTION OF DOUBLE DIPOLE DENSITIES IN THE CGC APPROACH.
In this section we discuss the evolution equation for the double gluon densities in the framework of the CGC approach. As mentioned this equation has been derived in Refs. [41][42][43], here we give a brief review both for completeness of the presentation, and to display the main features of the multi gluon densities, which are not the main subject of these papers. The partonic wave function can be expanded as the sum of Fock states with fixed multiplicity of partons (colorless dipoles): The colorless dipole is characterized by two variables: its size r i and its impact parameter b i . However, in this paper we will sometimes use a different set of variables: x i for the position of the quark and y i for the position of the anti-quark in the dipole. One can see that r i = x i − y i and b i = 1 2 (x i + y i ). α 2 n is the probability to find n dipoles with the same value of rapidity Y : The QCD cascade can be written as the linear functional equation for the following functional [41,63]: where u(r i , b i ) ≡ u i is an arbitrary function of r i and b i and y = y 1 = y 2 . It follows immediately from Eq. (5) that the functional obeys the condition: at u i = 1 The physical meaning of Eq. (7) is that the sum over all probabilities is one.

A. Balitsky-Kovchegov parton cascade
To write the evolution equation, we need to specify the QCD processes with color dipoles. The parton cascade that leads to the Balitsky-Kovchegov non-linear equation for the scattering amplitude, stems from the process of the decay of one dipole to two dipoles, which gives the main contribution at the leading order of perturbative QCD at large N c [63]. The probability of this decay is equal to Bearing Eq. (8) in mind, we can write the linear equation for Z: with the definitions and .
The functional derivative with respect to u(r, b), plays the role of an annihilation operator for a dipole of size r, at impact parameter b. The multiplication by u(r, b) corresponds to a creation operator for this dipole. Therefore,Eq. (9) is a typical cascade equation in which the first term describes the depletion of the probability due to splitting into n + 1 dipoles, while the second term is responsible for the growth due to splitting of (n − 1) dipoles into n dipoles. From Eq. (6), one can see that the multi dipole density ρ (n) (Y − Y 0 ; {r i , b i }) can be found as follows which gives the probability of finding n-dipoles with the given kinematics. From Eq. (9) we obtain For ρ (2) we have +2 There are two main features of the equation: (i) there are no non-linear corrections, and (ii) we have two contributions: BFKL evolution of ρ (2) and the contribution to ρ (2) from single parton showers. This structure is the same as in the DGLAP evolution (see Refs. [1-3, 5? , 6]).

B. 1/N 2 c corrections to the Balitsky-Kovchegov cascade
x ( y ) 2 1 x ( y ) 4 2 x 3 r 1 r 2 The lines of the same colors indicate the colorless dipole which decays in two dipoles due to emission of a gluon with coordinate z.
In Ref. [43] it is suggested to add the following term to Eq. (9): The process that is described by Eq. (15) is shown in Fig. 1. Indeed, if two dipoles consist of quarks and antiquarks with the same colors, the dipole, in which one quark from one dipole and the antiquark from another dipole, can create a different dipole which decays into q,q and a gluon. The term of Eq. (15) generates the1/(N 2 c − 1) corrections to Eq. (14) which have the following form: One can see that the structure of Eq. (16) is similar to that of Eq. (14): i.e. the production of two dipoles from the single parton cascade, and the evolution of the double density with a kernel which is different from Eq. (14), although it consists of the same elements. We need to re-write Eq. (14) with additional term of Eq. (16) in the momentum representation, to obtain the more familiar form of the double gluon density evolution equation. However, we chose a different strategy: we derive the equation directly in the momentum representation. First, we believe that we can relax the assumption that y 1 = y 2 and second, we hope that this derivation will clarify the physical meaning of the additional term of Eq. (16).
We would like to stress that the derivation, which we have discussed, is very instructive for understanding the different contributions to the evolution due to the clear physics interpretation of the dipole approach in perturbative QCD.
In this section we re-write Eq. (14) in the momentum representation using two lessons from the derivation of the previous section. The double parton density is not affected by the shadowing (screening) corrections and obeys the linear BFKL equations, and the equation should match Eq. (14) at y 1 = y 2 .
In the BFKL region we consider thatᾱ S ln (1/ 1, andᾱ S 1 and the evolution equations sum the contributions of the order of (ᾱ S ln(1/x)) n (leading log(1/x) approximation (LLA)). In the region of small x i only a gluon can be produced [44], and for the double gluon density we expect to have two equations of the following forms(see Fig. 2) where K (p T , k T ) is the BFKL kernel which is equal to [46] K (p T , k T ; q T ) = 1 The non -homogenous term takes into account the possibility to produce two gluons from a single gluon cascade. The expression for these terms is written directly from the second diagram of Fig. 2. Function Γ (p 1T , p 2T ; y 1 , y 2 ) has to be found. It is clear that forᾱ S |y 2 − y 1 | 1 two emitted gluons in Fig. 2 with rapidities y 1 and y 2 can emit gluons, and the observed gluons will be amongst them. Therefore, the general diagram, which determines the non-homogenous term is the triple BFKL Pomeron diagram of Fig. 3

-a.
This diagram can be written in the following general form where Γ 3I P denotes the triple BFKL Pomeron vertex which we will discuss below. The single gluon densities in Eq. (22) φ pr and φ are different, since φ pr is the density of the gluons in the projectile (hadron), while φ is the density Y y' Fig. 3-a:The graphical representation of the triple BFKL Pomeron diagrams. For simplicity we show the evolution at qT = 0. Fig. 3-b: the structure of the partonic wave function in the light-cone perturbative approach for the production of two gluons from the single parton cascade at x1 = x2 (y1 = y2). We denote by Dn the dominator of the propagator for the state that has n-gluons.
of the gluons in the cascade of the single gluon with rapidity Y − y . This term in the evolution was considered for the first time in Ref. [39]. We refer our readers to this paper or more detail. The contributions of this diagram to the evolution equation have the following form for y 2 > y 1 : From Eq. (17) and Eq. (18) one can see that the first term in the both equations is included in the first terms of the evolution equations since We also see that Eq. (23) does not generate the non-homogenous term in Eq. (17) if y 1 < y 2 , in the case which we consider here: y 2 > y 1 . We need to specify the triple BFKL Pomeron vertex and φ 0, p T , p 2,T , q T in Eq. (23). Actually In the diagram of Fig. 3-a the triple Pomeron vertex enters with the momentum transferred along the upper Pomeron being equal to zero. We believe that we can find this vertex directly from the non-linear Balitsky-Kovchegov evolution equation [65] for the scattering dipole amplitude: where x ik = x i − x k and b denotes the impact factor. Eq. (26) in the momentum representation, which we define as has the following form: We can build the diagram of Fig. 3-a by iterating Eq. (28), using the non-linear term as the first iteration. Returning to Eq. (23), we see the the non-homogeneous term has the form: Finally, the set of evolution equations can be re-written in the following formf orᾱ S |y 1 − y 2 | gg 1 : where ϑ (y 12 ) is the step function. Comparing Eq. (30) and Eq. (31) with Eq. (17) and Eq. (18), we see that we have found the exact form of the function Γ (p 1T , p 2T ; y 1 , y 2 ) forᾱ S |y 1 − y 2 | 1. Recall, that φ y 2 − y 1 , p 2,T , p 1,T , q T denotes the multiplicity of gluons with rapidities y 1 and transverse momenta p 1,T in the gluon with rapidity y 2 and transverse momentum p 2,T .
These equations are written for y 2 y 1 1. For y 2 ∼ y 1 we can replace φ y 2 − y 1 , p 2,T , p 1,T , q T by the DGLAP single parton density. However, we discuss the Bose-Einstein correlation which are essential at y 1 = y 2 . In the kinematic regionᾱ S |y 2 − y 1 | 1 we need to re-write the non-homogeneous term.
In the framework of the LLA, the gluon with the fraction of energy x can produce two gluons with x 1 ≈ x 2 x in the subsequent decay g(x, p T ) → g(x , p T ) + g(x 1 , p 1,T ) and g(x , p T ) → g(x , p T ) + g(x 2 , p 2,T ) as is shown in Fig. 3 We calculate the contributions of these decays to the partonic wave function using light-cone perturbative theory (see Refs. [44,64,67]). The wave function of Fig. 3-b can be written in the following form where polarization vectors λ µ (p) are defined as λ µ (p) = 0, The light-cone denominators are defined as where P − is the light-cone energy of the incoming hadron. In Eq. (32) we have omitted the color indices.
The particular solution to Eq. (40) has a simple form: whereφ in can be found from the initial conditions for the single gluon density.
where ξ 1 = ln p 2 1,T and ξ 2 = ln p 2 2,T . The general solution will be a sum of the particular solution and the solution to the homogeneous equation, and it has the following form The integrals over γ 1 and γ 2 in the first term of Eq. (46) can be evaluated using the method of steepest descend with the saddle point for both γ's close to 1 2 , where we can use the diffusion approximation (see Eq. (41)) for χ(γ). The values of γ's at the saddle point are the following: After integration over γ 1 and γ 2 in the vicinities of these saddle points we obtain the contribution: This contribution is proportional to φ (Y − y, ξ 1 ) φ (Y − y, ξ 2 ) and in agreement with Eq. (1). In the second term the integration over γ 1 + γ 2 can be taken using the method of steepest descend, leading to This integration generates the contribution which is proportional to exp ᾱ S ω 0 (Y − y) − (ξ 1 + ξ 2 ) 2 / (4ᾱ S D (Y − y)) .
Comparing this contribution with Eq. (95), one can see that it is suppressed at large values of ( Y − y) .

IV. THE INTERFERENCE DIAGRAM IN THE BFKL EVOLUTION
The interference diagram is shown in Fig. 5-a. In this diagram the t-channel gluons with the same color and with rapidities larger than y , are in colorless states. For rapidities, that are less than y , t channel gluons with rapidities y 1 and with y 2 are in a colorless state. The arguments for such a color structure of this diagram stem from the first diagram with the exchange of two identical gluons, shown in Fig. 5-b. In this diagram all emitted gluons with rapidities larger than y , can be absorbed in the solution of the evolution equation without the Bose-Einstein enhancement, and can be used as a solution of Eq. (1). In this solution the double gluon density can be viewed as the exchange of two BFKL Pomerons, shown in Fig. 5 These Pomerons carry transferred momenta Q T and −Q T , respectively, where Q T = k − p 2,T . Therefore, we first need to deal with the BFKL Pomeron with non zero transfer momentum . However, before discussing this problem we calculate the diagram of Fig. 5-c, which is the diagram of Fig. 5-b in the Born approximation. Denoting the wave function of the colorless dipole ( the onium state of a heavy quark and antiquark) by Ψ (q T , z), where q T is the transverse momentum of the quark and z its fraction of the energy. We obtain that the component of the gluonic wave function with one emitted gluon with transverse momentum p 1,T and rapidity y 1 is equal to [44,63] Ψ (1) q T , z; p 1,T , y 1 = gλ a p 1,T · λ where λ a denotes the Gell-Mann matrix and λ is the polarization vector of the gluon with the helicity λ. The single gluon density has the form In the integral over z, the lower limit is e −Y +y , but we assumed that this integral is convergent and we can safely take this limit equal to zero.
The emission of the second gluon with y 2 and p 2,T leads to the wave function Ψ (1,1) q T , z; p 1,T , y 1 ; p 2,T , y 2 = (53) Using the expression for the Lipatov vertex Γ µ which has the form [44,46] where l µ = k T,µ − k T,µ is the momentum of the emitted gluon, as well as Eq. (53) we obtain for Fig. 5-c the following contribution where ∆p 12,T = p 1,T − p 2,T and p 2,T , k T is equal to (see Ref. [51]) In Eq. (53) for simplicity, we have omitted all color coefficients and coupling constants. We see that the largest contribution stems from the region: |p 2,T − k T | ∝ 1/r, where r is the size of the dipoles. Assuming that p 1,T and p 2,T are larger that 1/r this region leads to the contribution which is equal with q = k T − p 2,T . Eq. (57) generates |q| ∼ 1/r p i,T and shows that the Bose-Einstein enhancement is essential for |∆p 12,T | p i,T and the value of |∆p 12,T | is determined by the scale of the initial conditions for the double gluon density.
Returning to the diagram of Fig. 5-b we see that generally the integration over k T enters the integration over momentum transfer of the BFKL Pomeron: Q T = p 2,T − k T , and integration over k T that characterize the size of the dipole in the Pomeron vertices (see Fig. 5-b). The typical transverse momentum in the BFKL Pomeron vertices are about p 1,T (p 2,T ) or about that of the saturation scale, at the rapidity of the vertex. On the other hand, the typical Q T ∼ 1/r where r is the size of the largest dipole of the two interacting dipoles, which constitute the exchange of the BFKL Pomeron. For the diagrams of Fig. 5-b this largest dipole has a size of the hadron, whose double parton density we discuss.
The Green function of the BFKL Pomeron G (r, R; Q T ; Y ) is known in the mixed representation, where r and R are the sizes of two interacting dipoles, Q T denotes the momentum transferred by the Pomeron, and Y the rapidity between the two dipoles. This Green function has the following form [69,70]: where ω (ν, n) = 2ᾱ S Re ψ 1 2 and ω (ν, 0) =ᾱ S χ 1 2 + iν ; with n = 0, 1, 3 . . . and χ (γ) from Eq. (41). Each term in Eq. (59) has a very simple structure, being the typical contribution of a Regge pole exchange: the product of two vertices, which depend on the size of the dipole and Q T , and the Regge-pole propagator e ω(ν,n) Y . From Eq. (60) one can see that at large Y the main contribution comes from the term with n = 0, and in what follows we will concentrate on this particular term.
The vertices with n = 0 have been determined in Refs. [69,70], and they have a simple form in the complex number representation for the point on the two dimensional plane: viz.
For r(x, y) : ρ = x + iy; ρ * = x − iy; For Using this notation the vertices have the following structure: At Q T → 0 this vertex takes the form: Using that Returning to Eq. (59), one can see that the exchange of the BFKL Pomeron turns out to be small for R Q T > 1, where R is the size of the larger of the two interacting dipoles. In the diagram of Fig. 5-b, the size of the smallest dipole is about r ∝ 1/p 2,T , while R is the dipole in the hadron which has a size of the order 1/µ soft , µ soft denotes the soft scale. In other words, we expect that Q T of the BFKL Pomerons is rather small, Q T ≤ µ soft . Since |Q T | = |k T − p 2,T | ≤ µ soft p 2,T we safely use for vertex V ν (r, Q T ) Eq. (63) which gives for the vertex in the momentum representation: the following expression: Actually, the second term does not contribute to the scattering amplitude at small Q T (see Refs. [69,71]), and therefore the diagram of Fig. 5-b gives the following contribution: where Φ 2I P denotes the double gluon density due to exchange of two BFKL Pomerons. All other notations are shown in Fig. 5-b. The second line of the equation is written assuming that ν 1 1 and ν 2 1 at high energies, in accord with the diffusion approximation (see Eq. (41)). It is easy to see that this contribution is the Fourier transform of the emission term with ρ (2) in Eq. (16) in momentum representation. We need to add the gluon reggeization term to Eq. (68) which is the Fourier transform of the second term with ρ (2) in Eq. (16). Therefore, we do not see any other contribution except the Bose-Einstein enhancement in the first diagrams. In the last line of the equation we consider Q T p 2T and replace k T by p 2,T . From Eq. (69) one can see that the double parton density with q T = 0 (see Eq. (3)) can be obtained only if we know the double parton density for q T = 0. For q T = 0 the diagram of Fig. 5-b can be re-written in the form:

V. BFKL EVOLUTION WITH BOSE-EINSTEIN ENHANCEMENT
We need to change Eq. (38) by adding the interference diagram. To do this we have to change Eq. (68), replacing Φ 2I P (Y − y, k T , p 2,T ; Q T ) by Φ (Y − y, k T , p 2,T ; Q T ) and taking into account the complete BFKL kernel of Eq. (19). Therefore, the contribution of the interference diagram to the evolution equation takes the form where K denotes the kernel of Eq. (19). In Eq. (72) we have taken into account that q T = 0. Substituting this kernel we reduce Eq. (72) to the form As we have discussed in the previous section, the typical Q T is determined by the soft scale from the initial condition, as in Fig. 5-b, or by the saturation scale at rapidity y > y. Both are much smaller than p iT , or the saturation momentum at rapidity y. Hence, we can neglect Q T = k T − p 1T , as it is much smaller than p i,T . Finally, Eq. (72) takes the form: The second term in {. . . } stems from the gluon reggeization, in which we neglect q T in comparison with p i,T . Bearing Eq. (73) in mind, Eq. (38) takes the form: We simplify the equation by first neglecting the q T dependence of the BFKL kernel in the first two terms of the R.H.S. of the equation, since as has been discussed, q T k T (p i,T ). As the second step we go to the impact parameter representation: Eq. (74) then has the form: We need to re-write the kernel 1 ∆p 12,T − Q T − q T 2 to regularize the infrared singularity and to take into account that |Q T | p i,T . The last constraint was used in deriving the equation. We now replace the kernel by the expression One can see that the term with µ 2 regularizes the infrared divergency, and the second term guarantees that only Q T and q T less than p 1,T , contribute to the integral. Calculating the Fourier transform of Eq. (77) we obtain Since at b → 0 (b p 1,T ) → ln p 2 1,T /µ 2 , the reggeization term with ω G (p 1,t ) = ln p 2 T /µ 2 cancels the infrared We first find the solution to the homogeneous equation re-writing it in the Mellin transform of Eq. (39). It has the form: K (γ 2 , γ 1 − γ 1 ; b) = p 1,T dp 1,T p 2,T dp 2,T p 2 Integrating over γ 1 (γ 2 ) we get the following equation We solve this equation using the iteration procedure with respect to the small parameter δ, assuming that the solution without the interference term is equal to Plugging this solution in Eq. (79), we obtain the following equation for the spectrum of the homogeneous equation: In general for the BFKL kernel (see Eq. (19)) we cannot integrate the integral analytically. Instead, we use the diffusion approximation to obtain the analytical result: Searching for the solution with the new intercept ω − 2ᾱ S ω 0 = ω (1) (γ 1 , γ 2 ), we obtain the solution for γ 1 = γ 2 → 1.
Therefore, we see that the value of the intercept for the double parton density (Y − y) dependence, is larger than that for the product of single parton densities (see Eq. (1)), which is equal to Therefore, the difference ∆ 2 − 2ᾱ S ω 0 turns out to be proportional to 1/ N 2 c − 1 2/3 , which is a small number at large N c . However, for N c = 3 (see Eq. (87)) ∆ 2 − 2ᾱ S ω 0 ≈ 0.2ᾱ S ≈ 0.07ᾱ S ω 0 . This value of the correction leads to an effect of the order of 1 at Y ≈ 20.
We can calculate the corrections of the order of δ 2 to the intercept, using for the iteration Φ (1) (ω, γ 1 , γ 2 ) the form However, we expect negligible values for this correction, and we proceed to find the solution of Eq. (74) using Eq. (88), as the solution for the homogeneous equation.
One can see that this contribution is not proportional to φ (Y − y, ξ 1 ) φ (Y − y, ξ 2 ), and contradicts Eq. (1). It should be stressed that Eq. (1) violates both the Y − y dependence due to the intercept ∆ 2 > 2ᾱ S ω 0 , and the ξ 1 (ξ 2 ) dependence (compare this equation with Eq. (95)). Note that the change in the ξ shape of the distribution, has no suppression of 1/ N 2 c − 1 .

VII. CONCLUSIONS
In this paper we found that in the BFKL evolution, the Bose-Einstein enhancement leads to a faster increase of the double parton densities than the product of two single parton distributions. This effect has been discussed by us for the DGLAP evolution in the region of low x [57]. On the qualitative level, the DGLAP and BFKL evolution lead to large correlations at high energies, due to the correlations of the identical gluons. It should be noted that all 1/(N 2 c − 1) corrections in the double gluon densities stem from the Bose-Einstein enhancement. The BFKL evolution generates the power dependence on x ( Φ ∝ (1/x) ∆2 with ∆ 2 − 2ω 0 > 0, where ω 0 is the intercept of the BFKL Pomeron, this difference turns out to be numerically small, since it is proportional tō α S / N 2 c − 1 2/3 .
In particular, these correlations clarify the physical meaning of the increase of the anomalous dimension of the twist four operator that has been discussed in Refs. [9,10,[53][54][55]66]. It should be stressed that we obtain the intercept for the double gluon density is proportion; to 1/ N 2 c − 1 2/3 which is quite different from the twist four intercept, which is proportional to 1/ N 2 c − 1 2 . However, in the case of the anomalous dimensions, corrections other than the Bose-Einstein enhancement, of the order of 1/(N 2 c − 1), contribute so making the calculation of the energy behavior of the twist four operator to be a different problem.
We view this paper as the next step in understanding the role of the identical parton correlations in the parton evolution. The next project that we plan to consider, is to include the identical parton correlations in the DGLAP evolution at finite x. We believe that the most interesting question in his part of the program is to include the Pauli blocking for quarks and antiquarks (see Ref. [72]).