Extended collinearly-improved Balitsky-Kovchegov evolution equation in target rapidity
EExtended collinearly-improved Balitsky-Kovchegov evolution equation intarget rapidity
Wenchang Xiang,
1, 2, ∗ Yanbing Cai, † Mengliang Wang, ‡ and Daicui Zhou § Guizhou Key Laboratory in Physics and Related Areas,and Guizhou Key Laboratory of Big Data Statistic Analysis,Guizhou University of Finance and Economics, Guiyang 550025, China Department of Physics, Guizhou University, Guiyang 550025, China Key Laboratory of Quark and Lepton Physics (MOE), and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China
An extended collinearly-improved Balitsky-Kovchegov evolution equation in the target rapidityrepresentation 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 evolutionequation, as one can see that there are extra contributions to the evolution kernel once the run-ning coupling is included. To identify the significance of the corrections, we numerically solve theevolution equation with and without the running coupling contributions during the S -matrix expan-sion. The numerical results show that the scattering amplitude is largely suppressed by the runningcoupling corrections, which indicate that one needs to consider the running coupling contributionsduring 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 todescribe 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–5] equationand its mean field version known as the Balitsky-Kovchegov (BK) equation[1, 6], have been able to qualita-tively 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–13].However, it has been found that the LO CGC theory is insufficient for direct applications to the phenomenol-ogy, since the evolution speed of the scattering amplitude resulting from the LO BK equation is too fast toquantitatively describe the experimental data[14, 15].There are tremendous developments in the calculations of the next-to-leading order (NLO) correctionsto the JIMWLK and BK equations in the literature over the past fifteen years[16–22]. The pioneer worktowards the NLO corrections to the BK equation was performed by including the running coupling effectin Refs.[16] and [17]. A running coupling Balitsky-Kovchegov (rcBK) equation was obtained. It was shownthat the growth of the dipole-hadron scattering amplitude resulting from the rcBK equation is significantlyslowed down as compared to the LO one. Furthermore, it was found that the rcBK equation gives a rathersuccessful description of the HERA data[14]. However, the running coupling effect is not the only largehigher order corrections to the LO BK or JIMWLK equations. Except the running coupling, the authors inRef.[18] derived the full NLO BK evolution equation by including the quark loops, gluon loops, as well as thetree gluon diagrams with quadratic and cubic non-linearities, they obtained other contributions which areenhanced by double transverse logarithms. Unfortunately, the numerical study of the full NLO BK equationfound that the equation is unstable, since the scattering amplitude can decrease as rapidity increasing andcan even turn to negative value for small dipoles[23, 24]. The instability was traced back to the radiativecorrections enhanced by double transverse logarithms.To cure the instability problem, one has to resum the radiative corrections enhanced by double transverselogarithms to all orders. Two methods related to this specific issue were proposed[25, 26], which use differentrecipes 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 LOBK equation, but is non-local in rapidity Y ; (ii) the resummation of the leading double logarithms has beenperformed in the evolution kernel in Ref.[26] resulting in a collinearly-improved Balitsky-Kovchegov (ciBK) ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] The JIMWLK is the abbreviation of Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov, Kovner. a r X i v : . [ h e p - ph ] F e b equation with a modified kernel, which is still local in rapidity Y . It has been shown that these two methodsare equivalent in the resummation of the leading double transverse logarithms. It is known that both ofthem lead to stable evolution equation and give good fits to the small- x HERA data[27–29]. However, theauthors in Ref.[30] found there are some inconsistencies with the original analysis in Refs[25, 26]. Theyfound that the instability of the full NLO BK is caused by the wrong choice of the rapidity variable whichplays the role of the evolution time. Moreover, it has been found that the growth of the saturation exponentwith the running coupling constant in the asymptotic region of rapidity has a strong scheme dependence onthe 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 isthat of the projectile rapidity Y . In terms of the previous experience with the NLO BFKL equation[31–33], where a similar issues were met and eventually cured, the evolution variable should be the rapidity ofthe hadronic target, η , rather than the rapidity of the projectile. Inspired by the experience on handlinginstability problem in NLO BFKL equation, a novel method was proposed to derive the collinearly-improvedBK equation in η -representation (ciBK- η ) from the corresponding evolution equation in Y -representationby the change of variable η = Y − ρ , where ρ is defined as ρ = ln( Q /Q ) with Q and Q to be the hardscale of the projectile and soft scale of the target, respectively[30]. As a consequence, the evolution variablein the ciBK- η equation is the physical rapidity η = ln(1 /x ), not the rapidity of the dipole projectile, and theciBK- η 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 foundthat 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 partof the ciBK- η equation, which is called “canonical” Balitsky-Kovchegov equation (caBK- η ), was used tostudy the data in order to avoid the cumbersomely numerical calculations, thus a set of NLO terms of order O ( α s ) were abandoned; (ii) the LO BK approximation was used in the replacement of the first derivativein the expansions of the “real” S -matrix 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 therunning coupling corrections during the expansions of the “real” S -matrix. To see the significance of thecorrections of the running coupling, we firstly recall the derivation of the ciBK- η equation with the LO BKapproximation in the replacement of the first derivative in the expansion of the “real” S -matrix. Second, wederive an extended collinearly-improved Balitsky-Kovchegov equation in η -representation by emphasizingthe running coupling corrections in the expansion of the “real” S -matrix. We obtain an extended “canon-ical” Balitsky-Kovchegov equation (exBK- η ) whose evolution kernel is modified by the running couplingcorrections as compared to the caBK- η equation. The exBK- η equation is analytically solved in saturationregion. Its analytic solution shows that the exponent of the S -matrix has a linear dependence on rapidityinstead 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 couplingcorrections. We compare our extended collinearly-improved Balitsky-Kovchegov equation in η with its orig-inal form to see how big difference is. It is easy to find that there are eight extra terms resulting fromthe running coupling corrections, which indicate that the precision of the expansion of the S -matrix has asignificant impact on the evolution equation.Finally, we numerically solve the evolution equations in η -representation to test the analytic outcomesmentioned above. The numerical results confirm our analytic findings. The saturation exponents ¯ λ areextracted from the solutions of the LO BK, caBK- η and exBK- η equations. As expected, the runningcoupling effect has a large influence on the evolution speed of the front, which largely suppresses the rapidityevolution 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 LOBK and rcBK equations in Y -representation. These two equations shall be used in the Taylor expansionof the “real” S -matrix in the derivation of the collinearly-improved BK equation in η -representation in thenext section. The BFKL is the abbreviation of Balitsky-Fadin-Kuraev-Lipatov. The “real” means a really measuring the scattering of the soft gluon.
A. Leading order BK equation
We consider the high energy scattering between a dipole which is consisted of a quark-antiquark pairmoving towards the positive direction of the longitudinal axis with momentum ( p + , p − , p ), and a hadronictarget moving along the negative direction with momentum ( p +0 , p − , p ). The scattering is treated inthe eikonal approximation, thus the transverse coordinates of the quark ( x ) and the antiquark ( y ) are notmodified by the collision. One can write the dipole scattering matrix as a correlator of two Wilson lines[18] S xy ( Y ) = 1 N c (cid:10) Tr { U ( x ) U † ( y ) } (cid:11) Y , (1)where the (cid:104)· · · (cid:105) Y means the average over target gluon field configurations at Y . Here, the Y is the rapiditydifference between the dipole and hadronic target Y = ln p + p +0 = ln 2 p + p − Q = ln sQ , (2)with the center of mass energy squared s = 2 p + p − and the typical momentum of the target Q . The U inEq.(1) is the time ordered Wilson line U ( x ) = P exp (cid:104) ig (cid:90) d x − A + ( x − , x ) (cid:105) , (3)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] ∂∂Y S xy ( Y ) = ¯ α s π (cid:90) d z ( x − y ) ( x − z ) ( z − y ) [ S xz ( Y ) S zy ( Y ) − S xy ( Y )] , (4)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 anon-linear term. The second term in the right hand side of Eq.(4) is referred as “virtual” term which depictsthe survival probability of the original dipole at the time of scattering. Note that the BK equation resumsonly the leading logarithmic α s ln(1 /x ) corrections in the fixed coupling case, thus it is a LO evolutionequation.In the later section of this paper, we shall be interested in the limiting form of the S -matrix in saturationregion where the dipole has very large size, such that r Q s (cid:29)
1. Here, the Q s is the typical transversemomentum of the saturated gluons and is rapidly increasing with Y . For later comparison, we analyticallysolve the Eq.(4) in the following. In saturation regime, the S -matrix approaches the black-disk limit, S → S in Eq.(4) can be neglected. Moreover, the saturation condition impliesthat the two daughter dipoles are larger than the typical transverse size r s ∼ /Q s . Therefore, the Eq.(4)can reduce to ∂S ( r, Y ) ∂Y (cid:39) − ¯ α s π (cid:90) r /Q s d r r r r S ( r, Y ) . (5)with r = | x − y | , r = | x − z | , r = | z − y | be the transverse size of the parent dipole, and two daughterdipoles, respectively. Further, one can find that the integral in Eq.(5) is governed by the case if one of thedaughter dipoles is much smaller than the parent dipole, r (cid:28) r and r ∼ r , or r (cid:28) r and r ∼ r [35, 36].We select to work in the case r ∼ r , the Eq.(5) can be rewritten as ∂S ( r, Y ) ∂Y (cid:39) − ¯ α s (cid:90) r /Q s d r r S ( r, Y ) , (6)where the right hand side includes a factor of 2 to take into account that the smaller dipole with size r canbe any one of the two daughter dipoles. Now we carry out the integrals in Eq.(6), which yield the analyticsolution of the LO BK equation as [35–38], S ( r, Y ) = exp (cid:20) − λ ¯ α s Y − Y ) (cid:21) S ( r, Y ) , (7)where we have used ln[ r Q s ( Y )] (cid:39) λ ¯ α s ( Y − Y ) with Y the rapidity scale at which Q s ( Y ) = 1 /r . TheEq.(7) is usually called as Levin-Tuchin formula, since it was firstly derived by them in Ref.[37]. From theabove 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 -matrixhas a quadratic dependence on rapidity, which renders the S -matrix too small at large rapidities, in theother words, the rapidity evolution speed of the dipole amplitude N is too fast, with N = 1 − S . It hasbeen recognized that the aforemention drawbacks are the reasons why the LO BK equation is insufficient todescribe the experimental data at HERA[14, 15, 39]. Therefore, one has to include the NLO corrections tothe LO BK evolution equation, such as the running coupling effect which can suppress the evolution speedof 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 loga-rithmic α s ln(1 /x ) corrections with a fixed coupling. Beyond the leading logarithmic approximation, therewas a significant progress in the BK evolution equation via the resummation of α s N f to all orders, whichis called as the running coupling corrections. The running coupling Balitsky-Kovchegov (rcBK) equationhave been derived independently by Balitsky in Ref.[16] and Kovchegov and Weigert in Ref.[17]. This twogroups obtained an analogous structure of the rcBK equation but with different evolution kernels. In thisstudy, we shall use the Balitsky version of kernel, since it is favored by the HERA data[14]. Note that wedon’t plan to give a details of the derivation of the rcBK equation, it is out of interesting of this paper. ThercBK equation reads ∂S ( r, Y ) ∂Y = (cid:90) d r K rc ( r, r , r ) [ S ( r , Y ) S ( r , Y ) − S ( r, Y )] , (8)where K rc ( r, r , r ) is the running coupling evolution kernel[16] K rc ( r, r , r ) = N c α s ( r )2 π (cid:20) r r r + 1 r (cid:18) α s ( r ) α s ( r ) − (cid:19) + 1 r (cid:18) α s ( r ) α s ( r ) − (cid:19)(cid:21) , (9)with r min = min { r, r , r } . There are several running coupling prescriptions in the literature[16, 27, 40–42].At the very beginning, the argument of the running coupling α s in the rcBK equation was interpreted to thetransverse size of the parent dipole. Recently, it was found that the size of the smallest dipole is a properargument of the coupling, and the smallest running coupling prescription is favored by the HERA data thanothers[27]. So, we shall use the smallest dipole running coupling prescription in this study. In addition, therunning coupling at one loop accuracy is used α s ( r ) = 1 b ln (cid:0) r Λ (cid:1) , (10)with b = (11 N c − N f ) / π .Now, let us turn to analytically solve the rcBK equation. As it was done in the previous section, we solvethe rcBK equation in saturation region where the dipoles have very large size, r , r , r ≥ /Q s . The integralover the r in Eq.(8) is governed by the case if one of the daughter dipoles is much smaller than the parentone, while the rest of dipole has a similar size as the parent dipole, that is r (cid:28) r and r ∼ r , or r (cid:28) r and r ∼ r . We choose to work in the first case, the Eq.(8) reduces to ∂S ( r, Y ) ∂Y = 1 π (cid:90) r /Q s d r ¯ α s ( r ) r [ S ( r , Y ) S ( r , Y ) − S ( r, Y )] , (11)where the right hand size includes a factor of 2 due to the fact that the smaller dipole can come from eitherof 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 ∂S ( r, Y ) ∂Y (cid:39) − (cid:90) r /Q s d r ¯ α s ( r ) r S ( r, Y ) . (12)Substituting Eq.(10) into Eq.(12), and performing the integrals over r and Y , we can get the analyticsolution of the rcBK equation as[36, 38] S ( r, Y ) = exp (cid:40) − N c bπ ( Y − Y ) (cid:34) ln (cid:32) (cid:112) λ (cid:48) ( Y − Y )ln r Λ (cid:33) − (cid:35)(cid:41) S ( r, Y ) , (13)where the NLO saturation momentum is usedln Q s Λ = (cid:112) λ (cid:48) ( Y − Y ) + O ( Y / ) . (14)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 replacedby the linear rapidity dependence due to the running coupling corrections. This outcome indicates that theevolution speed of the dipole amplitude is slowed down by the running coupling corrections. This finding isconsistent with the theoretical expectations[16, 43]. Furthermore, the phenomenological applications of thercBK equation at HERA energies show that the rcBK equation gives a more reasonable description of theexperimental 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 physicalrapidity 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 wrongchoice of the “evolution time”, this refers to the choice of the rapidity variable[30]. However, if one simplytransforms the ciBK- Y equation to η -representation by change of variable η = Y − ρ , the instability problemis still existence due to NLO corrections enhanced by double collinear logarithms, but not as severe as forthe 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 onthe successive gluon emissions, one can know that the ordering of the successive emissions in longitudinalmomentum should be enforced in order to solve instability problem in η -representation. By doing this, theciBK- η was obtained in Ref.[30], which can directly apply to phenomenology and supposes to give a betterdescription of the HERA data than the other evolution equations in the literature. However, a very recentstudy 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 thefirst 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 isconcerned. The reason why the caBK- η does not give a superior description of the data, could possiblybe 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 levelof running coupling when the “real” S -matrix is expanded. An extended collinearly-improved BK equationin η is obtained, which has the same structure as the ciBK- η equation but with a running coupling modifiedkernel. In the next section, we use the numerical method to solve the caBK- η and exBK- η equations. Thenumerical results show that the running coupling corrections play a significant role in the suppression of theevolution 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] wherethe 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 xy ( Y ) ∂Y = ¯ α s π (cid:90) d z · ( K + K q + K g ) · (cid:16) S xz ( Y ) S zy ( Y ) − S xy ( Y ) (cid:17) + ¯ α s π (cid:90) d u d z · K · (cid:16) S xu ( Y ) S uz ( Y ) S zy ( Y ) − S xu ( Y ) S uy ( Y ) (cid:17) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:16) S xz ( Y ) S uy ( Y ) − S xu ( Y ) S uy ( Y ) (cid:17) (15)with K = ( x − y ) ( x − z ) ( z − y ) , (16) K q = ( x − y ) ( x − z ) ( z − y ) ¯ α s (cid:20) b ln( x − y ) µ − b ( x − z ) − ( y − z ) ( x − y ) ln ( x − z ) ( y − z ) (cid:21) , (17) K g = ( x − y ) ( x − z ) ( z − y ) ¯ α s (cid:20) − π − N f N c −
12 ln ( x − z ) ( x − y ) ln ( y − z ) ( x − y ) (cid:21) , (18) K = 1( u − z ) (cid:26) − x − u ) ( y − z ) + ( x − z ) ( y − u ) − x − y ) ( u − z ) ( x − u ) ( y − z ) − ( x − z ) ( y − u ) ln ( x − u ) ( y − z ) ( x − z ) ( y − u ) + ( x − y ) ( u − z ) ( x − u ) ( y − z ) (cid:20) x − y ) ( u − z ) ( x − u ) ( y − z ) − ( x − z ) ( y − u ) (cid:21) ln ( x − u ) ( y − z ) ( x − z ) ( y − u ) (cid:27) , (19)and K f = 1( u − z ) (cid:20) − ( x − u ) ( y − z ) + ( x − z ) ( y − u ) − ( x − y ) ( u − z ) ( x − u ) ( y − z ) − ( x − z ) ( y − u ) ln ( x − u ) ( y − z ) ( x − z ) ( y − u ) (cid:21) . (20)From Eq.(15), one can see that the full NLO BK equation has two main changes in the structure as comparedto 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 twoadditional partons except the original parent partons (quark and antiquark). In large N c limit, they arecorresponding to two consecutive emissions, the original parent dipole ( x , y ) emits a gluon at transversecoordinate u , which is equivalent to two daughter dipoles ( x , u ) and ( u , y ). Then the dipole ( u , y ) emits agluon at transverse coordinate z , which yields the dipoles ( u , z ) and ( z , y ).
1. The “canonical” BK equation
In order to transform the Eq.(15) from Y -representation into η -representation, we need to change variablesin terms of Y = η + ρ. (21)Now, we can rewrite the S -matrices in η representation as S xy ( Y ) = S xy ( η + ρ ) ≡ ¯ S xy ( η ) , (22) S xz ( Y ) = S xz ( η + ρ ) = S xz (cid:18) η + ln ( x − z ) ( x − y ) + ρ xz (cid:19) = ¯ S xz (cid:18) η + ln ( x − z ) ( x − y ) (cid:19) , (23) S zy ( Y ) = S zy ( η + ρ ) = S zy (cid:18) η + ln ( y − z ) ( x − y ) + ρ zy (cid:19) = ¯ S zy (cid:18) η + ln ( y − z ) ( x − y ) (cid:19) , (24)where we have used ρ = ln (cid:16) Q Q (cid:17) = ln (cid:18) x − y ) Q (cid:19) = ln (cid:18) ( x − z ) ( x − z ) ( x − y ) Q (cid:19) = ln (cid:18) ( x − z ) ( x − y ) (cid:19) + ρ xz , (25)and ρ = ln (cid:16) Q Q (cid:17) = ln (cid:18) x − y ) Q (cid:19) = ln (cid:18) ( z − y ) ( z − y ) ( x − y ) Q (cid:19) = ln (cid:18) ( z − y ) ( x − y ) (cid:19) + ρ zy . (26)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 expansionsof the S xz and S zy can be expressed as S xz ( Y ) = ¯ S xz (cid:18) η + ln ( x − z ) ( x − y ) (cid:19) (cid:39) ¯ S xz ( η ) + ln ( x − z ) ( x − y ) ∂ ¯ S xz ( η ) ∂η , (27)and S zy ( Y ) = ¯ S zy (cid:18) η + ln ( y − z ) ( x − y ) (cid:19) (cid:39) ¯ S zy ( η ) + ln ( y − z ) ( x − y ) ∂ ¯ S zy ( η ) ∂η . (28)To the order of interesting, we need only to keep the first non-trivial term in the above expansions, sinceeach ∂S/∂η is formally suppressed by a power of ¯ α s . For the derivative terms in Eqs.(27) and (28), in thissubsection we use the LO BK equation to evaluate them as what was done in Ref.[30], S xz ( Y ) = ¯ S xz (cid:18) η + ln ( x − z ) ( x − y ) (cid:19) (cid:39) ¯ S xz ( η )+ ¯ α s π (cid:90) d u ( x − z ) ( x − u ) ( u − z ) ln ( x − z ) ( x − y ) (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) , (29)and S zy ( Y ) = ¯ S zy (cid:18) η + ln ( y − z ) ( x − y ) (cid:19) (cid:39) ¯ S zy ( η ) + ¯ α s π (cid:90) d u ( y − z ) ( y − u ) ( u − z ) ln ( y − z ) ( x − y ) (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) . (30)From the above equations, one can see that the rapidity shift in the argument of S -matrices is equivalent toadding a term of order O (¯ α s ). We would like to point out that the LO BK equation is a rough approximationto estimate the derivative terms in Eqs.(27) and (28). Actually, the LO BK equation is insufficient due toits lower precision. To reach the interested order of accuracy, one should use at least the level of rcBKequation 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 BKequation in η [30] ∂S xy ( η ) ∂η = ¯ α s π (cid:90) d · K · (cid:16) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z · ( K q + K g ) · (cid:16) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z d u · K l · ¯ S xu ( η ) (cid:16) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:17) + ¯ α s π (cid:90) d u d z · K · (cid:16) ¯ S xu ( η ) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:16) ¯ S xz ( η ) ¯ S uy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) , (31)with K l = ( x − y ) ( x − u ) ( u − z ) ( z − y ) ln ( u − y ) ( x − y ) . (32)The third term in the right hand side of Eq.(31) is resulting from the expansion of the “real” S -matrices interms of rapidity shift with the LO BK approximation in the evaluating of the derivative term. Note thatthe rapidity shift is neglected in the all NLO terms when Eq.(31) is derived. All the terms of order O (¯ α 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 (¯ α s ). While we are only interested in the terms up tothe order of O (¯ α s ), thus all the terms beyond ¯ α s are abandoned in this paper. Moreover, two mathematicaltricks 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 ofEq.(31) is relabelled in terms of u ↔ z in order to keep consistence with the physics picture mentionedabove.The Eq.(31) is a NLO evolution equation in η -representation, which is a local equation in rapidity. Bycomparing 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, oneanticipates that the change of variable from Y to η can eliminate the instabilities occurred in Eq.(15). Asexpected that the instabilities caused by the violations of time-ordering (double anti-collinear logarithms) The reason why we call it as semi-finished equation is that it is still a unstable equation, which needs to do the resummationsof the double collinear logarithms to totally get rid of the instability. are disappeared in Eq.(31), since the time-ordering property is automatically guaranteed in the η evolu-tion. Unfortunately, it has been shown that the change of variable triggers off another type of instabilitiesassociated with double collinear logarithms[30].To cure the instabilities mentioned above, it is known that the successive gluon emissions during therapidity evolution have to be simultaneously ordered in lifetime and longitudinal momentum, τ p (cid:29) τ k (cid:29) τ , (33)and p + (cid:29) k + (cid:29) p +0 ⇒ p + Q Q (cid:29) k + k k (cid:29) p +0 Q Q ⇒ τ p Q (cid:29) τ k k (cid:29) τ Q , (34)where the ( p + , p − , p ), ( p +0 , p − , p ), 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 (cid:28) Q ) ortoo hard ( k (cid:29) Q ). So, one has to put the constraint on the evolution equation. Before doing that, weneed to rewrite the constraint, Eq.(34), in a proper form. As we know ρ = ln (cid:16) Q Q (cid:17) , Y = ln (cid:16) p + p +0 (cid:17) , (35)and ρ = ln (cid:16) k Q (cid:17) , Y = ln (cid:16) k + p +0 (cid:17) . (36)Thus, the target rapidities can be expressed as η = Y − ρ = ln (cid:16) p + p +0 (cid:17) − ln (cid:16) Q Q (cid:17) = ln (cid:16) τ p τ (cid:17) , (37)and η = Y − ρ = ln (cid:16) k + p +0 (cid:17) − ln (cid:16) k Q (cid:17) = ln (cid:16) τ k τ (cid:17) . (38)Using Eqs.(37) and (38), one can rewrite the constraint, Eq.(34), as η − ln( k Q ) (cid:29) η (cid:29) ln( Q k ) . (39)Moreover, according to the lifetime constraint in Eq.(33), one can get η (cid:29) η (cid:29) . (40)Combining the two constraints, Eqs.(39) and (40), we obtain the final rapidity constraint asmin (cid:26) η, η − ln k Q (cid:27) > η > max (cid:26) , ln Q k (cid:27) . (41)By using Eq.(35) and (36), the above equation can be written in another form asΘ( − ρ ) | ρ | (cid:28) η (cid:28) η − Θ( ρ − ρ )( ρ − ρ ) , (42)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¯ S xy ( η ) = S (0) xy + ¯ α s π (cid:90) d z ( x − y ) ( x − z ) ( z − y ) η − Θ( ρ − ρ )( ρ − ρ ) (cid:90) Θ( − ρ ) | ρ | d η (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) (43)which turns to differential format as ∂ ¯ S xy ( η ) ∂η = ¯ α s π (cid:90) d z ( x − y ) ( x − z ) ( z − y ) Θ (cid:0) η − δ xyz (cid:1) Θ (cid:0) η − Θ( − ρ ) | ρ | (cid:1) × (cid:2) ¯ S xz ( η − δ xyz ) ¯ S zy ( η − δ xyz ) − ¯ S xy ( η − δ xyz ) (cid:3) , (44)with the rapidity shift δ xyz as δ xyz = max (cid:26) , ln ( x − y ) min { ( x − z ) , ( z − y ) } (cid:27) . (45)It has been checked that the Eq.(44) has little scheme dependence on the prescription of rapidity shift. Inthis 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 ( η ) , (46)with δ xz ; r = max (cid:26) , ln r ( x − z ) (cid:27) , (47)and δ zy ; r = max (cid:26) , ln r ( z − y ) (cid:27) . (48)Substituting Eq.(46) into Eq.(44), one gets the caBK- η equation as[30] ∂ ¯ S xy ( η ) ∂η = ¯ α s π (cid:90) d z · K · Θ (cid:0) η − δ xyz (cid:1)(cid:2) ¯ S xz ( η − δ xz ; r ) ¯ S zy ( η − δ zy ; r ) − ¯ S xy ( η ) (cid:3) , (49)which is a non-local evolution equation in rapidity η .
2. Full collinear improved BK equation at O (¯ α s ) The caBK- η equation can be generalized to full NLO accuracy by adding all the NLO corrections at O (¯ α s )in Eq.(31). Before adding those terms, we need to subtract the the O (¯ α s ) piece in Eq.(49). The O (¯ α s )piece in Eq.(49) can be identified by expanding the ¯ S xz ( η − δ xz ; r ) ¯ S zy ( η − δ zy ; r ) term as what we have donein Eqs.(29) and (30). One finds the O (¯ α s ) piece as − ¯ α s π (cid:90) d z d u ( x − y ) ( x − u ) ( u − z ) ( z − y ) δ uy ; r ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) . (50)Adding the O (¯ α s ) pieces from Eq.(31) to Eq.(49) and subtracting the above O (¯ α s ) piece, one obtains theciBK- η equation as[30] ∂S xy ( η ) ∂η = ¯ α s π (cid:90) d z · K · Θ (cid:0) η − δ xyz (cid:1)(cid:16) ¯ S xz ( η − δ xz ; r ) ¯ S zy ( η − δ zy ; r ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z · ( K q + K g ) (cid:16) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z d u · K · ¯ S xu ( η ) (cid:16) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:17) + ¯ α s π (cid:90) d u d z · K · (cid:16) ¯ S xu ( η ) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:16) ¯ S xz ( η ) ¯ S uy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) , (51)with K = K l + δ uy ; r . (52)0One can see that Eq.(51) does not include double collinear logarithms now. These double logarithms areincluded in the first line in the right hand side of Eq.(51) through the rapidity shift. In addition, the doubleanti-collinear logarithm term in kernel K g in the second line in the right hand side of Eq.(51) is canceled bya relative piece generated by the integral over u in the third line. All the unstable factors are under controlin 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 HERAdata than the kcBK and ciBK equations as the theoretical expectations[34]. The reason why the outcomesresulting from the caBK- η equation are not as desired, possibly comes from the insufficient accuracy of theexpansion of the “real” S -matrices (Eqs.(29) and (30)) and the integral LO BK equation Eq.(43).
3. 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 thedaughter is much smaller than the other one, while the larger daughter dipole has comparable size as theparent dipole. As we know that the non-locality is only important for the S -matrix which is associated withsmaller dipole. Thus, the Eq.(49) can simplify to ∂ ¯ S ( r, η ) ∂η (cid:39) α s π ¯ S ( r, η ) (cid:90) r /Q s d zz Θ (cid:16) η − ln r z (cid:17) (cid:20) ¯ S (cid:16) z, η − ln r z (cid:17) − (cid:21) , (53)where a factor 2 is taken into account due to the fact that the smaller dipole can come from either of the twodaughter dipoles, Q s is the saturation momentum which is associated with ¯ Q s as r Q s = ( r ¯ Q s ) / (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 conditionrequires the dipole size z larger than the typical size 1 /Q s (lower integral bound in Eq.(53)), which leads tothe “real” S -matrix S ( z, η − ln r z ) is negligibly small;(ii) the integral over z become logarithmic when z ismuch smaller than r . Applying the two rules, Eq.(53) reduces to ∂ ¯ S ( r, η ) ∂η (cid:39) − ¯ α s ¯ S ( r, η ) (cid:90) r /Q s d z z , (54)whose solution is ¯ S ( r, η ) = exp (cid:20) − ¯ α s λ α s ¯ λ (cid:0) η − η (cid:1) (cid:21) ¯ S ( r, η ) , (55)where we have assumed the saturation momentum in η -representation as ¯ Q s = Q exp(¯ λη ). By comparingEq.(55) with Eq.(7), one can see that they have similar form, but the solution of the caBK- η equation hasan extra suppression factor in the exponent, which leads to the evolution speed of the dipole amplitude isslowed down. The numerical solutions of these two equations shall be done in the next section, where thenumerical 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) areapproximately replaced by the LO BK equation, which are insufficient. To achieve the interested orderof accuracy, we shall derive the collinearly-improved BK equation in η by using the rcBK equation (8) toreplace LO BK equation (4) in the expression of the first derivative terms.
1. Extended caBK- η equation Using the rcBK equation, one can re-expand the S -matrices in Eqs.(27) and (28) as S xz ( Y ) = ¯ S xz (cid:18) η + ln ( x − z ) ( x − y ) (cid:19) (cid:39) ¯ S xz ( η ) + ln ( x − z ) ( x − y ) ∂ ¯ S xz ( η ) ∂η (cid:39) ¯ S xz ( η ) + (cid:90) d u K rc ( x , z , u ) ln ( x − z ) ( x − y ) (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) , (56)1and S zy ( Y ) = ¯ S zy (cid:18) η + ln ( y − z ) ( x − y ) (cid:19) (cid:39) ¯ S zy ( η ) + ln ( y − z ) ( x − y ) ∂ ¯ S xz ( η ) ∂η (cid:39) ¯ S zy ( η ) + (cid:90) d u K rc ( z , y , u ) ln ( z − y ) ( x − y ) (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) , (57)with the running coupling evolution kernel K rc ( x , y , z ) = ¯ α s π (cid:20) ( x − y ) ( x − z ) ( y − z ) + 1( x − z ) (cid:18) α xzs α yzs − (cid:19) + 1( y − z ) (cid:18) α yzs α xzs − (cid:19)(cid:21) , (58)where we use the shorthand notation α xzs = α s (( x − z ) ) and similarly for others. Note that Eq.(58)is another form of Eq.(9), expressed in terms of transverse coordinates of the dipoles. Substituting theEqs.(56) and (57) into Eq.(15), we obtain a semi-finished local collinearly-improved BK euqation in η aftersome complicated algebra calculations (for the detailed derivation, see Appendix A), ∂S xy ( η ) ∂η = ¯ α s π (cid:90) d z · K · (cid:16) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z · ( K q + K g ) (cid:16) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z d u · K rc · ¯ S xu ( η ) (cid:16) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:17) + ¯ α s π (cid:90) d u d z · K · (cid:16) ¯ S xu ( η ) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:16) ¯ S xz ( η ) ¯ S uy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) , (59)with K rc = ln ( u − y ) ( x − y ) (cid:20) ( x − y ) ( x − u ) ( u − z ) ( y − z ) + ( x − y ) ( x − u ) ( y − u ) u − z ) (cid:18) α uzs α yzs − (cid:19) + ( x − y ) ( x − u ) ( y − u ) ( y − z ) (cid:18) α yzs α uzs − (cid:19) + ( u − y ) ( x − u ) ( u − z ) ( y − z ) (cid:18) α xus α yus − (cid:19) + 1( x − u ) ( u − z ) (cid:18) α xus α yus − (cid:19) (cid:18) α uzs α yzs − (cid:19) + 1( x − u ) ( y − z ) (cid:18) α xus α yus − (cid:19) (cid:18) α yzs α uzs − (cid:19) + 1( u − z ) ( y − z ) (cid:18) α yus α xus − (cid:19) + 1( y − u ) ( u − z ) (cid:18) α yus α xus − (cid:19) (cid:18) α uzs α yzs − (cid:19) + 1( y − u ) ( y − z ) (cid:18) α yus α xus − (cid:19) (cid:18) α yzs α uzs − (cid:19) (cid:21) . (60)By comparing Eq.(59) with Eq.(31), there are extra eight terms resulting from the running coupling correc-tions, which have a significant impact on the evolution speed of the dipole amplitude. However, these extraterms do not cure the unstable issue of the evolution equation. Based on the discussion after Eq.(31), weknow 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 beordered 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,¯ S xy ( η ) = S (0) xy + (cid:90) d z K rc ( x , y , z ) η − Θ( ρ − ρ )( ρ − ρ ) (cid:90) Θ( − ρ ) | ρ | d η (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) . (61)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 rcBK2equation instead of the LO BK equation as the start point to derive the ciBK- η equation. Performingderivatives over η in Eq.(61), one can get the exBK- η equation ∂ ¯ S xy ( η ) ∂η = (cid:90) d z K rc ( x , y , z )Θ (cid:0) η − δ xyz (cid:1)(cid:2) ¯ S xz ( η − δ xz ; r ) ¯ S zy ( η − δ zy ; r ) − ¯ S xy ( η ) (cid:3) , (62)which has the same structure as the caBK- η equation (49), but with a running coupling modified kernel. Interms of the experience from Y -representation, we deduce that the rapidity evolution of the dipole amplitudeis also suppressed by the modified kernel, which shall be approved by the numerical calculations in the nextsection.
2. Extended full collinearly improved BK equation at O (¯ α s ) Based on the discussion in previous subsection, one can extend the exBK- η equation to full collinearly-improved BK equation in η by adding the NLO terms from Eq.(59). Before writing down the extended fullcollinearly-improved BK equation in η -representation, we need to identify the O (¯ α s ) piece in the right handside of Eq.(62). First, we expand the ¯ S xz and ¯ S zy to linear order in the rapidity shift, and then use thercBK equation to replace the derivative terms,¯ S xz ( η − δ xz ; r ) (cid:39) ¯ S xz ( η ) − δ xz ; r ∂ ¯ S xz ( η ) ∂η (cid:39) ¯ S xz ( η ) − (cid:90) d u K rc ( x , z , u ) δ xz ; r (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) , (63)and ¯ S zy ( η − δ zy ; r ) (cid:39) ¯ S zy ( η ) − δ zy ; r ∂ ¯ S zy ( η ) ∂η (cid:39) ¯ S zy ( η ) − (cid:90) d u K rc ( z , y , u ) δ zy ; r (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) . (64)Substituting Eqs.(63) and (64) into the non-linear term in Eq.(62), we get the O (¯ α s ) piece in the right handside of Eq.(62), − (cid:90) d z d u K rc ( x , y , u ) K rc ( u , y , z ) δ uy ; r ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) , (65)where we have used the property that the running coupling terms are invariant under x − z → z − y (forthe detailed derivation, see Appendix B).Now, we are arriving the final stage to get the extended evolution equation in η -representation. Subtractingthe O (¯ α s ) piece from Eq.(62), and adding the O (¯ α s ) pieces from Eq.(59), we obtain the extended fullcollinearly-improved BK equation in η as ∂S xy ( η ) ∂η = ¯ α s π (cid:90) d z · K · Θ( η − δ xyz ) (cid:16) ¯ S xz ( η − δ xz ; r ) ¯ S zy ( η − δ zy ; r ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z · ( K q + K g ) (cid:16) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:17) + ¯ α s π (cid:90) d z d u · K · ¯ S xu ( η ) (cid:16) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:17) + ¯ α s π (cid:90) d u d z · K · (cid:16) ¯ S xu ( η ) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:16) ¯ S xz ( η ) ¯ S uy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:17) , (66)with K = K rc + δ uy ; r , (67)which is non-local in η , and a stable equation. The double anti-collinear logarithmic term in the second linein the right hand side of Eq.(66) is canceled by the relevant piece generated via the integral over u in the3third term. The Eq.(66) is free of double collinear logarithms, since all these logarithms are fully included inthe first term. By comparing Eq.(66) with Eq.(51), one can find that in the right hand side of Eq.(66) thereare eight extra terms which are resulting from the running coupling corrections in the Taylor expansion ofthe S -matrix. These extra terms play a significant role in the suppression of the evolution speed of thedipole amplitude.
3. 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, weknow that one of the two daughter dipoles has similar size as the parent dipole, but the size of the rest oneis much smaller than the parent dipole. Moreover, it is known that the non-locality is only important forthe S -matrix which is associated with small size. Therefore, we can reduce Eq.(62) to ∂ ¯ S ( r, η ) ∂η (cid:39) π ¯ S ( r, η ) (cid:90) r /Q s d zz ¯ α s ( z )Θ (cid:16) η − ln r z (cid:17) (cid:20) ¯ S (cid:16) z, η − ln r z (cid:17) − (cid:21) , (68)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 inSec.III A 3. The one loop running coupling, Eq.(10), is used in the calculations. The Eq.(68) becomes ∂ ¯ S ( r, η ) ∂η (cid:39) − N c π ¯ S ( r, η ) (cid:90) r /Q s d z z b ln z Λ . (69)Performing the integral over z , we get ∂ ¯ S ( r, η ) ∂η (cid:39) − N c bπ (cid:20) ln (cid:18) ln Q s Λ (cid:19) − ln (cid:18) ln 1 r Λ (cid:19)(cid:21) , (70)whose solution is¯ S ( r, η ) = exp (cid:40) − N c bπ ( η − η ) (cid:34) ln (cid:32) √ ¯ λ (cid:48) ( η − η ) + √ ¯ λ (cid:48) ln r Λ ( √ η − η + √ ¯ λ (cid:48) ) ln r Λ (cid:33) − (cid:35)(cid:41) ¯ S ( r, η ) , (71)with the saturation momentum in NLO case,ln ¯ Q s Λ = (cid:113) ¯ λ (cid:48) ( η − η ) + O ( η / ) , (72)and r Q s (cid:39) (cid:2) r ¯ Q s (cid:3) (cid:114) ¯ λ (cid:48) η . (73)By comparing Eq.(71) with Eq.(55), one can see that the quadratic rapidity dependence in the exponent ofthe S -matrix is replaced by the linear rapidity dependence once the running coupling corrections are takeninto account. This change implies that the evolution speed of the dipole amplitude is suppressed by theNLO corrections. IV. NUMERICAL ANALYSIS
To test the analytic results obtained in the above section, we shall numerically solve the evolution equa-tions in this section. The Eqs.(49) and (62) are integro-differential equations which can be numericallystraightforward solved on a lattice. To simplify the computation, we neglect impact parameter dependenceof the dipole amplitude throughout this numerical calculations, which imply that the dipole amplitude doesnot depend on angle, N ( r , Y ) = 1 − S ( r , Y ) = 1 − S ( | r | , Y ). Thus, we can view the evolution equations as aset 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 spacebetween r min = 10 − GeV − and r max = 50GeV − . The GNU Scientific Library (GSL) is a good candi-date to solve them, since the GSL contains almost all the routines required by our purpose, such as the4Runge-Kutta method for solving differential equations, adaptive integral routines for performing numericalintegrals, 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], N GBW ( r, η = 0) = (cid:26) − exp (cid:20) − (cid:18) r Q (cid:19) p (cid:21)(cid:27) /p , (74)with p = 4, and Q s0 = 0 . 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 thesmallest dipole running coupling prescription which means the argument of coupling is the smallest dipoleamong the parent and daughter dipoles, α s ( r ) = α s (cid:0) min { ( x − y ) , ( x − z ) , ( z − y ) } (cid:1) . (75)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 exponentis numerically calculated ¯ λ = d ln ¯ Q s ( η )d η , (76)where the saturation moment ¯ Q s ( η ) is determined by N ( r = 1 / ¯ Q s , η ) = κ with κ to be a constant of order1. . . . . . − N ( R , η ) R = ln (1 /r ) I.C. GBWLO BK η = 6 caBK- ηη = 12 caBK- ηη = 18 caBK- η . . . . . − N ( R , η ) R = ln (1 /r ) I.C. GBWLO BK η = 6 exBK- ηη = 12 exBK- ηη = 18 exBK- η FIG. 1: The numerical solutions to LO BK, caBK- η , and exBK- η equations for 4 different rapidities. The lefthand panel shows the comparisons of the evolution speed between the LO BK and caBK- η dipole amplitudes. Theright hand panel gives the comparisons of the evolution speed between LO BK and exBK- η dipole amplitudes. Thezooming diagrams show the relevant results in the saturation region. The left-hand panel of Fig.1 gives the solutions of the LO BK and caBK- η equations for 4 differentrapidities. By comparing the solutions for each respective rapidities, one can see that the values of thecaBK- η dipole amplitude are smaller than the LO BK ones, which indicate that the NLO correctionsenhanced by the double transverse logarithms suppress the evolution speed of the dipole amplitude. Azooming in diagram is provided to clearly show the numerical results in the saturation region. One cansee that the evolution is also slowed down in the saturation region. This numerical outcome is consistentwith the analytic results in Eqs.(7) and (55), where the analytic solutions of the LO BK and caBK- η haveanalogous expression, but with different coefficients in the exponent. The coefficient in the caBK- η caseis 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 for4 different rapidities. We plot a inner zooming in diagram for a clear comparison between the LO BKand exBK- η dipole amplitudes in saturation region. One can see, from the zooming in diagram, that therespective solution of the exBK- η equation in each rapidity is much smaller than the LO BK one, whichimplies that the evolution speed of the dipole amplitude is significantly suppressed. The suppression is5 . . . . . − N ( R , η ) R = ln (1 /r ) I.C. GBW η = 6 exBK- ηη = 12 exBK- ηη = 18 exBK- ηη = 6 caBK- ηη = 12 caBK- ηη = 18 caBK- η . . . . . . ¯ λ η LO BKcaBK- η exBK- η FIG. 2: The left hand panel gives the comparisons of the evolution speed between the caBK- η and exBK- η dipoleamplitudes. The inner diagrams are the zooming in amplitudes in the saturation region. The right hand panel showsthe saturation exponent as a function of η in the LO BK, caBK- η , and exBK- η cases. much more than the caBK- η case. This outcome confirms the analytic result in Sec.III B 3, where the dipoleamplitude 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- η equationsfor 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 inthe saturation region from the inner zooming in diagram. This outcome supports the analytic findings inEq.(71), where the dipole amplitude is further suppressed by the running coupling corrections on top of thedouble logarithm resummation. Finally, we present the η dependence of the saturation exponent predictedby the LO BK, caBK- η , and exBK- η in the right-hand panel of Fig.2. As expected, the ¯ λ resulting fromexBK- η 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. Acknowledgments
This work is supported by the National Natural Science Foundation of China under Grant Nos.11765005,11305040, 11947119 and 11847152; the Fund of Science and Technology Department of Guizhou Provinceunder Grant Nos.[2018]1023, and [2019]5653; the Education Department of Guizhou Province underGrant No.KY[2017]004; the National Key Research and Development Program of China under GrantNo.2018YFE0104700, and Grant No.CCNU18ZDPY04.
Appendix A: Local collinearly-improved BK equation in η with running coupling corrections inTaylor expansion We give a detailed derivation of the collinearly-improved BK equation in η -representation by using thercBK equation to expand out S xz ( Y ) in Eq.(56) and S zy ( Y ) in Eq.(57). To simplify the calculations, werewrite the running coupling evolution kernel, Eq.(58), as K rc ( x , y , z ) = ¯ α s π ( K + K q (cid:48) ) , (A1)with K q (cid:48) = 1( x − z ) (cid:18) α xzs α yzs − (cid:19) + 1( y − z ) (cid:18) α yzs α xzs − (cid:19) . (A2)6The full NLO BK equation, Eq.(15), is our starting point of this derivation. Using Eq.(A1), Eq.(15) isrewritten as ∂S xy ( Y ) ∂Y = (cid:90) d z K rc ( x , y , z ) (cid:2) S xz ( Y ) S zy ( Y ) − S xy ( Y ) (cid:3) + ¯ α s π (cid:90) d z · K g · [ S xz ( Y ) S zy ( Y ) − S xy ( Y )]+ ¯ α s π (cid:90) d u d z · K · [ S xu ( Y ) S uz ( Y ) S zy ( Y ) − S xu ( Y ) S uy ( Y )]+ ¯ α s π N f N c (cid:90) d u d z · K f · [ S xz ( Y ) S uy ( Y ) − S xu ( Y ) S uy ( Y )] . (A3)Note that the kernel K rc is in the order of O (¯ α s ), which means the expansions of S xz ( Y ) and S zy ( Y ) inEqs.(56) and (57) are equivalent to adding a term of order O (¯ α s ). So we can deduce the rapidity shift byusing a similar scheme in LO expansion case. For the first term (to be denoted as T rc ) in the right handside of Eq.(A3), we have T rc = (cid:90) d z K rc ( x , y , z ) (cid:2) S xz ( Y ) S zy ( Y ) − S xy ( Y ) (cid:3) = (cid:90) d z K rc ( x , y , z ) (cid:20) ¯ S xz (cid:18) η + ln ( x − z ) ( x − y ) (cid:19) ¯ S zy (cid:18) η + ln ( y − z ) ( x − y ) (cid:19) − ¯ S xy ( η ) (cid:21) . (A4)Substituting Eqs.(56) and (57) into Eq.(A4), one can get T rc = (cid:90) d z K rc ( x , y , z ) (cid:26)(cid:20) ¯ S xz ( η ) + (cid:90) d u K rc ( x , z , u ) ln ( x − z ) ( x − y ) (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3)(cid:21) × (cid:20) ¯ S zy ( η ) + (cid:90) d u K rc ( z , y , u ) ln ( z − y ) ( x − y ) (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3)(cid:21) − ¯ S xy ( η ) (cid:27) . (A5)After some algebra calculations, Eq.(A5) becomes T rc = (cid:90) d z K rc ( x , y , z ) (cid:110) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + ¯ S xz ( η ) (cid:90) d u K rc ( z , y , u ) ln ( z − y ) ( x − y ) × (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) + ¯ S zy ( η ) (cid:90) d u K rc ( x , z , u ) ln ( x − z ) ( x − y ) (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) + (cid:90) d u d v K rc ( x , z , u ) K rc ( z , y , v ) ln ( x − z ) ( x − y ) ln ( z − y ) ( x − y ) (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) × (cid:2) ¯ S zv ( η ) ¯ S vy ( η ) − ¯ S zy ( η ) (cid:3) (cid:111) . (A6)Using the property that the running coupling terms are invariant under x − z → z − y , we can see that thesecond term is equal to the third term in the brace of Eq.(A6). In addition, the fourth term in the brace ofEq.(A6) is in the order of O (¯ α s ). It becomes O (¯ α s ) order because of an extra K rc factor in front the braceand can be discarded. Therefore, Eq.(A6) can be reduced to T rc = (cid:90) d z K rc ( x , y , z ) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + 2 (cid:90) d z d u K rc ( x , y , z ) K rc ( z , y , u ) ln ( z − y ) ( x − y ) ¯ S xz ( η ) (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) . (A7)In order to keep consistence with the physics picture, we relabel the integral variables in the second termin the right hand side of Eq.(A7) through variable transformation u ↔ z . Then Eq.(A7) becomes T rc = (cid:90) d z K rc ( x , y , z ) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + 2 (cid:90) d z d u K rc ( x , y , u ) K rc ( u , y , z ) ln ( u − y ) ( x − y ) ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) . (A8)7The last three terms in the right hand side of Eq.(A3) are in order of O (¯ α s ). For these three terms, wecan simply replace the rapidity shift, like S xu ( Y ) → S xu ( η ), since the rapidity shift makes them to be oforder O (¯ α s ) which can be safely neglected in our study. Combining these terms with T rc in Eq.(A8), onecan obtain ∂ ¯ S xy ( η ) ∂η = (cid:90) d z K rc ( x , y , z ) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + 2 (cid:90) d z d u K rc ( x , y , u ) K rc ( u , y , z ) ln ( u − y ) ( x − y ) ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) + ¯ α s π (cid:90) d z · K g · (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + ¯ α s π (cid:90) d u d z · K · (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:3) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:2) ¯ S xz ( η ) ¯ S uy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:3) . (A9)The first integral term (to be denoted as T LK ) and second integral term (to be denoted as T SK ) inEq.(A9) contain linear K rc factor and quadratic K rc factor, respectively. For a clear comparison with theLO expansion case, we expend out these two terms by using Eq.(A1) T LK = (cid:90) d z K rc ( x , y , z ) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) = ¯ α s π (cid:90) d z · K · (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + ¯ α s π (cid:90) d z · K q · (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) , (A10)and T SK = 2 (cid:90) d z d u K rc ( x , y , u ) K rc ( u , y , z ) ln ( u − y ) ( x − y ) ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) =2 (cid:90) d z d u ¯ α s π (cid:20) ( x − y ) ( x − u ) ( y − u ) + 1( x − u ) (cid:18) α xus α yus − (cid:19) + 1( y − u ) (cid:18) α yus α xus − (cid:19)(cid:21) × ¯ α s π (cid:20) ( u − y ) ( u − z ) ( y − z ) + 1( u − z ) (cid:18) α uzs α yzs − (cid:19) + 1( y − z ) (cid:18) α yzs α uzs − (cid:19)(cid:21) × ln ( u − y ) ( x − y ) ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) . (A11)After some algebra calculations, Eq.(A11) becomes T SK = ¯ α s π (cid:90) d z d u · K rc · ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) , (A12)where the K rc is defined in Eq.(60).Substituting the Eqs.(A10) and (A12) into Eq.(A9), one can obtain a semi-finished collinearly-improvedBK equation in η∂ ¯ S xy ( η ) ∂η = ¯ α s π (cid:90) d z · K · (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + ¯ α s π (cid:90) d z · ( K q + K g ) · (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) + ¯ α s π (cid:90) d z d u · K rc · ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) + ¯ α s π (cid:90) d u d z · K · (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:3) + ¯ α s π N f N c (cid:90) d u d z · K f · (cid:2) ¯ S xz ( η ) ¯ S uy ( η ) − ¯ S xu ( η ) ¯ S uy ( η ) (cid:3) . (A13)8 Appendix B: Identify an O (¯ α s ) piece in Eq.(62) In this appendix, we give the details of the derivation to identify an O (¯ α s ) piece in the right hand side ofEq.(62). Substituting Eqs.(63) and (64) into Eq.(62) and neglecting the step function in it, one can get ∂ ¯ S xy ( η ) ∂η = (cid:90) d z K rc ( x , y , z ) (cid:26)(cid:20) ¯ S xz ( η ) − (cid:90) d u K rc ( x , z , u ) δ xz ; r (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3)(cid:21) × (cid:20) ¯ S zy ( η ) − (cid:90) d u K rc ( z , y , u ) δ zy ; r (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3)(cid:21) − ¯ S xy ( η ) (cid:27) . (B1)Expanding the terms in the brace in the right hand side of Eq.(B1), we obtain ∂ ¯ S xy ( η ) ∂η = (cid:90) d z K rc ( x , y , z ) (cid:110) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) − ¯ S xz ( η ) (cid:90) d u K rc ( z , y , u ) δ zy ; r × (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) − ¯ S zy ( η ) (cid:90) d u K rc ( x , z , u ) δ xz ; r (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) + (cid:90) d u d v K rc ( x , z , u ) K rc ( z , y , v ) δ xz ; r δ zy ; r (cid:2) ¯ S xu ( η ) ¯ S uz ( η ) − ¯ S xz ( η ) (cid:3) (cid:2) ¯ S zv ( η ) ¯ S vy ( η ) − ¯ S zy ( η ) (cid:3) (cid:111) , (B2)which has a similar structure as Eq.(A6). So we use the same scheme to simplify it. Using the propertythat the running coupling terms are invariant under x − z → z − y and discarding the last term in the orderof O (¯ α s ), Eq.(B2) can be reduced to ∂ ¯ S xy ( η ) ∂η = (cid:90) d z K rc ( x , y , z ) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) − (cid:90) d z d u K rc ( x , y , z ) K rc ( z , y , u ) δ zy ; r ¯ S xz ( η ) (cid:2) ¯ S zu ( η ) ¯ S uy ( η ) − ¯ S zy ( η ) (cid:3) . (B3)In order to keep consistence with the physics picture mentioned in the previous sections, we relabel theintegral variables according to u ↔ z for the second term in the right hand side Eq.(A12), it becomes ∂ ¯ S xy ( η ) ∂η = (cid:90) d z K rc ( x , y , z ) (cid:2) ¯ S xz ( η ) ¯ S zy ( η ) − ¯ S xy ( η ) (cid:3) − (cid:90) d z d u K rc ( x , y , u ) K rc ( u , y , z ) δ uy ; r ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) . (B4)It is clear that the O (¯ α s ) piece in the above equation is − (cid:90) d z d u K rc ( x , y , u ) K rc ( u , y , z ) δ uy ; r ¯ S xu ( η ) (cid:2) ¯ S uz ( η ) ¯ S zy ( η ) − ¯ S uy ( η ) (cid:3) . (B5) [1] I. Balitsky, Nucl. Phys.
B463 (1996) 99.[2] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert,
Nucl. Phys.
B504 (1997) 415.[3] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert,
Phys. Rev.
D59 (1998) 014014.[4] E. Iancu, A. Leonidov, and L. McLerran,
Nucl. Phys.
A692 (2001) 583.[5] E. Ferreiro, E. Iancu, A. Leonidov, and L. McLerran,
Nucl. Phys.
A703 (2002) 489.[6] Y. Kovchegov,
Phys. Rev.
D60 (1999) 034008; ibid.
D61 (1999) 074018.[7] E Iancu, K. Itakura, and S. Munier,
Phys. Lett.
B590 (2004) 199.[8] M. Kozlov, A. Shoshi, and W. Xiang,
JHEP (2007) 020.[9] D. Kkarzeev, and E. Levin,
Phys. Lett.
B523 (2001) 79.[10] A. Dumitru, A. Hayashigaki, and J. Jalilian-Marian,
Nucl. Phys.
A765 (2006) 464.[11] E. Levin, and A. Rezaeian,
Phys. Rev.
D82 (2010) 074016.[12] G. Chirilli, B. Xiao, and F. Yuan,
Phys. Rev. Lett. (2012) 122301.[13] T. Lappi, and H. M¨antysaari,
Phys. Rev.
D88 (2013) 114020.[14] J. Albacete, N. Armesto, J. Milhano, and C. Salgado,
Phys. Rev.
D80 (2009) 034031.[15] J. Albacete, N. Armesto, J. Milhano, P. Quiroga-Arias, and C. Salgado,
Eur. Phys. J.
C71 (2011) 1705. [16] I. Balitsky, Phys. Rev.
D75 (2007) 014001.[17] Y. Kovchegov and H. Weigert,
Nucl. Phys.
A784 (2007) 188.[18] I. Balitsky, and G. Chirilli,
Phys. Rev.
D77 (2008) 014019.[19] A. Kovner, M. Lublinsky, and Y. Mulian,
Phys. Rev.
D89 (2014) 061704.[20] A. Kovner, M. Lublinsky, and Y. Mulian,
JHEP (2014) 030.[21] A. Kovner, M. Lublinsky, and Y. Mulian,
JHEP (2014) 114.[22] D. Zheng, and J. Zhou,
JHEP (2019) 177.[23] T. Lappi, and H. M¨antysaari,
Phys. Rev.
D91 (2015) 074016.[24] T. Lappi, and H. M¨antysaari,
Phys. Rev.
D93 (2016) 094004.[25] G. Beuf,
Phys. Rev.
D89 (2014) 074039.[26] E. Iancu, J. Madrigal, A. Mueller, G. Soyez, and D. Triantafyllopoulos,
Phys. Lett.
B744 (2015) 293.[27] E. Iancu, J. Madrigal, A. Mueller, G. Soyez, and D. Triantafyllopoulos,
Phys. Lett.
B750 (2015) 643.[28] Y. Cai, W. Xiang, M. Wang, and D. Zhou,
Chin. Phys.
C44 (2020) 074110.[29] W. Xiang, M. Wang, Y. Cai, and D. Zhou,
Chin. Phys.
C45 (2021) 014103.[30] B. Ducloue, E. Iancu, A. Mueller, G. Soyez, and D. Triantafyllopoulos,
JHEP (2019) 081.[31] V. Fadin, M. Kotsky, and R. Fiore,
Phys. Lett.
B359 (1995) 181.[32] G. Camici, and M. Ciafaloni,
Phys. Lett.
B412 (1997) 396.[33] M. Ciafaloni, and G. Camici,
Phys. Lett.
B430 (1998) 349.[34] G. Beuf, H. Hanninen, and T. Lappi,
Phys. Rev.
D102 (2020) 074028.[35] A. Mueller, hep-ph/0111244.[36] W. Xiang,
Phys. Rev.
D79 (2009) 014012.[37] E. Levin, and K. Tuchin,
Nucl. Phys.
B573 (2000) 83.[38] W. Xiang, S. Cai, and D. Zhou,
Phys. Rev.
D95 (2017) 116009.[39] C. Contreras, E. Levin, R. Meneses, and I. Potashnikova,
Phys. Rev.
D94 (2016) 114028.[40] W. Xiang, Y. Cai, M. Wang, and D. Zhou,
Phys. Rev.
D101 (2020) 076005.[41] J. Albacete,
Nucl. Phys.
A957 (2017) 71.[42] J. Cepila, J. Contreras, and M. Matas,
Phys. Rev.
D100 (2019) 054015.[43] J. Albacete, and Yu.V. Kovchegov,
Phys. Rev.
D75 (2007) 125021.[44] B. Ducloue, E. Iancu, G. Soyez, and D. Triantafyllopoulos,
Phys. Lett.
B803 (2020) 135305.[45] K. Golec-Biernat, and M. Wusthoff,
Phys. Rev.