Extended collinearly-improved Balitsky-Kovchegov evolution equation in target rapidity

An extended collinearly-improved Balitsky-Kovchegov evolution equation in the target rapidity representation is derived by including the running coupling corrections during the expansion of the"real"$S$-matrix. We find that the running coupling brings important corrections to the evolution equation, as one can see that there are extra contributions to the evolution kernel once the running coupling is included. To identify the significance of the corrections, we numerically solve the evolution equation with and without the running coupling contributions during the $S$-matrix expansion. The numerical results show that the scattering amplitude is largely suppressed by the running coupling corrections, which indicate that one needs to consider the running coupling contributions during the derivation of the non-linear evolution equation in the target rapidity representation.


I. INTRODUCTION
The Color Glass Condensate (CGC) effective theory has been approved to be as a powerful theory to describe the strong interactions associated with high energy and density environments. The leading order (LO) CGC calculations which refer to the derivation of the non-linear Balitsky-JIMWLK 1 [1][2][3][4][5] equation and its mean field version known as the Balitsky-Kovchegov (BK) equation [1,6], have been able to qualitatively describe many phenomenological results, such as the reduced cross-section in deep inelastic scattering (DIS) [7,8], and single and double inclusive particle production in high energy heavy ion collisions [9][10][11][12][13]. However, it has been found that the LO CGC theory is insufficient for direct applications to the phenomenology, since the evolution speed of the scattering amplitude resulting from the LO BK equation is too fast to quantitatively describe the experimental data [14,15].
There are tremendous developments in the calculations of the next-to-leading order (NLO) corrections to the JIMWLK and BK equations in the literature over the past fifteen years [16][17][18][19][20][21][22]. The pioneer work towards the NLO corrections to the BK equation was performed by including the running coupling effect in Refs. [16] and [17]. A running coupling Balitsky-Kovchegov (rcBK) equation was obtained. It was shown that the growth of the dipole-hadron scattering amplitude resulting from the rcBK equation is significantly slowed down as compared to the LO one. Furthermore, it was found that the rcBK equation gives a rather successful description of the HERA data [14]. However, the running coupling effect is not the only large higher order corrections to the LO BK or JIMWLK equations. Except the running coupling, the authors in Ref. [18] derived the full NLO BK evolution equation by including the quark loops, gluon loops, as well as the tree gluon diagrams with quadratic and cubic non-linearities, they obtained other contributions which are enhanced by double transverse logarithms. Unfortunately, the numerical study of the full NLO BK equation found that the equation is unstable, since the scattering amplitude can decrease as rapidity increasing and can even turn to negative value for small dipoles [23,24]. The instability was traced back to the radiative corrections enhanced by double transverse logarithms.
To cure the instability problem, one has to resum the radiative corrections enhanced by double transverse logarithms to all orders. Two methods related to this specific issue were proposed [25,26], which use different recipes to impose a kinematical constraint to the successive gluon emissions during the rapidity evolution, (i) the kinematical constrain introduced in Ref. [25] leads to an evolution equation which is similar to the LO BK equation, but is non-local in rapidity Y ; (ii) the resummation of the leading double logarithms has been performed in the evolution kernel in Ref. [26] resulting in a collinearly-improved Balitsky-Kovchegov (ciBK) equation with a modified kernel, which is still local in rapidity Y . It has been shown that these two methods are equivalent in the resummation of the leading double transverse logarithms. It is known that both of them lead to stable evolution equation and give good fits to the small-x HERA data [27][28][29]. However, the authors in Ref. [30] found there are some inconsistencies with the original analysis in Refs [25,26]. They found that the instability of the full NLO BK is caused by the wrong choice of the rapidity variable which plays the role of the evolution time. Moreover, it has been found that the growth of the saturation exponent with the running coupling constant in the asymptotic region of rapidity has a strong scheme dependence on the resummation method, which should not occur in practice.
It is improtant to point out that the rapidity generally used in all above mentioned BK equations is that of the projectile rapidity Y . In terms of the previous experience with the NLO BFKL 2 equation [31][32][33], where a similar issues were met and eventually cured, the evolution variable should be the rapidity of the hadronic target, η, rather than the rapidity of the projectile. Inspired by the experience on handling instability problem in NLO BFKL equation, a novel method was proposed to derive the collinearly-improved BK equation in η-representation (ciBK-η) from the corresponding evolution equation in Y -representation by the change of variable η = Y − ρ, where ρ is defined as ρ = ln(Q 2 /Q 2 0 ) with Q 2 and Q 2 0 to be the hard scale of the projectile and soft scale of the target, respectively [30]. As a consequence, the evolution variable in the ciBK-η equation is the physical rapidity η = ln(1/x), not the rapidity of the dipole projectile, and the ciBK-η equation shows a little scheme dependence on the resummation prescriptions. Although the ciBK-η equation has a significant advantage over the relevant one in Y -representation, it has recently been found that the evolution equation in η or Y -representation gives a very similar description of the HERA data [34]. The reason for this unexpected result could be attributed to several factors: (i) only the dominant part of the ciBK-η equation, which is called "canonical" Balitsky-Kovchegov equation (caBK-η), was used to study the data in order to avoid the cumbersomely numerical calculations, thus a set of NLO terms of order O(α 2 s ) were abandoned; (ii) the LO BK approximation was used in the replacement of the first derivative in the expansions of the "real" S-matrix 3 when the ciBK equation in Y -representation was transformed to η-representation [30], which leads to lose precision at the level of NLO accuracy.
In this paper, we shall derive the next-to-leading order BK equation in η-representation by including the running coupling corrections during the expansions of the "real" S-matrix. To see the significance of the corrections of the running coupling, we firstly recall the derivation of the ciBK-η equation with the LO BK approximation in the replacement of the first derivative in the expansion of the "real" S-matrix. Second, we derive an extended collinearly-improved Balitsky-Kovchegov equation in η-representation by emphasizing the running coupling corrections in the expansion of the "real" S-matrix. We obtain an extended "canonical" Balitsky-Kovchegov equation (exBK-η) whose evolution kernel is modified by the running coupling corrections as compared to the caBK-η equation. The exBK-η equation is analytically solved in saturation region. Its analytic solution shows that the exponent of the S-matrix has a linear dependence on rapidity instead of a quadratic rapidity dependence in the caBK-η case, which obeys a similar law as results in Y -representation that the evolution speed of the dipole amplitude is suppressed by the running coupling corrections. We compare our extended collinearly-improved Balitsky-Kovchegov equation in η with its original form to see how big difference is. It is easy to find that there are eight extra terms resulting from the running coupling corrections, which indicate that the precision of the expansion of the S-matrix has a significant impact on the evolution equation.
Finally, we numerically solve the evolution equations in η-representation to test the analytic outcomes mentioned above. The numerical results confirm our analytic findings. The saturation exponentsλ are extracted from the solutions of the LO BK, caBK-η and exBK-η equations. As expected, the running coupling effect has a large influence on the evolution speed of the front, which largely suppresses the rapidity evolution of the dipole amplitude.

II. LEADING ORDER, AND RUNNING COUPLING BK EQUATIONS IN Y -REPRESENTATION
In order to collect the basic elements of the dipole evolution equations, we give a brief recall of the LO BK and rcBK equations in Y -representation. These two equations shall be used in the Taylor expansion of the "real" S-matrix in the derivation of the collinearly-improved BK equation in η-representation in the next section.

A. Leading order BK equation
We consider the high energy scattering between a dipole which is consisted of a quark-antiquark pair moving towards the positive direction of the longitudinal axis with momentum (p + , p − , p), and a hadronic target moving along the negative direction with momentum (p + 0 , p − 0 , p 0 ). The scattering is treated in the eikonal approximation, thus the transverse coordinates of the quark (x) and the antiquark (y) are not modified by the collision. One can write the dipole scattering matrix as a correlator of two Wilson lines [18] where the · · · Y means the average over target gluon field configurations at Y . Here, the Y is the rapidity difference between the dipole and hadronic target with the center of mass energy squared s = 2p + p − 0 and the typical momentum of the target Q 0 . The U in Eq.(1) is the time ordered Wilson line with A + (x − , x) as the gluon field of the hadronic target.
In the mean field approximation, the rapidity evolution of the S-matrix satisfies the BK equation [1,6] withᾱ s = α s N c /π, and x, y, z as the transverse coordinates of the quark, antiquark, and emitted gluon, respectively. In large N c limit, the Eq.(4) depicts the parent dipole (x, y) evolved into two daughter dipoles (x, z) and (z, y) as the rapidity increasing. The first term in the right hand side of Eq.(4) is called as "real" term which describes the two daughter dipoles scattering with the target simultaneously, thus it is a non-linear term. The second term in the right hand side of Eq.(4) is referred as "virtual" term which depicts the survival probability of the original dipole at the time of scattering. Note that the BK equation resums only the leading logarithmic α s ln(1/x) corrections in the fixed coupling case, thus it is a LO evolution equation.
In the later section of this paper, we shall be interested in the limiting form of the S-matrix in saturation region where the dipole has very large size, such that r 2 Q 2 s 1. Here, the Q s is the typical transverse momentum of the saturated gluons and is rapidly increasing with Y . For later comparison, we analytically solve the Eq.(4) in the following. In saturation regime, the S-matrix approaches the black-disk limit, S → 0. Therefore, the quadratic term in S in Eq.(4) can be neglected. Moreover, the saturation condition implies that the two daughter dipoles are larger than the typical transverse size r s ∼ 1/Q s . Therefore, the Eq.(4) can reduce to with r = |x − y|, r 1 = |x − z|, r 2 = |z − y| be the transverse size of the parent dipole, and two daughter dipoles, respectively. Further, one can find that the integral in Eq.(5) is governed by the case if one of the daughter dipoles is much smaller than the parent dipole, r 1 r and r 2 ∼ r, or r 2 r and r 1 ∼ r [35,36]. We select to work in the case r 2 ∼ r, the Eq.(5) can be rewritten as where the right hand side includes a factor of 2 to take into account that the smaller dipole with size r 1 can be any one of the two daughter dipoles. Now we carry out the integrals in Eq. (6), which yield the analytic solution of the LO BK equation as [35][36][37][38], where we have used ln[r 2 Q 2 s (Y )] λᾱ s (Y − Y 0 ) with Y 0 the rapidity scale at which Q 2 s (Y 0 ) = 1/r 2 . The Eq.(7) is usually called as Levin-Tuchin formula, since it was firstly derived by them in Ref. [37]. From the above derivation, we know that the exponent in Eq. (7) is known only to leading double logarithmic accuracy, which means that the sub-leading terms are not under control. In addition, the exponent of the S-matrix has a quadratic dependence on rapidity, which renders the S-matrix too small at large rapidities, in the other words, the rapidity evolution speed of the dipole amplitude N is too fast, with N = 1 − S. It has been recognized that the aforemention drawbacks are the reasons why the LO BK equation is insufficient to describe the experimental data at HERA [14,15,39]. Therefore, one has to include the NLO corrections to the LO BK evolution equation, such as the running coupling effect which can suppress the evolution speed of the dipole amplitude by modifying the evolution kernel of the BK equation [16,17].

B. Running coupling BK equation
It is known that the LO BK equation discussed above considers only the resummation of leading logarithmic α s ln(1/x) corrections with a fixed coupling. Beyond the leading logarithmic approximation, there was a significant progress in the BK evolution equation via the resummation of α s N f to all orders, which is called as the running coupling corrections. The running coupling Balitsky-Kovchegov (rcBK) equation have been derived independently by Balitsky in Ref. [16] and Kovchegov and Weigert in Ref. [17]. This two groups obtained an analogous structure of the rcBK equation but with different evolution kernels. In this study, we shall use the Balitsky version of kernel, since it is favored by the HERA data [14]. Note that we don't plan to give a details of the derivation of the rcBK equation, it is out of interesting of this paper. The rcBK equation reads where K rc (r, r 1 , r 2 ) is the running coupling evolution kernel [16] K rc (r, r 1 , with r min = min{r, r 1 , r 2 }. There are several running coupling prescriptions in the literature [16,27,[40][41][42].
At the very beginning, the argument of the running coupling α s in the rcBK equation was interpreted to the transverse size of the parent dipole. Recently, it was found that the size of the smallest dipole is a proper argument of the coupling, and the smallest running coupling prescription is favored by the HERA data than others [27]. So, we shall use the smallest dipole running coupling prescription in this study. In addition, the running coupling at one loop accuracy is used . Now, let us turn to analytically solve the rcBK equation. As it was done in the previous section, we solve the rcBK equation in saturation region where the dipoles have very large size, r, r 1 , r 2 ≥ 1/Q s . The integral over the r 1 in Eq. (8) is governed by the case if one of the daughter dipoles is much smaller than the parent one, while the rest of dipole has a similar size as the parent dipole, that is r 1 r and r 2 ∼ r, or r 2 r and r 1 ∼ r. We choose to work in the first case, the Eq.(8) reduces to where the right hand size includes a factor of 2 due to the fact that the smaller dipole can come from either of the two daughter dipoles. In saturation region, the S-matrix approaches the black-disk limit, one has S → 0, thus the quadratic term in S in Eq.(11) can be discarded. The Eq.(11) simplifies to Substituting Eq.(10) into Eq. (12), and performing the integrals over r 1 and Y , we can get the analytic solution of the rcBK equation as [36,38] where the NLO saturation momentum is used If one compares the solution of the rcBK equation, Eq. (13), with the solution of the LO BK equation, Eq. (7). It is easy to find that the quadratic rapidity dependence in the exponent of the S-matrix is replaced by the linear rapidity dependence due to the running coupling corrections. This outcome indicates that the evolution speed of the dipole amplitude is slowed down by the running coupling corrections. This finding is consistent with the theoretical expectations [16,43]. Furthermore, the phenomenological applications of the rcBK equation at HERA energies show that the rcBK equation gives a more reasonable description of the experimental data than the LO BK equation [14,15].

III. COLLINEARLY-IMPROVED BK EQUATION IN η-REPRESENTATION
In previous section, all the dipole evolution equations are studied in the Y -representation. In this section, we shall discuss the dipole evolution equations in η-representation due to two key reasons. On one hand, a recent study in Ref. [44] realized that the rapidity η = ln(1/x) of the hadronic target is the physical rapidity used in the DIS experiments at HERA rather than the projectile rapidity Y . On the other hand, it was found that the reason for the instability of the full NLO BK equation is a consequence of the wrong choice of the "evolution time", this refers to the choice of the rapidity variable [30]. However, if one simply transforms the ciBK-Y equation to η-representation by change of variable η = Y − ρ, the instability problem is still existence due to NLO corrections enhanced by double collinear logarithms, but not as severe as for the corresponding issue in Y . Based on the previous experience with the full NLO BK equation in Y , where a similar problem was identified and eventually cured by enforcing the time-ordering constrains on the successive gluon emissions, one can know that the ordering of the successive emissions in longitudinal momentum should be enforced in order to solve instability problem in η-representation. By doing this, the ciBK-η was obtained in Ref. [30], which can directly apply to phenomenology and supposes to give a better description of the HERA data than the other evolution equations in the literature. However, a very recent study in Ref. [34] showed that three different evolution equations (kinematical constraint BK (kcBK) [25], ciBK [27], and caBK-η [30] equations) result in a very similar description of the HERA data. Note that the first two equations are presented in Y -representation, while the last equation is given in η-representation. It is easy to understand that the kcBK and ciBK equations give a equally good depiction of the data, since they are equivalent in the sense that the resummation of the leading double transverse logarithms is concerned. The reason why the caBK-η does not give a superior description of the data, could possibly be attributed to the insufficient accuracy used in the expansion of the "real" S-matrix when the caBK-η equation was derived. In this section, we shall derive the collinearly-improved BK equation in η at the level of running coupling when the "real" S-matrix is expanded. An extended collinearly-improved BK equation in η is obtained, which has the same structure as the ciBK-η equation but with a running coupling modified kernel. In the next section, we use the numerical method to solve the caBK-η and exBK-η equations. The numerical results show that the running coupling corrections play a significant role in the suppression of the evolution of the dipole amplitude.

A. Collinearly-improved BK equation in η: LO BK approximation in expansion of S-matrix
To obtain the collinearly-improved BK equation in η, we follow the same strategy used in Ref. [30] where the change of variable is employed to transform the ciBK equation from the Y -representation to the ηrepresentation.
We start with the full NLO BK equation in Y -representation [16] ∂S with and From Eq. (15), one can see that the full NLO BK equation has two main changes in the structure as compared to the LO BK equation. The first term in the right hand side receives corrections from quark loops (K q ) and gluon loop (K g ). The last two terms in the right hand side refer to partonic fluctuations involving two additional partons except the original parent partons (quark and antiquark). In large N c limit, they are corresponding to two consecutive emissions, the original parent dipole (x, y) emits a gluon at transverse coordinate u, which is equivalent to two daughter dipoles (x, u) and (u, y). Then the dipole (u, y) emits a gluon at transverse coordinate z, which yields the dipoles (u, z) and (z, y).

The "canonical" BK equation
In order to transform the Eq.(15) from Y -representation into η-representation, we need to change variables in terms of Now, we can rewrite the S-matrices in η representation as S xy (Y ) = S xy (η + ρ) ≡S xy (η), where we have used and ρ = ln When one works at NLO inᾱ s , one can expand out the "real" S-matrices (S xz and S zy ) in Taylor series, since the rapidity shift in the argument of S-matrices is typically much smaller than η itself. The expansions of the S xz and S zy can be expressed as and To the order of interesting, we need only to keep the first non-trivial term in the above expansions, since each ∂S/∂η is formally suppressed by a power ofᾱ s . For the derivative terms in Eqs. (27) and (28), in this subsection we use the LO BK equation to evaluate them as what was done in Ref. [30], and (30) From the above equations, one can see that the rapidity shift in the argument of S-matrices is equivalent to adding a term of order O(ᾱ s ). We would like to point out that the LO BK equation is a rough approximation to estimate the derivative terms in Eqs. (27) and (28). Actually, the LO BK equation is insufficient due to its lower precision. To reach the interested order of accuracy, one should use at least the level of rcBK equation to evaluate the derivative terms, which shall study in the next subsection. Now, substituting Eqs. (29) and (30) into Eq. (15), one can get a semi-finished collinearly-improved BK equation 4 with The third term in the right hand side of Eq.(31) is resulting from the expansion of the "real" S-matrices in terms of rapidity shift with the LO BK approximation in the evaluating of the derivative term. Note that the rapidity shift is neglected in the all NLO terms when Eq.(31) is derived. All the terms of order O(ᾱ 2 s ) in Eq.(15) are simply replaced like S xz (Y ) →S xz (η), since the rapidity shift take a contribution of order O(ᾱ s ) which renders all the NLO terms of order O(ᾱ 3 s ). While we are only interested in the terms up to the order of O(ᾱ 2 s ), thus all the terms beyondᾱ 2 s are abandoned in this paper. Moreover, two mathematical tricks are used when Eq.(31) is derived, (i) the property that the LO term is invariant under x − z → z − y, is exploited to combine some terms; (ii) the integral variables in the third term in the right hand side of Eq.(31) is relabelled in terms of u ↔ z in order to keep consistence with the physics picture mentioned above.
The Eq.(31) is a NLO evolution equation in η-representation, which is a local equation in rapidity. By comparing Eq.(31) with Eq.(15), one can see that the difference between them is only by an extra term (resulting from the change of variable) in the third line in the right hand side of Eq.(31). Originally, one anticipates that the change of variable from Y to η can eliminate the instabilities occurred in Eq. (15). As expected that the instabilities caused by the violations of time-ordering (double anti-collinear logarithms) are disappeared in Eq.(31), since the time-ordering property is automatically guaranteed in the η evolution. Unfortunately, it has been shown that the change of variable triggers off another type of instabilities associated with double collinear logarithms [30].
To cure the instabilities mentioned above, it is known that the successive gluon emissions during the rapidity evolution have to be simultaneously ordered in lifetime and longitudinal momentum, and where the (p + , p − , p), (p + 0 , p − 0 , p 0 ), and (k + , k − , k) denote the light cone momenta of the projectile, target and emitted gluon, respectively. The first constraint is automatically satisfied as mentioned above. However, the second constraint may be violated when the radiated gluon is either too soft (k 2 Q 2 0 ) or too hard (k 2 Q 2 ). So, one has to put the constraint on the evolution equation. Before doing that, we need to rewrite the constraint, Eq.(34), in a proper form. As we know and Thus, the target rapidities can be expressed as and Using Eqs. (37) and (38), one can rewrite the constraint, Eq.(34), as Moreover, according to the lifetime constraint in Eq.(33), one can get Combining the two constraints, Eqs. (39) and (40), we obtain the final rapidity constraint as By using Eq.(35) and (36), the above equation can be written in another form as which is a proper form, and can be directly used. Now, we apply the constraint, Eq.(42), to the integral form of BK equation and get which turns to differential format as with the rapidity shift δ xyz as It has been checked that the Eq.(44) has little scheme dependence on the prescription of rapidity shift. In this study, we choose to work with the "canonical" one [30] S xz (η−δ xyz )S zy (η−δ xyz )−S xy (η−δ xyz ) −→S xz (η−δ xz;r )S zy (η−δ zy;r )−S xy (η), with δ xz;r = max 0, ln and δ zy;r = max 0, ln r 2 (z−y) 2 . (48) Substituting Eq.(46) into Eq. (44), one gets the caBK-η equation as [30] ∂S xy (η) ∂η =ᾱ which is a non-local evolution equation in rapidity η.
Adding the O(ᾱ 2 s ) pieces from Eq.(31) to Eq.(49) and subtracting the above O(ᾱ 2 s ) piece, one obtains the ciBK-η equation as [30] ∂S xy (η) ∂η =ᾱ with One can see that Eq.(51) does not include double collinear logarithms now. These double logarithms are included in the first line in the right hand side of Eq.(51) through the rapidity shift. In addition, the double anti-collinear logarithm term in kernel K g in the second line in the right hand side of Eq.(51) is canceled by a relative piece generated by the integral over u in the third line. All the unstable factors are under control in Eq.(51). So, it is a stable equation which can directly apply to the phenomenological studies. However, a very recent study showed that the caBK-η equation does not give a superior description of the HERA data than the kcBK and ciBK equations as the theoretical expectations [34]. The reason why the outcomes resulting from the caBK-η equation are not as desired, possibly comes from the insufficient accuracy of the expansion of the "real" S-matrices (Eqs. (29) and (30)) and the integral LO BK equation Eq.(43).

Analytic solution to the caBK-η equation in saturation region
Let us turn to analytically solve the caBK-η equation in saturation region. In this regime, one of the daughter is much smaller than the other one, while the larger daughter dipole has comparable size as the parent dipole. As we know that the non-locality is only important for the S-matrix which is associated with smaller dipole. Thus, the Eq.(49) can simplify to where a factor 2 is taken into account due to the fact that the smaller dipole can come from either of the two daughter dipoles, Q s is the saturation momentum which is associated withQ s as r 2 Q 2 s = (r 2Q2 s ) 1/(1+λ) [30]. For simplicity, we denote z as the size of the smaller dipole in Eq.(53).
To solve Eq.(53) in saturation region, we follow two rules to do the calculations, (i) the saturation condition requires the dipole size z larger than the typical size 1/Q s (lower integral bound in Eq.(53)), which leads to the "real" S-matrix S(z, η − ln r 2 z 2 ) is negligibly small;(ii) the integral over z become logarithmic when z is much smaller than r. Applying the two rules, Eq.(53) reduces to whose solution isS where we have assumed the saturation momentum in η-representation asQ 2 s = Q 2 0 exp(λη). By comparing Eq.(55) with Eq. (7), one can see that they have similar form, but the solution of the caBK-η equation has an extra suppression factor in the exponent, which leads to the evolution speed of the dipole amplitude is slowed down. The numerical solutions of these two equations shall be done in the next section, where the numerical calculations support the aforementioned analytic result.

B. Collinearly-improved BK equation in η: rcBK approximation in expansion of S-matrix
In the previous subsection, the first derivative terms in the Taylor expansions in Eqs. (27) and (28) are approximately replaced by the LO BK equation, which are insufficient. To achieve the interested order of accuracy, we shall derive the collinearly-improved BK equation in η by using the rcBK equation (8) to replace LO BK equation (4) in the expression of the first derivative terms.

Extended caBK-η equation
Using the rcBK equation, one can re-expand the S-matrices in Eqs. (27) and (28) as and with the running coupling evolution kernel where we use the shorthand notation α xz s = α s ((x − z) 2 ) and similarly for others. Note that Eq.(58) is another form of Eq.(9), expressed in terms of transverse coordinates of the dipoles. Substituting the Eqs.(56) and (57) into Eq.(15), we obtain a semi-finished local collinearly-improved BK euqation in η after some complicated algebra calculations (for the detailed derivation, see Appendix A), with By comparing Eq.(59) with Eq. (31), there are extra eight terms resulting from the running coupling corrections, which have a significant impact on the evolution speed of the dipole amplitude. However, these extra terms do not cure the unstable issue of the evolution equation. Based on the discussion after Eq.(31), we know that the Eq.(59) still has the instabilities caused by double collinear logarithms.
To curve the instability problem, the successive gluon emissions during the rapidity evolution must be ordered in lifetime and longitudinal momentum simultaneously. So, we need to put the constraints, Eqs. (33) and (34), on the successive gluon emissions as what we have done in the previous section, Note that originally, the integral form of the LO BK equation was used in the derivation of the ciBK-η equation [30]. In order to achieve the interested order of accuracy, we use the integral form of the rcBK equation instead of the LO BK equation as the start point to derive the ciBK-η equation. Performing derivatives over η in Eq.(61), one can get the exBK-η equation which has the same structure as the caBK-η equation (49), but with a running coupling modified kernel. In terms of the experience from Y -representation, we deduce that the rapidity evolution of the dipole amplitude is also suppressed by the modified kernel, which shall be approved by the numerical calculations in the next section.

Extended full collinearly improved BK equation at O(ᾱ 2 s )
Based on the discussion in previous subsection, one can extend the exBK-η equation to full collinearlyimproved BK equation in η by adding the NLO terms from Eq.(59). Before writing down the extended full collinearly-improved BK equation in η-representation, we need to identify the O(ᾱ 2 s ) piece in the right hand side of Eq.(62). First, we expand theS xz andS zy to linear order in the rapidity shift, and then use the rcBK equation to replace the derivative terms, andS zy (η − δ zy;r ) S zy (η) − δ zy;r with

Analytic solution to the exBK-η equation in saturation region
Let us move to analytically solve the exBK-η equation, Eq.(62), in saturation region. In this regime, we know that one of the two daughter dipoles has similar size as the parent dipole, but the size of the rest one is much smaller than the parent dipole. Moreover, it is known that the non-locality is only important for the S-matrix which is associated with small size. Therefore, we can reduce Eq.(62) to where the factor 2 accounts for the smaller size dipole coming from any one of the two daughter dipoles, and the smallest size of the dipoles (z) is used to be as the argument of the QCD coupling.
To solve Eq.(68) analytically in saturation region, the same strategy is used as what we have done in Sec.III A 3. The one loop running coupling, Eq.(10), is used in the calculations. The Eq.(68) becomes Performing the integral over z, we get whose solution is with the saturation momentum in NLO case, and By comparing Eq.(71) with Eq.(55), one can see that the quadratic rapidity dependence in the exponent of the S-matrix is replaced by the linear rapidity dependence once the running coupling corrections are taken into account. This change implies that the evolution speed of the dipole amplitude is suppressed by the NLO corrections.

IV. NUMERICAL ANALYSIS
To test the analytic results obtained in the above section, we shall numerically solve the evolution equations in this section. The Eqs.(49) and (62) are integro-differential equations which can be numerically straightforward solved on a lattice. To simplify the computation, we neglect impact parameter dependence of the dipole amplitude throughout this numerical calculations, which imply that the dipole amplitude does not depend on angle, N (r, Y ) = 1 − S(r, Y ) = 1 − S(|r|, Y ). Thus, we can view the evolution equations as a set of differential equations and solve them at discrete values of transverse separation. To be more specific, we discretize the dipole transverse size r into 800 points which are equally located in the logarithmic space between r min = 10 −8 GeV −1 and r max = 50GeV −1 . The GNU Scientific Library (GSL) is a good candidate to solve them, since the GSL contains almost all the routines required by our purpose, such as the Runge-Kutta method for solving differential equations, adaptive integral routines for performing numerical integrals, and the cubic spline interpolation codes for interpolating the data points not located on the lattice.
The initial condition for the evolution equations is parameterized at rapidity η = 0. We use the Golec-Biernat and Wusthoff (GBW) parametrization as the initial condition [45], with p = 4, and Q s0 = 0.362GeV [44]. For the strong coupling constant, we use the one-loop running coupling, Eq. (10), with N f = 3 and N c = 3. According to the performance of the running coupling in the fit to HERA data [27,29], we choose to use the smallest dipole running coupling prescription which means the argument of coupling is the smallest dipole among the parent and daughter dipoles, We freeze the running coupling α s (r fr ) = 0.75 when r > r fr in order to regularize the infrared behavior.
To show the impact of the running coupling corrections on the speed of the fronts, the saturation exponent is numerically calculatedλ where the saturation momentQ s (η) is determined by N (r = 1/Q s , η) = κ with κ to be a constant of order 1. The left-hand panel of Fig.1 gives the solutions of the LO BK and caBK-η equations for 4 different rapidities. By comparing the solutions for each respective rapidities, one can see that the values of the caBK-η dipole amplitude are smaller than the LO BK ones, which indicate that the NLO corrections enhanced by the double transverse logarithms suppress the evolution speed of the dipole amplitude. A zooming in diagram is provided to clearly show the numerical results in the saturation region. One can see that the evolution is also slowed down in the saturation region. This numerical outcome is consistent with the analytic results in Eqs. (7) and (55), where the analytic solutions of the LO BK and caBK-η have analogous expression, but with different coefficients in the exponent. The coefficient in the caBK-η case is smaller than the LO BK one, which leads to the caBK-η dipole amplitude smaller than LO BK one. The right-hand panel of Fig.1 shows the comparison of the solutions of LO BK and exBK-η equation for 4 different rapidities. We plot a inner zooming in diagram for a clear comparison between the LO BK and exBK-η dipole amplitudes in saturation region. One can see, from the zooming in diagram, that the respective solution of the exBK-η equation in each rapidity is much smaller than the LO BK one, which implies that the evolution speed of the dipole amplitude is significantly suppressed. The suppression is much more than the caBK-η case. This outcome confirms the analytic result in Sec.III B 3, where the dipole amplitude is suppressed by the running coupling corrections. The left-hand panel of Fig.2 gives the comparisons of the solutions of the caBK-η and exBK-η equations for 4 different rapidities. We can see that the respective dipole amplitudes resulting from the exBK-η equation are smaller than the ones from caBK-η equation. Especially, it clearly show the suppression in the saturation region from the inner zooming in diagram. This outcome supports the analytic findings in Eq.(71), where the dipole amplitude is further suppressed by the running coupling corrections on top of the double logarithm resummation. Finally, we present the η dependence of the saturation exponent predicted by the LO BK, caBK-η, and exBK-η in the right-hand panel of Fig.2. As expected, theλ resulting from exBK-η equation is the smallest one among them, which is consistent with the theoretical expectations, since the exBK-η equation includes the running coupling corrections on top of the collinear resummations.