Energy evolution of J/$\mathbf{\psi}$ production in DIS on nuclei

In this paper we show, that the $J/\psi$ production in DIS, is the main source of information about the events with two parton shower production. We attempt to develop our theoretical acumen of this process, to a level compatible with the theoretical description of inclusive DIS. We revisit the problem of the linear evolution equation for the double gluon densities, and include Bose-Einstein enhancement to these equations. We find that the Bose-Einstein correlations lead to an increase of the anomalous dimension, which turns out to be suppressed as $1/(N^2_c -1)^2$, in agreement with the estimates for the twist four anomalous dimension. We believe that understanding what happens to these contributions at ultra high energies, is a key question for an effective theory, based on high energy QCD. We derive the evolution equation for the scattering amplitude of two dipoles with a nucleus, taking into account the shadowing corrections, and investigate the analytical solutions in two distinct kinematic regions: deep in the saturation region, and in the vicinity of the saturation scale. The suggested non-linear evolution equation is a direct generalization of the Balitsky-Kovchegov equation, which has to be solved with the initial condition that depends on the saturation scale $Q_s(Y=Y_0,b)$. With the goal of finding a new small parameter, it is instructive to compare the solution of the non-linear equation with the qusi-classical approximation, in which in the initial condition we replace $Q_s(Y=Y_0,b)$ by $Q_s(Y,b)$. Our final result is that the shadowing corrections in the elastic amplitude generate the survival probability, which suppresses the growth of the amplitude with energy, caused by the Bose-Einstein enhancement.


I. INTRODUCTION
We believe that the process for the production of the J/ψ meson in the fragmentation region in DIS, is the simplest process aside from inclusive DIS, which allows us to investigate the scattering amplitudes in the wide kinematic region, from short distances which are the subject of perturbative QCD, to long distances of about 1/Q s (Q s is the saturation momentum), which can be described by the Colour Glass Condensate (CGC) approach(see Ref. [1] for a review). For distances > 1/Q s , the non-perturbative corrections have to be included. The nucleus target has additional advantages: a new parameterᾱ S A 1/3 , which results in a more rigorous theoretical approach; and the fact that the dominant mechanism is inelastic J/ψ production, accompanied by the production of two parton showers, (see Fig. 1), which shows the process at the initial energy in the Born approximation of perturbative QCD, for dipolenucleon scattering [2,3]). Therefore, we believe that this process is the main source of the experimental information about the event with two parton showers, in the same way as inclusive DIS is a source for the structure of the one parton shower event. The goal of this paper is to develop our theoretical understanding of J/ψ production, to a level compatible with the theoretical description of inclusive DIS. At high energy this description includes the BFKL evolution equation [4] and the CGC/saturation approach for the scattering amplitude [21][22][23][24][25][26]. As we have eluded, the process of J/ψ production, is intimately related to the double parton densities (double parton distribution functions). The DGLAP [5] evolution equation is known for the double parton distribution functions (see Ref. [6]), also the BFKL evolution is known for these processes [7][8][9]. In this paper we write the BFKL equation for the particular double parton distribution that enters the J/ψ production cross sections, and we build the generating functional for all multi-parton distributions. We also wish to amend these equations, to include the Bose-Einstein enhancement, resulting from the correlations of identical gluons. This effect has a suppression of the order of 1/(N 2 c − 1), where N c is the number of colours, and is closely related to the 1/(N 2 1 − 1) corrections to the anomalous dimension of the twist four operator, which is larger than the sum of the anomalous dimensions of two twist 2 operators, as has been shown in Refs. [10][11][12]. Bose-Einstein correlations have drawn considerable attention recently, since they give essential contributions to the azimuthal angle correlations [14][15][16][17][18][19][20]. It has been shown (see Ref. [17] for example), that the Bose-Einstein enhancement provides a significant contribution to the measured angle correlations. We believe that this fact calls for a generalization of the evolution equation by also including this enhancement.
Unfortunately, only in the linear approximation, is the cross section for J/ψ production, determined by the double gluon density. The total cross section of DIS, is also affected by strong shadowing corrections [21][22][23][24][25][26]. We propose a non-linear evolution for this cross section, which is a direct generalization of the Balitsky-Kovchegov evolution. We need to solve this non-linear equation, using the results of Refs. [2,3] as the initial condition for this equation. It should be stressed that all formulae in these papers were derived using the Born approximation of perturbative QCD for the scattering amplitude of a colourless dipole. In making estimates, the so called quasi-classical approach is widely used, replacing the dipole-proton cross section by σ = r 2 Q 2 s /4, where Q s denotes the saturation scale at arbitrary Y , but not only at initial Y = Y 0 , r is the size of the dipole. This is certainly not the correct procedure. For the total DIS cross section, this approximation means that we can use the McLerran-Venugopalan [23] formula with the above replacement, instead of the correct BK non-linear equation. We can view our equation as a simplification of the general approach in the framework of the CGC, given in Ref. [27], neglecting the 1/N c correction, and using the approximations related to physical nuclei.
In section VI, we show that the saturation for the elastic amplitude leads to the survival probability, which is so large that it suppresses the power-like increase with energy, which is generated by the Bose-Einstein enhancement. We summarize all our results in the Conclusions.

II. THE FIRST DIAGRAM: TWO BFKL POMERON EXCHANGES
At present, the effective theory for QCD at high energies, exists in two forms: the CGC/saturation approach [21][22][23][24][25][26] and the BFKL Pomeron calculus [4,7,21,[28][29][30][31][32][33] . It has been proven that in general, these two approaches are equivalent in a limited region of the rapidities [34]. However, for the Balitsky-Kovchegov cascade, which we discuss here for the interaction with the nucleus target, the equivalence holds in the entire kinematic region in rapidity. The interpretation of processes at high energy appear quite different in each approach, since they have different structural elements. The CGC/saturation approach , being more microscopic, describes the high energy interactions in terms of colourless dipoles, their density, distribution over impact parameters, evolution in energy and so on. This approach can be easily applied to inclusive processes at high energies, generating a new typical scale: i.e. saturation momentum Q s .
The BFKL Pomeron calculus which deals with BFKL Pomerons and their interactions, is similar to the old Reggeon theory [35], and is suitable for describing diffractive physics and correlations in multi-particle production, as we can use the Mueller diagram technique [36]. The relation between different processes at high energy are very often more transparent in this approach, since in addition to the Mueller diagram technique we can use the AGK cutting rules [37], which are useful in spite of the restricted region of their application [38].
In this paper we use the BFKL Pomeron calculus for the process of interest, since we wish to find the relationship between the J/ψ production and inclusive DIS. This is the main difference to Ref. [27], in which the same process has been treated using the CGC approach.
In the framework of the BFKL Pomeron calculus, the first diagram that describes J/ψ production, is the exchange of two BFKL Pomerons shown in Fig. 2. Since we are interested in inelastic J/ψ production, the Pomerons in Fig. 2 are cut Pomerons in which all gluons are produced. From the unitarity constraints for the elastic amplitude of the dipole of size r, rapidity Y and at the impact parameter b, we have Its contribution to the total cross section for J/ψ production is equal to and Q T = k T − k ′ T . In Eq. (3) Ψ γ * (Q, r) denotes the wave function of the dipole of size r in the photon with virtuality Q, while Ψ J/ψ (r) is the wave function of the same dipole in the J/ψ meson. The amplitude of the BFKL Pomeron N BFKL cut k T + 1 2 q T , Q T can be re-written in the case of the interaction of the Pomeron with the nucleus as follows, taking into account Eq. (1).
where N BFKL N k T + 1 2 q T , Q T = 0 is the scattering amplitude of the dipole with the nucleon.
where ρ (b, z) denotes the nucleon density in the nucleus. S A (b) is the number of the nucleons at fixed impact parameter b. The last equation in Eq. (4) follows from the fact that the typical value Q T in the nucleus form factor is about 1/R A , where R A denotes the radius of the nucleus. On the other hand Q T in N FKL N is of the order of the soft scale µ sof t , or the saturation scale Q s and, therefore, turns out to be much larger than 1/R A , and can be neglected.
Using Eq. (4) we can re-write Eq. (2) in the form where l T = k T + 1 2 q T and l ′ T = −k t + 1 2 q T . Introducing the double transverse momentum densities which characterize how many partons with (x 1 , p 1,T ) and (x 2 , p 2,T ) are in the parton cascade, and which is equal to where Ψ {x i , k i,T }, x 1 , p 1,T − 1 2 Q T ; x 2 , p 2,T + 1 2 Q T denotes the partonic wave function, we can re-write Eq. (6) in the form In Eq. (7) a + and a denote the creation and annihilation operators for gluons that have longitudinal momentum x i and transverse momentum p i,T .
For the double parton density in the coordinate representation which is equal to Eq. (8) takes the form where

III. GENERATING FUNCTIONAL AND BFKL EVOLUTION FOR THE MULTI-GLUON DENSITIES
The evolution equations for the multi-parton densities has been given in Ref. [9]. In this section we briefly review this proof, and derive the non-linear evolution equation for the generating functional with two colourless dipoles as the initial condition. This will give us the linear BFKL evolution for all partonic densities ρ (n) ({r i , b i }). Following Ref. [24] we consider a colourless dipole of size r i at impact parameter b i , as the parton in the fast projectile which can be described by the partonic wave function Ψ ({r i , b i }). For a system of n-partons with coordinates {r i , b i } with i = 1, . . . , n and rapidity Y , we can introduce the probability (probability density) to find them in the projectile wave function: P n (Y ; {r i , b i }). For P n we write the classical cascade equation [9,24,32]: The two terms of Eq.(11) have a simple meaning: the first one describes the decrease in probability to find n dipoles, due to a decay of one dipole into two dipoles of arbitrary size. This probability is equal toᾱ where ρ is an infrared cutoff. The second term, reflects the increase in probability to find n dipoles due to a creation of a new dipole from n − 1 dipoles, with probabilitȳ To find the probabilities, we need to add the initial condition at Y = Y 0 to Eq. (11). We wish to stress that we are only dealing with the decay of the dipoles, neglecting the possibilities for the dipoles to recombine. In terms of the Pomeron calculus, it means that we only take into account 'fan' Pomeron diagrams, which give the dominant contribution for DIS with the nuclear target [21,24,25]. Eq. (11) can be re-written in more elegant form introducing a generating functional Z where u(r i , b i ) ≡ u i is an arbitrary function of r i and b i . It follows immediately from Eq. (14) that the functional obeys the condition: at u i = 1 The physical meaning of (15) is that the sum over all probabilities is one. Multiplying Eq. (11) by the product n i=1 u i and integrating over all r i and b i , we obtain the following linear equation for the generating functional: with the definitions and . 6 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. Having this in mind, we can see that the n-dipole densities in the projectile are defined as follows: Using Eq. (19) we obtain the following linear evolution equations for (20) can be re-written in the following form . In the momentum representation Eq. (21) takes the form which is similar to the DGLAP evolution ( see Ref. [6] and the references therein): where K (p i , k) denotes the BFKL kernel that has the form Eq. (20) can be resolved if we know the initial conditions. Assuming that at we see that at Y = Y 0 the initial condition for the generating functional implies The linear equation Eq. (16) has a general solution of the form and using Eq. (23) we obtain the non-linear equation of Ref. [24]: As we have demonstrated for the J/ψ production, the initial conditions at and all ρ (n) = 0 for n = 2. In other words, for Z we have the initial condition Using this equation we obtain

IV. BOSE-EINSTEIN ENHANCEMENT IN THE EVOLUTION OF DOUBLE PARTON DENSITIES
Comparing Eq. (7),Eq. (8) and Eq. (8), one can see that the double gluon density has a general form given by Eq. (27) indicates that two gluons are produced from two different parton showers, and have no correlations between them, at least at large values of rapidities. However, we can see from Fig The origin of interference diagrams which violate Eq. (27), and should be included in the evolution.

A. The interference diagram in the DLA
In this section for a sake of simplicity, we will treat the interference diagrams in double log approximation(DLA) of perturbative QCD, in which we considerᾱ S Y ≥ 1,ᾱ S ln p 2 T ≥ 1, butᾱ S ≪ 1. This approximation appears more credible not for the gluon densities, but for double gluon distribution function which is defined as At large Q, the cross section of J/ψ production depends on D (2) (x, x, Q), as can be seen from Eq. (8). First, we consider the diagram of Fig. 4-a which shows the emission of an extra gluon in one parton shower. This emission leads to the DLA increase of ρ (1) 2.
we see that For For D (2) the sum of the diagrams of Fig. 4-a give where ξ =ᾱ S ln Q 2 . We wish to stress that D (2) (Y, Y, ξ; extra gluon) can be calculated only using The interference diagram of Fig. 4-b can be written in the following form: where the colour coefficient reflects the fact that only identical gluons contribute to this diagram. The calculation of the colour coefficients are shown in Fig. 4-c.
One can see that for We would like to draw the reader's attention to the integration over Q T and Q ′ T , since, as has been noticed in Refs. [10][11][12], without this integration the interference diagrams will not acquire the term proportional to ln(Q 2 ) * . * We thank Jochen Bartels for discussion on the importance of Q T integration.

B. DGLAP evolution with Bose-Einstein enhancement: anomalous dimension
The evolution equation for the double parton distribution function can be derived from Eq. (22),and in the DLA and it has the following form: Note that, the last term in Eq. (36) depends on x 1 + x 2 = 2x, but in the leading log(1/x) approximation, we do not take into account the difference betweenỸ = ln (1/(2x)) and Y = ln(1/x), on which the first two terms depend.
To include the Bose-Einstein enhancement, we need a more general evolution equation, which is shown in Fig. 5. In this equation we add to Eq. (36) the term, which is given by Eq. (35), and which describes the contributions of the diagrams of Calculating the colour coefficient (see Fig. 4-c), we assumed that the main contribution in the region of small x, comes from the colourless states in the t-channel for the gluons with the same values of x i , as is shown in Fig. 4. In Eq. (37) we use a new notation Y = 2 N c ln(1/x). e) x x x We first solve the homogeneous equation: We introduce the following Mellin transform to solve this equation: where ω Σ = 1 2 (ω 1 + ω 2 ) and ω D = 1 2 (ω 1 − ω 2 ). Plugging Eq. (39) into Eq. (38) we obtain where δ = 2/(N 2 c − 1). Looking for a solution of the form we obtain for the spectrum of the linear equation: Utilizing the smallness of δ, we see that we have an anomalous dimension, whose value is larger than at δ = 0. This fact has been noted and discussed, for the twist four anomalous dimension, in Refs. [10][11][12][39][40][41][42]. The solution to Eq. (42) at ω D = 0 has the following form: One can check that the general form of the anomalous dimension is We wish to stress, that the corrections to the anomalous dimension of the double parton distributions, are of the . It has the same dependence on N c , as the correction to the anomalous dimension of the twist four operator [10][11][12].
The numerical value forδ in Eq. (44) is about 0.016 for N c = 3, which indicates that this correction is rather small. However, these corrections generate the term which increases both with Y and ξ, and could be important in the region of small x. We will examine them in the next section.
The procedure just described, is a method to find the first corrections with respect to δ, to the value of the anomalous dimension. Replacing Eq. (41) by the following expression we obtain an equation for the next order correction to the value of the anomalous dimension. It has the form Eq. (46) leads to the next order corrections for the value of the anomalous dimension: C. Evolution with Bose-Einstein enhancement: solution The solution to the non-homogenous equation has a general form:

general solution of homogenous equation
Taking the integral over γ, by closing the contours around the poles: γ = γ G (ω 1 + ω 2 ) and γ = γ an (ω Σ , ω D ), we obtain in (ω Σ , ω D ) e γan(ωΣ,ωD)ξ + d (1) in (ω 1 + ω 2 ) δ (ω D ) e γan(ωΣ,ωD) ξ − e γG(ω1+ω2) ξ In the DLA both Y and ξ are large, and we can use the method of steepest descent to evaluate the integral. The equations for the saddle point values for ω 1 = ω SP 1 and ω 2 = ω SP 2 have the form : From Eq. (50), one can see that ω SP D = 0. Therefore Eq. (50) reduces to Integrating over ω D expanding γ an (ω Σ , ω D ) = γ an (ω Σ , ω D = 0) + ω 2 D /ω 3 Σ and neglecting the contributions that are proportional to δ everywhere, except for the anomalous dimension in the exponent, we have the solution in the form: Recall, that we have changed our notation, and Y = ln (1/x) in Eq. (52) . The principle difference of the behaviour of the double parton density without Bose-Einstein enhancement, is that we did not have a term which is proportional to the product of two single parton densities. On the other hand, the violation of Eq. (27) is small, and we reach a value larger than 10%, only at Y ξ ∼ > 4.
At the LHC energies, Y ∼ 16, and ξ ≥ 0.25 appears to be a reasonable region for measurements.

A. The equation
We have discussed the double parton densities, but the cross section for J/ψ production depends on the scattering amplitude for which the two BFKL Pomerons contribution of the diagram of Fig. 2 and the double gluon density in general, only give a contribution in the linear approach. The scattering amplitude depends crucially on the shadowing corrections (see [1] for the review). In this section we write the generalization of the BK equation [9,25] for the scattering amplitude for two dipoles with a target. Following Ref. [25] the generating functional for the scattering amplitude N is defined where γ n denotes the amplitude for simultaneous scattering of n dipoles off the target, and ρ (n) is the n-density of the fast projectile, i.e. the two dipoles. As we saw from The equation (55) can be solved using the ansatz: N (Y, [γ]) = N (γ 1 (Y ), γ 2 (Y ) . . .) and the initial condition that at Y = Y 0 : (56) Using Eq. (56) we can reduce Eq. (55) to the following form where N (1) denotes the solution to the Balitsky-Kovchegov equation. Note, that in Eq. (57) r 3 ≡ r 1 .

B. Solution deep in the saturation region
In this section, using the approach developed in Ref. [43], we find the solution to Eq. (57), deep in the saturation region. In this region both N (2) and N (1) are close to unity, and can be written as with ∆ (2) ≪ 1 and ∆ (1) ≪ 1. In this kinematic region all τ i = r 2 i Q 2 s (Y, b i ) ≫ 1. Plugging Eq. (58) into Eq. (57), and neglecting the contribution of order ∆ (2) ∆ (1) , we obtain the following equation for ∆ (2) .
First we solve the homogenous equation: Inside the saturation region [43] ω (r i ) = ln (τ i ) = ln r 2 i Q 2 s (Y, b) = z i and therefore the last term of the l.h.s. of the equation has a contribution: 1 2 (z 1 + z 2 )∆ (2) . Bearing in mind that we are searching for the solution of Eq. (60) in the form ∆ (2) (z 1 + z 2 ). To find the solution we go to the Mellin transform given by The term in {. . . } in Eq. (60) gives a contribution which is equal to 1 2 χ (γ) ∆ (2) (γ), where χ is the BFKL kernel [1,4]: where ψ(x) = d ln Γ(x)/dx and Γ is the Euler gamma function [44]. Taking into account that ln Q 2 1−γcr , where γ cr in the leading log approximation, which we use in this paper, stems from the solution of the equation: and γ cr ≈ 0.37. Finally, Eq. (60) takes the form: Solving this equation we obtain where C is the arbitrary constant. From Eq. (61) we obtain that Taking the integral over γ using the method of steepest descent, we obtain the following equations for the saddle point γ SP at large z 1 + z 2 we see that where h.eq. denotes homogenous equation.
As has been discussed in Refs. [21,45], one does not need to know the exact form of the shadowing corrections to find the behaviour of the scattering amplitude (N (1) (Y, r, b)), in the vicinity of the saturation scale. One should find the solution to the linear equation and the equation for the saturation momentum has the form: where N 0 is an arbitrary constant, which is numerically small ( say N 0 = 1/3). We can find the solution to the linear equation using the Mellin transform: where ξ i = ln r 2 i µ 2 where µ is the soft scale. For d (2) (ω, γ 1 , γ 2 ) the equation takes the form: Therefore, the amplitude N (2) (Y, ξ 1 , Y, ξ 2 ; b) is equal to Taking the integral using the method of steepest descent, we obtain the following equations for the saddle points γ SP 1 and γ SP 2 : The line, on which N (2) (Y, ξ 1 , Y, ξ 2 ; b) is constant, is given by the equation: One can see that the difference ∆γ SP = γ SP 1 − γ SP 2 turns out to be small. Indeed, In Fig. 6-a we show that d 2 χ(γ)/dγ 2 turns out to be large, as well as the value of Y . Therefore, for reasonable values of ξ 1 − ξ 1 ≤ Y , Eq. (80) leads to a large value of ∆γ SP . Bearing this in mind we see that the solution of Eq. (79) coincides with the solution of Eq. (63). Fig. 6-b gives the solution of Eq. (79) when γ SP 1 = γ SP 2 . One can see that a significant difference for ∆γ SP can be obtained for (ξ 1 − ξ 2 )/Y ≈ 10.
Generally, the solution to Eq. (77) has the following dependence on r 1 and r 2 : where γ cr 1 and γ cr 2 are solutions to Eq. (79) (see Fig. 6-b) and

D. Initial conditions
Taking the shadowing corrections into account for the scattering amplitude, we can re-write Eq. (9) for the cross section of J/ψ production in the form: From Eq. (3) and Eq. (10) we find that Eq. (83) takes the following form: To solve Eq. (57), we need to put in the initial conditions. First, we need to find the solution to the BK equation, and the initial conditions for N (1) (Y, r, b) for dipole scattering with nuclei, this is given by the McLerran-Venugopalan formula [23] where the dipole-nucleon cross section (σ dn ) in the Born approximation of perturbative QCD is equal to σ dn = 4πα 2 S r 2 ln(1/(r 2 µ 2 ) where µ is the soft scale. In the quasi-classical approach (QCA) we replace this cross section by σ dn = 1 8 r 2 Q 2 s (Y = Y 0 , b). The initial condition for . . . in Eq. (84) can be found using the approach of Refs. [2,3]. For J/ψ production at low energy, the qq pair propagate from first to the last inelastic interaction (see Fig. 1-a points z 1 and z 2 , respectively), undergoing inelastic interactions with the cross section 1 16 (r 1 − r 2 ) 2 Q 2 s (Y 0 , b). Before the first inelastic collisions and after the last one, the quark-antiquark pair has only elastic rescatterings. Plugging in the expression for the first and the last inelastic interaction, given in Refs. [2,3], we obtain In Ref. [27]( see also Refs. [46,47]) it is shown how Eq. (86) stems from the scattering amplitude of two dipoles with the nucleus target in the large N c limit.

E. Wave functions for virtual photon and J/ψ
For the completeness of presentation we need to specify what expression we use for the wave functions in Eq. (57). We adopt the wave functions from Refs. [48,49] which have the following form, Denoting Ψ γ * (Q, where e f denotes the charge of the c-quark and m c is its mass.Q 2 f = m 2 c z(1 − z) + Q 2 . L and T denote the polarizations of the photon. For φ L,T T (r, z) we suggest to use the simplest Gaussian parameterization.

VI. BOSE-EINSTEIN ENHANCEMENT AND SHADOWING
At first sight, the increase of the double parton density, which is stronger than the product of two single parton densities, is potentially harmful, and could lead to a violation of the unitarity constraints in the framework of the BFKL Pomeron calculus. On the other hand the JIMWLK equation [26] which includes all corrections of the order of 1/(N 2 c − 1) does not exhibit any problem with unitarity. The goal of this section is to show that it is sufficient that to ensure unitarity at high energies. The reason for this feature can be illustrated using Eq. (57) in which we insert N (1) = 1. Such an equation determines the asymptotic solution for N (2) . The resulting equation takes the form: The first line of this equation shows that N (2) = Const is the solution of the equation. The second line illustrates that N (2) approaches the asymptotic solution with the function ∆ (2) (Y, r 1 , (2) ). It should be stressed that the fact that N (2) approaches a constant, does not depend on the value ofᾱ S which determines the intercept of the energy (Y ) dependence of the solution to the linear equation. This screening is provided by the elastic amplitude, which is equal to unity at high energy. Hence, the elastic rescattering of patrons creates a survival probability [50,51], which suppresses the power-like increase of the double parton amplitudes. We wish to show that 1/(N 2 c − 1) corrections to the evolution equation have the same structure as Eq. (89). For the sake of simplicity, we discuss the shadowing corrections in the framework of the DLA of perturbation QCD (see Eq. (37) and Fig. 5). Fig. 7 illustrates that the emission of an extra gluon can be screened by the exchange of the BFKL Pomeron, both for normal DGLAP evolution, as well as for the interference diagram.
The exchange of the extra BFKL Pomeron for DGLAP ( Fig. 7-a) and interference ( Fig. 7-b) diagrams in DLA of perturbative QCD. Fig. 7-a and Fig. 7-b show the emission an extra gluon in one of the parton showers and its interaction with the target due to the BFKL Pomeron exchange.
The advantage of the evolution equation (see Eq. (57)), is that we can take into account the shadowing corrections in the same way for both terms of the equation for the double gluon structure function. The homogenous part of the resulting equation takes the form: N (2) (Y, Y, ξ) denotes the scattering amplitude, which in the linear approximation coincides with the double parton distribution function. However, this amplitude, as can be seen from Eq. (90) and Fig. 7, contains shadowing corrections, while the double parton distribution function cannot be affected by the shadowing, as we have seen above. This screening is provided by the elastic amplitude which is equal to unity at high energy. Hence, the elastic rescattering of patrons creates a survival probability [50,51], which suppress the power-like increase of the double parton ampli-tudesEq. (90), and has a typical form for the shadowing correction in DLA of perturbative QCD [1,21,22]. However, to see all numerical coefficients and to discuss the influence of shadowing on the asymptotic behavior of the scattering amplitude, it is desirable to re-write Eq. (90) in the coordinate representation. The amplitude N (2) (Y, Y, ξ) in the coordinate representation is intimately related to the scattering amplitude of the dipole of size x 01 (see Fig. 8) with the nucleus which interacts with the target two or more times. In the Born approximation such an amplitude is equal tõ whereξ = ln 1/r 2 . In DLA it is sufficient to change ξ toξ for transforming from momentum to coordinate representation. The energy evolution of this amplitude is given by the emission of the extra gluon (see Fig. 8-c) and replacesξ in Eq. (95) by the function ofξ. The contribution of the interference diagram is shown in Fig. 8-a. One can see that gluons with transverse momenta p 2T and −p 2T interacts with the dipole x 12 . However, gluons with momenta p 1T and −p 1T interact with the dipoles of the different sizes: x 12 and x 01 . In the DLA x 12 ≫ x 01 . The interaction with these two gluons can be written as follows: Using the result of Eq. (96), we have the following expression for the diagram of Fig. 8-a: From Eq. (94) we can re-write Eq. (90) in the coordinate representation in a more transparent form. It has the structure where for the terms without Bose-Einstein enhancement). This general equation has the following form forÑ (2) Y, Y,ξ; b .
Each term of this equation for N (1) = 1 has the same form as Eq. (89), including the term proportional toᾱ S /(N 2 −1) and, therefore, the asymptotic solution givesÑ (2) In other words, the fast increase of the double gluon densities, which is generated by the Bose-Eintstein enhancement, turned out to be damped by the screening. It happens in spite of the power -like increase of this amplitude with energy in the linear evolution, which is faster, than the increase for the elastic amplitude of dipole-target interaction.

VII. CONCLUSIONS
In this paper we discussed several facets of the energy evolution of the J/ψ inelastic production cross section in deep inelastic processes. We showed that this cross section in the linear approximation, can be written by means of the double gluon density, for which we can use the evolution equation in the BFKL kinematic region. This equation has been proposed previously in Refs. [7][8][9] (see Ref. [6] for a review of the equations in the DGLAP kinematic region). We completed this study by suggesting an equation for the generating functional which allows us to calculate multi-gluon density, that can be generated in the process of J/ψ production.
We include in the evolution equation the Bose-Einstein enhancement which occurs at all values of energy (for example Ref. [17]). We did this in the simplest case: DGLAP evolution in DLA approximation of perturbative QCD. Indeed, the particular form of the double gluon distribution function on which the cross section of J/ψ production depends: d 2 Q T D (2) x 1 , p 1,T + 1 2 Q T ; x 2 , p 2,T − 1 2 Q T , has a DLA limit. We found that the Bose-Einstein enhancement leads to the double gluon distribution function which increases faster that the product of two single gluon distribution function. We discussed the value of the anomalous dimension as well as the solution of the resulting evolution equations. It turns out that the contribution of the Bose-Einstein enhancement is rather small, and proportional toᾱ S / N 2 c − ) 2 , in accord with the estimates for twist four anomalous dimension [10][11][12]. However, we believe that understanding what happens to these contributions at ultra high energies, is a principal question for the effective theory, based on high energy QCD.
Bearing this in mind, we derived the evolution equation for the scattering amplitude of two dipoles with a nucleus, taking into account the shadowing corrections. We investigated the analytical solutions in two distinct kinematic regions: deep in the saturation region, and in the vicinity of the saturation scale. Therefore, we have prepared the ground for building saturation models to describe the J/ψ production cross section. Most data exists for the J/ψ production, not in DIS but for the proton-nucleus interaction. As has been demonstrated in Refs. [2,3,27] this process is closely related to the production of J/ψ in DIS. Hence, we need to re-write our formulae, since the distribution over transverse momentum of the produced J/ψ was measured, which has not been discussed in this paper. We plan to discuss this process elsewhere.
The suggested non-linear evolution equation is a direct generalization of the BK equation, and has to be solved using the formulae of Refs. [2,3], at the initial energy (Y = Y 0 ). The initial conditions depend on the saturation scale Q s (Y = Y 0 , b). We need to compare the solution to the non-linear equation in quasi-classical approach, in which we replace Q S (Y = Y 0 , b) by Q s (Y, b) at arbitrary values of rapidities Y . Such substitution looks attractive from the physical point of view, since it takes into account the fact that the main properties of the dense system of patrons depend on the saturation scale. However, we do not observe a small parameter that allows us to do this. Consequently, it is instructive to compare the solution of the non-linear equation with this approximation, with the goal to find a new small parameter.
Our last result is that the shadowing corrections in the elastic amplitude, generate the survival probability that suppressed the energy behavior of the double gluon densities with the Bose-Einstein enhancement. We found that the cross section reaches a constant at ultra high energies. More precisely, we found that a constant is the asymptotic solution for the non-linear equation. Nevertheless, this constant could be equal to zero. In this case the cross section decreases as exp − z1+z2) 2 8κ being the solution of Eq. (89). If quasi-classical approach reflects the main property of the solution to the non-linear equation, we can expect from Eq. (86) that, indeed, the cross section vanishes at high energies.
In the paper we used the BFKL Pomeron calculus which is not as general as the CGC approach, being equivalent to this approach only, in the limited range of rapidity. On the other hand, the BFKL Pomeron calculus has two advantages: the possibility to treat on the same footing elastic (diffractive) and inelastic processes, and to consider processes with different multiplicities.