High energy asymptotic behavior of the $S$-matrix in the saturation region with the smallest dipole running coupling prescription

We present results from analytic solutions to the running coupling, full next-to-leading order, and collinearly improved next-to-leading order Balitsky-Kovchegov equations in the saturation region with the smallest dipole size QCD running coupling prescription. The analytic results of the $S$-matrix of the latter two equations show that the $\exp(-\mathcal{O}(Y^{3/2}))$ rapidity dependence of the solutions are replaced by $\exp(-\mathcal{O}(Y))$ dependence once the running coupling prescription is switched from parent dipole to the smallest dipole prescription, which indicate that the $S$-matrix has a strong dependence on the choice of running coupling prescription. We compute the numerical solutions of these Balitsky-Kovchegov equations with the smallest and parent dipole running coupling prescriptions, the numerical results confirm the analytic outcomes. The rare fluctuations of the $S$-matrix on top of next-to-leading order corrections are also studied under the smallest dipole running coupling prescription in the center of mass frame. It shows that the rare fluctuations are strongly suppressed and less important in the smallest dipole running coupling prescription case as compared to the parent dipole running coupling prescription case.


I. INTRODUCTION
Perturbative QCD predicts that the gluon density in high energy hadronic collisions is rapid growth with decreasing Bjorken x (or increasing energy) due to each emitted gluon itself as a source of further emission, which leads to fill up the available phase space for gluon radiation, and forms a state of saturated gluons called the Color Glass Condensate (CGC) [1]. The rapidity evolution of the CGC is described by the Balitsky-JIMWLK 1 [2][3][4][5][6] infinite hierarchy of renormalization group equations for the multi-point correlators of Wilson lines which describe a high energy quark or a gluon traveling through a dense target color field. In the mean field approximation, the hierarchy Balitsky-JIMWLK equations decouple and result in a single non-linear integro-differential equation known as the Balitsky-Kovchegov (BK) equation [2,7]. The BK equation is a closed equation, which significantly simplifies the direct applications to phenomenologies, such as proton structure function in deep inelastic scattering at HERA, and particle production in heavy ion collisions at LHC. However, the BK equation resums only large logarithms ∼ α s ln(1/x) to all orders with a fixed coupling constant α s , thus it is a leading order (LO) equation. It has been found that although the models (i.e. Iancu-Itakura-Munier model [9]) inspired by the LO BK equation can give a successfully qualitative fits to the small-x HERA data [8,[10][11][12], there is some tensions when one uses the LO BK equation to quantitatively compare with the experimental data, since the higher order corrections can be very large [13][14][15][16][17][18].
Over the past decade, the understanding of the non-linear evolution in QCD beyond leading order accuracy has received important developments, which refers to the first calculations of the next-to-leading order (NLO) corrections to the Balitsky-JIMWLK and BK equations [19,20]. The NLO corrections to the BK equation were calculated by the resummation of the α s N f contributions to all orders with N f to be as the number of flavors, which allow one to estimate running coupling corrections to the evolution kernel and result in the running coupling Balitsky-Kovchegov (rcBK) equation. In the language of Feynman diagram, this type of NLO corrections refer to the quark loop contributions. Soon after the rcBK equation was derived [19,20], Balitsky and Chirilli in Ref. [21] found that the gluon loops and the tree gluon diagrams with quadratic and cubic nonlinearities have also significant contributions to the BK evolution equation. By combining the quark part and gluon part contributions, they obtained a full NLO BK evolution equation. However, the first numerical studies to the full NLO BK equation found that the solution strongly depends on the details of the initial condition, and the scattering amplitude decreases with increasing rapidity and can even turn to negative for small dipoles [16,22], which indicate that the full NLO BK equation is unstable. It was found that the instability comes from a large double transverse logarithm in the BK evolution kernel. Fortunately, the authors in Ref. [23] devised a novel method to resum those single and double collinear logarithms to all orders and got a stablized evolution equation known as collinearly improved BK equation.
The stablized full NLO BK equation was solved numerically and its numerical solution gives a rather good description to the small x HERA data [17,18,24]. However, this numerical process is very cumbersome to use in practice, due to the intricate programming for solving an integro-differential equation and time-costing for running the program. In order to facilitate the use of the BK equation in direct phenomenological application, a lot of efforts were devoted to establish an analytical dipole scattering amplitude in the literature [25][26][27][28][29][30]. Among them, we analytically solved the full NLO BK equation in the saturation region [29]. We find that the exp(−O(Y )) rapidity dependence of the scattering S-matrix in the case of running coupling corrections (including only quark loop contributions) is replaced by exp(−O(Y 3/2 )) in the full NLO case (with both quark and gluon loop contributions). In addition, the authors in Ref. [26] used another approach to independently solve the full NLO BK equation in the saturation region, and also got exp(−O(Y 3/2 )) rapidity dependence of the S-matrix. They established a piecewise dipole scattering amplitude based on the analytic solution. It was found that the piecewise amplitude gives a rather successful fits to the small x HERA data. We would like to note that all the calculations mentioned above are used the size of the parent dipole as the argument of the running coupling constant α s . For several years, as we know that the argument of the coupling constant α s in the NLO BK equation was interpreted to the size of the parent dipole, especially in the applications to phenomenology [14,15,31], although the authors in Ref. [21] has pointed out that the proper interpretation of the argument of the coupling constant is the size of the smallest dipole rather than the size of the parent dipole. Recently, It was found that a very good description of the small x HERA data was obtained by using the collinearly improved NLO BK equation with the smallest dipole running coupling (SDRC) prescription among several different prescriptions of QCD running coupling. The significance of the argument of the coupling constant α s was aroused [17,18].
In this paper, we shall solve analytically the NLO BK equations with the SDRC prescription in the saturation region. To see the significance of the prescription of the running coupling, we firstly recall the analytic solutions of the rcBK, full NLO BK, and collinearly improved NLO BK equations with the parent dipole running coupling (PDRC) prescription in the saturation region, and we shall use these solutions for the latter comparisons. Secondly, we analytically solve the rcBK, full NLO BK, and collinearly improved NLO BK equations again with emphasizing on the SDRC prescription, and we compare the solutions resulting from two different running coupling prescriptions to see how big difference is.
Interestingly, we get that the solutions of the rcBK equation are exactly same under the PDRC and SDRC two different prescriptions, which means that its solution is independent on the choice of the running coupling prescription. We find that the exp(−O(Y 3/2 )) rapidity dependence of the S-matrix (solution of the full NLO BK equation) obtained by our previous studies with the PDRC prescription in Ref. [29] is replaced by exp(−O(Y )) rapidity dependence once the running coupling prescription is switched to the SDRC prescription, which indicate that the SDRC prescription suppresses the evolution of the dipole scattering amplitude and renders the rapidity dependence of the exponent of the S-matrix keeping a linear dependence as the one in rcBK case, see Fig.1 for a diagrammatic depiction. With the SDRC prescription, we now get that all the solutions of the NLO BK (rcBK, full NLO BK and collinearly improved NLO BK) equations have the same rapidity dependence exp(−O(Y )).
In order to test these analytic outcomes, we numerically solve the above mentioned NLO BK equations with the focusing on physics in the saturation region. The numerical results support the analytic findings. In addition, we also study the rare fluctuations of the S-matrix on top of the full NLO corrections with the SDRC prescription. We find that the rare fluctuation effects take a negligible contribution to the suppression of the evolution of the dipole amplitude, which is unlike the one obtained with the PDRC prescription in our previous publication [32] where a factor of √ 2 suppression of the exponential factor of the S-matrix is occurred when the rare fluctuation effects are included.

II. SOLUTIONS TO THE EVOLUTION EQUATIONS WITH PARENT DIPOLE RUNNING COUPLING PRESCRIPTION
To introduce notations and explain the kinematics, we review the rc, full NLO and collinearly improved NLO BK equations and their solutions in the case of PDRC prescription. These solutions will be used for the latter comparisons with the results obtained under the SDRC prescription in the next section.

A. Leading order Balitsky-Kovchegov equation and its solution
We consider the high energy scattering between a dilute projectile, consisting of a quark-antiquark color dipole with quark leg at transverse coordinate x ⊥ and an antiquark leg at transverse coordinate y ⊥ , and a dense target which may be another dipole, a hadron or a nucleus. The evolution of the S-matrix in the fixed coupling case is described by the BK equation, which we write as [2,7] where the evolution kernel is given by withᾱ s = α s N c /π. Here we use the notation r = x ⊥ − y ⊥ as the transverse size of parent dipole and r 1 = x ⊥ − z ⊥ and r 2 = z ⊥ − y ⊥ as the transverse sizes of the two emitted daughter dipoles, respectively. The BK equation resums only the leading logarithmic α s ln(1/x) in the fixed coupling case, therefore it is a leading-order equation. We would like to note that the BK equation is a mean field version of the Balitsky-JIMWLK hierarchy [3][4][5][6] equations, in which higher order correlations are neglected. Due to the simple structure of the BK equation, one can analytically solve it in the saturation region, in which the parton density is so high that the dipole scattering amplitude approaches unit, N ∼ 1 (N = 1−S), thus S ∼ 0. So, we can neglect the quadratic term of the S-matrix in the BK equation in the saturation region, the Eq.(1) becomes where Q s is the saturation momentum. The Q s is an intrinsic momentum scale which provides a separation between dense and dilute parton system. The lower bound of the integral in Eq.(3) is set to 1/Q s since the saturation condition requires that the transverse dipole size should be larger than the typical transverse dipole size, r s ∼ 1/Q s . We would like to note that although there are few radiated daughter dipoles having transverse size larger than parent dipole, we still set the upper bound of the integral to r because the kernel has a rapid decay when the size of the daughter dipole is larger than the parent dipole, therefore the contribution from the region, r 1 > r, is negligible. It is not hard to find that the integral is governed by the region either from the radiated dipole closing to quark leg of the parent dipole, 1/Q s |r 1 | |r| and |r 2 | ∼ |r|, or the radiated dipole closing to antiquark leg of the parent dipole, 1/Q s |r 2 | |r| and |r 1 | ∼ |r|. In this study, we choose to work in the region |r 2 | ∼ |r|, the Eq.(3) becomes Note that the factor 2 on the right hand side of Eq.(4) comes from the symmetry of the two regions mentioned above. Now we can easily obtain the analytic solution of the LO BK equation in the saturation region by integrating over r 1 and Y [25,33], where we have used Q 2 Eq. (5) shows that in the saturation region the S-matrix has a quadratic rapidity dependence in its exponent. It has been found that this S-matrix cannot precisely describe the experimental data from HERA [15], since the evolution speed of the scattering amplitude N is too fast, in other words the S-matrix is too small. So, one needs to take into account the NLO corrections which can make the S-matrix becoming larger, or slow down the evolution speed of the dipole scattering amplitude.

B. Running coupling Balitsky-Kovchegov equation and its solution
The first NLO corrections to the LO BK equation were derived by Refs. [19,20], which consider the quark loop contributions to the BK evolution kernel and resum α s N f to all order. We usually call this modifications as running coupling corrections. The evolution equation including the quark loop corrections reads [14] where we divide the right hand side of Eq.(6) into two parts, 'running coupling part' and 'subtraction part'. The 'running coupling part', resumes all power of α s N f corrections to evolution kernel, which has the same structure as the LO BK equation but with a modified evolution kernel. The modified kernel has Balitsky type [20] K rcBal (r, r 1 , and Kovchegov-Weigert type [19] K rcKW (r, r 1 , with which depend on the decomposition scheme of the two parts, see Ref. [14] for the details of separation schemes. It is easy to see that the couplings in Eqs. (8) and (9) are a function of the transverse size of a dipole, which do not like the one in leading order case where the coupling is fixed. We choose the running coupling at one loop accuracy in this study . We would like to point out although we have found that the Balitsky and Kovchegov-Weigert kernels reduce to the same one, in the saturation region in our previous studies [28], it has been found by numerical studying the rcBK equation that the Balitsky kernel gives a more reasonable solution as suggested in Ref. [20]. Therefore, we will only use the Balitsky kernel in the all following studies. The second term on the right hand side of Eq.(6) is the 'subtraction part', which has form as [14] S where the α µ and w ⊥ are the bare coupling and subtraction point, respectively. The subtraction point can be chosen to be the transverse coordinate of the emitted gluon z ⊥ [19] or the transverse coordinate of either the quark z ⊥1 or the antiquark z ⊥2 [20]. The evolution kernel in Eq.(13) is written as where the K g 1 (x ⊥m , y ⊥n ; z ⊥1 , z ⊥2 ) is the resummed JIMWLK kernel. When one substitutes Eq. (14) into Eq.(13), it seems that the Eq.(13) is so complicated that the evolution equation, Eq.(6), cannot be solved analytically. Fortunately, in this study we only focus on the physics in the high-energy region, in which the unitarity corrections become important and the S-matrix is small, thus the quadratic terms of the S-matrix can be neglected. Keeping the linear term of the S-matrix only, the Eq.(6) simplifies to where the factor 2 on the right hand side of Eq.(15) comes from considering of the symmetry of the two integral regions as the case of fixed coupling. Performing the integral over the dipole size r 1 in Eq.(15), one gets whose solution is [28] with the saturation momentum in the NLO case We would like to note that in the exponent of the S-matrix the rapidity is changed from quadratic dependence in the fixed coupling case, Eq. (5), to linear dependence in the running coupling case, Eq. (17). This change significantly slows down the evolution speed of the dipole scattering amplitude, which is supported by the HERA data [15,34].

C. Full next-to-leading order Balitsky-Kovchegov equation and its solution
As we know that the rcBK equation only considers part of NLO corrections from quark loop contributions. It has been shown in Ref. [21] that a comprehensive corrections should include both the contributions from the quark and gluon loops as well as from the tree gluon diagrams with quadratic and cubic nonlinearities. Including all the above mentioned NLO corrections, one obtains the full NLO BK evolution equation as [21] where the kernels are We use the notation which are the transverse size of dipoles. In Eq. (20), the b is the first coefficient of the β function, and µ is the renormalization scale. To simplify the calculation, we firstly absorb the term involving the renormalization scale µ into the running coupling α s , then the terms involving b are absorbed into α s by using the Balitsky running coupling scheme which was developed in Ref. [20]. The kernel K 1 can be rewritten as Note that the Balitsky running coupling prescription resums all α s β contributions, especially the term ∼ b ln We focus on dipole scattering in the saturation region in which the unitarity corrections are very important or S approaches to zero. Thus, one can neglect the non-linear terms in Eq. (19). The full NLO BK evolution equation simplifies to Let's turn to analytically solve Eq.(24) in either 1/Q s |r 1 | |r|, |r 2 | ∼ |r| or 1/Q s |r 2 | |r|, |r 1 | ∼ |r| region as mentioned in the LO case. If we choose to work in the first regime, the NLO kernel K 1 reduces to Substituting the simplified kernel into Eq. (24), the evolution equation becomes and has a solution as following [29] S where the NLO saturation momentum has the same definition as Eq. (18) and . Let's look at the solutions in LO Eq. (5), running coupling Eq. (17), and full NLO Eq. (27), and compare the variation of these solutions. We can see that the NLO corrections slow down the evolution speed of the dipole scattering amplitude, the running coupling effect (quark loop contribution) makes the exponent of the S-matrix changing from quadratic rapidity dependence to linear rapidity dependence. However, the full NLO effects, which include quark loop and gluon loop contributions, force the linear rapidity dependence back to the rapidity raised to the power of 3/2 dependence due to gluon loops binging part of compensation to offset the decrease, see Fig.1 for a diagrammatic depiction.

D. Collinearly improved next-to-leading order Balitsky-Kovchegov equation and its solution
We would like to point out that although we get a analytic solution of the full NLO BK equation in the saturation region, it has been found that the numerical solutions of the full NLO BK equation in Ref. [16] become unstable for many values of initial conditions. The dipole amplitude can decrease with growing energy and can switch to a negative value, which is in disagreement with the theoretical expectations.
It has been shown that the main source for such an instability comes from a large double-logarithmic contribution [23]. To solve this unstable problem, one needs to resum double transverse logarithms to all orders under double logarithmic approximation (DLA). There is also a large single transverse logarithms (STL) which appear in the NLO corrections to the BK equation. Such single logarithms must be kept under control to ensure a good convergence of the α s expansion. The single and double logarithmic resummations were done in Refs. [17,23], they gave an collinearly improved evolution equation as where the collinearly-improved kernel is [22] K CI with the DLA kernel and the STL kernel The constant A 1 in above equations comes from the DGLAP anomalous dimension, the J 1 is the Bessel function with ρ = ln r 2 1 /r 2 ln r 2 2 /r 2 . Note that when ln r 2 1 /r 2 ln r 2 2 /r 2 < 0, then an absolute value is used and the Bessel function J 1 turns into I 1 [17]. Here, we would like to note several points, (i) we focus on the physics in saturation region in this paper, thus we neglect the uninteresting terms in Eq.(28) as compared to Eq. (19); (ii) the double logarithmic resummation is included in the DLA kernel, which is to remove the double logarithmic term from the last line in Eq. (29) in order to avoid double counting; (iii) the second (subtraction) term in Eq. (29) is to subtract the α 2 s part of the single transverse logarithm included in K 2 to avoid double counting, see Ref. [22] for more detailed discussions about the subtraction term. Now let's turn to analytically solve Eq. (28) in the saturation region in which the quadratic term of Smatrix can be neglected due to very small S. As it was done in previous subsections, we choose to work in the region, 1/Q s |r 1 | |r|, |r 2 | ∼ |r|. In this regime, the ρ is equal to zero leading to the DLA kernel K DLA 1, which corroborates the statement that the double logarithm only important in phase-space where the scattering is weak [23]. So, we neglect the quadratic term of S-matrix and double logarithmic corrections, then the Eq.(28) becomes with a simplified kernel whose solution is If one compares the two solutions Eq. (27) and Eq.(34), one can find that the dominant terms in the exponent of the S-matrix are the same, which indicate that in the saturation region the resummations of double and single transverse logarithms is negligible under the parent dipole running coupling prescription. However, in the next section we shall show that the situation has a dramatic change once the parent dipole running coupling prescription is replaced by the smallest dipole prescription.

III. SOLUTIONS TO THE EVOLUTION EQUATIONS WITH THE SMALLEST DIPOLE RUNNING COUPLING PRESCRIPTION
In the last section, we discuss the BK evolution equations with the PDRC prescription. However, the recent studies in Refs. [17,18,35] found that the evolution equation with the SDRC prescription has been advocated to be the correct QCD running coupling prescription, since it is favored by the HERA data at a phenomenological level. It was also pointed out in Ref. [21] that the proper interpretation of the argument of the QCD coupling is running according to the size of the smallest dipole for the BK equation at NLO. In this section, we start with our discussion on the rcBK evolution equation, since the coupling of the LO BK evolution equation is fixed, thus the argument of the coupling does not affect the solution of the LO BK equation.

A. Balitsky-Kovchegov equation and its solution in the running coupling case
As it is shown in Section II B that the rcBK equation can be written as with the Balitsky kernel Here, in order to clearly see the effect of the coupling argument we keep only the running coupling part of the evolution equation and neglect the subtraction terms in Eq. (35). Indeed, the subtraction terms are dropped eventually, since they are proportion to the square of the S-matrix, while the S-matrix is small in this paper's interesting region (saturation regime). We choose to work in the region, 1/Q s << |r 1 | << |r|, |r 2 |, as done in previous section, the Balitsky kernel reduces to [21,28] K rc (r, r 1 , From the above equation, one can see that the argument of the QCD coupling is the smallest dipole size r 1 . Interestingly, the kernel in Eq.(37) is exactly same as the one in Eq. (12), which implies that the rcBK kernel is independent of the choice of the running coupling prescription in the saturation region. So, the solution of the Eq. (35) is which is exactly the same as Eq. (17).

B. Balitsky-Kovchegov equation and its solution in the full next-to-leading order case
In the full NLO case, the BK evolution equation is given by Eq. (19). Here, we focus on the physics in the saturation region in which the S-matrix is small, therefore we neglect the quadratic terms of the S-matrix, the full NLO evolution equation reduces to, We need to emphasize that the argument of the QCD coupling in the second term on the right hand side of the above equation is the smallest dipole size r 1 rather than the parent dipole size as in Eq. (26). One can see that this change shall significantly affect the behavior of the dipole scattering amplitude.
To solve the Eq.(39), we use the running coupling at one loop accuracy which is given by Eq.(11), the equation becomes whose solution is where the saturation momentum is given by Eq. (18) and B = 67/36 − π 2 /12 − 5N f /18N c . By comparing Eq.(41) with Eq. (27), one can see that in the exponent of S-matrix the rapidity raised to power of 3/2 dependence is replaced by the linear rapidity dependence once the SDRC prescription is used. The numerical studies in the next section shall prove this finding. We would like to point out that the solution of the full NLO BK equation now has a similar linear rapidity dependence as the one obtained by solving to the rcBK equation with the SDRC prescription, although the evolution speed of the scattering amplitude is slightly rebounded due to the coefficient of the first term in the exponent becoming larger than the running coupling one. The rebound is caused by the gluon loop contributions as we discussed in our previous studies [29].

C. Collinearly improved Balitsky-Kovchegov equation and its solution in the full next-to-leading order case
In the saturation region, the NLO BK equation with resummation, Eq.(28), can be rewritten as under the SDRC description. We need to point out that the argument of the QCD coupling in the above equation is the smallest dipole size r 1 instead of parent dipole size r, which gives rise to a great impact on the evolution speed of the dipole scattering amplitude. Substituting the one loop running coupling, Eq.(11), into Eq.(42), one gets the evolution equation as whose solution is This result deserves several important comments, • One can see that in the exponent the rapidity raised to power of 3/2 dependence with the PDRC prescription, Eq.(34), is now replaced by linear rapidity dependence under the SDRC prescription, which indicates that the argument of the coupling plays an important role on the rapidity dependence of the S-matrix.
• By comparing the coefficients of the dominant terms between Eq.(41) and Eq.(44), one can find that in the saturation region the resummation of single transverse logarithms takes a significant effect on the suppression of the evolution of the scattering amplitude, although the resummation of the double logarithmic corrections has a negligible effect on it. This finding is supported by the numerical results performed in the next section. To conclude, we can see from the above discussion that the prescription of the QCD running coupling has a strong impact on the rapidity dependence of the S-matrix. The SDRC prescription takes a dramatic suppression on the evolution speed of the dipole scattering amplitude, which is favored by HERA data [18,35] and satisfies the theoretical expectations [21].

IV. NUMERICAL ANALYSIS
In this section, we perform the numerical studies to the NLO BK evolution equations in order to test the analytic solutions obtained in the above sections. The translational invariant approximation is used in our numerical studies, we suppose that the scattering matrix is independent of the impact parameter of the collision, S = S(|r|, Y ). The evolution equations are complicated integro-differential equations, but they can numerically straightforward solve on a lattice. The variable r is discretized into 800 points which are equally separated in the logarithmic space between r min = 10 −8 GeV −1 and r max = 50GeV −1 . We solve these equations by using the GNU Scientific Library (GSL). To be more specific, the Runge-Kutta method with a step size in rapidity ∆Y = 0.1 is used to solve the differential equations, all the integrals are performed using adaptive integration routines, and the cubic spline interpolation method is employed to interpolate the data points.
To solve the evolution equations, we use the McLerran-Venugopalan (MV) model as the initial condition [36], where we set Q s (Y ) = 1 at Y = 0, and put Λ = 0.2GeV. For the running coupling α s (r 2 ) in the evolution equations, we use the one-loop running coupling, Eq. (11), with N f = 3. We freeze the coupling to a fixed value α s (r fr ) = 0.75 for larger dipole size, r > r fr , to regularize the infrared behavior. The left hand panel of Fig.2 shows the solutions of the full NLO BK equation as a function of the dipole transverse size in the cases of PDRC and SDRC prescriptions for 5 different rapidities. To clearly show the numerical solutions in the saturation region, we plot a zooming in diagram in Fig.2. By comparing the full NLO dipole scattering amplitudes with these two different prescriptions, we can see that the amplitudes with the SDRC prescription are smaller than the ones with the PDRC prescription, respectively. This numerical result is in agreement with the analytic finding, Eq.(41) in Section III B, in which the rapidity raised to power of 3/2 dependence in the exponent of the S-matrix is replaced by the linear rapidity dependence once the SDRC prescription is used. The right hand panel of Fig.2 gives the comparison between the solutions of full NLO BK equation and collinearly improved NLO BK equation (including a resummation of large single and double transverse logarithms) with the SDRC prescription for 5 different rapidities. We would like to note that the inner zooming in diagram is to clearly show the corresponding solutions in the saturation region. One can see that the evolution of the dipole scattering amplitude is significantly slowed down by the resummation corrections. This outcome supports the analytic results in Section III B, where the amplitude, Eq.(44), is suppressed by the resummation corrections as compared to the one, Eq.(41), without resummations.
The left hand panel of Fig.3 gives the solutions of the rc and full NLO BK equations as a function of the dipole size with the SDRC prescription for 5 different rapidities. By comparing the corresponding dipole scattering amplitudes between the rc and full NLO cases for each rapidity, we can see that the full NLO scattering amplitude is larger than the running coupling one, which indicates that the gluon loop effect compensates part of decrease made by quark loop effect (running coupling effect). The outcomes is similar as the findings in our previous paper [29] in which the solution was calculated with the PDRC prescription. One can see that the gluon loop has a rebound effect regardless of which running coupling prescription is used. However, we would like to point out that the compensation to the decrease of the dipole scattering amplitude in the SDRC prescription is less than the parent dipole case. In other words, the SDRC prescription slows down the evolution speed of the dipole scattering amplitude more than the PDRC prescription, which is consistent with the phenomenological desire [17]. The right hand panel of Fig.3   the dipole scattering amplitudes in the rcBK case and the collinearly improved NLO BK (resummations of single and double transverse logarithms) for 5 different rapidities in the case of the SDRC prescription. If one compares the corresponding amplitudes between two cases respectively, one finds that the amplitude is dramatically suppressed by the resummation corrections, especially by the resummation of single logarithm. The resummation of the single logarithm takes effect not only in the non-saturation region but also in the saturation region (see the zooming in diagram on the right hand panel of Fig.3), it is unlike the resummation of double logarithm which takes a negligible effect in the saturation region [29].

V. RARE FLUCTUATIONS OF THE S-MATRIX WITH THE SMALLEST DIPOLE RUNNING COUPLING PRESCRIPTION
As we have studied in Refs. [28,32] that the rare fluctuations can play an important effect in the suppression of the evolution speed of the dipole scattering amplitude in the saturation region. It was shown that in the case of fixed coupling (LO) the exponential factor of the S-matrix without the rare fluctuation effect is twice as large as the result when the rare fluctuations are taken into account [37]. We also demonstrated that in the saturation regime the exponential factor of the full NLO S-matrix is √ 2 as large as the result which emerges when the rare fluctuation effects are considered. Although it seems that the rare fluctuations are important in the LO and full NLO BK cases, our studies showed that in the rcBK case the rare fluctuations take a tiny effect in the suppression of the evolution speed of the dipole scattering amplitude, thus they can be ignored. We need to point out that all just aforementioned running couplings are used the PDRC prescription. In this study, we want to see if one switches the QCD running coupling prescription of the NLO BK equations to the SDRC prescription, whether the rare fluctuations are still important or not in the suppression of the evolution speed of the dipole amplitude. In this section, we only present the rare fluctuations of the S-matrix on top of the full NLO effect in the center of mass (CM) frame, since the relevant results are exactly the same when one works in a general frame, please see the Appendix A.
Following the framework of Ref. [37], we consider a high energy scattering of a left-moving dipole on a right-moving dipole at zero impact parameter in the center of mass frame. We let the right-moving and left-moving dipoles have normal Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution [38,39] in the rapidity intervals 0 < y < Y 0 /2 and −Y 0 /2 < y < 0, respectively. In order to produce the rare configuration, we need to constrain the rapidity evolution of the system in such a way that the wavefunctions of the right-moving and left-moving dipoles consist of only parent dipole with size r 0 in the rapidity intervals Y 0 /2 < y < Y /2 and −Y /2 < y < Y 0 /2, respectively. However, the aforemention constraint of the evolution is a optimal case which cannot be obtained in a real dipole scattering, since one cannot require all the evolutions are absent. What we can do is to only allow the evolutions to produce dipoles with size much smaller than r 0 in order to avoid the daughter dipoles evolving into ones with similar size as r 0 in intermediate rapidities, see Fig.4. We require that the gluon emission from the parent dipoles, as part of the evolution which forms the left and right moving states which scatter on each other, is forbidden if the gluon has k ⊥ and y lying in the shaded triangles of Fig.4. The upper line in Fig.4 is determined by the requirement that gluons on the right hand side of that line cannot evolve by normal BFKL evolution into shaded triangle, and the same requirement is applied to the lower line. According to the above description, one knows that the probability of rare configuration S(r, Y − Y 0 ) satisfies the evolution equation [28,32,37] with the kernel whose solution is with B = 67/36 − π 2 /12 − 5N f /18N c . Note that Eq.(49) is derived by using the SDRC prescription instead of the PDRC prescription, since we focus on studying the rare fluctuation effects of the S-matrix in the case of the SDRC prescription in this study.
To obtain the S-matrix including the rare fluctuations in the center of mass frame, we need to compute then the S-matrix can be calculated as If one compares Eq.(51) with Eq.(41), one can find that the dominant terms in the exponent of the Smatrix are almost same, which indicate that the rare fluctuations are not important in the case of the SDRC prescription. Thus, the rare fluctuation effects can be neglected when one works with the SDRC prescription. It is simple to find that the rare fluctuations are also not important in the collinearly improved NLO BK case if one uses the SDRC prescription instead of the PDRC prescription. However, it is necessary to remind that we found a 1/ √ 2 suppression to the exponential factor of the S-matrix when the rare fluctuations are taken into account in our previous studies [32], where the PDRC prescription was used. To ensure the relevant results of the S-matrix are independent of the frame choice, we study the rare fluctuations of the S-matrix on top of the full NLO corrections with the SDRC description in a general frame. We let a right-moving dipole with size r 0 and rapidity Y − Y 2 scattering off a left-moving dipole with size r 1 and rapidity −Y 2 , see Fig.5 for the scattering picture. For later convenience, we assume Y 2 ≤ (Y − Y 0 )/2, where Y 0 is the rapidity difference between the two scattering dipoles at the 'time' of the onset of unitarity corrections. Since we suppose that Y 2 ≤ (Y − Y 0 )/2, which indicates that the left-moving dipole takes the smaller rapidity. Thus, it is easy to suppress its evolution. We do not allow any additional daughter dipoles created by gluon emission, which means that we have to suppress the emission of those dipoles which, after a normal BFKL evolution during the intermediate rapidities −Y 2 < y < 0, could become of size 1/Q s or larger. Therefore, we require that the gluon emission is forbidden if the gluon has k ⊥ and y lying in the lower shaded triangle of Fig.5. For the right-moving dipole, we suppress the evolution in such a way that the gluons laying on the right hand side of the line cannot evolve into the upper shaded triangle in Fig.5 through normal BFKL evolution in the rapidity interval Y 1 + Y 0 < y < Y − Y 2 . We let the parent dipole undergo a normal evolution during Y 1 < y < Y 0 + Y 1 , where Y 1 is determined by maximizing the S-matrix later. The unshaded triangle, 0 < Y < Y 1 , denotes the saturation region where the right-moving dipole has already evolved into a CGC state. In terms of the scattering picture depicted above, we can calculate the the rare fluctuations of the S-matrix on top of the NLO corrections as [32] S(r 0 , r 1 where the S is the S-matrix of a elementary dipole scattering off a CGC. The Now we turn to compute the suppression factors, S R , and S L , which denote no gluon emission from the right-moving and left-moving dipoles, and can be estimated in terms of the area of the upper and lower shaded triangles in Fig.5 [32,37], and