DGLAP evolution for DIS diffraction production of high masses

Abstract

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.

1 Introduction

In this paper we continue to re-visit the process of diffractive production in the deep inelastic scattering in the framework of CGC/saturation approach (see Ref. [1] for review). In this approach the diffraction production is characterized by two new scales (two saturation momenta): $Q_s\left( Y - Y_0,Y_0\right)$, which describes the dense system of produced gluons (see Fig. 1a) and $Q_s\left( Y_0\right)$, which corresponds to the dense system of gluons that is responsible for $N_{el}$ (see Fig. 1b). The equations, that govern the emission of gluons in this process, were proven long ago [2] (see also Refs. [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,14]). 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 [15] 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 $\beta$ (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 ($x_{01}$¹ for which $x^2_{01} \,Q^2_s\left( Y - Y_0, Y_0\right) \,\,<\,\,1$. On the other hand, we will show that $x^2_{01}\,Q^2_s\left( Y_0\right) \,\sim \,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. 1b. Bearing this in mind we develop the DGLAP approach in the coordinate representation directly from the equations of Ref. [2].

Fig. 1

The graphic representation of the processes of diffraction production in the region of high mass. The wavy lines denote the BFKL Pomerons. $Y = \ln \left( 1/x_{Bj}\right) , Y_0 = \ln \left( 1/x_{I\!\!P}\right)$ and $\beta = Q^2/\left( Q^2 + M^2_X\right)$, $x_{I\!\!P}=\left( Q^2 + M^2_X\right) /s, x_{bj} = Q^2/s$ where Q is the photon virtuality, s the squared of energy and $M_X$ is the produced mass. $r \equiv x_{01}$

The DGLAP evolution has been discussed (see for example, Refs. [16,17,18,19] and reference therein) but mostly, using the Ingelman-Schlein factorization [20] 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 Sect. 3 we solve these equation in double log approximation (DLA) in the kinematic region where $x^2 _{01} Q^2_s\left( Y_0,b\right) \ll \,1$ and show that this solution can describe the experimental data. In Sect. 4 we continue to discuss the DLA and expand it to the region $x^2_{01} Q^2_s\left( Y_0,b\right) \, \le \,1$. We put our main attention on fixing the initial conditions for the DLA equation. In Sect. 5 we present the DGLAP evolution equation for the diffractive reduced cross sections. In the conclusions we summarize our results.

2 The CGC/saturation equations for DIS diffractive productions

A sketch of the process of diffraction production in DIS is shown in Fig. 1a, from this figure one can see that the main formula he form

$$\begin{aligned}&\sigma ^\mathrm{diff}\left( Y, Y_0, Q^2\right) \,\,\,\nonumber \\&\quad =\,\,\int \,\,d^2 x_{01} \int \,d z\,\, |\Psi ^{\gamma ^*}\left( Q^2; r_{\perp }, z\right) |^2 \,\,\sigma _\mathrm{dipole}^{diff}\left( x_{01}, Y, Y_0\right) ,\nonumber \\ \end{aligned}$$

(1)

where $Y = \ln \left( 1/x_{Bj}\right)$ and $Y_0$ is the minimum rapidity gap for the diffraction process (see Fig. 1b). In other words, we consider diffraction production, in which all produced hadrons have rapidities larger than $Y_0$. For $\sigma _\mathrm{dipole}^{diff} (x_{01}, Y, Y_0)$ we have a general expression

$$\begin{aligned} \sigma _\mathrm{dipole}^{diff} (x_{01},Y, Y_0) \,\,=\,\,\,\,\int \,d^2 b\,N^D(x_{01},Y, Y_0;\varvec{b}), \end{aligned}$$

(2)

where the structure of the amplitude $N^D$ is shown in Fig. 1a.

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 $\beta$ and $x_{I\!\!P}$ are very small ($Y - Y_0 \,=\,\ln (1/\beta ) \,\gg \,1$ and $Y_0 = \ln (1/x_{{I\!\!P}}) \,\gg \,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

$$\begin{aligned} \mathcal{N}\left( Y, Y_0; x_{01}, b \right) \,\,\equiv \,\,2 N_\mathrm{el}\left( Y; x_{01}, b \right) \,\,-\,\,N^D(Y, Y_0; x_{01}, b) \end{aligned}$$

(3)

the equation has the same form as Balitsky-Kovchegov equation [21,22,23]: viz.

$$\begin{aligned}&\displaystyle {\frac{\partial \mathcal{N}\left( Y, Y_0; \varvec{x}_{01}, \varvec{b} \right) }{\partial Y}}\,\nonumber \\&\quad =\,\displaystyle {\bar{\alpha }_S\int \frac{d^2 {\mathbf {x_{2}}}}{2 \pi }\frac{{\mathbf {x^2_{01}}}}{{\mathbf {x^2_{02}}}\, {\mathbf {x^2_{12}}}}\Bigg \{ \mathcal{N}\left( Y,Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12}\right) } \nonumber \\&\qquad \displaystyle {+ \, \mathcal{N}\left( Y,Y_0; \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \,-\,\mathcal{N}\left( Y,Y_0; \varvec{x}_{01},\varvec{b} \right) }\, \nonumber \\&\qquad -\,\displaystyle {\mathcal{N}\left( Y,Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12}\right) \mathcal{N}\left( Y,Y_0; \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \Bigg \}}\nonumber \\ \end{aligned}$$

(4)

with $\bar{\alpha }_S\,=\,N_c\alpha _S/\pi$ where $\alpha _S$ is QCD coupling and $N_c$ is the number of colours.

Note, that $\varvec{r} = \varvec{x}_{01}$, and the kernel of the equation describes the decay of a dipole to two dipoles: $\varvec{x}_{01}\,\rightarrow \,\varvec{x}_{02}\,+\,\varvec{x}_{12}$. The initial condition for Eq. (4) has the following form:

$$\begin{aligned}&\mathcal{N}\left( Y = Y_0, Y_0; \varvec{x}_{01}, \varvec{b} \right) \,\,=\,\,2\,N_\mathrm{el}\left( Y = Y_0; \varvec{x}_{01}, \varvec{b} \right) \,\nonumber \\&\quad -\,N^2_\mathrm{el}\left( Y = Y_0; \varvec{x}_{01}, \varvec{b} \right) \end{aligned}$$

(5)

Re-writing Eq. (4) as the equation for $N^D$ we obtain the second form of the set of the equations:

$$\begin{aligned}&\frac{ \partial N^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) }{ \partial Y} \,\nonumber \\&\quad =\bar{\alpha }_S\int \,\frac{ d^2 \varvec{x}_2}{2 \pi } \frac{{\mathbf {x^2_{01}}}}{{\mathbf {x^2_{02}}}\, {\mathbf {x^2_{12}}}}\Bigg \{\, N^D \left( Y,Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12}\right) \nonumber \\&\qquad +N^D\left( Y,Y_0; \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \,-\,N^D \Big (Y,Y_0; \varvec{x}_{01},\varvec{b} \Big )\nonumber \\&\qquad +\,\, N^D\left( Y, Y_0; \varvec{x}_{02},\varvec{ b}- \frac{1}{2}\varvec{x}_{12} \right) N^D\left( Y,Y_0; \varvec{x}_{12}, \varvec{ b}- \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&\qquad - 2 \,N^D\left( Y, Y_0; \varvec{x}_{02},\varvec{ b}- \frac{1}{2}\varvec{x}_{12} \right) N_\mathrm{el}\left( Y; \varvec{x}_{12}, \varvec{ b}- \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&\qquad - 2\,N^D\left( Y, Y_0; \varvec{x}_{12},\varvec{ b}- \frac{1}{2}\varvec{x}_{02} \right) N_\mathrm{el} \Big (Y; \varvec{x}_{02}, \varvec{ b}- \frac{1}{2}\varvec{x}_{12}\Big )\, \nonumber \\&\qquad +2\, N_\mathrm{el}\left( Y; \varvec{x}_{02},\varvec{ b}- \frac{1}{2}\varvec{x}_{12} \right) N_\mathrm{el} \left( Y; \varvec{x}_{12}, \varvec{ b}- \frac{1}{2}\varvec{x}_{02}\right) \,\Bigg \}. \end{aligned}$$

(6)

The initial conditions are

$$\begin{aligned} N^D\left( Y = Y_0 ,Y_0; \varvec{x}_{01}, \varvec{b}\right) \,\,=\,\,N^2_\mathrm{el} (Y_0; \varvec{x}_{01}, \varvec{ b}) \end{aligned}$$

(7)

Equation (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

$$\begin{aligned}&n^D\left( Y , \text{ rapidity } \text{ gap }= Y_0; \varvec{x}_{01},\varvec{b}\right) =- \frac{\partial N^D\left( Y ,Y_0; \varvec{x}_{01},\varvec{b}\right) }{\partial Y_0}\,\nonumber \\&\quad =\, \frac{\partial \mathcal{N}\left( Y ,Y_0; \varvec{x}_{01},\varvec{b}\right) }{\partial Y_0}\, \end{aligned}$$

(8)

It should be noted that the initial condition for $n^D\left( Y,\delta Y = Y - Y_0; \varvec{x}_{01},\varvec{b}\right) \,=\, \,N^2_\mathrm{el} (Y_0; \varvec{x}_{01}, \varvec{ b})\,\delta \left( Y - Y_0\right)$. The appearance of $\delta$-function is the artifact of the LLA, in which we only sum contributions with large $\delta Y$. However, it has been shown [7, 11, 13] that more careful estimates of the $q {\bar{q}}$ production leads to smearing of the $\delta$ function and, in the first approximation, the contribution of this production to $n^D$ can be written as

$$\begin{aligned} n^D\left( Y ,\delta Y = Y - Y_0; \varvec{x}_{01},\varvec{b}\right) \,\,=\,\,N^2_\mathrm{el} (Y_0; \varvec{x}_{01}, \varvec{ b})\,e^{ - \delta Y} \end{aligned}$$

(9)

Equation (9) shows that the production of $q {\bar{q}}$ pairs decreases as $d \sigma /d M^2_X\,\propto \,1/M^4_x$, which corresponds to the high energy behaviour of the amplitude due to $q {\bar{q}}$ exchange (see for example Ref. [1], section 3.2).

3 Double log approximation for the produced gluons for $\tau _0\,=\,r^2\,Q^2_s\left( Y_0, b\right) \, \ll \,1$

In this section we consider the kinematic region in which $N_\mathrm{el}\left( Y_0; \varvec{x}_{01}, \varvec{b}\right)$ is in the vicinity of the saturation scale $Q_s\left( Y_0,b \right)$, but at $x^2_{01}\,Q^2_s\left( Y_0,b\right) \, \ll \, 1$. As it was found in Refs. [24, 25] we have geometric scaling behaviour in this region, and the amplitude behaves as

$$\begin{aligned} \displaystyle {N_\mathrm{el}\left( Y_0; \varvec{x}_{01}, \varvec{b}\right) \,\,=\,\,\mathrm{c}\,\left( x^2_{10}\,Q^2_s\left( Y_0,b\right) \right) ^{1 - \gamma _{cr}} }\,\,\ll \,\,1 \end{aligned}$$

(10)

in the leading order approximation of perturbative QCD (LOA) with $\gamma _{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 $\left( N^D\right) ^2$ at $Y-Y_0 \ll Y_0$, since it is of the order of $N^4_\mathrm{el}$, and the term which is proportional to $N^D N_\mathrm{el} \propto N^3_\mathrm{el}$, while we have to keep all other terms.

Therefore, the equation takes the form:

$$\begin{aligned}&\frac{ \partial N^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) }{ \partial Y} \, \nonumber \\&\quad =\bar{\alpha }_S\int \, \frac{d^2 \varvec{x}_2}{ 2 \pi } \frac{{\mathbf {x^2_{01}}}}{{\mathbf {x^2_{02}}}\, {\mathbf {x^2_{12}}}} \Bigg \{\, N^D \left( Y,Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12}\right) \nonumber \\&\qquad +N^D\left( Y,Y_0; \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \,-\,N^D \Big (Y,Y_0; \varvec{x}_{01},\varvec{b} \Big ) \nonumber \\&\qquad +2\, N_\mathrm{el}\left( Y; \varvec{x}_{02},\varvec{ b}- \frac{1}{2}\varvec{x}_{12} \right) N_\mathrm{el}\Big (Y; \varvec{x}_{12}, \varvec{ b} \nonumber \\&\qquad -\frac{1}{2}\varvec{x}_{02}\Big )\,\Bigg \}. \end{aligned}$$

(11)

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} \approx 0.4 - 0.5$.

Taking derivatives with respect to $Y_0$, we re-write Eq. (11) for the amplitude $n^D\left( Y,Y_0,\varvec{x}_{01},\varvec{b}\right)$ that has been introduced in Eq. (8). It has the form of the linear equation:

$$\begin{aligned}&\frac{ \partial n^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) }{ \partial Y} \,\nonumber \\&\quad =\bar{\alpha }_S\int \,\frac{ d^2 \varvec{x}_2}{ 2 \pi } \frac{{\mathbf {x^2_{01}}}}{{\mathbf {x^2_{02}}}\, {\mathbf {x^2_{12}}}} \Bigg \{\, n^D \left( Y,Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12}\right) \nonumber \\&\qquad +\,n^D\left( Y,Y_0; \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \,-\,n^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b} \right) \Bigg \}\nonumber \\ \end{aligned}$$

(12)

The initial condition for this equation should be taken from Eq. (9) and it has the form:

$$\begin{aligned} \displaystyle {n^D\left( Y_0, Y_0; \varvec{x}_{01}, \varvec{b}\right) \,\,=\,\,N^2_\mathrm{el}\left( Y_0; \varvec{x}_{01}, \varvec{b}\right) }\nonumber \\ \end{aligned}$$

(13)

The elastic amplitude is:

$$\begin{aligned} N_\mathrm{el} (Y; x_{10}, b) =\mathrm{c}\left( x^2_{10}\,Q^2_s\left( Y,b\right) \right) ^{{\bar{\gamma }}} \equiv \mathrm{c}\,\,e^{{\bar{\gamma }} \left( \bar{\alpha }_S\kappa \left( Y - Y_0\right) \,-\, \xi \right) }\end{aligned}$$

(14)

where ${\bar{\gamma }} = 1 - \gamma _{cr}$ and $\xi \,\equiv \,\ln \left( 1/\left( x^2_{10}\,Q^2_s\left( Y_0, b\right) \right) \right)$. All other parameters of Eq. (14) will be defined in Eq. (19) below.

Taking the double Mellin transform,

$$\begin{aligned}&\displaystyle n^D\left( Y, Y_0, \xi , b\right) \nonumber \\&\quad = \int _{C_1}\frac{ d \gamma d \omega }{ (2 \pi i)^2}\,\phi (\omega , \gamma )\,e^{\bar{\alpha }_S\, \omega \, \left( Y- Y_0\right) \,\,+\,\,\left( \gamma - 1\right) \xi } \end{aligned}$$

(15)

we obtain the solution to the equation of Eq. (12) in the following form:

$$\begin{aligned}&\displaystyle n^D\left( Y, Y_0, \xi , b\right) \,\,\nonumber \\&\quad = \int _{C_1}\frac{ d \gamma d \omega }{ (2 \pi i)^2}\,\frac{\phi _{in}\left( \gamma ,Y_0\right) }{\omega - \chi \left( \gamma \right) }\,e^{\bar{\alpha }_S\, \omega \left( Y - Y_0\right) \,\,+\,\,\left( \gamma - 1\right) \xi } \end{aligned}$$

(16)

where $\phi _{in}$ has to be determined from the initial condition of Eq. (24), and it has the form

$$\begin{aligned} \phi _{in}\left( \gamma ,Y_0\right) \,\,=\,\,\,\frac{c^2}{\gamma - {\tilde{\gamma }}} \end{aligned}$$

(17)

In this paper we will use the value of $\gamma _{cr}$ which comes from the leading order estimates: $\gamma _{cr}\,\approx \,\,0.37$ which gives ${\bar{\gamma }}\,=\,\,1 - \gamma _{cr}=0.63$ ,${\tilde{\gamma }} = -1+ 2 \gamma _{cr} = - 0.26$.

Therefore, the solution has the form (see Fig. 1 for notations):

$$\begin{aligned}&n^D\left( Y, Y_0, \xi , b\right) \,\,\nonumber \\&\quad =\mathrm{c}^2 \int _{C_1}\frac{ d \gamma }{ 2 \pi i}\, \frac{1}{ \gamma - {\tilde{\gamma }}}\,e^{\bar{\alpha }_S\, \chi \left( \gamma \right) \left( Y - Y_0\right) \,\,+\,\,\left( \gamma - 1\right) \xi } \end{aligned}$$

(18)

with $\chi \left( \gamma \right)$ and $\kappa$ are equal to

$$\begin{aligned} \kappa \,\,=\,\,\frac{\chi \left( 1 - \gamma _{cr}\right) }{1 - \gamma _{cr}}; \quad \chi \left( \gamma \right) \,=\,2 \psi \left( 1\right) - \psi \left( \gamma \right) - \psi \left( 1 - \gamma \right) ; \end{aligned}$$

(19)

where $\psi (x) = d \ln \Gamma (x)/d x$ and $\Gamma$ is the Euler gamma function [26]. The choice of the contour of integration over $\gamma$ (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\left( Y_0,b\right) \, <\, 1$ and $\xi > 0$, we can safely move this contour, and for large values of $\delta Y = Y - Y_0$ and $\xi$, we can take the integral using the method of steepest decent. For $\alpha _S\delta Y \,\gg \,\xi$ we evaluate the integral by this method, integrating along the contour $C_2$ which crosses the real axis at $\gamma$ close to $\frac{1}{2}$. Here, we develop the double log(1/x) approach in which we replace the kernel of Eq. (19) by the $\chi \left( \gamma \right) \,=\,1/\gamma$. At $\xi \,\gg \,\bar{\alpha }_S\,\delta Y$, we can integrate by the method of steepest descend, but moving contour $C_2$ closer to y-axis in Fig. 2. For $\delta Y = 0$ we can close the counter over pole $\gamma = {\bar{\gamma }}$. However, for $\delta Y \sim 1$ we cannot use the same method, since at $\gamma = 0$ we have singularities in the kernel $\chi (\gamma )$.

Fig. 2

The contours of integration over $\gamma$

For developing the DLA , let us analyze the solution iterating the equation keeping $\bar{\alpha }_S\,\delta \, Y \,\,<\,1$. To obtain the solution as a sum of $\left( \delta Y\right) ^n$ contributions we need to expand

$$\begin{aligned}&e^{ \bar{\alpha }_S\chi \left( \gamma \right) \left( Y - Y_0\right) }\nonumber \\&\quad =\sum ^\infty _{n=0} \,\frac{\left( \bar{\alpha }_S\chi (\gamma )\,\delta Y\right) ^n}{n!} \,\,\xrightarrow {~~ \gamma \rightarrow 0~~~}\,\,\sum ^\infty _{n=0} \,\frac{1}{n!} \left( \frac{ \bar{\alpha }_S\,\delta Y}{\gamma }\right) ^n \end{aligned}$$

(20)

For each term of this series, we need to plug in our solution, and integrate over $\gamma$. This integral takes the following form for the third term in Eq. (18) for $n\,\ge \,1$:

$$\begin{aligned}&\oint _{C_3} \,d \gamma \,\left( \frac{ \bar{\alpha }_S\,\delta Y}{\gamma }\right) ^n\frac{e^{ \left( \gamma - 1\right) \xi } }{ \gamma - {\tilde{\gamma }}}\nonumber \\&\quad = \left( \bar{\alpha }_S\,\delta Y\right) ^n\Bigg \{\frac{1}{{\tilde{\gamma }}}\frac{1}{(n - 1)!} \left( \frac{e^{ \left( \gamma - 1\right) \xi } }{ \gamma - {\tilde{\gamma }}}\right) ^{(n)}_{\gamma ,\gamma \rightarrow 0}+ \frac{1}{{\tilde{\gamma }}^n}\,e^{ \left( \tilde{ \gamma } - 1\right) \xi }\Bigg \}\nonumber \\ \end{aligned}$$

(21)

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 $\gamma =0$, and over the pole $\gamma = {\tilde{\gamma }}$. The BFKL kernel also has poles at $\gamma \,=\,- n,\, n = 1,2,3 \dots$, but their contributions are exponentially suppressed with $\xi$ leading to the next twists contributions. Eq. (21) can be re-written as follows

$$\begin{aligned}&\oint _{C_3} \,d \gamma \,\left( \frac{ \bar{\alpha }_S\,\delta Y}{\gamma }\right) ^n\frac{e^{ \left( \gamma - 1\right) \xi } }{ \gamma - {\tilde{\gamma }}}\nonumber \\&\quad = \left( \frac{\bar{\alpha }_S\,\delta Y}{{\tilde{\gamma }}}\right) ^n\, e^{- \xi }\,\Bigg \{\frac{1}{{\tilde{\gamma }}}\frac{1}{(n - 1)!}\,\sum _{k=0}^{n - 1}\frac{\left( {\tilde{\gamma }} \, \xi \right) ^k}{k!} \,\, +\,\,e^{{\tilde{\gamma }} \,\xi }\Bigg \}\nonumber \\ \end{aligned}$$

(22)

The last term in Eq. (22) is the contribution at the pole $\gamma = {\tilde{\gamma }}$, 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

$$\begin{aligned}&n^D\left( Y, Y_0,\xi , b\right) \,\,\nonumber \\&\quad =\mathrm{c}^2\,\Bigg \{\int _{C_1 - C_4}\frac{d \gamma }{ 2 \pi i}\, \frac{1}{ \gamma - {\tilde{\gamma }}}\,e^{\bar{\alpha }_S\, \chi \left( \gamma \right) \left( Y - Y_0\right) + \left( \gamma - 1\right) \xi }\nonumber \\&\qquad +\, e^{ \left( \tilde{ \gamma } - 1\right) \xi } \,e^{\bar{\alpha }_S\chi \left( {\tilde{\gamma }}\right) \, \left( Y - Y_0\right) }\Bigg \} \end{aligned}$$

(23)

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 $\gamma$, 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 $\delta Y = 0$, and the first term generates all perturbative logs with respect to the dipole sizes.

In the situation when $\bar{\alpha }_S\,\delta Y\, \xi \,\ge \,1$ while $\bar{\alpha }_S\,\ll \,1$, the first term reduces to the double log approximation generating the contribution

$$\begin{aligned}&n^D_\mathrm{1-st\, term \,of\, Eq.~(23)}\left( Y,Y_0, \xi \right) \,\,\nonumber \\&\quad =\,\, \mathrm{c}^2\,\frac{\bar{\alpha }_S\delta Y}{|{\bar{\gamma }}|}\,e^{- \xi }\,\sum ^\infty _{n=1}\frac{1}{n! (n - 1)!} \left( \bar{\alpha }_S\delta Y \xi \right) ^{n - 1}\nonumber \\&\quad =\mathrm{c}^2 \,\frac{1}{|{\bar{\gamma }}|} \, \sqrt{\frac{\bar{\alpha }_S\delta Y}{\xi }}\,e^{- \xi } I_1\left( 2 \sqrt{\bar{\alpha }_S\delta Y\,\xi }\right) \end{aligned}$$

(24)

which stems from the term with $\xi ^{n - 1}$ in Eq. (22).

Finally the double log contribution takes the form:

$$\begin{aligned} n^D\left( Y,Y_0, \xi \right)= & {} \mathrm{c}^2\Bigg \{\frac{\bar{\alpha }_S}{|{\tilde{\gamma }}|} \sqrt{\frac{\bar{\alpha }_S\delta Y}{\xi }}\,e^{- \xi } I_1\left( 2 \sqrt{\bar{\alpha }_S\delta Y\,\xi }\right) \,\nonumber \\&+e^{ \left( \tilde{ \gamma } - 1\right) \xi } \,e^{- \frac{\bar{\alpha }_S}{{\tilde{\gamma }}}\, \left( Y - Y_0\right) }\Bigg \} \end{aligned}$$

(25a)

$$\begin{aligned}= & {} \mathrm{c}^2\frac{\bar{\alpha }_S}{|{\tilde{\gamma }}|} \sqrt{\frac{\bar{\alpha }_S\delta Y}{\xi }}\,e^{- \xi } I_1\left( 2 \sqrt{\bar{\alpha }_S\delta Y\,\xi }\right) \,\,\nonumber \\&+N^2_{el}\left( Y; \xi \right) \,e^{- \lambda _1\,\delta Y} \end{aligned}$$

(25b)

In the leading order estimates the value of $\lambda _1 = \bar{\alpha }_S/{\tilde{\gamma }} + 2\, \bar{\alpha }_S/\gamma _{cr}$. However, this term describes the quark-antiquark production whose $\delta 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 $\lambda _1$ as a free parameter in the description of the experimental data. We expect $\lambda _1 \approx 1$ from Eq. (9), but we will discuss this term below in more detail.

In appendix A we remove the assumption that $\bar{\alpha }_S\,\delta Y\, \xi {\gg }~1$.

Fig. 3

$\sigma ^\mathrm{diff}\left( Y, Y_0, Q^2\right) = x_{I\!\!P}\,\sigma _r$ versus $Q^2$ at fixed $\beta$ and $x_{I\!\!P}$ (a) and versus $x_{I\!\!P}$ at fixed $\beta$ and $Q^2$ (b). The data are taken from Ref. [27]

Using Eq. (25b) we attempted to describe the combined set of the inclusive diffractive cross sections measured by H1 and ZEUS collaboration at HERA [27]. The measured cross sections were expressed in terms of reduced cross sections, $\sigma _r^{D(4)}$, which is related to the measured ep cross section by

$$\begin{aligned}&\frac{\mathrm{d} \sigma ^{ep \rightarrow eXp}}{\mathrm{d} \beta \mathrm{d} Q^2 \mathrm{d} x_{I\!\!P}\mathrm{d} t} \nonumber \\&\quad = \frac{4\pi \alpha ^2}{\beta Q^4} \ \ \left[ 1-y+\frac{y^2}{2}\right] \ \ \sigma _r^{D(4)} (\beta ,Q^2,x_{I\!\!P},t) \ . \end{aligned}$$

(26)

In the paper, the table of $x_{I\!\!P}\sigma _r^{D(3)} (\beta ,Q^2,x_{I\!\!P}) = x_{I\!\!P}\int d t~ \sigma _r^{D(4)} (\beta ,Q^2,x_{I\!\!P},t)$ are presented at different values of $Q, \beta$ and $x_{I\!\!P}$. This cross section is equal to $\,\frac{Q^2}{4 \pi ^2}\,\sigma ^\mathrm{diff}\left( Y, Y_0, Q^2\right)$ where $\sigma ^\mathrm{diff}\left( Y, Y_0, Q^2\right)$ is given by Eq. (1). In the experimental data from Ref.   [27] the integral in t was performed in the region $0.09 \,\le \,t \,\le 0.55\, GeV^2$, while our formulae are derived for the integration in t from 0 to $\infty$. Assuming that the b-dependence of the saturation scale $Q_s\left( Y_0,b\right)$ is the same as was suggested in Ref. [28] we estimate that the experimentally measured region in t leads to factor 0.52 in $x_{I\!\!P}\sigma _r^{D(3)} (\beta ,Q^2,x_{I\!\!P})$.

For the fit 61 experimental points were selected which satisfy the following criteria: $Q^2\, \le \,26.5\,GeV^2$, $\beta \,\le \,0.18$ and $x_{I\!\!P}\,\le \,0.025$. The region of $Q^2$ and $x_{I\!\!P}$ was chosen from the description of the HERA data for inclusive DIS [28], 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 $\beta$ can be considered as outcome of the fit. For larger $\beta$, our approach cannot describe the data. For this sample we obtain the fit with $\bar{\alpha }_S= 0.119$ and $\lambda _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 $\chi ^2$ is equal 62 leading to $\chi ^2/\mathrm{d.o.f.}\,= 1.02$. In Fig. 3 we show examples of the comparison of Eq. (25b) with the experimental data As expected $\lambda _1$ turns out to be close to 1. $N_{el}$ has been taken from Ref. [28], 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 [15], one can see that we obtain a much better agreement.

Fig. 4

The average multiplicities of the emitted gluon at different values of $\beta$ as function of $Q^2$ using the parameters of the fit

In Fig. 4 we show the average multiplicity of the emitted gluons $n = \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }I_0\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi } \right) \big / I_1\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi } \right)$. 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,14]. 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 $\bar{\alpha }_S\,\ll \,1,\bar{\alpha }_S\,\xi \ll \,1, \,\bar{\alpha }_S\, \delta Y \,\ll \,1$ but $\bar{\alpha }_S\,\delta Y \,\xi \,\sim \,1$. We checked that describing the experimental data we have $\bar{\alpha }_S\,\xi \le 0.35$ and $\bar{\alpha }_S\delta Y \,\le \,0.6$ while $0.3 \,<\,\bar{\alpha }_S\delta Y \xi \,<\,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 $\chi ^2 = 66$ for 61 experimental points. The values of parameters $\lambda _1=0.608$ and $\bar{\alpha }_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 $\bar{\alpha }_S$. We firmly believe that the small values of the fitted $\bar{\alpha }_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 $\tau _0 \,=\,x^2_{01}\,Q^2_s\left( Y_0,b\right) \,\le \,1$. However, the region of $\tau _0 \,\ge 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 $\chi ^2 = 100$ for 61 points but the values of parameters $\bar{\alpha }_S$ and $\lambda _1$ turn out to be quite different for this fit : $\bar{\alpha }_S= 0.128$ and $\lambda _1= 0.95$. This shows that we need to consider the kinematic region $\tau _0 \,\le \,1$. We do this in the next section.

4 Double log approximation for the produced gluons for $\tau _0\,=\,r^2\,Q^2_s\left( Y_0, b\right) \,\le \,1$

In this section we will consider Eq. (6), assuming that the elastic amplitude $N_{el}$ is in the saturation region with $\tau _0 \,\approx \,1$. We consider that in this kinematic region $N_{el}\left( \tau _0=1\right)$ 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,14] and we briefly review the results for the completeness of the presentation (Fig. 5).

Fig. 5

The diffraction production of the quark-antiquark pair with mass $M_X$. $\lambda , \sigma$ and $\sigma '$ are photon, quark and antiquark helicities

4.1 $q {\bar{q}}$ production.

The total cross section of the quark-antiquark pair production is determined by $N^2_\mathrm{el}\left( r_{\perp }, Y; \varvec{b}\right)$ and it can be found from Eq. (1) with

$$\begin{aligned} \sigma _\mathrm{dipole}^{diff} (x_{01}, Y)\,\,=\,\,\int d^2 b\, N^2_\mathrm{el}\left( Y; x_{01}, \varvec{b}\right) \end{aligned}$$

(27)

We need to re-write Eq. (1) in the following form to obtain the contribution of the $q {\bar{q}}$ production to the cross section with fixed produced mass $M_X$:

$$\begin{aligned}&\sigma ^\mathrm{diff} (Y, Q^2)\nonumber \\&\quad =\int d M^2_X \,\int \frac{d^2 k_T}{(2 \pi )^2}\,\delta \left( M^2_X \,- \frac{k^2_T}{z (1 - z)}\right) \nonumber \\&\qquad \times \left( \int d^2 x_{01} \Psi ^{\gamma ^*}\left( Q; x_{01}, z\right) N_\mathrm{el}\left( Y; x_{01}, \varvec{b}\right) \Psi ^{q {\bar{q}}}\left( \varvec{k}_T; \varvec{x}_{01}\right) \right) ^2\nonumber \\ \end{aligned}$$

(28)

From Eq. (28) one can see that

$$\begin{aligned}&\left( M^2_X + Q^2\right) \frac{d \sigma ^\mathrm{diff} (Y, Y_M \,=\,\delta Y; Q^2)}{d M^2_X}\nonumber \\&\quad =x_{{I\!\!P}} \frac{d \sigma ^\mathrm{diff} (Y, Y_M \,=\,\delta Y; Q^2)}{d x_{{I\!\!P}} }\nonumber \\&\quad =\frac{d \sigma ^\mathrm{diff} (Y, Y_M \,=\,\delta Y ; Q^2)}{ d Y_0}\,\,\nonumber \\&\quad =\left( M^2_X + Q^2\right) \int ^1_0 d z \int d^2 b \int \frac{d^2 k_T}{(2 \pi )^2}\,\delta \left( M^2_X \,-\, \frac{k^2_T}{z (1 - z)}\right) \nonumber \\&\qquad \times \left( \int d^2 x_{01} \Psi ^{\gamma ^*}\left( Q;\varvec{x}_{01}, z\right) N_\mathrm{el}\left( Y; x_{01}, \varvec{b}\right) \Psi ^{q {\bar{q}}}\left( \varvec{k}_T; \varvec{x}_{01}\right) \right) ^2 \nonumber \\&\quad = \frac{Q^2}{\beta }\,\int d^2 b \int ^1_0\!\!\!\!\! z (1-z) \,d z \left( \int d^2x_{01}\Psi ^{\gamma ^*}\left( Q;\varvec{x}_{01}, z\right) N_\mathrm{el}\right. \nonumber \\&\qquad \times \left( Y; x_{01} \varvec{b}\right) \Psi ^{q {\bar{q}}}\left( \varvec{k}_T; \varvec{x}_{01}\right) \bigg )^2\bigg |_{k^2_T = M^2_X\,z \,(1-z)} \end{aligned}$$

(29)

The wave functions of photon are known (see Ref. [1] Eqs. 4.18 and 4.20) and have the following form

$$\begin{aligned}&\Psi _T\left( Q; r_\perp ,z\right) \, \nonumber \\&\quad =\,\frac{e Z_f}{2 \pi }\Bigg \{ (1 - \delta _{\sigma , \sigma })\,\left( 1 - 2\,z\,-\,\sigma \,\lambda \right) \, i a_f \frac{\varvec{\epsilon }^\lambda _\perp \cdot \varvec{r}_{\perp }}{ r_{\perp }} K_1\left( r_\perp \,a_f\right) \,\nonumber \\&\qquad +\,\delta _{\sigma \sigma '}\frac{m_f}{\sqrt{2}}\left( 1 + \sigma \lambda \right) K_0\left( r_\perp \,a_f\right) \Bigg \}; \end{aligned}$$

(30a)

$$\begin{aligned}&N_T =\,\,Q^2 z (1-z)\,\,+\,\,m^2_f ~~~~\text{ where }~~m_f ~~ \text{ mass } \text{ of } \text{ the } \text{ quark }; \nonumber \\\end{aligned}$$

(30b)

$$\begin{aligned}&\Psi _L\left( Q; r_\perp ,z\right) \, =\,\frac{e Z_f}{2 \pi } 2 Q z (1-z) \left( 1 - \delta _{\sigma , \sigma '}\right) K_0\left( r_\perp \,a_f\right) \nonumber \\ \end{aligned}$$

(30c)

Plugging Eq. (30a) and Eq. (30c) into Eq. (29) and taking into account that $\Psi ^{q{\bar{q}}}$ are plane waves we obtain

$$\begin{aligned}&\frac{d \sigma ^\mathrm{diff}_T\left( Y, Y_M \,=\,\delta Y ; Q^2\right) }{ d Y_0}\nonumber \\&\quad =\,\frac{ 2 N_c\,\alpha _\mathrm{e.m.} }{ \pi ^2}\sum Z^2_f \frac{Q^2}{\beta }\int d^2 b \int ^1_z z(1-z) d z\left( N_T F^2_T\,\right. \nonumber \\&\qquad \left. +\,m^2_f \,F^2_L\right) ; \end{aligned}$$

(31a)

$$\begin{aligned}&\frac{d \sigma ^\mathrm{diff}_L\left( Y, Y_M \,=\,\delta Y ; Q^2\right) }{ d Y_0}\nonumber \\&\quad =\,\frac{ 2 N_c\,\alpha _\mathrm{e.m.} }{ \pi ^2}\sum Z^2_f\frac{4 Q^4}{\beta }\int d^2 b \int ^1_z \left( z(1\right. \nonumber \\&\qquad \left. -z)\right) ^3 d z\,F^2_L; \end{aligned}$$

(31b)

where $N_T = z + (1-z)^2$ and

$$\begin{aligned} F_T \,= & {} \,\int \frac{r_\perp d r_\perp }{2} \,a_f\,K_1\left( a_f\,r_{\perp }\right) \,J_1\nonumber \\&\times \left( Q \sqrt{z (1-z))} \sqrt{\frac{1 - \beta }{\beta }}\,r_\perp \right) \,N_\mathrm{el}\left( Y; r_\perp , b\right) ; \nonumber \\\end{aligned}$$

(32a)

$$\begin{aligned} F_L \,= & {} \,\int \frac{r_\perp d r_\perp }{2} \,\,K_0\left( a_f\,r_{\perp }\right) \,J_0\nonumber \\&\times \left( Q \sqrt{z (1-z)) }\sqrt{\frac{1 - \beta }{\beta }}\,r_\perp \right) \,N_\mathrm{el}\left( Y; r_\perp ,\, b\right) \nonumber \\ \end{aligned}$$

(32b)

4.2 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:

$$\begin{aligned}&\frac{ \partial }{ \partial \,\delta Y} N^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) \nonumber \\&\quad =\frac{\bar{\alpha }_S}{ 2 \,\pi }\, x^2_{10}\, \int \, \frac{d^2 \varvec{x}_{02}}{\varvec{x}^2_{02}\,\varvec{x}^2_{12}} \Bigg \{ -\, N^2_\mathrm{el} \left( Y_0; \varvec{x}_{01},\varvec{b} \right) \nonumber \\&\qquad +\Bigg (\, N_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12}\right) +N_\mathrm{el} \left( Y_0; \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&\qquad - N_{\mathrm{el}} \Bigg (Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{12} \Bigg ) \, N_{\mathrm{el}} \Bigg ( Y_0; \varvec{x}_{12},\varvec{b} \nonumber \\&\qquad - \frac{1}{2}\varvec{x}_{02}\Bigg )\Bigg )^2 \Bigg \} \end{aligned}$$

(33)

In Eq. (33) we have two region of integrations $x_{01}\,\gg x_{02}$ (or $x_{01}\,\gg x_{12}$ ) and $x_{20}\,\gg \,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}\,\gg x_{02}$. In this region Eq. (33) takes the form:

$$\begin{aligned}&\frac{ \partial }{ \partial \,\delta Y} N^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) \,\,\equiv \,\,I_1\nonumber \\&\quad =\frac{1}{2}\,\bar{\alpha }_S\, \int ^{x^2_{01}}\, \frac{d \varvec{x}^2_{02}}{\varvec{x}^2_{02}} \Bigg \{ -\, N^2_\mathrm{el} \left( Y_0; \varvec{x}_{01},\varvec{b} \right) \nonumber \\&\qquad + \Bigg (\, N_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{01}\right) \nonumber \\&\qquad +\,N_\mathrm{el} \left( Y_0; \varvec{x}_{01},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&\qquad -\, N_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{01}\right) \, N_\mathrm{el} \left( Y_0; \varvec{x}_{01},\varvec{b} \right. \nonumber \\&\qquad \left. - \frac{1}{2}\varvec{x}_{02}\right) \Bigg )^2 \Bigg \} \end{aligned}$$

(34)

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 $(\dots )^2$ and the resulting integral over $x_{02}$ leads to the contribution which is proportional to $N^2_\mathrm{el}\left( Y_0, \varvec{x}_{01}, \varvec{b} \right) \,\sim \,\tau _0^{2{\bar{\gamma }}}$ (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}\,\gg \,x_{10}$ which generates $\ln \left( x^2_{10}\,Q^2_s\left( Y_0,b\right) \right)$. The reason for this is that this term does not generate the log contribution ($\propto \xi ^n$) and can be neglected in the DLA.

One can see that the the first iteration in this kinematic region is

$$\begin{aligned}&\frac{ \partial }{ \partial \,\delta Y} \frac{N^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) }{ x^2_{10}}\nonumber \\&\quad =\,\frac{\bar{\alpha }_S}{2\,\pi } \int _{x_{10}} \, \frac{d^2 \varvec{x}_{02}}{\varvec{x}^4_{20}} \Bigg (\,2\, N_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&\qquad - N^2_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \Bigg )^2 \end{aligned}$$

(35)

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 $q {\bar{q}}$ 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

$$\begin{aligned} N_\mathrm{el} \left( Y; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) =\left\{ \begin{array}{l} \propto \,\,\left( x_{02}^2 \,Q^2_s\left( Y_0, b\right) \right) ^{1 - \gamma _{cr}}\quad \text{ for } \quad x_{02}^2 \,Q^2_s\left( Y_0, b\right) \,\le \,1\,\,\text{ with }\,\,\gamma _{cr} = 0.37;\\ \\ 1 - \mathrm{Const} \exp \left( - z^2/(2 \kappa )\right) \,\,\,\,\,\text{ for }\,\,\, x_{02}^2 \,Q^2_s\left( Y_0, b\right) \,\ge \,1; \end{array} \right. \end{aligned}$$

(36)

where $z \,\,=\,\,\ln \left( x_{02}^2 \,Q^2_s\left( Y_0, b\right) \right)$.

Therefore, we can integrate in Eq. (35) from $x_{02}=0$. The RHS of Eq. (35) is $n^D\left( Y=Y_0,Y_0,x_{01};b\right)$ (see Eq. (8)). Finally, the initial condition for $n^D$ takes the form

$$\begin{aligned}&n^D\left( Y=Y_0,Y_0,x_{01}; b\right) \nonumber \\&\quad =\frac{\bar{\alpha }_S}{2\,\pi } x^2_{01}\,\int _0 \, \frac{d^2 \varvec{x}_{02}}{\varvec{x}^4_{20}} \Bigg (\,2\, N_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&\qquad - N^2_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \Bigg )^2\nonumber \\&\quad = \frac{\bar{\alpha }_S}{2\,\pi }\,\,x^2_{01}\,{\bar{Q}}^2\left( Y_0,b\right) \end{aligned}$$

(37)

with

$$\begin{aligned} {\bar{Q}}^2\left( Y_0,b\right) \,\,= & {} \,\,\int \frac{d^2 \varvec{x}_{02}}{\varvec{x}^4_{20}} \Bigg (\,2\, N_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \nonumber \\&- N^2_\mathrm{el} \left( Y_0; \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \Bigg )^2\,\,\nonumber \\= & {} Q^2_s\left( Y_0,b\right) \int \frac{ \pi d \tau _0}{\tau ^2_0}\Bigg (\,2\, N_\mathrm{el} \left( \tau _0\right) \,-\, N^2_\mathrm{el} \left( \tau _0 \right) \Bigg )^2\nonumber \\= & {} \,\,\varepsilon \,\,Q^2_s\left( Y_0,b\right) \, \end{aligned}$$

(38)

where $\varepsilon$ is a constant. Since the integral over $\tau _0$ is convergent, the typical value of $\tau _0$ cannot be found analytically and we used our model [28] to estimate it. It turns out that $\varepsilon \,\approx \,\,1.2$.

It should be stressed that this contribution to $n^D$ is proportional to $x^2_{01}\,Q^2_s\left( Y_0,b\right)$ and it is much larger than the contribution of the order of $\tau _0^{2 {\bar{\gamma }}}$ that stems from the region $x^2_{01} \,\gg \,x^2_{02}$.

Finally, collecting all terms for the initial condition of $n^D$ we obtain

$$\begin{aligned} n^D\left( Y=Y_0,Y_0; x_{01}, b\right)= & {} \frac{\bar{\alpha }_S}{2\,\pi }\,\,x^2_{01}\,{\bar{Q}}^2\left( Y_0,b\right) \nonumber \\&-5.62\,\,\bar{\alpha }_SN_\mathrm{el}^{2}(Y_{0};x_{10},b)\end{aligned}$$

(39)

where the last term we discuss in the appendix B.

The first term in Eq. (39) was firstly derived in Ref. [13], considering the extra gluon emission in perturbative QCD. Here, we derived it directly from Eq. (6) and and Eq. (8).

4.3 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

$$\begin{aligned}&\frac{ \partial }{ \partial \,\delta Y} \Bigg (\frac{n^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) }{x^2_{01}}\Bigg )\nonumber \\&\quad =\bar{\alpha }_S\, \int _{x^2_{01}}\!\!\!\frac{d \varvec{x}^2_{02}}{\varvec{x}^2_{02}} \,\Bigg ( \frac{n^D\left( Y,Y_0; \varvec{x}_{02},\varvec{b}\right) }{x^2_{02}} \Bigg )\end{aligned}$$

(40)

Equation (40) can be solved by the iteration with the answer:

$$\begin{aligned} n^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right)= & {} \frac{\bar{\alpha }_S}{2\,\pi }\,\,x^2_{01}\,{\bar{Q}}^2\left( Y_0,b\right) \,\,I_0\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }\right) \,\,\nonumber \\= & {} \,\,\varepsilon \,\frac{\bar{\alpha }_S}{2\,\pi }\,\,x^2_{01}\,Q_s^2\left( Y_0,b\right) \,I_0\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }\right) \nonumber \\ \end{aligned}$$

(41)

The cross section has the following form

$$\begin{aligned}&\sigma ^\mathrm{diff} (Y, Y_0, Q^2)\nonumber \\&\quad =\,\,\int \,\,d^2 r_{\perp } \int \,d z\,\, |\Psi ^{\gamma ^*} (Q^2; r_{\perp }, z)|^2\,\nonumber \\&\qquad \times \int d^2 b\,n^D_{Eq.~(41)}\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) \,\,\nonumber \\&\qquad +\left( 1 - \,5.62\,\bar{\alpha }_S\right) \left( \frac{d \sigma ^\mathrm{diff}_T\left( Y, Y_M \,=\,\delta Y ; Q^2\right) }{ d Y_0}\bigg |_{Eq.~(31a)}\,\right. \nonumber \\&\qquad \left. +\,\,\,\frac{d \sigma ^\mathrm{diff}_L\left( Y, Y_M \,=\,\delta Y ; Q^2\right) }{ d Y_0}\bigg |_{Eq.~(31b)}\right) \end{aligned}$$

(42)

We compare Eq. (42) with the experimental data, choosing 42 experimental points (see Ref. [27]) with $\beta \,\le \,0.056$. All other kinematic variables were the same, as in Sect. 3, in this comparison. We obtain the description with $\chi ^2 = 82$ and $\bar{\alpha }_S= 0.063$. One can see that the value of $\bar{\alpha }_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 $\xi$ 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.

4.4 $n^D$ in the vicinity of the saturation scale $Q_s\left( Y - Y_0, b\right)$

Equation (41) can be re-written at large $\delta Y$ as

$$\begin{aligned}&n^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) \nonumber \\&\quad =\varepsilon \,\frac{\bar{\alpha }_S}{2\,\pi }\,\,x^2_{01}\,Q_s^2\left( Y_0,b\right) \,\,I_0\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }\right) \nonumber \\&\qquad \times \xrightarrow {\bar{\alpha }_S\,\delta Y\,\xi \,\gg \,1}\, \varepsilon \,\frac{\bar{\alpha }_S}{2\,\pi } \frac{1}{\sqrt{4 \pi \sqrt{\bar{\alpha }_S\delta Y\,\xi }}}\exp \nonumber \\&\qquad \times \left( 2 \sqrt{\bar{\alpha }_S\delta Y \,\xi } \,-\,\xi \right) \end{aligned}$$

(43)

From equation $2 \sqrt{\bar{\alpha }_S\delta Y \,\xi } \,-\,\xi \,\,=\,0$ we find

$$\begin{aligned} \xi _s \,\,=\,\,\ln \left( \frac{1}{x^2_s\,Q_s\left( Y_0\right) }\right) \,\,=\,\,4\,\bar{\alpha }_S\,\delta Y \end{aligned}$$

(44)

which is the equation for the saturation scale: $\xi _s = \bar{\alpha }_S\kappa \,\delta Y$, in the case of the DLA. However, we see that the solution is not a constant, but smoothly (logarithmically) depends on $\delta Y$. Therefore, the simple equation for $Q_s = Q_0 \exp \left( \bar{\alpha }_S\kappa \,\delta Y\right)$ has to be corrected [1, 25, 29]. It turns out [25] that the corrected formula for energy dependence of the saturation momentum has the following form:

$$\begin{aligned} Q^2_s\left( \delta Y\right) \,=\,Q^2_0 \left( \frac{1}{\bar{\alpha }_S\delta Y}\right) ^{\frac{3}{ 2\,{\bar{\gamma }}}}\,e^{\bar{\alpha }_S\,\kappa \,\delta Y} \end{aligned}$$

(45)

From Eq. (44) one can see that $Q_0^2 = Q^2_s\left( Y_0\right)$. Eq. (41) in the vicinity of the saturation scale, where $x^2_{01} Q^2_s\left( \delta Y\right) \sim 1$, gives

$$\begin{aligned} n^D= & {} \mathrm{Const}\left( x^2_{01}\,Q^2_s\left( \delta Y\right) \right) ^{{\bar{\gamma }}}\nonumber \\= & {} \mathrm{Const}\left( x^2_{01} Q^2_s\left( Y = 0\right) \right) ^{{\bar{\gamma }}}\,\frac{e^{ \bar{\alpha }_S\,\kappa \,{\bar{\gamma }} Y}}{\left( \bar{\alpha }_SY_0\,\bar{\alpha }_S\delta Y\right) ^{3/2}} \end{aligned}$$

(46)

Equation (46) has simple meaning, that the scale of the dense system of emitted gluons is determined by the saturation scale $Q_s\left( \delta Y\right)$, which is equal for $\delta Y = 0$ to the saturation scale of the parton system which leads to $N_{el}\left( x_{01},Y_0\right)$. This equation has been derived recently [35]² from different and more microscopic insight in the evolution of the emitted dipoles.

5 DGLAP evolution for the emitted gluons and quarks at small $\beta$

Equation (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\left( Y,Y_0,x_{01}\right) \,\,\equiv \,\,n^D\left( Y,Y_0; \varvec{x}_{01},\varvec{b}\right) \big /x^2_{01}$ for which the equation can be re-written in the form

$$\begin{aligned} \frac{\partial g^D\left( Y,Y_0,x_{01}\right) }{\partial \xi }= & {} \bar{\alpha }_S\int ^Y_{Y_0} d Y' \, g^D\left( Y',Y_0,x_{01}\right) \nonumber \\= & {} \bar{\alpha }_S\int ^{\delta Y}_{0} d\, \delta Y' \, g^D\left( \delta Y',Y_0,x_{01}\right) \ \end{aligned}$$

(47)

Equation (47) is written for small $\beta$ but can be easily generalized for any values of $\beta$ replacing $d Y' = - d \beta '/\beta '$ integration by the DGLAP kernel: viz.

$$\begin{aligned}&\frac{\partial \, \beta G^D\left( \beta ,Y_0,\xi \right) }{\partial \xi }\nonumber \\&\quad =\bar{\alpha }_S\int ^1_\beta d \beta ' P_{gg}\left( \beta '\right) \,\left( \frac{\beta }{\beta '}\, G^D\left( \frac{\beta }{\beta '}, \xi \right) \right) \end{aligned}$$

(48)

where $\beta \,G\left( \beta ,Y_0,\xi \right) \,\equiv \,n^D\left( \beta , Y_0,\xi \right) \,e^{ - \xi }$, which is the gluon structure function for diffractively produced gluons in coordinate representation.

$$\begin{aligned} \beta \,G^D\left( \beta ,Y_0,\xi = 0 \right) \,\,=\,\,\frac{\bar{\alpha }_S}{2 \pi }\varepsilon \end{aligned}$$

(49)

The general form of the DGLAP [30,31,32] equation takes the form:

$$\begin{aligned}&\frac{\partial }{\partial \,\xi } \left( \begin{array}{c} \beta \,G^D\left( \beta , Y_0,\xi \right) \nonumber \\ \beta \,\Sigma ^D\left( \beta , Y_0,\xi \right) \end{array} \right) \\&\quad =\bar{\alpha }_S\,\int ^1_\beta d \beta '\, \left( \begin{array}{cc} P_{gg}\left( \beta '\right) &{} P_{gf}\left( \beta '\right) \nonumber \\ P_{fg}\left( \beta '\right) &{} P_{ff}\left( \beta '\right) \end{array} \right) \nonumber \\&\qquad \times \left( \begin{array}{c} \frac{\beta }{\beta '}\,G^D\left( \frac{\beta }{\beta '}, Y_0,\xi \right) \\ \frac{\beta }{\beta '} \,\Sigma ^D\left( \frac{\beta }{\beta '}, Y_0,\xi \right) \end{array} \right) \end{aligned}$$

(50)

In Eq. (50) $\Sigma \left( \beta ,Y_0,\xi \right)$ is the structure function of sea quarks and antiquarks in the coordinate representation. The initial condition in our approach has the form

$$\begin{aligned} \Sigma ^D\left( \beta , Y_0, \xi = 0\right) \,\,=\,\,0 \end{aligned}$$

(51)

The splitting functions are well known and can be found in any text book (see Refs. [1, 33] for example).

We need to use the $\omega$-representation, viz.

$$\begin{aligned} f\left( \omega , \xi \right)= & {} \int ^1_0 d \beta \beta ^\omega \,f\left( \beta , \xi \right) \nonumber \\= & {} \int ^\infty _0\,d Y \,e^{\omega \,Y}\left( \beta f\left( \beta , \xi \right) \right) ;\quad f\left( \beta , \xi \right) \nonumber \\= & {} \int ^{\epsilon + i \infty }_{\epsilon - i \infty }\frac{d \omega }{2 \,\pi \,i} \beta ^{- \omega - 1 }\,f\left( \omega , \xi \right) , \end{aligned}$$

(52)

to solve Eq. (50). In this representation the equation takes the form:

$$\begin{aligned}&\frac{\partial }{\partial \,\xi } \left( \begin{array}{c} \,G^D\left( \omega , Y_0,\xi \right) \\ \,\Sigma ^D\left( \omega , Y_0,\xi \right) \end{array} \right) \nonumber \\&\quad =\,\,\bar{\alpha }_S\, \left( \begin{array}{cc} \gamma _{gg}\left( \omega \right) &{} \gamma _{gf}\left( \omega \right) \\ \gamma _{fg}\left( \omega \right) &{} \gamma _{ff}\left( \omega \right) \end{array} \right) \left( \begin{array}{c} \,G^D\left( \omega , Y_0,\xi \right) \\ \,\Sigma ^D\left( \omega , Y_0,\xi \right) \end{array} \right) \end{aligned}$$

(53)

where

$$\begin{aligned} \gamma _{i,j}\,\,=\,\,\int ^1_0 d z z^\omega \,P_{i,j}\left( z \right) \end{aligned}$$

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 $\beta$, has been discussed in details in Ref. [33]. The solution to the secular(characteristic) equation that corresponds to Eq. (53) has the following form:

$$\begin{aligned} \lambda _{\pm }= & {} \frac{1}{2}\Bigg \{ \gamma _{ff}\left( \omega \right) \,+\,\gamma _{gg}\left( \omega \right) \nonumber \\&\pm \sqrt{\left( \gamma _{ff}\left( \omega \right) \,-\,\gamma _{gg}\left( \omega \right) \right) ^2 \,\,+\,\,4\,\gamma _{fg}\left( \omega \right) \,\gamma _{gf}\left( \omega \right) }\Bigg \}\nonumber \\ \end{aligned}$$

(54)

In Fig. 6 we show the dependence on $\omega$ for $\lambda _{\pm }$. One can see that only $\lambda _{+}$ is large, at the small values of $\omega$ [33]. Indeed, the expansion of $\lambda _{\pm }$ [33] at $\omega = 0$ gives (for $N_c = N_f =3$)

$$\begin{aligned} \lambda _{+}\left( \omega \right) \,\,=\,\,\frac{1}{\omega }\,\,-\,\,\frac{101}{108};\quad \lambda _{-}\,\,=\,\,- \frac{4}{27}; \end{aligned}$$

(55)

In Ref. [33] it was noted that $\lambda _{+}\left( \omega \right) \,\,= \,\,\lambda _\mathrm{appr}\left( \omega \right) \,\,=\,\,1/\omega -1$ within 5% accuracy for the values of $\omega \le 3$. It should be stressed that both $\lambda _{+}$ and $\lambda _\mathrm{approx}$ have zeros at $\omega =1$ which follow from the energy conservation.

Fig. 6

$\lambda _{\pm }$ versus $\omega$ together with $\lambda _\mathrm{approx}\left( \omega \right) \,\,=\,\,\frac{1}{\omega } - 1$

The general form of the solution has the following form [33]:

$$\begin{aligned} \Bigg \{\mathbf {P}_{+}\,e^{\bar{\alpha }_S\, \lambda _{+}\left( \omega \right) \,\xi }\,\,+\,\,\mathbf {P}_{-}\,e^{ \bar{\alpha }_S\,\lambda _{-}\left( \omega \right) \,\xi }\Bigg \} \left( \begin{array}{c} \,G\left( \omega , Y_0,\xi = 0\right) \\ \,\Sigma \left( \omega , Y_0,\xi = 0\right) \end{array} \right) \end{aligned}$$

(56)

where operators $\mathbf {P}_{\pm }$ are the projectors on the eigenfunction of Eq. (53): $\mathbf {P}_{+}\,+\, \mathbf {P}_{-}\,\,\equiv \,\,\mathbf {I}$ where $\mathbf {I}$ is the unit operator. In Ref. [33] the operators $\mathbf {P}_{\pm }$ are found in the region of small $\omega$ (small $\beta$) and they have the form:

$$\begin{aligned} \mathbf {P}_{+}= & {} \left( \begin{array}{cc} 1 &{} \frac{C_F}{C_A} \\ 0 &{} 0 \end{array} \right) \,\,+\,\,\omega \, \left( \begin{array}{cc} \rho &{} -\,C_F \\ C_F &{} C_A \end{array} \right) \,\,=\,\, \left( \begin{array}{cc} 1 &{} 4/9 \\ 0 &{} 0 \end{array} \right) \nonumber \\&+\omega \, \left( \begin{array}{cc} 85/54 &{} -\,4/3 \\ 4/3 &{} 3 \end{array} \right) \end{aligned}$$

(57)

where $C_A = N_c$ , $\rho \,=\frac{C_F\,C_A}{4\,T_R\,N_f} + \frac{C_F}{2} - \frac{2 C^2_F}{C_A}$ and $T_R = \frac{1}{2}$. The second equation shows the projector operator for $N_c = N_f = 3$.

We need to find the initial conditions in $\omega$-representation. For $G\left( \omega , Y_0,\xi \right)$ we need to re-write in $\omega$-representation Eq. (49) which takes the form

$$\begin{aligned} G^D\left( \omega , Y_0,\xi \right) \,\,=\,\,\frac{\bar{\alpha }_S\,\varepsilon }{ 2 \,\pi }\frac{1}{ \omega } \end{aligned}$$

(58)

The contribution of G to Eq. (56) for $n^D$ is equal to (see Eq. (58))

$$\begin{aligned} n^D_{G}= & {} e^{-\xi }\int ^{\epsilon + i \infty }_{\epsilon - i \infty }\frac{ d \,\omega }{2\,\pi \,i} \,\frac{\bar{\alpha }_S\,\varepsilon }{ 2 \,\pi }\frac{1}{ \omega }\,\Bigg \{e^{ \omega \,\delta Y\,\,\,+\,\,\bar{\alpha }_S\, \lambda _{+} (\omega )\,\xi }\,\,\nonumber \\&+\,\,\rho \,\omega \,\left( e^{ \omega \,\delta Y\,\,\,+\,\,\bar{\alpha }_S\, \lambda _{+} (\omega )\,\xi }\,\,-\,\, e^{ \omega \,\delta Y\,\,\,+\,\,\bar{\alpha }_S\, \lambda _{-} (\omega )\,\xi } \right) \Bigg \}\nonumber \\ \end{aligned}$$

(59)

For $\lambda _{+}\left( \omega \right)$ we use $\lambda _\mathrm{approx} \,=\,1/\omega - 1$. This substitution simplifies the first term which takes the following form:

$$\begin{aligned} n^D_1\,=\,\,\frac{\bar{\alpha }_S}{2 \pi } \left( x^2_{01}\,{\bar{Q}}^2\left( Y_0,b\right) \right) \,e^{- \bar{\alpha }_S\,\xi }\,I_0\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }\right) \end{aligned}$$

(60)

The term with $\lambda _{-}$ decreases at large $\xi$ but it has singularities at $\omega = - k$ with $k=1,2, \dots$. Having this in mind we can find the corrections which could be essential at large $\beta$. In vicinity of $\omega = -1$, $\lambda _{-}$ takes a form (for $N_c=N_f = 3$)

$$\begin{aligned} \lambda _{-}\left( \omega \right) \,\,=\,\,-\, \frac{7}{9 (\omega + 1)} \end{aligned}$$

(61)

Replacing $\lambda _{-}$ by Eq. (61) and $\lambda _{+}$ by $\lambda _\mathrm{approx}$ we obtain for $n^D$

$$\begin{aligned} n^D\,= & {} \,\,\frac{\bar{\alpha }_S}{2 \pi } \left( x^2_{01}\,{\bar{Q}}^2\left( Y_0,b\right) \right) \Bigg \{e^{- \bar{\alpha }_S\,\xi }\left( \,I_0\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }\right) \right. \nonumber \\&+\left. \rho \sqrt{\frac{\bar{\alpha }_S\,\xi }{\delta Y}}\,I_1\left( 2 \sqrt{\bar{\alpha }_S\,\delta Y\,\xi }\right) \right) \,\,\nonumber \\&+\,\,\rho e^{- \delta Y} \sqrt{\frac{7 \,\bar{\alpha }_S\,\xi }{9\,\delta Y}}\,J_1\left( 2\sqrt{ \frac{7}{9} \bar{\alpha }_S\,\delta Y\,\xi }\right) \Bigg \} \end{aligned}$$

(62)

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. [27]) with $\beta \,\le \,0.056$, having the same other kinematic variables as in Sects. 3 and 4. We obtain the description with $\chi ^2 = 97$ and $\bar{\alpha }_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.  

6 Conclusions

In this paper we developed the DGLAP evolution in the region of small $x_{I\!\!P}$ and small $\beta$ 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_\mathrm{el}$. Both approaches can describe the available experimental data but the price for this description is a small value of the QCD coupling ($\bar{\alpha }_S\sim 0.1$). The conclusion from these attempts is that the system of produced gluons are dilute, even at small values of $\beta$ as $\beta = 0.0056$.

At first sight, the origin of this dilute system of gluons, stems from the small value of the elastic amplitude at $\tau _0 = r^2 Q^2_s\left( Y_0,b\right) \,=\,1$. In Refs. [15, 28] we found that $N_\mathrm{el}\left( \tau _0 =1\right) = 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: $\frac{1}{2}\bar{\alpha }_S\,x^2_{01} {\bar{Q}}^2_s\left( Y_0,b\right)$³, where ${\bar{Q}}^2_s\left( Y_0,b\right)$ is given by the integral over all $\tau$ (distances) in Eq. (38). This integral cannot be estimated only from $N_\mathrm{el}$ in the vicinity of the saturation scale, and we used the saturation models of Refs. [15, 28] to its estimate it. The value, which we obtain ${\bar{Q}}^2_s\left( Y_0,b\right) = 1.2 \,Q^2_S\left( Y_0,b\right)$, shows that the smallness of $N_\mathrm{el}$ at $\tau _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 $\bar{\alpha }_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.

Appendices

Appendix A: Solution for $\bar{\alpha }_S\delta Y\,\xi \sim 1$

In this appendix we obtain the cross section for the diffractive production in the kinematic region, where $r^{2}Q^{2}_{s} (Y_{0},b)\approx 1$ and $\delta Y\ll 1$, but do not use the assumption, that $\bar{\alpha }_S\xi \gg 1$, which we used in section IIC-3. As in this section we are dealing with the integral

$$\begin{aligned} \displaystyle {\sum _{n=0}^{\infty }\frac{(\bar{\alpha }_S\delta Y)^{n} }{n!} \int _{C_{3}} \frac{d\gamma }{2\pi i} \dfrac{1}{\gamma ^{n}}\dfrac{e^{(\gamma -1)\xi }}{\gamma -\tilde{\gamma }}. } \end{aligned}$$

(A1)

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 $\gamma = {\tilde{\gamma }}$. For each $n\ge 1$, we can write the contribution as the convolution integral

$$\begin{aligned}&\displaystyle { \int _{C_{3}} \frac{d\gamma }{2\pi i} \dfrac{1}{\gamma ^{n}}\dfrac{e^{(\gamma -1)\xi }}{\gamma -\tilde{\gamma }}}\nonumber \\&\quad = \displaystyle {e^{({\tilde{\gamma }}-1)\,\xi } \int _{x}^{1}\dfrac{1}{\Gamma (n)} (-\ln (t))^{n-1} t^{\tilde{\gamma }} \dfrac{dt}{t} }, \end{aligned}$$

(A2)

with $x=e^{-\xi }$.

Plugging (A2) into (A1) we obtain the following expression

$$\begin{aligned}&n^{D} (Y,Y_{0}; \xi ,b)\nonumber \\&\quad = \,\,\,e^{(\tilde{\gamma }-1)\xi }\left( 1\right. \nonumber \\&\qquad \left. +\displaystyle { \int _{x}^{1}\dfrac{\bar{\alpha }_S\delta Y}{\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)}}I_{1} (2\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)}) t^{\tilde{\gamma }}\dfrac{dt}{t}}\right) .\nonumber \\ \end{aligned}$$

(A3)

The next step is to express the integral on the r.h.s. in term of Lommel’s function. Using that

$$\begin{aligned}&\dfrac{d}{dt}I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)})\nonumber \\&\quad =-\dfrac{\bar{\alpha }_S\delta Y}{\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)}}I_{1} (2\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)})\dfrac{1}{t}. \end{aligned}$$

(A4)

the expression in (A3) can be written as follows

$$\begin{aligned}&n^{D} (Y,Y_{0}; \xi , b)\,\nonumber \\&\quad =\,\,e^{(\tilde{\gamma }-1)\xi }\left( 1- \displaystyle { \int _{x}^{1} \dfrac{d}{dt}I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)}) t^{\tilde{\gamma }} dt }\right) .\quad \end{aligned}$$

(A5)

Taking integration by parts into (A5) we obtain

$$\begin{aligned}&n^{D} (Y,Y_{0};\xi , b)\nonumber \\&\quad = e^{-\xi } I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\xi })\nonumber \\&\qquad +\tilde{\gamma }e^{(\tilde{\gamma }-1)\xi } \int _{x}^{1} I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\ln (1/t)})\, t^{\tilde{\gamma }} \dfrac{dt}{t}. \end{aligned}$$

(A6)

Using $u=\sqrt{\ln (1/t)/\ln (1/x)}$, the integral in the above expression can be rewritten as

$$\begin{aligned}&n^{D} (Y,Y_{0}; \xi ,b)\nonumber \\&\quad = e^{-\xi } I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\xi })\nonumber \\&\qquad +2\tilde{\gamma }\xi e^{-\xi } \int _{0}^{1} e^{\tilde{\gamma }\xi (1- u^{2})}I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\,\xi }\, u) u du . \end{aligned}$$

(A7)

Introducing $z=2i\sqrt{\bar{\alpha }_S\delta Y\xi }$, $w=-2i\tilde{\gamma }\xi$ into (A7) yields the following representation

$$\begin{aligned} n^{D} (Y,Y_{0}; \xi ,b)= & {} \,\,\left( e^{-\xi } I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\xi })\,-\, e^{-\xi }\dfrac{1}{i}\big ( U_{1} (w,z)\right. \nonumber \\&\left. +i\,U_{2} (w,z) \right) \Bigg ), \end{aligned}$$

(A8)

where $U_{\nu } (w,z)$ denotes the Lommel function of two variables. The series representation

$$\begin{aligned} \begin{array}{rcl} U_{1} (w, z)&{}=&{}i\displaystyle { \sum _{m=0}^{\infty } \left( - \dfrac{\tilde{\gamma }\xi }{\sqrt{\bar{\alpha }_S\delta Y\xi }} \right) ^{2m+1}I_{2m+1} (2\sqrt{\bar{\alpha }_S\delta Y\xi }) }\\ U_{2} (w, z)&{}=&{}-\displaystyle { \sum _{m=0}^{\infty } \left( - \dfrac{\tilde{\gamma }\xi }{\sqrt{\bar{\alpha }_S\delta Y\xi }} \right) ^{2m+2}I_{2m+2} (2\sqrt{\bar{\alpha }_S\delta Y\xi }) } \end{array} \end{aligned}$$

(A9)

becomes

$$\begin{aligned}&n^{D} (Y,Y_{0}; \xi ,b)\nonumber \\&\quad =\left( e^{-\xi } I_{0} (2\sqrt{\bar{\alpha }_S\delta Y\xi })\right. \nonumber \\&\qquad \left. +\, e^{-\xi }\displaystyle {\sum _{m=1}^{\infty }\left( \dfrac{\tilde{\gamma }\xi }{\sqrt{\bar{\alpha }_S\delta Y\xi }} \right) ^{m} I_{m} (2\sqrt{\bar{\alpha }_S\delta Y\xi }) }\right) . \end{aligned}$$

(A10)

As $\xi \gg 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

$$\begin{aligned} e^{(z/2)\left( t+t^{-1} \right) }=\displaystyle { \sum _{m=-\infty }^{\infty }t^{m}I_{m} (z) } \end{aligned}$$

(A11)

Introducing Eq. (A11) into Eq. (A10), and using that $I_{m} (z)=I_{-m} (z)$, we obtain

$$\begin{aligned}&n^{D} (Y,Y_{0}; \xi ,b)\nonumber \\&\quad =\left( e^{(\tilde{\gamma }-1) \xi +\dfrac{\bar{\alpha }_S\delta Y}{\tilde{\gamma }}}\,\,\right. \nonumber \\&\qquad \left. -\,\,\displaystyle {e^{-\xi } \sum _{m=1}^{\infty }\left( \dfrac{1}{\tilde{\gamma }}\sqrt{\dfrac{\bar{\alpha }_S\delta Y}{\xi }} \right) ^{m} I_{m} (2\sqrt{\bar{\alpha }_S\delta Y\xi }).}\right) \end{aligned}$$

(A12)

At large $\bar{\alpha }_S\delta Y\xi$, Eq. (A12) takes the form

$$\begin{aligned}&n^{D} (Y,Y_{0}; \xi ,b)\nonumber \\&\quad =\,\,\left( e^{(\tilde{\gamma }-1) \xi +\dfrac{\bar{\alpha }_S\delta Y}{\tilde{\gamma }}}\,\right. \nonumber \\&\qquad \left. -\,\frac{1}{ 2\sqrt{\pi }}\displaystyle { \left( \dfrac{(\bar{\alpha }_S\delta Y)}{\xi ^3} \right) ^{1/4}\,e^{-\xi +2\sqrt{\bar{\alpha }_S\delta Y\xi }}}\right) \end{aligned}$$

(A13)

If we consider $\sqrt{\bar{\alpha }_S\delta Y/\xi }$ being much larger than $\tilde{\gamma }$, 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 $\Vert \varvec{x}_{ij} \Vert ^{2}$, and $\varvec{x}_{20}$ as the variable in the integration., we write

$$\begin{aligned} x_{12}^{2}= & {} (x_{10}-x_{20})\cdot (x_{10}-x_{20})\nonumber \\= & {} x^{2}_{10}-2 \Vert x_{10}\Vert \Vert x_{20}\Vert \cos (\theta )+x^{2}_{20}, \end{aligned}$$

(B1)

where $\theta$ is the angle between $\varvec{x}_{10}$ and $\varvec{x}_{20}$. In the integral of Eq. (34) $\Omega _{1}$ is assumed to be $N_{el} (Y_{0}; x_{ij},b)\sim (Q^{2}_{s}(Y_{0},b)^{2}x^{2}_{ij})^{\overline{\gamma }}$ since we are calculating this integral in the kinematic region $x^2_{01}\,Q^2_s\left( Y_0,b\right) \,<\,1$. From this form of $N_{el}$ we obtain

$$\begin{aligned} \dfrac{N_{el}(Y_{0};x_{20},b)}{N_{el}(Y_{0};x_{10},b)}&= \left( x^{2}_{20}/x^{2}_{10}\right) ^{\overline{\gamma }}, \end{aligned}$$

(B2a)

$$\begin{aligned} \dfrac{N_{el}(Y_{0};x_{12},b)}{N_{el}(Y_{0};x_{10},b)}&= \left( 1-2\dfrac{\Vert x_{20} \Vert }{\Vert x_{10}\Vert } \cos (\theta ) +\dfrac{x^{2}_{20}}{x^{2}_{10}} \right) ^{\overline{\gamma }}. \end{aligned}$$

(B2b)

Using polar coordinate for integrating over $\varvec{x}_{02}$ and introducing a new radial variable considering $r^{2}=x^{2}_{20}/x^{2}_{10}$, we obtain Eq. (34) in the form

$$\begin{aligned} I_{1}= & {} \alpha _Sf_{1}N_\mathrm{el}^{2}(Y_{0};x_{10},b)+\alpha _Sf_{2}N_\mathrm{el}^{3}(Y_{0};x_{10},b)\nonumber \\&+\alpha _Sf_{3}N_\mathrm{el}^{4}(Y_{0};x_{10},b), \end{aligned}$$

(B3)

with

$$\begin{aligned} f_{1}&=\dfrac{1}{2\pi } \int _{0}^{1}\dfrac{dr^{2}}{r^{2}} \int _{0}^{\pi }\dfrac{d\theta }{1-2r\cos (\theta )+r^{2}}\Bigg [ (r^{2})^{2\overline{\gamma }} \nonumber \\&\quad \left. +\left( 1-2r\cos (\theta ) +r^{2}\right) ^{2\overline{\gamma }}-1\right. \nonumber \\&\quad +2(r^{2})^{\overline{\gamma }} (1-2r\cos (\theta ) +r^{2})^{\overline{\gamma }} \Bigg ] , \end{aligned}$$

(B4a)

$$\begin{aligned} f_{2}&=\dfrac{1}{2\pi } \int _{0}^{1}\dfrac{dr^{2}}{r^{2}} \int _{0}^{\pi }\dfrac{d\theta }{1-2r\cos (\theta )+r^{2}}\Bigg [(r^{2}) ^{2\overline{\gamma }}\left( 1 \right. \nonumber \\&\quad \left. \left. -2r\cos (\theta ) +r^{2}\right) ^{\overline{\gamma }}+(r^{2})^{\overline{\gamma }}\left( 1\right. \right. \nonumber \\&\quad \left. -2r\cos (\theta )+r^{2} \right) ^{2\overline{\gamma }} \Bigg ], \end{aligned}$$

(B4b)

$$\begin{aligned} f_{3}&=\dfrac{1}{2\pi } \int _{0}^{1}\dfrac{dr^{2}}{r^{2}} \int _{0}^{\pi }\dfrac{d\theta }{1-2r\cos (\theta )+r^{2}}\Bigg [ (r^{2})^{2\overline{\gamma }} (1 \nonumber \\&\quad -2r\cos (\theta ) + r^{2})^{2\overline{\gamma }} \Bigg ], \end{aligned}$$

(B4c)

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:

$$\begin{aligned}&\displaystyle {\int _{0}^{\pi }\dfrac{\sin ^{(2\mu -1)} (\theta )}{(1-2a\cos (\theta ) +a^{2})^{\nu } }}\nonumber \\&\quad =B(\nu ,1/2)\ _{2}F_{1} (\nu ,\nu -\mu +(1/2);a^{2}),\quad (a^{2}<1) \end{aligned}$$

(B5a)

$$\begin{aligned}&\displaystyle { \int \dfrac{x^{b}}{1-x} dx }=\dfrac{x^{b+1}}{b+1}\ _{2}F_{1} (1,b+1,b+2;x), \end{aligned}$$

(B5b)

$$\begin{aligned}&\displaystyle { \int x^{b} _{2}F_{1} (a,a,1,x) }\nonumber \\&\quad =\dfrac{x^{b+1}}{b+1}\ _{3}F_{2} (\{a,a,b+1\},\{1,b+2\};x), \end{aligned}$$

(B5c)

$$\begin{aligned}&\displaystyle { _{2}F_{1} (1,1,1,x) }=\dfrac{1}{1-x}, \end{aligned}$$

(B5d)

$$\begin{aligned}&\displaystyle { \lim _{x\rightarrow 1^{-}} \dfrac{_{2}F_{1} (a,b,c,x)}{-\ln (1-x)} }=\displaystyle { \dfrac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)} }\quad (c=a+b). \end{aligned}$$

(B5e)

Here B denotes the Beta function whereas $_{p}F_{q}$ correspond to the (generalized) hypergeometric function represented by the following series:

$$\begin{aligned}&_{p}F_{q}\left( \{p_{1},\dots ,p_{m}\},\{q_{1},\dots ,q_{n}\};x\right) \nonumber \\&\quad =\displaystyle {\sum _{k=0}^{\infty }\dfrac{\prod _{j=1}^{m} (p_{j})_{k}}{\prod _{j=1}^{n} (q_{j})_{k} }\dfrac{x^{k}}{k!} }\quad (\vert x\vert <1), \end{aligned}$$

(B6)

with $(a)_{k}=\Gamma (a+k)/\Gamma (a)$ (see Ref. [26] for more details for integral representations of (B5) and for the special functions, respectively). It should be stressed that the equation of (B5e) shows that we do not face any logarithmic infrared divergency in our procedure of the estimates. Considering $\mu =1/2$ in (B5a) (and knowing $B(1/2,1/2)=\pi$), the closed expression for angular integrals in (B4a) is obtained. Finally, $f_1$ takes the form

$$\begin{aligned} f_{1}= & {} \dfrac{1}{2}\int _{0}^{1}dz\left[ \dfrac{z^{2\overline{\gamma }-1}-1}{1-z}+\dfrac{1}{z} \left( _{2}F_{1}\left( 1-2\overline{\gamma },1\right. \right. \right. \nonumber \\&\left. \left. \left. -2\overline{\gamma },1;z\right) -1\right) +2 z^{\overline{\gamma }-1} _{2}F_{1}\left( 1-\overline{\gamma },1\right. \right. \nonumber \\&\left. -\overline{\gamma },1;z\right) \Bigg ],\ \text {(with }z=r^{2}) \end{aligned}$$

(B7)

where

$$\begin{aligned}&\displaystyle {\int _{0}^{1}dz\dfrac{z^{2\overline{\gamma }-1}-1}{1-z}}= \psi (1)-\psi (2\overline{\gamma }),\nonumber \\&\displaystyle {\int _{0}^{1}dz \dfrac{1}{z}\ \left( _{2}F_{1}\left( 1-2\overline{\gamma },1-2\overline{\gamma },1;z\right) -1\right) }\nonumber \\&\quad =\displaystyle { \sum _{k=1}^{\infty } \dfrac{\left( \left( 1-2\overline{\gamma }\right) _{k}\right) ^{2}}{(k!)^{2}k} }, \nonumber \\&\displaystyle {\int _{0}^{1}dz \ _{2}F_{1}\left( 1-\overline{\gamma },1-\overline{\gamma },1;z\right) z^{\overline{\gamma }-1}}\nonumber \\&\quad =\displaystyle { \sum _{k=0}^{\infty } \dfrac{\left( \left( 1-\overline{\gamma }\right) _{k}\right) ^{2}}{(k!)^{2} (\overline{\gamma }+k)}}, \end{aligned}$$

(B8)

and therefore

$$\begin{aligned} \begin{array}{rcl} f_{1}&{}=&{}\dfrac{ \psi (1)-\psi (2\overline{\gamma })}{2}+ \dfrac{1}{2}\sum _{k=1}^{\infty } \dfrac{((1-2\overline{\gamma })_{k})^{2}}{(k!)^{2}k}\\ &{}&{}+ \sum _{k=0}^{\infty } \dfrac{((1-\overline{\gamma })_{k})^{2}}{(k!)^{2} (\overline{\gamma }+k)} , \\ &{}=&{}1.59 \end{array} \end{aligned}$$

(B9)

Similarly, using (B5) the terms $f_{2}$ and $f_{3}$ can be reduced to the following forms:

$$\begin{aligned} \begin{array}{rcl} f_{2}&{}=&{}-\int _{0}^{1}dz\left[ z^{2\overline{\gamma }-1}\ _{2}F_{1}\left( 1-\overline{\gamma },1-\overline{\gamma },1;z\right) \right. \\ &{}&{}\left. +z^{\overline{\gamma }-1}\ _{2}F_{1}\left( 1-2\overline{\gamma },1-2\overline{\gamma },1;z\right) \right] \\ &{}=&{}-\displaystyle { \sum _{k=0}^{\infty }\dfrac{\left( \left( 1-\overline{\gamma }\right) _{k}\right) ^{2}}{(k!)^{2} (2\overline{\gamma }+k)} - \sum _{k=0}^{\infty }\dfrac{\left( \left( 1-2\overline{\gamma }\right) _{k}\right) ^{2}}{(k!)^{2} (\overline{\gamma }+k)}},\\ &{}=&{}-\,2.53.\\ f_{3}&{}=&{}\displaystyle {\dfrac{1}{2}\int _{0}^{1}dz \ _{2}F_{1}\left( 1-2\overline{\gamma },1-2\overline{\gamma },1;z\right) z^{2\overline{\gamma }-1} }\\ &{}=&{}\displaystyle { \dfrac{1}{2} \sum _{k=0}^{\infty }\dfrac{\left( \left( 1-2\overline{\gamma }\right) _{k}\right) ^{2}}{(k!)^{2} (2\overline{\gamma }+k)}}\\ &{}=&{}0.413876 \end{array} \end{aligned}$$

(B10)

Hence, the integral $I_1$ has the value

$$\begin{aligned} I_{1}= & {} 1.59\,\bar{\alpha }_S\,N_\mathrm{el}^{2}(Y_{0};x_{10},b)\,- \,2.53\bar{\alpha }_S\, N_\mathrm{el}^{3}(Y_{0};x_{10},b)\,\nonumber \\&+\,0.414\,\bar{\alpha }_S\,N_\mathrm{el}^{4}(Y_{0};x_{10},b), \end{aligned}$$

(B11)

and since $\tau _{0}\lessapprox 1$, we can safely assume that the value of $I_1$ is given by

$$\begin{aligned} I_{1}=1.59\,\bar{\alpha }_S\,N_\mathrm{el}^{2}(Y_{0};x_{10},b)\end{aligned}$$

(B12)

We now discuss the contribution of the following integral

$$\begin{aligned} I_2 \,\,= & {} \,\,\,\frac{\bar{\alpha }_S}{ 2 \,\pi }\, x^2_{10}\, \int _{x^2_{01}}\, \frac{d^2 \varvec{x}_{02}}{\varvec{x}^2_{02}\,\varvec{x}^2_{12}}\Bigg [ \Bigg \{ -\, N^2_\mathrm{el} \left( Y_0, \varvec{x}_{01},\varvec{b} \right) \nonumber \\&+\Bigg (\, N_\mathrm{el} \left( Y_0, \varvec{x}_{02},\varvec{b}\right) \,+\,N_\mathrm{el} \left( Y_0, \varvec{x}_{12},\varvec{b} \right) \nonumber \\&-\, N_\mathrm{el} \left( Y_0, \varvec{x}_{02},\varvec{b} \right) \, N_\mathrm{el} \left( Y_0, \varvec{x}_{12},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \Bigg )^2 \Bigg \}\nonumber \\&-\Bigg \{ 2\, N_\mathrm{el} \left( Y_0, \varvec{x}_{02},\varvec{b} \right) \,-\, N^2_\mathrm{el} \left( Y_0, \varvec{x}_{02},\varvec{b} \right) \Bigg \}^2\Bigg ] \end{aligned}$$

(B13)

It is clear that for the integral $I_2$ we have the same expansion as for $I_1$: viz.

$$\begin{aligned} I_{2}\,\,= & {} \,\,\bar{\alpha }_S\, {\tilde{f}}_{1}N_\mathrm{el}^{2}(Y_{0};x_{10},b)\,\nonumber \\&+\,\bar{\alpha }_S\,{\tilde{f}}_{2}N_\mathrm{el}^{3}(Y_{0};x_{10},b)\,+\,\bar{\alpha }_S\, {\tilde{f}}_{3}N_\mathrm{el}^{4}(Y_{0};x_{10},b),\nonumber \\ \end{aligned}$$

(B14)

where each ${\tilde{f}}_{i}$ is described similarly as (B4) but for r being in the range $1\,\le \,r\,\le \,\,\infty$ . Changing the variable $r \rightarrow z$ with $z\,=\,1/r^2$ we obtain

$$\begin{aligned} {\tilde{f}}_{1}= & {} \dfrac{1}{2}\int _{0}^{1}dz\left[ \dfrac{z^{1-2\overline{\gamma }}-1}{1-z} +z^{-2\overline{\gamma }}\left( {}_{2}F_{1}\left( 1-2\overline{\gamma },1\right. \right. \right. \nonumber \\&\left. \left. \left. -2\overline{\gamma },1,z\right) \,-\,1\right) + 2 z^{-2\overline{\gamma }}\left( {} _{2}F_{1}\left( 1-\overline{\gamma },1\right. \right. \right. \nonumber \\&\left. \left. -\overline{\gamma },1,z\right) \,-\,1\right) \Bigg ], \end{aligned}$$

(B15)

Taking integrals we obtain

$$\begin{aligned} {\tilde{f}}_{1}= & {} -\,\frac{1}{2}\,H_{1-2 \overline{\gamma }} \nonumber \\&+\frac{{}_3F_2\left( 1-2 \overline{\gamma },1-2 \overline{\gamma },1-2 \overline{\gamma };1,2-2 \overline{\gamma };1\right) -1}{2\left( 1-2 \overline{\gamma }\right) } \,\nonumber \\&+\frac{{} _3F_2\left( 1-2 \overline{\gamma },1-\overline{\gamma },1-\overline{\gamma };1,2-2 \overline{\gamma };1\right) -1}{1-2 \overline{\gamma }}\nonumber \\ \end{aligned}$$

(B16)

where $H_{1 - 2 \overline{\gamma }}\,\,=\,\,\psi (1) - \psi ( 2(1 - \overline{\gamma }))$ is harmonic number.

For the numerical value for $\overline{\gamma }$ we obtain that

$$\begin{aligned} {\tilde{f}}_{1}\,\,=\,\,0.448\,\bar{\alpha }_S\,\,\,N_\mathrm{el}^{2}(Y_{0};x_{10},b)\end{aligned}$$

(B17)

Calculating ${\tilde{f}}_{1}$ we assumed that the integral of Eq. (B15) is diverged in the region of $\tau _0 \sim 1$ where we can replace $N_{el} (Y_{0},x_{ij},b)$ by $(Q^{2}_{s}(Y_{0},b)^{2}x^{2}_{ij})^{\overline{\gamma }}$ . However, for estimates of the term $\,\bar{\alpha }_S\, {\tilde{f}}_{3}N_\mathrm{el}^{4}(Y_{0};x_{10},b)$ in Eq. (B14) we cannot use this replacement and have to calculate it using the saturation model of Ref. [28].

The last contribution, which proportional to $N_\mathrm{el}^2$, stems from Eqs. (37) and (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} \ge x^2_{01}$. Hence we need to subtract the following integral

$$\begin{aligned} I_4\,= & {} \,\frac{\bar{\alpha }_S}{2} x^2_{01}\,\int ^{x^2_{01}}_0 \, \frac{d x^2_{02}}{\varvec{x}^4_{20}} \Bigg (\,2\, N_\mathrm{el} \left( Y_0, \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \,\nonumber \\&- N^2_\mathrm{el} \left( Y_0, \varvec{x}_{02},\varvec{b} - \frac{1}{2}\varvec{x}_{02}\right) \Bigg )^2 \nonumber \\= & {} \bar{\alpha }_S\left( \frac{2}{ 2 \overline{\gamma }-1} N_\mathrm{el}^{2}(Y_{0};x_{10},b)-\frac{2}{3 \overline{\gamma }{-}1} N_\mathrm{el}^{3}(Y_{0};x_{10},b)\,\right. \nonumber \\&+\left. \frac{1}{2 (4 \overline{\gamma }-1)}N_\mathrm{el}^{4}(Y_{0};x_{10},b)\right) \end{aligned}$$

(B18)

Collecting all contributions we can see that we need to add

$$\begin{aligned} f^{\Sigma }_1\,\,=\,\,-\, 5.62\,\bar{\alpha }_S\,N_\mathrm{el}^{2}(Y_{0};x_{10},b)\end{aligned}$$

(B19)

Appendix C: Anomalous dimensions of the DGLAP equations

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, 33]).

$$\begin{aligned} \gamma _{ff} (\omega ) \,&= \,\frac{ C_F }{2 \,N_c}\, \left\{ \frac{3}{2} + \frac{1}{(1+\omega ) \, (2 + \omega )} \right. \nonumber \\&\quad \left. - 2 \, \psi (\omega + 2) + 2 \, \psi (1) \right\} \end{aligned}$$

(C1a)

$$\begin{aligned} \gamma _{gf} (\omega ) \,&= \,\,\frac{ C_F }{2 \,N_c}\, \left\{ \frac{1}{2+\omega } + \frac{2}{\omega \, (1+\omega )} \right\} \end{aligned}$$

(C1b)

$$\begin{aligned} \gamma _{fg} (\omega ) \,&= \, \frac{ N_f }{2 \,N_c}\, \left\{ \frac{1}{1+\omega } - \frac{2}{(2+ \omega ) \, (3+\omega )} \right\} \end{aligned}$$

(C1c)

$$\begin{aligned} \gamma _{gg} (\omega ) \,&= \, \frac{11 \, N_c - 2 \, N_f}{12\,N_c} + \left\{ \frac{1}{\omega \, (1+\omega )} {+}\frac{1}{(2{+} \omega ) \, (3{+}\omega )} \right. \nonumber \\&\quad - \psi (\omega + 2) + \psi (1) \Bigg \} \end{aligned}$$

(C1d)

where $\psi (w) = \Gamma '(w)/\Gamma (w)$ is the digamma function. Note that $\psi (1) = - \gamma _E$ where $\gamma _E$ is Euler’s constant. $N_f$ is the number of the fermions (quarks) and $C_F = (N^2_c - 1)/(2\,N_c)$.

Appendix D: Energy evolution of the saturation scale

In this appendix we discuss Eq. (45) as well as the fact that $\mathrm{Const}$ in Eq. (46) is proportional to $\ln \left( x^2_{01} Q^2_s\left( Y\right) \right)$. Eq. (45) can be written as

$$\begin{aligned} \ln \left( Q^2_s\left( \delta Y\right) \big / Q^2_s\left( Y_0\right) \right) \,\,=\,\,\bar{\alpha }_S\,\kappa \,\delta Y\,-\,\frac{3}{2\,{\bar{\gamma }}}\,\ln \left( \delta Y\right) \end{aligned}$$

(D1)

First term was derived in Ref. [34] and the second in Ref. [25]. Let us introduce a new variable $z_0 = \ln \left( x^2_{01}\,Q^2_s\left( Y_0\right) \right) \,+\,\bar{\alpha }_S\,\kappa \,\delta Y$. Eq. (D1) stems from two observations. First, that for $z_0 \,<\,0$ we can use for $n^D$ the solution of the linear equation:

$$\begin{aligned} n^D\left( z_0 \,<\,0\right) \,= & {} \,\int \frac{d \gamma }{2 \,\pi \,i} n^D_\mathrm{in}\left( \gamma ; \{\dots \}\right) \exp \Bigg \{\left( \chi \left( {\bar{\gamma }}\right) \right. \nonumber \\&\left. +\,\chi '\left( {\bar{\gamma }}\right) \,\left( \gamma - \overline{\gamma }\right) \,\right. \nonumber \\&\left. +\,\frac{1}{2}\chi ''\left( {\bar{\gamma }}\right) \,\left( \gamma - {\bar{\gamma }}\right) ^2\right) \,\bar{\alpha }_S\, \delta Y\,-\,\left( 1 - \gamma \right) \xi \Bigg \}\nonumber \\= & {} \int \frac{d \gamma }{2 \,\pi \,i} n^D_\mathrm{in}\left( \gamma ; \{\dots \}\right) \exp \Bigg \{{\bar{\gamma }}\,z_0 \, \,\nonumber \\&+\,\frac{1}{2}\chi ''\left( {\bar{\gamma }}\right) \,\left( \gamma - {\bar{\gamma }}\right) ^2\,\,\bar{\alpha }_S\, \delta Y\,-\,\left( {\bar{\gamma }} - \gamma \right) \,z_0\Bigg \}\nonumber \\ \end{aligned}$$

(D2)

The second observation is that $n^D$ for $z_0\,=\,0$ is not a constant by $n^D\left( z_0 = 0\right) \,\xrightarrow {\delta Y \gg 1} 1/\left( \bar{\alpha }_S\delta Y\right)$. This follows directly from semi-classical equations (see, for example, Ref. [1], section 4.5.3, Eq. 4.187). Since on the critical line $z_0=0$ we can use the solution of Eq. (D2) we see that $n^D\left( \gamma \{\dots \}\right) \,\,=\,\,\left( 2/\bar{\alpha }_S\right) \,\left( {\bar{\gamma }} - \gamma \right)$. For $\xi \,\ne \,\bar{\alpha }_S\,\kappa \,\delta Y$, the integral of Eq. (D2) taken by the method of steepest descent gives:

$$\begin{aligned} n^D\left( z_0 \,<\,0\right) \,=\,2\sqrt{\frac{2\,\pi }{\chi ''\left( {\bar{\gamma }}\right) }}\,\frac{z_0}{\left( \bar{\alpha }_S\,\delta Y\right) ^{3/2}}\,e^{{\bar{\gamma }} \,z_0} \end{aligned}$$

(D3)

In Eq. (D3) we consider $\delta Y \,\gg \,z_0$ and neglect the terms which rate proportional to $z^2_0/\delta Y$.

We need to take Eq. (D1) for the saturation momentum to compensate this behaviour. Plugging the new saturation momentum in Eq. (D2) we obtain that in the vicinity of the new saturation momentum, $n^D$ takes the following form

$$\begin{aligned} n^D\,\,=\,\,\mathrm{Const}\,\ln \left( x^2_{01}\,Q^2\left( Y\right) \right) \,\left( x^2_{01}\,Q^2_s\left( \delta Y\right) \right) ^{{\bar{\gamma }}} \end{aligned}$$

(D4)

Equation (D4) reproduces $n^D$ in the form obtained in Ref. [35].