Solution to the Sudakov suppressed Balitsky-Kovchegov equation and its application to the HERA data

We analytically solve the Sudakov suppressed Balitsky-Kovchegov evolution equation with the fixed and running coupling constants in the saturation region. The analytic solution of the $S$-matrix shows the $\exp(\mathcal{O}(\eta^2))$ rapidity dependence of the solution with the fixed coupling constant is replaced by $\exp(\mathcal{O}(\eta^{3/2}))$ dependence in the smallest dipole running coupling case rather than obeying the law found in our previous publication, in which all the solutions of the next-to-leading order evolution equations comply with $\exp(\mathcal{O}(\eta))$ rapidity dependence once the QCD coupling is switched from the fixed coupling to the smallest dipole running coupling prescription. This finding indicates that the corrections of the sub-leading double logarithms in the Sudakov suppressed evolution equation are significant, which compensate part of the evolution decrease of the dipole amplitude made by running coupling effect. To test the analytic findings, we calculate the numerical solutions of the Sudakov suppressed evolution equation, the numerical results confirm the analytic outcomes. Moreover, we use the numerical solutions of the evolution equation to fit the HERA data. It shows that the Sudakove suppressed evolution equation can give good quality fit to the data.

is unstable [22]. Mathematically, the reason for this difficulty was traced back to a large double transverse logarithmic correction in the evolution kernel of the full NLO BK equation. The physics behind the double logarithms (also called anti-collinear logarithms) is the time-ordering of the successive gluon emissions.
To cure the instability issues, one has to resum the radiative corrections enhanced by the double transverse logarithms to all orders. There are two approaches proposed to perform the resummations [17,21]. One of the strategies for enforcing the time-ordering in the evolution with rapidity Y of dipole projectile is to put kinematical constraints in the evolution kernel leading to a non-local equation in Y [17]. The other is to resum the double logarithmic corrections to all orders giving rise to a local collinearly improved Balitsky-Kovchegov (ciBK) equation in Y [21]. These two approaches are equivalent to each other in the leading double logarithmic level, although they bring the modifications to the structure of the evolution equation. Moreover, the condition of the time-ordering also takes the modifications to the corresponding initial conditions which are required for solving the evolution equations. However, the modifications of the initial conditions have been not properly implemented in the both approaches aforementioned, which lead to a matter of fact that the modifications in the initial condition impact not only on the evolution of the dipole amplitude in the region of low Y but also on the asymptotic behavior in the region of large Y [26]. This is an unexpected result, in fact that the high energy asymptotic behavior of the dipole amplitude should be inappreciably affected by the formulation of the initial condition.
To overcome the instability problems and solve it fundamentally, an effective method was proposed by the authors in Ref. [26] inspired from previous experience on handling similar issues of the NLO BFKL 2 equation [30][31][32]. In Ref. [26], they used a new rapidity variable (rapidity of the target, η) instead of the rapidity of the projectile (Y ) to be as the "evolution time", and reorganized the perturbative QCD theory for the evolution of the dipole amplitude. The advantage of this method is that the time-ordering condition is automatically satisfied, and then the anti-collinear contributions are absent in the target rapidity evolution. Furthermore, this choice of evolution variable is more reasonable, since the rapidity of the target is the one indeed used in the DIS rather than Y . A new version of non-local collinearly improved Balitsky-Kovchegov (non-local ciBK) equation was obtained [26], which was shown to give a rather good fits to the HERA data [33]. Soon after the non-local ciBK equation was established, it was found that there are important corrections to the evolution kernel from the sub-leading double logarithms located beyond the strong time-ordering region, it was shown that the sub-leading double logarithms arise from the incomplete cancellation between real and virtual corrections and are the typical Sudakov type ones [27]. When these double logarithms are resummed to all orders, a Sudakov suppressed Balitsky-Kovchegov (SSBK) equation in the evolution of the target rapidity is obtained [27]. The kernel of the SSBK equation is modified significantly by the sub-leading double logarithms.
In this paper, we shall solve analytically the SSBK and non-local ciBK equations in the saturation region with the smallest dipole running coupling prescription (SDRCP). Note that the reason why we use the SDRCP is that it was shown to be an effective QCD coupling in our previous publication [34]. To see the variances in the solutions of two types of evolution equations (based on rapidity of projectile Y and on rapidity of target η) due to the change of evolution variable from Y to η, we firstly recall the analytic solutions of the LO BK, rcBK and ciBK equations in Y , and we shall use these solutions for the comparisons with the ones of the respective evolution in η. We find that the analytic solutions of the non-local ciBK and SSBK in η with the fixed coupling Eqs. (37) and (44) are similar to the one gained at LO BK in Y Eq. (9), except that the coefficients in the exponent are different. We also find that the solution of the non-local ciBK in η with the SDRCP Eq.(39) is analogous to the one obtained at rcBK in Y Eq. (20). Surprisingly, the analytic solution of the SSBK equation with SDRCP shows that the rapidity in the exponent of the S-matrix has rapidity raised to the power of 3/2 dependence, exp(O(η 3/2 )), instead of linear rapidity dependence exp(O(η)), which do not obey the law found in our previous studies [34], in which we showed that the solutions of all kinds of NLO BK equations in Y with SDRCP have linear rapidity dependence in the exponent of the S-matrix. Coincidentally, the rapidity dependence of the solution to the SSBK in η with SDRCP is similar to the one obtained at the full NLO BK in Y , exp(O(Y 3/2 )), with the parent dipole running coupling prescription (PDRCP) [25,28,34].
To test the analytic findings mentioned above, we numerically solve the SSBK and non-local ciBK equations with the fixed and running coupling constants, via focusing on the physics in the saturation region. The numerical results confirm our analytic findings, see Fig.2. Furthermore, the SSBK equation is used to fit the HERA data. It shows that the theoretical calculations almost overlap with all the data points, see Fig.3, and a reasonable value of χ 2 /d.o.f = 1.128 is obtained from the fit, which indicate that the SSBK equation can give a rather good description to the data.

II. THE LEADING ORDER, RUNNING COUPLING AND COLLINEARLY IMPROVED EVOLUTION EQUATIONS OF COLOR DIPOLES IN Y -REPRESENTATION
In this section, we give a brief review of the LO BK, rcBK, ciBK equations in Y -representation in order to collect the basic elements of the BK equations, which shall be useful for the latter discussion. The review shall also give us a chance to introduce notations and explain the kinematics of color dipoles.
A. The Balitsky-Kovchegov equation and its analytic solution at leading order Consider a high energy dipole which is consisted of a quark-antiquark pair, scattering off a hadronic target, in the eikonal approximation the dipole scattering matrix can be written as a correlator of two Wilson lines [16] S(x, y; where x and y are the transverse coordinates of the quark and antiquark, and U is the time ordered Wilson line with A + (x − , x) to be as the gluon field of the target hadron. Note that the average in Eq.(1) is given by the average over the target gluon field configurations at a fixed rapidity. The rapidity evolution of the dipole scattering matrix can be described by the BK evolution equation [1,6] where the K LO (x, y, z) is LO evolution kernel describing the dipole splitting probability density, and has form with coupling being redefined asᾱ s = α s N c /π. The r = x − y, r 1 = x − z and r 2 = z − y in Eq.(4) are the transverse sizes of the parent dipole and two emitted daughter dipoles, respectively. In the large N c limit, the Eq.(3) describes the evolution of the original dipole (x, y) splitting into two daughter dipoles, (x, z) and (z, y) sharing a common transverse coordinate of the emitted gluon z. The non-linear term in S on the right hand side (r.h.s) of Eq.(3) depicts the two daughter dipoles interaction with the target simultaneously, which is usually called as "real" term due to its really measurement of the scattering of the soft gluon. While the linear term in Eq.(3) is referred as "virtual" term since it gives the survival probability for the original dipole at the time of scattering. We would like to note that the BK equation is a LO evolution equation since it resums only the leading logarithmic α s ln(1/x) corrections in the fixed coupling case. Meanwhile, the BK equation is a mean field version of the Balitsky-JIMWLK hierarchy [2][3][4][5] equations, in which higher order correlations are neglected. Now, let us solve the BK equation analytically in the saturation region. In this region we know that the parton density in the target is so high that the interaction between the dipole and target is very strong, which leads to the scattering amplitude close to unit, namely T ∼ 1. Via using the relation between the scattering amplitude and scattering matrix, T = 1 − S, one can get S ∼ 0. Therefore, one can neglect the contribution from the non-linear term, the Eq.(3) becomes To get the solution of the Eq.(5), we need to know the upper and lower integral bounds. The lower integral bound can be set to 1/Q s since the saturation condition requires that the size of the dipole should be larger than the typical transverse size r s ∼ 1/Q s . Here the Q s is the saturation scale which is an intrinsic momentum scale playing a role of the separation of dilute region from saturation region. The upper integral bound can be set to the original dipole size r although few daughter dipoles may have size larger than r, however the evolution kernel is rapid decrease when r 1 (r 2 ) > r. Hence, the contribution from those dipoles can be negligible. With the upper and lower bounds, the Eq.(3) can be rewritten as In the saturation region, one can find that the integral in Eq.(6) is governed by the region either from the transverse coordinate of the emitted gluon approaching to quark leg of the parent dipole, 1/Q s |r 1 | |r| and |r 2 | ∼ |r|, or the transverse coordinate of the emitted gluon approaching to antiquark leg of the parent dipole, 1/Q s |r 2 | |r| and |r 1 | ∼ |r|, see Fig.1. In this study, we work in the region |r 2 | ∼ |r|, the evolution kernel simplifies to and the Eq.(6) can be rewritten as where the factor 2 on the r.h.s comes from the symmetry of the aforementioned two integral regions. If one performs the integrals over variables r 1 and Y in Eq.(8), one can get the analytic solution of the LO BK equation in the saturation region as [35,36], where the c is a constant coming from the saturation momentum Q 2 The solution in Eq.(9) has been firstly derived in Ref. [35], which is called as Levin-Tuchin formula. One can see that the exponent of the scattering matrix S has a quadratic rapidity dependence, which leads to the scattering matrix too small when the rapidity is large. In terms of the relation between the scattering matrix S and scattering amplitude T , T = 1 − S, we know that the evolution speed of scattering amplitude is too fast, which renders the LO BK equation insufficiently to describe the experimental data from HERA [8,9]. So, one has to take into account the NLO order corrections to the LO BK equation, such as running coupling effect which can bring modifications to the evolution kernel leading to reduce the evolution speed of the dipole amplitude.

B. The Balitsky-Kovchegov equation and its analytic solution in the case of running coupling
The strong coupling constant α s in the LO BK equation (3) was assumed to be constant when the LO BK equation was derived, which makes the LO BK equation to be as a leading logarithm (LL) accuracy evolution equation. A naive way to promote it to include NLO correction is to replace α s → α s (r 2 ) in the LO BK equation [37] K PDRCP (r, r 1 , where the argument of the coupling constant is the transverse size of the parent dipole. For the α s , we shall use the running coupling at one loop accuracy There is another running coupling prescription proposed recently in Refs. [33,34,38,39] where they found that using the size of the smallest dipole to be as the argument of the coupling constant is favored by the HERA data at a phenomenological level, which is referred to the SDRCP. In the case of the SDRCP, the kernel can be written as We would like to point out that all the following studies referred to running coupling in this paper shall use the SDRCP unless otherwise specified prescription. The Eq.(10) is a naive way to include the running coupling effect into the LO BK equation, which is insufficient. The real running coupling corrections to the LO BK equation were calculated by including the contribution from the quark bubbles in the gluon lines. The calculations have been performed with two prescriptions by Balitsky in Ref. [14] and Kovchegov and Weigert in Ref. [15]. We don't go details of the derivations of the rcBK equation, and just quote the results here, since it is out of interest of present paper.
The rcBK equation reads where K rc (r, r 1 , r 2 ) is the running coupling evolution kernel. One can find that the rcBK equation Eq. (13) has the same structure as the LO BK equation Eq.(3) but the evolution kernel modified by the running coupling effect. There are two different (Balitsky and Kovchegov-Weigert) prescriptions for the kernel of rcBK equation as mentioned above. In the case of Balitsky prescription [14], the running coupling kernel can be written as Under the Kovchegov-Weigert prescription [15], the running coupling kernel reads with Note that we have found an interesting result in our previous studies in Ref. [40], in which the Balitsky and Kovchegov-Weigert kernels reduce to the same one under the saturation condition which indicates that the running coupling kernel is independent of the choice of prescription in the saturation region.
To analytically solve the rcBK equation in the saturation region, we substitute the simplified kernel Eq.(17) into Eq.(13), and get Under the saturation condition, one knows that the scattering matrix is small, therefore one can neglect the quadratic term of the S-matrix in Eq. (18). The Eq.(18) reduces to where the upper and lower bounds are determined using the same way as the LO BK equation in section II A, and the factor 2 on the r.h.s of Eq. (19) results from the symmetry of the two integral regions, see Fig.1. Doing the integral over the variables r 1 and Y in Eq. (19), one obtains the analytic solution of the rcBK equation [40] S where the NLO saturation momentum is used In Eq. (20), one can see that the exponent of the S-matrix has a linear rapidity dependence in the running coupling case, while the exponent of the S-matrix quadratically depends on the rapidity in the fixed coupling case, see Eq. (9). The change of rapidity of the S-matrix from quadratic to linear dependence implies that the evolution speed of the dipole amplitude is slowed down by the running coupling effect, which is in agreement with the theoretical expectations [14,29]. In addition, the phenomenological studies of the HERA experimental data in Refs. [8,9] showed that the rcBK equation takes a significant improvement of the description of the HERA data, which indicate that the NLO corrections are important.

C. The collinearly improved Balitsky-Kovchegov equation and its analytic solution
It is known that the rcBK equation includes only the contributions from the quark bubbles, while a full NLO evolution equation should consider contributions from the quark and gluon bubbles as well as from the tree gluon diagrams with quadratic and cubic nonlinearities. The authors in Ref. [16] have performed a comprehensive derivation, and gotten a full NLO BK equation which includes all the NLO corrections just mentioned. However, it has been found that the full NLO BK equation is unstable, since the dipole amplitude resulting from the full NLO BK equation can decrease with increasing rapidity, it can even turn to a negative value [22]. The reason for this instability is traced to a large contribution from a double-logarithm in the kernel of the full NLO BK equation [22]. To solve the unstable problem, a novel method was devised in Ref. [21] to resum double transverse logarithms to all orders, they obtained a resummed BK equation which governs the evolution in the double logarithmic approximation (DLA). Soon after the DLA BK equation released, they found that the single transverse logarithms (STL) also have a large corrections to the BK equation. Combining the single and double logarithmic corrections together, they got a collinearly improved (ci) BK equation, which reads [38] ∂S(r, Y ) where the collinearly improved evolution kernel is [23] K ci (r, r 1 , r 2 ) =ᾱ s 2π with and The constant A 1 = 11/12 in Eq. (24) We would like to point out when ln r 2 1 /r 2 ln r 2 2 /r 2 < 0, then an absolute value is used and the Bessel function of the first kind J 1 changes to the modified Bessel function of the first kind I 1 [38].
To analytically solve Eq. (22) in the saturation region, one needs to employ the saturation condition which means the scattering matrix is very small. Thus, one can neglect the quadratic term in Eq. (22) and keeps only the linear term, then the Eq. (22) becomes where the upper and lower bounds of the integral are determined by using the same way as the ones in section II A. As it was done in previous subsections, we work in the regime, 1/Q s |r 1 | |r|, |r 2 | ∼ |r|. So, we have ln(r 2 2 /r 2 ) 0 leading to ρ = 0, which implies the DLA kernel K DLA 1. This outcome confirms the statement that the double logarithm only plays a significant role in the week scattering phase-space [21]. Under the approximation just mentioned above, one has whose solution is Note that the dominant terms in the exponents of Eqs. (29) and (20) have linear rapidity dependence once the SDRCP is applied, which is a law found in Ref. [34]. This outcome implies that the running coupling correction is a dominant effect over all the aforementioned NLO corrections, like resummations of double and single transverse logarithms, in the suppression of dipole evolution.

III. THE NON-LOCAL NEXT-TO-LEADING ORDER AND SUDAKOV SUPPRESSED EVOLUTION EQUATIONS OF COLOR DIPOLES IN η-REPRESENTATION
In the previous section, all the BK equations are derived by following the evolution in terms of the rapidity of the projectile (Y ). However, the recent studies in Ref. [26] found that the NLO BK equation in Y must to be re-established according to the rapidity of the dense target (η), since the evolution of the rapidity Y could cause the instability to the NLO BK equation. In this section, we shall discuss the non-local ciBK equation in η and its extended version which includes sub-leading double logarithms from the region beyond the strong time-ordering [27]. The evolution equations in η-representation shall be solved analytically in the fixed and running coupling cases in the saturation region, respectively. The results are compared to the solutions obtained in the Y -representation.

A. The non-local collinearly improved BK equation in η and its analytic solution
The non-local ciBK equation was reformulated via change of variables η ≡ Y − ρ to transform the BK equation from the Y -representation to the η-representation, where ρ = ln 1/(r 2 Q 2 0 ). The non-local ciBK equation can be written as [26] ∂S(r, η) ∂η = d 2 r 1ᾱ s 2π r 2 r 2 1 r 2 2 Θ(η − δ r,r1,r2 ) S (r 1 , η − δ r1,r )S(r 2 , η − δ r2,r ) −S(r, η) r 2 + δ r 2 ,r S (r 1 , η) S (r , η)S(r 2 , η) −S(r 2 , η) where the rapidity shifts δ r,r1,r2 and δ r1,r are defined as δ r,r1,r2 = max 0, ln and δ r1,r = max 0, ln and similarly for δ r2,r . Note that on the r.h.s of Eq.(31) there are only two NLO terms explicit display. All the other NLO terms are collectively denoted as "regular". The roles of the rapidity shift δ r,r1,r2 as well as the step function in Eq.(31) are to introduce constraints on the soft-to-hard evolution [26]. Now let's turn to solve the non-local ciBK equation analytically in the saturation region. As we know from the previous integrals in Eqs. (6) and (19), that the integral in Eq.(31) is governed by the region either 1/Q s |r 1 | |r| and |r 2 | ∼ |r|, or 1/Q s |r 2 | |r| and |r 1 | ∼ |r|, whereQ s takes the role of saturation scaleQ We choose to work in the former region, where the S-matrix is negligibly small, thus we can neglect the quadratic term in Eq. (31). In addition, ln(r 2 2 /r 2 ) ∼ 0 leads to the second term on the r.h.s of Eq.

Solution with the fixed coupling constant
In the fixed coupling case, we know that the QCD coupling can be viewed as a constant. Thus, theᾱ s in Eq.(35) can be factorized out of the integral, whose solution is [26] S(r, η) = exp −cᾱ If one compares Eq.(37) with Eq.(9), it is easy to find that both of them have quadratic rapidity dependence except for different coefficients in the exponent, although they are in different rapidity representation. It is known that the value ofc is smaller than c [26], so the S-matrix in the η-representation is larger than the one in the Y -representation for the same value of rapidity.

solution with the running coupling constant
In the running coupling case, the QCD coupling is a function of the smallest dipole size. The coupling constant in Eq.(35) cannot be factorized out of the integral, the evolution equation (35) becomes whose solution isS where the NLO saturation momentum is used Note that the Eq.(39) has similar rapidity dependence as the Eqs. (20) and (29) except for different coefficients in the exponent, although they are in different rapidity representation. This result means that the solution of the non-local ciBK equation also complies in the law what we found in Ref. [34]. However, one shall find in the next subsection that the solution of the evolution equation including the corrections of the sub-leading double logarithms violates the law mentioned above.

B. The Sudakov suppressed BK equation and its analytic solution
It has been shown in Ref. [27] that there are significant corrections coming from the regime beyond the strong ordering region, where the sub-leading double logarithms are induced due to the incomplete cancellation between the real corrections and virtual corrections. These double logarithms have typical Sudakov feature and can be resummed into an exponential type resulting in a Sudakov suppressed Balitsky-Kovchegov equation [27]: where the Sudakov suppressed evolution kernel is

Solution with the fixed coupling constant
In the fixed coupling case, one can set the QCD coupling in Eq.(42) to be as a constant. Thus, one can factorize it out of the integral during the analytic solving Eq.(41) in the saturation region. As it was done in the previous section, we need to use the saturation condition which indicates that the S-matrix is very small, the non-linear term on the r.h.s of Eq.
where the upper and lower bounds are determined by the same way as the ones in the previous section. The factor 2 on the r.h.s of Eq.(43) comes from the symmetry of the two integral regions, see Fig.1. We perform the integral over r 1 and Y in Eq.(43), and get its solution as where we have used ln r 2Q2 s =cᾱ s (η − η 0 ) withc to be as a constant. If one compares the solution of the SSBK equation in η (44) with the solution of the LO BK equation in Y (9), it is easy to find that the S-matrix in the Sudakov case is larger than the LO one, since the exponential factor in the second term of Eq.(9) on the r.h.s is almost twice as large as the one in Eq.(44). In other words, the scattering amplitude T in the Sudakov case becomes smaller than the LO one, which indicates that the evolution speed of the dipole amplitude is slowed down by the Sudakov effect as compared to the LO one.

Solution with the running coupling constant
In the running coupling case, one cannot factorize the QCD coupling out of the integral in Eq.(41), since the QCD coupling could be a function of integral variable. As it was done in the fixed coupling case, the Eq.(41) is solved in the saturation region, therefore we can neglect the non-linear term and just keep the linear term in S, the Eq.(41) simplifies to Note that the kernel in Eq.(45) is simplified by the fact that the integral is governed by the region either from the transverse coordinate of the emitted gluon approaching to quark leg of the parent dipole, 1/Q s |r 1 | |r| and |r 2 | ∼ |r|, or the transverse coordinate of the emitted gluon approaching to antiquark leg of the parent dipole, 1/Q s |r 2 | |r| and |r 1 | ∼ |r|, see Fig.1. The factor 2 on the r.h.s of Eq.(45) accounts for the symmetry of the integral regions just mentioned. We choose to work in the region, 1/Q s |r 1 | |r| and |r 2 | ∼ |r|, the evolution kernel becomes K SS (r, r 1 , r 2 ) =ᾱ s (r 2 min ) 2π Performing the integral over variables r 1 and η in Eq.(45), one can get its solution as where the saturation momentum in the NLO case is used, see Eq. (40). The solution in Eq.(47) deserves several important comments which are as follows: • If one compares the running coupling solution of the SSBK Eq.(47) with the fixed coupling solution of the SSBK Eq.(44), one can see that the rapidity dependence of the dominant term in the exponent changes from the quadratic rapidity dependence Eq.(44) to the rapidity raised to power of 3/2 dependence Eq.(47), rather than linear dependence. This result does not coincide with the law what we found in Ref. [34] in which the rapidity dependence of the solutions of all the NLO evolution equations have linear dependence once the SDRCP is applied. This outcome indicates that the sub-leading double logarithms compensate part of decrease made by the running coupling effect.
• One can see that the solutions of the fixed coupling Eq.(44) and running coupling Eq.(47) SSBK equations have two terms due to the kernel having two terms, see Eq. (42). The second terms in Eqs.(44) and (47) are similar as the respective ones in the non-local ciBK cases except an additional factor 1/2 difference in the exponents, which implies that the first terms in Eqs.(44) and (47) could come from the incomplete cancellation of the real and virtual corrections.
• Interestingly, we find that the dominant term (the first term on the r.h.s of Eq.(47)) of the solution of the SSBK equation with SDRCP exp(O(η 3/2 )) has similar rapidity dependence as the one obtained by solving the full NLO BK equation in Y with the PDRCP exp(O(Y 3/2 )).
To conclude, one can see that the sub-leading double logarithms resulting from the incomplete cancellation between the real and virtual corrections beyond the strong time-ordering region make significant change to the rapidity dependence of the dipole amplitude not only in fixed coupling case but also in the running coupling case.

IV. NUMERICAL SOLUTIONS OF THE SUDAKOV SUPPRESSED BALITSKY-KOVCHEGOV EQUATION
In this section, we shall use the numerical method to solve the SSBK equation in order to test their analytic solutions obtained in the above section, especially focusing on the comparison of the numerical results located in the saturation region with the analytic solutions. As one knows that the dipole evolution equations are a set of complicated integro-differential equations, the Runge-Kutta method is needed to solve them on the lattice. The integrals in these equations are performed by the adaptive integration routines. In addition, the interpolation is needed during the numerical calculations, since some data points not locating on the lattice should be estimated. So, the cubic spline interpolation method is employed in this study. To simplify the computation, we employ the translational invariant approximation and assume the S-matrix independent of the impact parameter of the collisions, S = S(|r|, Y ). In terms of the above discussion, we choose to use the GNU Scientific Library (GSL) to perform the numerical computation, since the GSL includes almost all the functions required by the numerical solution to the evolution equations.
In order to solve the SSBK evolution equations, the McLerran-Venugopalan (MV) model is used to be as the initial condition [41], where we setQ 2 s0 = 0.15 GeV 2 at η = 0, γ = 1, and Λ = 0.2 GeV for simplicity, but theQ 2 s0 and γ shall be free parameters when they use to fit the HERA data in the next section. The one-loop running coupling with N f = 3 and N c = 3, Eq. (11), is used to be as the QCD coupling in this numerical simulation. In order to regularize the infrared behavior, the coupling value is freezed to α s (r fr ) = 0.75 when r > r fr .
The left-hand panel of Fig.2 gives the solutions of the non-local ciBK and SSBK equations as a function of the dipole size in η-representation for 4 different rapidities. Note that we plot a zooming in diagram to clearly show the numerical results in the saturation region. One can see that the values of the SSBK dipole amplitudes are larger than the non-local ciBK ones for each corresponding rapidities in the saturation region. This outcome is consistent with the analytic findings, Eqs.(47) and (39) in Sec.III, in which the linear rapidity dependence in the exponent of the S-matrix in non-local ciBK case is replaced by the rapidity raised to power of 3/2 dependence due to the contribution of the sub-leading double logarithms. The righthand panel of Fig.2 shows the solutions of the SSBK equation with the fixed coupling (ᾱ s = 0.3) and running coupling constant for 4 different rapidities. The zooming in diagram is used to assist in a clear view of the numerical results in the saturation region. We can see that the respective dipole amplitude with the the running coupling constant is smaller than the corresponding one with fixed coupling constant. This numerical result agrees with the analytic results, Eqs.(44) and (47), in which the quadratic rapidity dependence in the exponent of the S-matrix is replaced by the linear rapidity dependence once the SDRCP is used. The result also indicates that the running coupling effect takes a significant role in suppression of the evolution of the dipole amplitude in both Y -representation and η-representation. In this section, we come to describe the HERA data [42] for the inclusive DIS cross-section with the SSBK equation. The actual quantity we shall fit is the reduced γ * p cross-section which can be expressed in terms of the transverse, σ γ * p T , and longitudinal, σ γ * p L , cross-sections: with y = Q 2 /(sx) to be as the inelasticity variable and s the squared center of mass collision energy. The transverse and longitudinal cross-sections in Eq.(49) can be written as [7,33,38] where |ψ (f ) T,L | 2 is the squared light-cone wave function representing the probability for a virtual photon splitting into a quark-antiquark pair with flavor f , and can be written as with K 0 and K 1 to be as the modified Bessel function of the second kind. TheQ f in Eqs. (51) and (52) is defined asQ 2 f = z(1 − z)Q 2 + m 2 f with m f to be as the quark mass. Note that we only use three light quarks with m u,d,s = 140 MeV in our fit.
The key ingredient in Eq.(50) is the dipole cross-section which includes the most important information about the scattering between dipole and target. Here, the S-matrix is calculated by numerical solving the SSBK equation, and the σ 0 is viewed as a free parameter whose value is determined by fitting to the HERA data. We would like to point out that we use the one loop QCD couplingᾱ s (r 2 ) =   The experimental data points are the combined measurement from H1 and ZEUS collaborations [42]. We only plot the reduced cross-section for some typical values of Q 2 . The others have the same good quality as well. to fit the HERA data in term of the rather successful running coupling scheme in Refs. [8,9]. We take Λ = 0.2 GeV and treat the C 2 as a free parameter in our fit. In summary, we totally have four free parameters, σ 0 ,Q 2 s0 , γ and C 2 . The first one comes from the dipole cross-section, Eq.(53), the second and third parameters are from the initial condition, Eq.(48), and the final parameter originates from the QCD running coupling, Eq.(54).
In the fit, we consider the combined data from HERA for the reduced cross-section in the kinematical range x ≤ 0.01 and 0.045 GeV 2 < Q 2 < 50 GeV 2 [42]. The low limit on Q 2 is chosen low enough to justify the use of the BK dynamics rather than DGLAP evolution, and the upper limit large enough to include a large amount of perturbative data points. By using the criteria, we have 252 data points in our fit.
With the set up described above, we obtain rather good fits to the HERA data on the reduced crosssection. In Table I, we show the values of the free parameters. The reasonable value of χ 2 /d.o.f indicates that the SSBK equation with the SDRCP gives a rather successful description of the data. In Fig.3, we show the reduced cross-section as a function of the Bjorken variable x. Note that the Fig.3 only plots the reduced cross-section for some typical values of Q 2 , we have check that the others have the same good quality description as well. By comparison the data points (solid points in Fig.3) with the theoretical calculations (triangles in Fig.3), one can see that the values of the reduced cross-section calculated by the SSBK equation are in agreement with the HERA data points. This outcome is satisfying the theoretical expectation [26]. To show the rapidity dependence of the saturation momentum, we extract the value of the Q s from the fit to the data via the definition T (r = 2/Q s , η) = 1/2. In Fig.4, we demonstrate the saturation momentum as a function of x. The Q s is shown on top of the data points that we use in the fit.