CGC/saturation approach: an impact-parameter dependent model for diffraction production in DIS

In the paper we discussed the evolution equations for diffractive production in the framework of CGC/saturation approach, and found the analytical solutions for several kinematic regions. The most impressive features of these solutions are, that diffractive production does not manifest geometric scaling behaviour i.e. being a function of one variable. Based on these solutions, we suggest an impact parameter dependent saturation model, which is suitable for describing diffraction production both deep in the saturation region, and in the vicinity of the saturation scale. Using the model we attempted to fit the combined data on diffraction production from H1 and ZEUS collaborations. We found that we are able describe both xIP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{I\!\!P}$$\end{document} and β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} dependence, as well as Q behavior of the measured cross sections. In spite of the sufficiently large χ2/d.o.f.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^2/d.o.f.$$\end{document} we believe that our description provides an initial impetus to find a fit of the experimental data, based on the solution of the CGC/saturation equation, rather than on describing the diffraction system in simplistic manner, assuming that only quark-antiquark pair and one extra gluon, are produced.


I. INTRODUCTION.
In this paper we discuss diffractive production in the deep inelastic scattering in the framework of CGC/saturation approach (see Ref. [1] for review).In spite of the fact that the equations for the diffractive production in this approach, were proven long ago [2](see also Ref. [3][4][5]) the intensive study, during the past two decades, has been concentrated on the simplified model 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]).Such models described the experimental data quite well, giving the impression that we do not need to search for the solution of the general equations.Indeed, we found only two attempts to solve the equations of Ref. [2] numerically (see Refs, [14,15]).
The main goal of this paper is to investigate the non-linear equations for diffractive production in DIS, and to find an analytical solution in different kinematic regions.Based on these analytical solutions we will suggest an impact parameter dependent model in the spirit of Refs.[16,17] which is based on Color Glass Condensate/saturation effective theory for high energy QCD.
The paper consists of two parts.In the first part, which is the most important contribution in this paper, we found the analytical solutions of the evolution equations for diffraction production [2] in different kinematic regions, mostly using the approach developed in Ref. [18].In the second part of the paper, we suggest an interpolation formula which satisfies two limits found analytically: deep in the saturation region, and in the vicinity of the saturation scale, putting into practice the key ideas of Ref. [19].This formula exhibits the main features of the DIS amplitude , given in Refs.[16,17,20], but it is different from the interpolation procedure that have been used in numerous attempts to build a such model in Refs.[8,19,[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39].We will attempt to describe the HERA data in the region of small x I P and small β. [40].
Unfortunately, we are still doomed to build models to introduce the main features of the CGC/saturation approach, since the CGC/saturation equations do not reproduce the correct behavior of the scattering amplitude at large impact parameter (see Ref. [21,41]).Real progress in theoretical understanding of the confinement of quarks and gluon has not yet been achieved and, as a result, we do not know how to formulate the CGC/saturation equations to incorporate the phenomenon of confinement.We have to build a model which includes both the theoretical knowledge that stems from the CGC/saturation equations, and the phenomenological large b behavior that does not contradict theoretical restrictions [42,43].In our modeling of the large b behaviour of the scattering amplitude, we follow the main ideas of all saturation models on the market (see for example Refs.[8,19,[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]): and only introduce the non-perturbative behaviour in the b-dependence of the saturation scale.
For the b behavior we use the procedure, suggested in Ref. [17] 1 : where S (b) is the Fourier image of S (Q T ) = 1/ 1 + , and we will discuss below the value of γ.Eq. ( 1) leads to the scattering amplitude which is proportional exp (−mb) at b ≫ 1/m in accord with the Froissart theorem [42].In addition, we reproduce the large Q T dependence of this amplitude, which is proportional to Q −4  T and follows from the perturbative QCD calculation [43].This impact parameter behaviour is the main phenomenological assumption that we used.

II. THEORETICAL INPUT
In this section we discuss our theoretical input that follows from the Colour Glass Condensate(CGC)/saturation effective theory of QCD at high energies(see Ref. [1] for the basic introduction).

A.
The evolution equation for diffraction production in the framework of CGC.
A sketch of the process of diffraction production in DIS is shown in Fig. 1-c, from this figure one can see that the main formula takes the form where Y = ln (1/x Bj ) and Y 0 is the minimum rapidity gap for the diffraction process (see Fig. 1-c).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.

FIG. 1: The graphic representation of the processes of diffraction production
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 details and general descriptions 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 pom gg 1 are large.We are aware, that it is sufficient to describe most of th experimental data by only taking into account q q and q qG final states in diffraction production.Because of this, our main goal is not to describe the current experiments, but to study the solution to the equations in the LLA, which introduces the screening corrections to all the channels of the diffraction production.By comparing with the experimental data, we wish to determine in which kinematic region the shadowing corrections will become important, both for the elastic amplitudes in Fig. 1, as well as for the diffractive production of the large number of gluons.
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 [44]: viz.
Note, that r = x 01 and the kernel of the equation describe the decay of a dipole to two dipoles: The initial condition for Eq. ( 5) has the following form: Re-writing Eq. ( 5) as the equation for N D we obtain the second form of the set of the equations: The initial conditions are A general feature, is that the amplitude with fixed rapidity gap can be calculate as follows From Fig. 1  is the conjugate variable to the momentum transfer in Fig. 1 for the amplitudes N el .
B. N D deep in the saturation region: , and a violation of the geometric scaling behavior First, we consider the kinematic region, where In this region where both Y and Y 0 as well as the difference between them are large, we can expect that both N D and N el are close to 1. Therefore, we can use the procedure suggested in Ref. [18].In this region we can replace and linearize Eq. ( 5), neglecting ∆ D 2 terms.Indeed, Eq. ( 5) takes the form where we use notation and considered in Eq. ( 11) the impact parameter |b| ≫ |x 02 | and |x 12 |.The initial condition to Eq. ( 11), given by Eq. ( 6), can be re-written in the form In Eq. ( 12) we use the solution given in Re. [18] for ∆ el which has the form In Eq. ( 12) and Eq. ( 13) Q s (Y, b) is the saturation scale.Both these equations show the geometric scaling behaviour of the scattering amplitude [18,45,46], which depends on a single variable where ψ(x) = d ln Γ(x)/dx and Γ is the Euler gamma function [48].In this paper we will use the value of γ cr which comes from the leading order estimates: γ cr ≈ 0.37.
Neglecting the term proportional to (∆ D ) 2 in Eq. ( 5) and integrating over x 2 from 1/Q s (Y, b) to x 2 10 .The linear equation can be multiplied by arbitrary function of Y 0 and x 10 .Bearing this in mind, the solution has the following form: Note that function G can be found from the initial condition of Eq. ( 12), leading to the final answer: The most impressive feature of the solution is that the function does not show the geometric scaling behaviour i.e. being a function of one variable.The solution is the product of two functions: one has a geometric scaling behaviour depending on one variable z, and the second depends on z 0 , showing geometric scaling behaviour in the same way, as elastic scattering amplitude at Y = Y 0 .In this subsection we consider the kinematic region in which N el (Y 0 , x 01 , b) is in the vicitity of the saturation scale As it was found in Ref. [46] we have geometric scaling behaviour in this region, and the amplitude behaves as If Y is close to Y 0 , we can neglect the non-linear terms in Eq. ( 7), and we can solve the linear equation for N D with the initial condition of Eq. ( 8) 1. Solution in the region where Nel < 1 Considering Eq. ( 18), one can see that in this kinematic region we can in general neglect the expression of Eq. ( 7) two terms: 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 Taking derivatives with respect to Y 0 , we re-write Eq. ( 19) for the amplitude n D (Y, Y 0 , x 01 , b) that has been introduced in Eq. ( 9).It takes the form of the linear equation: The initial condition for this equation is the following: The elastic amplitude has the form: where ξ ≡ ln 1/ x 2 10 Q 2 0 .Taking the double Mellin transform, and defining we obtain the solution to the equation of Eq. ( 20) in the following form: where φ in has to be determined from the initial condition of Eq. (36), and it has the form Therefore, the solution takes the form (see Fig. 1 for notations): with γ cr = 0.37 and γ = 1 − γ cr = 0.63 ,γ = −1 + 2γ cr = −0.26.χ (γ) is given by Eq. ( 14).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. Since ξ > 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 .At ξ ≫ ᾱS δY , we can integrate by the same method, 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 χ(γ).We cannot use the method of steepest decent for such small values of δY .

Solution in the saddle point approximation
Using solution of Eq. ( 26) we can find the saturation momentum for the process of diffraction production, calculating the integral over γ by the method of steepest descent.Q s can be determined from the following two equations: Eq. ( 27) is the equation for the saddle point while Eq. ( 28) is the condition that the solution is a constant on the critical line x 2 10 = 1/Q 2 s .The solution of these two equation is well known γ = γ cr = 0.37 and In the vicinity of the saturation scale Eq. ( 26) behaves as One can see that there is no geometric scaling behavior of the scattering amplitude n D , even at large Y − Y 0 .We can also see that the solution does not satisfy the initial condition.It stems from Eq. ( 27) and Eq. ( 28), which both are correct only if δY ≫ 1, assuming that the saddle point value of γ is not close to γ.We need to re-write these equations to take into account the possibility that γ SP → γ taking into account the factor 1/(γ − γ).The equations for the saddle point take the form: The contribution of the additional term is essential, only if γ SP → γ, but even γ = 0 which corresponds the integration with the contour C 3 (see Fig. 2), is still not close to γ. Bearing this in mind, we prefer to treat δY ≪ 1 without using the method of steepest descent.In this section we will return to the discussion of Eq. ( 26) at δY = Y − Y 0 ≤ 1.In this kinematic region, we cannot use the method of steepest descent, and have to look for a different approach.First, let us analyze the solution iterating the equation keeping δ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 it over γ.This integral takes the following form for the third term in Eq. ( 26) for n ≥ 1: For n = 0, we have the contribution only of the second term in Eq. ( 33).In Eq. ( 32) 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. ( 33) can be re-written as follows The last term in Eq. ( 34) 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. ( 26) in the following form The first term in Eq. ( 35) 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 ξ ≫ 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. (34).
4. Solution in the region where Nel ∼ 1 In this kinematic region we need to keep the term which is proportional to N el N D in Eq. ( 6) and solve the equation which takes the form Therefore, we took into account terms n D N el in comparison with the previous sections.We consider that n D N D are sufficiently small, so that we can neglect the contributions n D N D ∼ N 4 el ≪ 1.This equation looks simpler in the momentum representation: In the momentum representation Eq. ( 7) takes the form: We solve this equation using the semi-classical approach.In this approach we are looking for the solution in the form with where are smooth functions of Y and ρ, and the following conditions are assumed: Plugging Eq. ( 41) and Eq. ( 42) in Eq. (36) we have For the equation in the form where S is given by Eq. ( 42), we can introduce the set of characteristic lines on which ρ(t), Y (t), S(t), ω(t), and γ(t) are functions of the variable t (which we call artificial time), that satisfy the following equations: In Eq. ( 48) we consider that with κ = χ (γ cr ) /(1 − γ cr ).
In Fig. 3 we plotted the numerical solutions of Eq. (48).One can see that the solution gives γ (Y ), which approach a constant at large Y − Y 0 .This feature stem from Eq. ( 48)-4 since N el → 0 at large ρ.
Bearing this in mind one can see that Eq. ( 46) and Eq. ( 48) degenerate to Eq. ( 19) in the semiclassical approach, at large Y .Indeed, Fig. 4 shows that the solution to Eq. ( 19) shown in dashed lines in Fig. 4, is close to the solution of Eq. ( 47) for Y ≥ 2.
Therefore, we need to consider the solution of Eq. ( 26) in the entire kinematic region where we can neglect the term proportional to N D 2 in Eq. (7).48).Constant c in Eq. ( 39) is chosen c=0.05 in accord with the description of the HERA data in Ref. [17].αS = 0.25.
D. Solution in the region where Nel ≪ 1 In this section we consider the solution, in the region where N el (Y 0 , x 01 , b) is a solution to the linear BFKL equation, which has the following general form where n in (γ, b) should be found from the initial condition at Y 0 = 0. Recall, that Y 0 ≡ ᾱS Y 0 , as we have discussed above and ξ = ln 1/(x 2 10 Q 2 0 .In the region of large Y 0 , we take only the leading term in χ (γ) γ≪1 − −− → 1/γ into account, and take the integral by the method of steepest descent.As the result we obtain a solution in the double log approximation (DLA) of perturbative QCD.At large Y 0 we can re-write N 2 el (Y 0 , x 01 , b) in the form where MOSD means the method of steepest descent.Bearing Eq. ( 51) in mind, we obtain the solution to the BFKL equation for N D in the form Calculating n D (Y 0 , γ) from Eq. ( 51) we obtain Taking the integral over γ in the DLA, the equation for the saddle point takes the form Eq. ( 54) has four solutions (see Fig. 5).Two of them have imaginary parts while other two are real.At . In this limit we see that our solution satisfies the initial conditions in the DLA.
For large Y − Y 0 the solution has then form From Eq. (55) one can see that the solution reaches the saturation bound at and behaves in the vicinity of this bound as Saddle point

III. THE MODEL
As we have discussed, the main idea of building the saturation model has been formulated in Ref. [19]: it is the matching of two analytical solutions in the vicinity of the saturation scale, and deep inside of the saturation domain.
A. The input: Nel (r, Y ) in our saturation model As we have mentioned, the initial condition for the equation for N D (see Eq. ( 8)) is determined by N el , which we have found from the HERA data for the deep inelastic structure function in Ref. [17].For completeness of presentation we describe the main formulae of this model which illustrates our procedure for the model building.
In the vicinity of the saturation scale or, in other words, for the dipole size r in the region: where N 0 is the phenomenological parameter that has been found in Ref. [17].Q s is the saturation momentum which we will discuss below.The values of γ cr can be found from the following equation: In Eq. ( 59) χ (γ) is the BFKL kernel that takes the form where ψ (z) is the digamma function.
Deep inside of the saturation domain, where τ ≡ r 2 Q 2 s (Y, b) ≫ 1, we use the analytical solution to the non-linear equation given in Ref. [18] where where ξ and λ are related to the behaviour of the saturation scale where Y in = ln(1/x in ) shows the initial value of x from which we start low x evolution.This is a phenomenological parameter of the model.The phenomenological dependence on b of the initial saturation scale Q 2 (Y = Y in , b) we have discussed in the introduction (see Eq. ( 1)).Finally, we use the following Using Eq. ( 14)-Eq.( 64) one can see that ξ = ln r 2 Q 2 (Y in ; b) in Eq. ( 14).The value of λ can be calculated and it is equal to λ = ᾱS χ (γ cr ) /(1 − γ cr ).However, in describing the experimental data in Ref. [17], we consider λ as the independent fitting parameter, since the next-to-leading correction turns out to be large.
Parameter A in Eq. ( 61) should be found from the matching procedure of two solution at z = z m : However, it turns out that Eq. ( 61) cannot satisfy Eq. ( 65).In Ref. [20] the correction to Eq. ( 61) has been found.The amplitude of Eq. ( 61) with the corrections takes the form: which has been used in the Eq. ( 65).
It should be stressed that all phenomenological parameters for the elastic amplitude has been extracted from the experimental data for F 2 at HERA (see Table 1).[17].λ, N0, m and Q 2 0 (Q 2 0 = m 2 x λ in ) are fitted parameters.The masses of quarks are not considered as fitted parameters and two sets of parameters, that are shown in the table, relate to two choices of the quark masses: the current masses and the masses of light quarks are equal to 140 M eV which is the typical infra-red cutoff in our approach.

B.
N D (r, Y, Y0) in the model: matching procedure In the kinematic region where 37) which we re-write as follows The first term in Eq. ( 67) corresponds to the first term in Eq. ( 37), which we simplify taking into account the experience in the description of the elastic amplitude (see [8,16,17,19,[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]).It was shown in these papers that we can describe the solution to the evolution equation taking the amplitude in the form of Eq. (58) replacing 1 − γ cr in Eq. ( 58) and in Eq. (67) by the following expression where λ = ᾱS (χ (γ cr ) / (1 − γ cr )) and κ = χ ′′ (γ cr ) /χ ′ (γ cr ).The factor Y − Y 0 is introduced to reflect the general features of the first term in Eq. ( 37) which is proportional to this factor.The saturation momentum where Q s (Y 0 , b ′ ) is the initial transverse momentum at Y = Y 0 .However, as we have discussed above, we introduce the non-perturbative corrections in the behaviour at large impact parameter in the b-dependence of the saturation scale.Bearing this in mind we use the following parameterization of where parameters Q 0 , λ and mass m in S (b) (see Eq. ( 64)) have been determined in our previous paper [17] from the fit of the elastic data.The parameter m 1 has to be extracted from the fit of the diffraction production as well as parameters N 1 and λ 1 .In the leading order of perturbative QCD λ 1 = χ (γ), but as well as for the energy behaviour of the saturation momentum, the higher order corrections are large, and we view these two parameters: λ and λ 1 as the phenomenological parameters which we have to extract from the experimental data.We need to use the matching procedure analogous to Eq. (65) for n D = − dN D dY 0 using Eq. ( 15) in the form The matching equations take the form: Eq. ( 72) allow us to specify function G (τ 0 , Y 0 ) in Eq. ( 71).We re-write Eq. (67) using τ (b) to simplify Eq. ( 72) in the form: From Eq. ( 73) we see that function G τ 0 , Y 0 , b, b ′ can be written as Eq. (72) degenerate to the following matching conditions for λ 1 (Y − Y 0 ) ≫ 1: One can see that Eq. ( 75) does not have a solutions for z m > 0. It has been found in Ref. [20] that the solution deep inside of the saturation scale has more general form than we used in Eq. (71): where Eq. (75) takes the form The solution to Eq. ( 78) is We can use Eq. ( 77) to determine the value of B. However, bearing in mind the many simplifications that we have assumed , we decided to view z m as a free parameter, which we will find from the fit of the experimental data.
In this region we use Eq. ( 16), which for n D takes the form: where z (b) is defined in Eq. ( 14).The elastic amplitude is equal to For practical purpose we define this region as z (b) ≥ z m , where z m is the matching point for DIS.

Kinematics and observables
The experimental data for diffractive production in DIS (see Ref. [40]) are presented using the following set of the kinematic variables: x Bj = β x I P ; where Q 2 is the virtuality of the photon and M X is the produced mass.The set of kinematic variables, that we used, has the following relation to Eq. ( 82): The main formulae that we use to calculate the experimental cross sections are given by Eq. ( 2) and Eq.(3).Eq. ( 3) can be re-written in the form, which includes the integration over impact parameter b ′ , in the form The expression for (Ψ * Ψ) γ * ≡ Ψ γ * (Q, r, z) Ψ γ * (Q, r, z) in Eq. ( 2) is well known (see Ref. [1] and references therein) where T(L) denotes the polarization of the photon and f is the flavours of the quarks.
C. Description of the HERA data Using Eq. (66) we attempted to describe the combined set of the inclusive diffractive cross sections measured by H1 and ZEUS collaboration at HERA [40].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 σ (β, Q 2 , x I P , t) are presented at different values of Q, β and x I P .This cross section is equal to Q 2 4π 2 σ diff Y, Y 0 , Q 2 where σ diff Y, Y 0 , Q 2 is given by Eq. ( 2).We view this paper as the next step in building the saturation model based on the CGC approach.The first step have been done in Refs.[17] where we build the saturation model for the DIS.The parameters that we found from this fit and which are shown in Table 1 we use for the diffractive production as given and we are not going to change them.The additional parameters that we used to parametrize the diffraction production cross section are N 1 , m 1 and λ 1 .N 1 is proportional to ᾱS which indicates that the typical values of N 1 is small.λ 1 = ᾱS χ (γ) ≈ 3.67 ᾱS in the leading order of perturbative QCD.However, we consider this as a fitting parameter since we expect that it will be heavily affected by the next order calculation.Recall, that the value of λ which is equal to λ = ᾱS χ (γ cr ) /(1−γ cr ) ≈ 4.88 ᾱS came out λ ≈ 0.2 from the fit and this value is in accord with the next to leading estimates.
First, we found a fit within parameters are equal N 1 = 7.7 10 −4 , λ 1 = 1.58 and m 1 = 2 GeV .The large value of m 1 which describes the non-perturbative behavior of Q s Y − Y 0 , b − b ′ led us to the idea that even the non-perturbative behavior of this saturation momentum stems from the CGC physics and determined by Q s (Y 0 , b).Therefore, we fitted the data fixing m 1 = Q s (Y 0 , b).It turns out that with the parameters of the fit: N 1 = 7. 10 −4 and λ 1 = 1.48, we found the description of the experimental data shown in Fig. 6.Actually, the first fit give the description of the data of the same quality as the second one; and the resulting curves for both fits look the same and cannot be differentiated in the figures.In spite of the fact that the quality of the fit is not good we see that Eq. (66) reproduces both x I P and Q dependence.6-a) and versus xIP at fixed β and Q 2 (Fig. 6-b).The data are taken from Ref. [40].The red curves show the results of the fit with m1 = Qs (Y0, b).The fit in which m1 was a fitted parameter turns out to be so close to this fit that cannot be clearly shown in the picture in spite of the fact that has higher value of χ 2 /d.o.f..
We do not expect a good description of the data as we have mentioned.As was expected the values of parameter λ 1 turns out to be quite different from the leading order estimates in perturbative QCD, which illustrate the need for the next to leading order corrections.We notice that in the experimental kinematic region Therefore, the most data are in the region which is outside of the saturation domain.Hence, Eq. ( 37) and Eq.(A14), which sums the emission of several gluons in perturbative QCD, has to be tried to describe the data.The success of the simple model for diffraction production: production of q q and q qG states[6-13], indicates, that taking into account the emission of several gluons, we will be able to describe the data.We are going to try in a separate further publications.On the other hand we see that the main contribution stems from the gluon emission.Indeed, in Fig. 7 we plot the two different terms of Eq. (66) writing it as n D = n D 1 + n D 2 where n D 2 ∝ N 2 el .One can see that the emission of gluons (the term n D 1 ) is certainly larger than the contribution of the diffractive production of quark-antiquark state (term n D 2 ).We view this fact as an argument that we have to take into account a large number of emitted gluons.
Β 0.0562, x P 0.0160 It should be mentioned that we have some hidden parameters which mostly specify the region of the applicability of the perturbative QCD estimates.For example, even in the case of deep inelastic processes, we can trust the wave function of perturbative QCD only, at rather large values of Q 2 ≥ Q 2 0 with Q 2 0 ≈ 0.7GeV 2 since for smaller Q it will be affected by the non-perturbative contributions.As we have mentioned that we consider the matching point z m of Eq. (78) as a fitting parameter due to a large uncertainty in the calculation of B given by Eq. (77).From the fit we specify them as β ≤ 0.056,x I P ≤ 0.025, 0.7 ≤ Q 2 ≤ 27 GeV 2 .We found that the fit does not depend on the value of the matching point z m .This confirms that the data are outside of the saturation region.

IV. CONCLUSIONS
In the paper we discussed the evolution equations for the diffractive production in the framework of CGC/saturation approach that have been proposed in Ref. [2] and found the analytical solutions in the several kinematic regions.The most impressive features of these solutions are that the diffractive production does not show the geometric scaling behaviour being a function of one variable.Even deep in the saturation regions for both diffractive and elastic amplitude the solution turns out to be the product of two functions: one has a geometric scaling behaviour depending on one variable z, and the second depends on z 0 showing the geometric scaling behaviour in the same way as elastic scattering amplitude at Y − Y 0 .
Based on these solutions we suggest a impact parameter dependent saturation model which is suited for the describing the diffraction production both deep in the saturation region and in the vicinity of the saturation scale.Since we are dealing in the diffraction production with two saturation scales: Q (Y, b) and Q Y − Y 0 , b − b ′ , where Y = ln(1/(x I P β) and Y 0 = ln(1/β), the model includes more information from the theoretical part of the paper than it has been needed for the inclusive DIS.However, the main key assumptions of the model are the same as for inclusive DIS: the non-perturbative impact parameter behavior is absorbed into two saturation scales.
Using the model we tried to fit the combined data on diffraction production from H1 and ZEUS collaborations [40].We fond that we are able describe both x I P and β dependence as well as Q behavior of the measured cross sections.In spite of the sufficiently large χ 2 /d.o.f.we believe that our description give the starting impetus to find a fit of the experimental data based on the solution of the CGC/saturation equation rather than on describing the diffraction system in smlistic way assuming that onle quark=antiquark pair and one extra gluons are produced.
Ae the result of the fit we found out that the experimental data are concentrated in the region outside the saturation domain for the produced diffractive system and we intend to try to sum the multi-gluon production in perturbative QCD approach using the formulae which we found in this paper.We are going to publish the result of this kind of approach elsewhere.
We believe that this paper will revive the interest to the process of the diffractive production which is a unique process which description needs the understanding both the multi-particle generation process and the elastic (diffractive) rescattering at high energy.

FIG. 6 :
Fig. 6-a Fig. 6-b FIG.6:σ diff Y, Y0, Q 2 = xIP σr versus Q 2 atfixed β and xIP (Fig.6-a) and versus xIP at fixed β and Q 2 (Fig.6-b).The data are taken from Ref.[40].The red curves show the results of the fit with m1 = Qs (Y0, b).The fit in which m1 was a fitted parameter turns out to be so close to this fit that cannot be clearly shown in the picture in spite of the fact that has higher value of χ 2 /d.o.f..

FIG. 7 :
FIG.7: Contributions of two different sources of the diffractive production: the production of quark-antiquark state (n D2 ) and the multi-gluon production (n D 1 ).

TABLE I :
Parameters of the model which has been extracted from DIS experiment in Ref.