DGLAP evolution for DIS diffraction production of high masses

In this paper we develop the DGLAP evolution for the system of produced gluons in the process of diffractive production in DIS, directly from the evolution equation in Color Glass Condensate approach. We are able to describe the available experimental data with small value of the QCD coupling ($\bar{\alpha_S} \approx 0.1$). We conclude that in diffractive production, we have a dilute system of emitted gluons and in the order to describe them, we need to develop the next-to-leading order approach in perturbative QCD.

II. The CGC/saturation equations for DIS diffractive productions 2 III. Double log approximation for the produced gluons for τ 0 = r 2 Q 2 s (Y 0 , b) < 1 4 IV. Double log approximation for the produced gluons for τ 0 = r 2 Q 2 s (Y 0 , b) ≤ 1 7 A. qq production. V. DGLAP evolution for the emitted gluons and quarks at small β 12

VI. Conclusions 14
VII. Acknowledgements 14 Fig. 1-b). The equations, that govern the emission of gluons in this process, were proven long ago [2](see also Ref. [3][4][5]) but, they have not been investigated carefully for almost two decades. The efforts of high energy community were concentrated on simple models in which the diffractive production of quark-antiquark pair and one additional gluon has been considered (see Refs. [6][7][8][9][10][11][12][13]). The successful description of the old experimental data led to an elusive impression, that it is not necessary to investigate the dense system of produced gluons in this process. In the previous paper [14] we developed the saturation model, in which we describe the dense system of produced gluons. However, the comparison with the experimental data shows, that we failed to describe the data in spite of a number of the fitting parameters in the model. Based on this experience, we wish to study the emission of gluons in the DGLAP evolution with the hope, that the experimental data at small β (see Fig. 1) can be interpreted as the production of a rather dense system of gluons, which is not in the saturation region, but approaching it. In other words, we wish to search the DGLAP evolution for dipole sizes for which r 2 Q 2 s (Y − Y 0 , Y 0 ) < 1. On the other hand, we will show that r 2 Q 2 s (Y 0 ) ∼ 1 contribute to the elastic amplitude, which means that we have to take into account the saturation of gluons in the structure of N el in Fig. 1-b. Bearing this in mind we develop the DGLAP approach in the coordinate representation directly from the equations of Ref. [2].
The DGLAP evolution has been discussed (see for example, Ref. [15] and reference therein) but mostly, using the Ingelman-Schlein factorization [16] and reducing the evolution of the diffractive structure function to the DGLAP equation for the Pomeron structure function. In this paper we take a completely different approach based on the evolution equation for the diffractive cross section of Ref. [2], in which we do not introduce the so called soft Pomeron and its structure.
The paper is organized as follows. In the next section we give a brief review of the energy evolution of the processes of the diffractive production in DIS, which have been derived in CGC/saturation approach (see Refs. [1,2]). In section 3 we solve these equation in double log approximation (DLA) in the kinematic region where r 2 Q 2 s (Y 0 , b) ≪ 1 and show that this solution can describe the experimental data. In section 4 we continue to discuss the DLA and expand it to the region r 2 Q 2 s (Y 0 , b) ≤ 1. We put our main attention on fixing the initial conditions for the DLA equation. In section 5 we present the DGLAP evolution equation for the diffractive reduced cross sections. In the conclusions we summarize our results.

II. THE CGC/SATURATION EQUATIONS FOR DIS DIFFRACTIVE PRODUCTIONS
A sketch of the process of diffraction production in DIS is shown in Fig. 1-a, from this figure one can see that the main formula he form where Y = ln (1/x Bj ) and Y 0 is the minimum rapidity gap for the diffraction process (see Fig. 1-b). In other words, we consider diffraction production, in which all produced hadrons have rapidities larger than Y 0 . For σ dif f dipole (r ⊥ , Y, Y 0 ) we have a general expression where the structure of the amplitude N D is shown in Fig. 1-a. For N D the evolution equation has been derived in Ref. [2] in the leading log(1/x) approximation (LLA) of perturbative QCD(see Ref. [1] for detail and general description of the LLA). Hence, we hope to describe the experimental data only in the kinematic region where both β and x I P are very small (Y − Y 0 = ln(1/β) ≫ 1 and Y 0 = ln(1/x I P ) ≫ 1 are large).
The equation as has been shown in Ref. [2], can be written in two forms. First, it turns out that for the new function the equation has the same form as Balitsky-Kovchegov equation [17]: viz.
The graphic representation of the processes of diffraction production in the region of high mass. The wavy lines denote the BFKL Pomerons. Y = ln (1/xBj ) , Y0 = ln (1/xIP ) and β = Q 2 / Q 2 + M 2 X , xIP = Q 2 + M 2 X /s, x bj = Q 2 /s where Q is the photon virtuality, s the squared of energy and MX is the produced mass. withᾱ S = N c α S /π where α S is QCD coupling and N c is the number of colours.
Note, that r = x 01 , and the kernel of the equation describes the decay of a dipole to two dipoles: x 01 → x 02 + x 12 . The initial condition for Eq. (4) has the following form: Re-writing Eq. (4) as the equation for N D we obtain the second form of the set of the equations: The initial conditions are Eq. (7) accounts for the production of quark-antiquark pair integrated over its mass. A general feature, is that the amplitude with fixed rapidity gap can be calculate as follows It should be noted that the initial condition for . The appearance of δ-function is the artifact of the LLA, in which we only sum contributions with large δY . However, it has been shown [7,10,12] that more careful estimates of the qq production leads to smearing of the δ function and, in the first approximation, the contribution of this production to n D can be written as Eq. (9) shows that the production of qq pairs decreases as dσ/dM 2 X ∝ 1/M 4 x , which corresponds to the high energy behaviour of the amplitude due to qq exchange (see for example Ref. [1], section 3.2).

III.
DOUBLE LOG APPROXIMATION FOR THE PRODUCED GLUONS FOR τ0 = r 2 Q 2 s (Y0, b) < 1 In this section we consider the kinematic region in which N el (Y 0 ; x 01 , b) is in the vicinitity of the saturation scale Q s (Y 0 , b), but at x 2 01 Q 2 s (Y 0 , b) < 1. As it was found in Ref. [18,19] we have geometric scaling behaviour in this region, and the amplitude behaves as in the leading order approximation of perturbative QCD (LOA) with γ cr = 0.37. Considering Eq. (10), one can see that in this kinematic region we can, in general, neglect two terms in Eq. (6): the term which is proportional to N D 2 at Y − Y 0 ≪ Y 0 , since it is of the order of N 4 el , and the term which is proportional to N D N el ∝ N 3 el , while we have to keep all other terms. Therefore, the equation takes the form: In this equation we take into account the corrections of the order N 2 el , but neglected the terms of the order of N 3 el and N 4 el , assuming they are small. We believe that this equation will allow us to take into account the correction for N el ≈ 0.4 − 0.5.
Taking derivatives with respect to Y 0 , we re-write Eq. (11) for the amplitude n D (Y, Y 0 , x 01 , b) that has been introduced in Eq. (8). It has the form of the linear equation: The initial condition for this equation sshould be taken from Eq. (9) and it has the form: The elastic amplitude is: whereγ = 1 − γ cr and ξ ≡ ln 1/ x 2 10 Q 2 s (Y 0 , b) . All other parameters of Eq. (14) will be defined in Eq. (19) below. Taking the double Mellin transform, we obtain the solution to the equation of Eq. (12) in the following form: where φ in has to be determined from the initial condition of Eq. (24), and it has the form whereγ = 1 − γ cr . In this paper we will use the value of γ cr which comes from the leading order estimates: γ cr ≈ 0.37 which givesγ = 1 − γ cr = 0.63 ,γ = −1 + 2γ cr = −0.26.
Therefore, the solutionhas the form (see Fig. 1 for notations): with χ (γ) and κ are equal to where ψ(x) = d ln Γ(x)/dx and Γ is the Euler gamma function [20].The choice of the contour of integration over γ (see Fig. 2) is standard for the solution of the BFKL Pomeron, and correctly reproduces the calculation of the gluon emission in perturbative QCD. The contour of integration (C 1 ) is shown in Fig. 2. Recalling that x 2 01 Q 2 s (Y 0 , b) < 1 and ξ > 0, we can safely move this contour, and for large values of δY = Y − Y 0 and ξ, we can take the integral using the method of steepest decent. For α S δY ≫ ξ we evaluate the integral by this method, integrating along the contour C 2 which crosses the real axis at γ close to 1 2 . Here, we develop the double log(1/x) approach in which we replace the kernel of Eq. (19) by the χ (γ) = 1/γ. At ξ ≫ᾱ S δY , we can integrate by the method of steepest descend, but moving contour C 2 closer to y-axis in Fig. 2. For δY = 0 we can close the counter over pole γ =γ. However, for δY ∼ 1 we cannot use the same method, since at γ = 0 we have singularities in the kernel χ(γ). For developing the DLA , let us analyze the solution iterating the equation keepingᾱ S δ Y < 1. To obtain the solution as a sum of (δY ) n contributions we need to expand For each term of this series, we need to plug in our solution, and integrate over γ. This integral takes the following form for the third term in Eq. (18) for n ≥ 1: For n = 0, we only have the contribution of the second term in Eq. (21). In Eq. (20) we evaluated the integral, closing the contour over the singularities of the BFKL kernel, which is the pole at γ = 0, and over the pole γ =γ. The BFKL kernel also has poles at γ = −n, n = 1, 2, 3 . . . , but their contributions are exponentially suppressed with ξ leading to the next twists contributions. Eq. (21) can be re-written as follows The last term in Eq. (22) is the contribution at the pole γ =γ, while the first term is the sum of logs term giving the leading twist perturbative series.
Bearing this in mind we can re-write Eq. (18) in the following form The first term in Eq. (23) is the difference between two integrals with contour C 1 and C 4 , while the last term is the result of integrating over γ, with the contour C 4 . The advantage of this form for the equation, is that it satisfies the initial condition, since the first term is equal to zero at δY = 0, and the first term generates all perturbative logs with respect to the dipole sizes.
In the situation whenᾱ S δY ξ ≥ 1 whileᾱ S ≪ 1, the first term reduces to the double log approximation generating the contribution which stems from the term with ξ n−1 in Eq. (22).
Finally the double log contribution takes the form: In the leading order estimates the value of λ 1 =ᾱ S /γ + 2ᾱ S /γ cr . However, this term describes the quark-antiquark production whose δY behaviour is given by Eq. (9). Based on this equation and, having in mind that the next to leading order corrections are large, we consider λ 1 as a free parameter in the description of the experimental data. We expect λ 1 ≈ 1 from Eq. (9), but we will discuss this term below in more detail.
In appendix A we remove the assumption thatᾱ S δY ξ ≫ 1.
Using Eq. (25b) we attempted to describe the combined set of the inclusive diffractive cross sections measured by H1 and ZEUS collaboration at HERA [21]. The measured cross sections were expressed in terms of reduced cross sections, σ D(4) r , which is related to the measured ep cross section by In the paper, the table of x I P σ (1). In the experimental data from Ref. [21] the integral in t was performed in the region 0.09 ≤ t ≤ 0.55 GeV 2 , while our formulae are derived for the integration in t from 0 to ∞. Assuming that the b-dependence of the saturation scale Q s (Y 0 , b) is the same as was suggested in Ref. [22] we estimate that the experimentally measured region in t leads to factor 0.52 in x I P σ For the fit 61 experimental points were selected which satisfy the following criteria: Q 2 ≤ 26.5 GeV 2 , β ≤ 0.18 and x I P ≤ 0.025. The region of Q 2 and x I P was chosen from the description of the HERA data for inclusive DIS [22], as the specification of the region of small x I P , where the saturation model is able to fit the experimental data. The maximal value of β can be considered as outcome of the fit. For larger β, our approach cannot describe the data. For this sample we obtain the fit withᾱ S = 0.119 and λ 1 = 0.6 for the massless light quarks and for m c = 1.4 GeV , where m c is the mass of c-quark. The value of χ 2 is equal 62 leading to χ 2 /d.o.f. = 1.02. In Fig. 3 we show examples of the comparison of Eq. (25b) with the experimental data As expected λ 1 turns out to be close to 1. N el has been taken from Ref. [22], where the HERA data on F 2 , have been described in the CGC/saturation approach. It is instructive to note, that the value of parameter c, which is needed to fit the data, turns out to be about 0.3 − 0.4. We believe, we can apply our approximation for such values of c. Comparing with the description of the same data in our previous paper [14], one can see that we obtain a much better agreement.  Fig. 3-a) and versus xIP at fixed β and Q 2 ( Fig. 3-b). The data are taken from Ref. [21].
In Fig. 4 we show the average multiplicity of the emitted gluons n = √ᾱ S δY ξI 0 2 √ᾱ S δY ξ I 1 2 √ᾱ S δY ξ . One can see, that we cannot restrict ourselves by calculating only one emitted gluon, as it has been discussed in Refs. [6][7][8][9][10][11][12][13]. On the other hand, the density is not large to use the CGC/saturation approach for this system of gluons. As we have discussed we developed the DLA, assuming thatᾱ S ≪ 1,ᾱ S ξ ≪ 1,ᾱ S δY ≪ 1 butᾱ S δY ξ ∼ 1. We checked that describing the experimental data we haveᾱ S ξ ≤ 0.35 andᾱ S δY ≤ 0.6 while 0.3 <ᾱ S δY ξ < 2. Based on this estimates we believe that DLA can produce a good first approximation for understanding the structure of the produced gluons. On the other hand, these estimates show that we need to go beyond the DLA. The first corrections of this type we discuss in the appendix. Using Eq. (A12) we tried to describe the data and obtain a good fit with χ 2 = 66 for 61 experimental points. The values of parameters λ 1 = 0.608 andᾱ S = 0.149 turns out to be close to the previous fit, showing that the corrections work in the right direction increasing the value ofᾱ S . We firmly believe that the small values of the fittedᾱ S stems from the higher order corrections, which should be taken into account.
As has been mentioned we integrated in Eq. (1) only over kinematic region where τ 0 = r 2 Q 2 s (Y 0 , b) ≤ 1. However, the region of τ 0 ≥ 1 does not give a negligible contribution. We checked this, integrating the second term in Eq. (25b) over all r. We got reasonable description of the data reaching χ 2 = 100 for 61 points but the values of parametersᾱ S and λ 1 turn out to be quite different for this fit :ᾱ S = 0.128 and λ 1 = 0.95. This shows that we need to consider the kinematic region τ 0 ≤ 1. We do this in the next section. In this section we will consider Eq. (6), assuming that the elastic amplitude N el is in the saturation region with τ 0 ≈ 1. We consider that in this kinematic region N el (τ 0 = 1) is not a small value, but is of the order of 1, and, therefore, we cannot use the small sizes of N 3 el and N 4 el as we did deriving Eq. (12). We start with specification of Eq. (9) for the quark-antiquark production. This process has been discussed in Refs. [6][7][8][9][10][11][12][13] and we briefly review the results for the completeness of the presentation. A.
The total cross section of the quark-antiquark pair production is determined by N 2 el (r ⊥ , Y ; b) and it can be found from Eq. (1) with We need to re-write Eq. (1) in the following form to obtain the contribution of the qq production to the cross section with fixed produced mass M X : From Eq. (28) one can see that where m f mass of the quark; (30b) Plugging Eq. (30a) and Eq. (30c) into Eq. (29) and taking into account that Ψ qq are plane waves we obtain Initial condition for evolution of n D We need to return to Eq. (6) to find the initial condition for the emission of the gluons in n D . First, we will make the first iteration of this equation using Eq. (7) as the initial condition. Plugging in the initial condition for N D from Eq. (7) we obtain: In Eq. (33) we have two region of integrations x 01 ≫ x 02 (or x 01 ≫ x 12 ) and x 20 ≫ x 10 , which can lead to large logs and correspond to the singularities of the BFKL kernel. Let us first consider the region x 01 ≫ x 02 . In this region Eq. (33) takes the form: From Eq. (34) we see that the log term which is originated from the gluon reggeization (the first term in the RHS of Eq. (33) and Eq. (34)), cancels with the term (. . . ) 2 and the resulting integral over x 02 leads to the contribution which is proportional to N 2 el (Y 0 , x 01 , b) ∼ τ 2γ 0 (see appendix B for more detailed discussion of this contribution). We show below that this contribution is much smaller than the contribution that stem from the region of integration x 20 ≫ x 10 which generates ln x 2 10 Q 2 s (Y 0 , b) . The reason for this is that this term does not generate the log contribution (∝ ξ n ) and can be neglected in the DLA.
One can see that the the first iteration in this kinematic region is The term of Eq. (35) stems from the emission of an extra gluon from quark -antiquark pair, while Eq. (34) describes the contribution of the virtual gluon to the qq production. In the general equation (see Eq. (6) ) this term corresponds to the reggeization of gluons.
The first observation is that the integral in the LHS of Eq. (35) is converged since where z = ln x 2 02 Q 2 s (Y 0 , b) . Therefore, we can integrate in Eq. (35) from x 02 = 0. The RHS of Eq. (35) is n D (Y = Y 0 , Y 0 , x 01 ; b) (see Eq. (8)). Finally, the initial condition for n D takes the form where ε is a constant. Since the integral over τ 0 is convergent, the typical value of τ 0 cannot be found analytically and we used our model [22] to estimate it. It turns out that ε ≈ 1.2.
It should be stressed that this contribution to n D is proportional to x 2 01 Q 2 s (Y 0 , b) and it is much larger than the contribution of the order of τ 2γ 0 that stems from the region x 2 01 ≫ x 2 02 . Finally, collecting all terms for the initial condition of n D we obtain where the last term we discuss in the appendix B. The first term in Eq. (39) was firstly derived in Ref. [12], considering the extra gluon emission in perturbative QCD. Here, we derived it directly from Eq. (6) and and Eq. (8).

C. DLA
Based on experience of the previous section we do not expect the number of emitted gluons will be large. Hence, we need to solve linear evolution equation (see Eq. (12)), with the initial condition of Eq. (39). This equation in the DLA takes the form Eq. (40) can be solved by the iteration with the answer: The cross section has the following form We compare Eq. (42) with the experimental data, choosing 42 experimental points (see Ref. [21]) with β ≤ 0.056. All other kinematic variables were the same, as in section III, in this comparison. We obtain the description with χ 2 = 82 andᾱ S = 0.063. One can see that the value ofᾱ S turns out to be smaller than in the previous section. This fact can be an indication that the higher order corrections are essential, that we need to take into account the full DGLAP kernel for ξ evolution. The high order correction we leave for the future publications, and consider the full DGLAP approach in the next section.
As far as comparison with the experimental data is concerned, we consider the comparison as a good, especially since we had only one parameter. Eq. (41) can be re-written at large δY as which is the equation for the saturation scale: ξ s =ᾱ S κ δY , in the case of the DLA. However, we see that the solution is not a constant, but smoothly (logarithmically) depends on δY . Therefore, the simple equation for Q s = Q 0 exp (ᾱ S κ δY ) has to be corrected [1,19,23]. It turns out [19] that the corrected formula for energy dependence of the saturation momentum has the following form: Eq. (46) has simple meaning, that the scale of the dense system of emitted gluons is determined by the saturation scale Q s (δY ), which is equal for δY = 0 to the saturation scale of the parton system which leads to N el (x 01 , Y 0 ). This equation has been derived recently [27] 1 from different and more microscopic insight in the evolution of the emitted dipoles.

V. DGLAP EVOLUTION FOR THE EMITTED GLUONS AND QUARKS AT SMALL β
Eq. (40), which sums the diagrams in the leading log approximation, guides us in writing the DLAP evolution equations. Indeed, it shows that the physics observable, for which we can write the DLA is g D (Y, Y 0 , x 01 ) ≡ n D (Y, Y 0 ; x 01 , b) x 2 01 for which the equation can be re-written in the form Eq. (47) is written for small β but can be easily generalized for any values of β replacing dY ′ = −dβ ′ /β ′ integration by the DGLAP kernel: viz.
The general form of the DGLAP [24] equation takes the form: In Eq. (50) Σ (β, Y 0 , ξ) is the structure function of sea quarks and antiquarks in the coordinate representation. The initial condition in our approach uhas the form The splitting functions are well known and can be found in any text book (see Refs. [1,25] for example). We need to use the ω-representation, viz.
to solve Eq. (50). In this representation the equation takes the form: where and for the completeness of presentation their explicit forms in the leading order are given in the appendix C. The solution to Eq. (53) in the region of small β, has been discussed in details in Ref. [25]. The solution to the secular(characteristic) equation that corresponds to Eq. (53) has the following form: In Fig. 6 we show the dependence on ω for λ ± . One can see that only λ + is large, at the small values of ω [25]. Indeed, the expansion of λ ± [25] at ω = 0 gives (for N c = N f = 3) In ref. [25] it was noted that λ + (ω) = λ appr (ω) = 1/ω − 1 within 5% accuracy for the values of ω ≤ 3. It should be stressed that both λ + and λ approx have zeros at ω = 1 which follow from the energy conservation. The general form of the solution has the following form [25]: where operators P ± are the projectors on the eigenfunction of Eq. (53): P + + P − ≡ I where I is the unit operator. In Ref. [25] the operators P ± are found in the region of small ω (small β) and they have the form: The second equation shows the projector operator for N c = N f = 3.
We need to find the initial conditions in ω-representation. For G (ω, Y 0 , ξ) we need to re-write in ω-representation Eq. (49) which takes the form The contribution of G to Eq. (56) for n D is equal to (see Eq. (58)) 1 ω e ω δY +ᾱS λ+(ω) ξ + ρ ω e ω δY +ᾱS λ+(ω) ξ − e ω δY +ᾱS λ−(ω) ξ For λ + (ω) we use λ approx = 1/ω − 1. This substitution simplifies the first term which takes the following form: The term with λ − decreases at large ξ but it has singularities at ω = −k with k = 1, 2, . . . . Having this in mind we can find the corrections which could be essential at large β. In vicinity of ω = −1, λ − takes a form (for N c = N f = 3) λ − (ω) = − 7 9 1 ω + 1 Replacing λ − by Eq. (61) and λ + by λ approx we obtain for n D Finally, we need to replace n D Eq. (41) in Eq. (42), by n D Eq. (62) . We compare the resulting formula we compare with the 42 experimental points (see Ref. [21]) with β ≤ 0.056, having the same other kinematic variables as in sections III and IV. We obtain the description with χ 2 = 97 andᾱ S = 0.0549. From this comparison we can conclude that the improvement which is related to the full kernels of the DGLAP equation, is not essential for the description of the experimental data.

VI. CONCLUSIONS
In this paper we developed the DGLAP evolution in the region of small x I P and small β for the diffractive production in DIS in two different cases: the elastic amplitude at Y = Y 0 is small and outside of the saturation region; and the saturation effects are essential in N el .Both approachers can describe the available experimental data but the price for this description is a small value of the QCD coupling (ᾱ S ∼ 0.1). The conclusion from these attempts is that the system of produced gluons are dilute, even at small values of β as β = 0.0056.
At first sight, the origin of this dilute system of gluons, stems from the small value of the elastic amplitude at In Refs. [14,22] we found that N el (τ 0 = 1) = 0.1 − 0.3. However, we showed that from the master equations (see Eq. (4)), the initial condition for the gluons density of the emitted gluons has the form: is given by the integral over all τ (distances) in Eq. (38). This integral cannot be estimated only from N el in the vicinity of the saturation scale, and we used the saturation models of Refs. [14,22] to its estimate it. The value, which we obtainQ 2 shows that the smallness of N el at τ 0 = 1, does not induce a small initial gluon density.
Hence, the only reason for a small density of emitted gluons we is the small probabilities of their emission, which is characterized by smallᾱ S . The only reasonable explanation why we obtain a small value of the coupling, stems from the importance of the next-to-leading corrections in DGLAP and BFKL evolution, which we are planning to approach in the future publications.
We wish to calculate this integral using the special function in the most compact and economic way. First, we note that for n = 0, we can take the integral closing the contour of the integration over the pole γ =γ. For each n ≥ 1, we can write the contribution as the convolution integral Plugging (A2) into (A1) we obtain the following expression The next step is to express the integral on the r.h.s. in term of Lommel's function. Using that the expression in (A3) can be written as follows Taking integration by parts into (A5) we obtain Using u = ln(1/t)/ ln(1/x), the integral in the above expression can be rewritten as Introducing z = 2i √ᾱ S δY ξ, w = −2iγξ into (A7) yields the following representation where U ν (w, z) denotes the Lommel function of two variables. The series representation As ξ ≫ 1, the expression (A10) it is not a suitable representation for the asymptotic analysis. This problem is resolved considering the generating function of the Bessel functions that can be written as follows Introducing Eq. (A11) into Eq. (A10), and using that I m (z) = I −m (z), we obtain At largeᾱ S δY ξ, Eq. (A12) takes the form If we consider ᾱ S δY /ξ being much larger thanγ, one can see that this equation coincides with Eq. (25a) Appendix B: Initial conditions for the gluon emissions: terms that are proportional to N 2 el . In this appendix we consider Eq. (34) in more details. Using x 2 ij ,which denotes x ij 2 , and x 20 as the variable in the integration., we write where θ is the angle between x 10 and x 20 . In the integral of Eq. (34) Ω 1 is assumed to be N el Using polar coordinate for integrating over x 02 and introducing a new radial variable considering r 2 = x 2 20 /x 2 10 , we obtain Eq. (34) in the form with f 1 = 1 2π 1 0 dr 2 r 2 π 0 dθ 1 − 2r cos(θ) + r 2 (r 2 ) 2γ + (1 − 2r cos(θ) + r 2 ) 2γ − 1 + 2(r 2 ) γ (1 − 2r cos(θ) + r 2 ) γ , (B4a) f 2 = 1 2π 1 0 dr 2 r 2 π 0 dθ 1 − 2r cos(θ) + r 2 (r 2 ) 2γ (1 − 2r cos(θ) + r 2 ) γ + (r 2 ) γ (1 − 2r cos(θ) + r 2 ) 2γ , (B4b) It should be noted that all integrals in the above equations are convergent ones. Their values can be obtained through the hypergeometric functions. Specifically, in the procedure of the calculation we use the following relations: π 0 sin (2µ−1) (θ) (1 − 2a cos(θ) + a 2 ) ν = B(ν, 1/2) 2 F 1 (ν, ν − µ + (1/2); a 2 ), (a 2 < 1) (B5a) x b 2 F 1 (a, a, 1, x) = It is clear that for the integral I 2 we have the same expansion as for I 1 : viz.
where eachf i is described similarly as (B4) but for r being in the range 1 ≤ r ≤ ∞ . Changing the variable r → z with z = 1/r 2 we obtaiñ Taking integrals we obtaiñ where H 1−2γ = ψ(1) − ψ(2(1 − γ)) is harmonic number. For the numerical value for γ we obtain that Calculatingf 1 we assumed that the integral of Eq. (B15) is diverged in the region of τ 0 ∼ 1 where we can replace N el (Y 0 , x ij , b) by (Q 2 s (Y 0 , b) 2 x 2 ij ) γ . However, for estimates of the termᾱ Sf3 N 4 el (Y 0 ; x 10 , b) in Eq. (B14) we cannot use this replacement and have to calculate it using the saturation model of Ref. [22].
The last contribution, which proportional to N 2 el , stems from Eq. (37) and Eq. (38).In both these equations we integrated over x 2 02 from x 2 02 = 0 while actually this integration should be for x 2 02 ≥ x 2 01 . Hence we need to subtract the following integral For the completeness of presentation we include the well know formulae for the anomalous dimensions in the leading order of perturbative QCD (see, for examples, Refs. [1,25]).