Solution to the evolution equation at NLL for high parton density QCD
Abstract
We analytically solve the full nexttoleading logarithmic BalitskyKovchegov equation in the saturation regime, which includes corrections from quark and gluon loops, and large double transverse logarithms. The analytic result for the matrix in the saturation regime shows that the linear rapidity decrease with rapidity of the exponent in running coupling case is replaced by rapidity raised to power of 3/2 decreasing. The collinearlyimproved BalitskyKovchegov equation are also analytically solved in the saturation regime. It shows that the double collinear logarithms do not contribute to the matrix and the solution is the same as the one obtained from the leading order BalitskyKovchegov equation. The numerical solutions to the leading order and full nexttoleading logarithmic BalitskyKovchegov equations are performed in order to test the analytic results derived in the saturation regime.
I Introduction
At high energy dipolehadron scattering perturbative QCD (pQCD) predicts a rapid growth of the gluon density with increasing rapidity (or energy), which leads to nonlinear phenomena such as gluon saturation and multiple scattering. The nonlinear evolution of amplitude between dipole and hadron scattering is well established by the BalitskyJIMWLK^{1}^{1}1The JIMWLK is the abbreviation of JalilianMarian, Iancu, McLerran, Weigert, Leonidov, Kovner. hierarchyBa (); JIMWLK1 (); JIMWLK2 (); JIMWLK3 (); JIMWLK4 () and its mean field approximation known as the BalitskyKovchegov (BK) equationBa (); K () at leading logarithmic accuracy. The BK equation resums leading logarithms to all orders, where is the Bjorken variable. The phenomenological studies of the experimental data from structure functions in deep inelastic scattering at HERAJS (); JM (), and particle production in heavy ion collisions at LHCLM3 (); CX () have been found that the BK equation at leading logarithmic accuracy in pQCD is insufficient to describe these high energy scattering processes.
Over the past decade, our understanding of the BK equation beyond leading logarithmic accuracy has been improved due to including the running coupling correctionsBnlo (); KW (). A running coupling BK (rcBK) equation is obtained, which resums all order corrections associated with the coupling and is an improved version of the leading order (LO) BK equation with the fixed coupling in BK kernel replaced by a ¡°triumvirate¡± of the running couplings. The numerical studies on the solutions of the rcBK have been found that these corrections are numerically large and these corrections significantly slow down the evolution speed of the dipole amplitude with increasing rapidity. However, the running coupling is not the only large corrections to the LO BK equation. In addition to the running coupling corrections (in the language of Feynman diagram, contribution from quark loops), the authors in Ref.BC () have been found that the gluon loops also take a large contribution to the rapidity evolution of the dipole amplitude. Via including contributions from quark and gluon loops, Balitsky and Chirilli got a full nexttoleading (NLO) order BK equation. It is not hard to find that the double transverse logarithmic term in the full NLO BK equation is large when the dipole size is small, which leads to the equation enhanced by double (or collinear) logarithms, see Eq.(III.1). However, the full NLO BK equation suffers from a severe lack of stable problem. The solution of the full NLO BK equation strongly depends on the details of the initial conditionLM (). The solution can be negative with increasing rapidity, which is not a physically meaningful evolution. The origin of this instability can be traced back to the a large double transverse logarithmic NLO correction () in the full NLO BK equationLM (). To solve this instability problem, a resummation of the double transverse logarithms has to be performed.
Recently, the radiative corrections enhanced by double collinear logarithms are resummed by the authors in Ref.IM (). They resummed these corrections to all orders by solving a nonlocal evolution equation, which can be reformulated to a local equation in rapidity with modified kernel and initial condition. The instability behavior is solved once the double transverse logarithms are resummed to all orders. We denote the collinearlyimproved evolution equation as the double logarithmic approximation (DLA) BK equation in this study. One can extend the DLA evolution equation to full nexttoleading logarithmic (NLL) accuracy by adding the NLO BK corrections from BC (). The numerical solutions to the full NLO, DLA and full NLL BK equations are studied in Refs.LM (); IM (); LM2 (). The numerical results show the essential role of the double transverse logarithmic resummation in stabilizing. In addition, the dipole amplitude significantly is slowed down by the double logarithmic resummation. Although the numerical results can be directly applied to phenomenologies, like deep inelastic scatteringIM2 (); JL (), single particle production and twoparticle correlations in heavy ion collisionsRA (); LM3 (); SX (), they would be very cumbersome to use in practice due to the intricate. If the analytic solutions to these equations are available, then one can establish an elegant analytic dipole amplitude, like IIM amplitudeIIM () in LO BK, which would be more convenient to use in phenomenology than the numerical one.
In our previous publicationWX (), we analytically solved the rcBK equation in the saturation region. We found that the running coupling corrections modify the matrix a lot as compared to the fixed coupling case. The exponent in matrix is linear decrease with rapidity in the case of running coupling while the exponent in matrix decreases quadratically with rapidity in LO BK case, which indicates that the running coupling slows down the rapidity evolution of the scattering amplitude. In this work, we analytically solve the full NLO and NLL BK equations in the saturation region and obtain their analytic results for the matrix at high energies. We find that the DLA and LO BK equations have the same analytic solutions in the saturation region. Moreover, the full NLO and NLL BK equations also have the same solutions in the saturation regime. Interestingly, the analytic solution of the full NLL BK equation shows that the linear decrease with rapidity of the exponent in the solution to the rcBK equation is replaced by rapidity raised to power of 3/2, which indicates the gluon loop contributions partially compensate the reduction from the quark loops in the saturation region. To test these outcomes, we numerically solve the above equations with focusing on the saturation region. The numerical results support the analytic findings.
Ii The LO, running coupling BK equations and their solutions
To motivate the higher order corrections and for comparison with more refined results which shall be obtained later, in this section we recall the low level (LO and rc) BK equations which describe the rapidity evolution of the matrix, = 1  , of a dipole scattering off a target which may be another dipole, a hadron or a nucleus. The scattering amplitude is small in the dilute target region, while at the dense target regime it approaches the unitarity limit (), which indicates the matrix approaching zero at the saturation region. For simplicity we analytically solve the BK equations in terms of the matrix instead of the scattering amplitude in the following studies.
ii.1 The LO BK equation and its analytic solution
The BK equation describes the rapidity evolution of the matrix of a quarkantiquark dipole (or a projectile) with a quark leg at transverse coordinate and an antiquark at transverse coordinate scattering off a hadronic target. To see how the matrix evolving with rapidity , we take boosting a small increment . If we put the evolution at dipole framework and let the target fixed, then the dipole has a probability to emit a gluon due to the increment of M (),
(1) 
where is transverse coordinate of the emitted gluon, and . For convenient later on calculations, we denote , and as the sizes of parent and of the new daughter dipoles produced by the evolution, respectively. In the large limit, the quarkantiquarkgluon state can be viewed as two daughter dipoles – one of the dipoles consists of the initial quark and the antiquark part of the gluon while the other dipole is constituted by the quark part of the gluon and the initial antiquark. In the dipole framework, the increment of matrix, , due to the boosting of can be written by multiplying the probability with the matrixM ()
(2) 
where expresses two daughter dipoles simultaneous scattering off the target. The last term in (2) describes the scattering of a single dipole on the target. Eq.(2) is an integrodifferential equation and gives the scattering amplitude at all rapidities . A few more words about the interpretation of the Eq.(2) are in order. The energy evolution follows from the gluon emission becoming possible when the dipole is boosted to higher rapidity. Integration over the rapidity interval is corresponding to multiple gluon emissions and thus we have large number of dipoles in the projectile wave function. Eq.(2) is an infinite hierarchy of coupled evolution equations due to lower number of dipoles scattering always coupling to higher number of dipoles, for instance two dipoles scattering coupling to three dipoles. Therefore, Eq.(2) is almost impossible for direct applications to phenomenology, since is not known. To get a closed equation, we need to employ the large limit, the can be simplified to
(3) 
then we obtain the BK equationBa (); K ()
(4) 
Now let’s analytically solve the BK equation in the saturation region where the target density is large and the scattering is strong with , thus . The nonlinear term in (4) is negligible. In the logarithmic regime of integration, the Eq.(4) becomes
(5) 
Note that we work in the saturation region in which the dipole size is much larger than characteristic size , the is saturation momentum which controls the separation between dilute and dense regimes. Thus, the lower and upper bounds of integration in (5) are and , respectively. The factor 2 in (5) comes from the symmetry of the two regions dominating the integral, either from or , see Figure 1. Via performing the integrations over transverse coordinate and rapidity, we obtain the analytic solution of the BK equationHM (); LT ()
(6) 
with and . One can see that the exponent has a quadratic rapidity decrease.
The analytic solution to the BK equation in the saturation region has been found by Levin and Tuchin LT (). The reason why we have gone through such a detailed “derivation” of (6) since one of the main purposes of this study is to show how the leading order matrix is modified by the higher order contributions like running coupling, full NLO and large double transverse logarithmic corrections.
ii.2 The running coupling BK equation and its analytic solution
Beyond the leading logarithmic approximation, the authors in Refs.KW (); Bnlo () calculated the higher order perturbative corrections by including fermion (quark) bubble diagrams which bring in a factor of and modify the BK evolution kernel. Once including the higher order corrections there are two impacts on the Feynman diagrams which describe the dipole structure generated under evolution. First, the propagator of the emitted gluon, which comes from the emission of the original parent dipole when the parent dipole boosts to higher rapidities, is now dressed with quark loops in contrast to leading order or fixed coupling case, which will modify the emission probability of gluon (or BK kernel), but keeping the leading order interaction terms untouched. Second, the quarkantiquark pair is added to the evolved wave function, which not only modifies the BK kernel but also changes the interaction structure of the evolution equation. In this case the evolution of the matrix is proportional to the product of two matrixes of the newly created dipoles. Concerning two above aspects and resumming to all orders, one can get the running coupling BK equationJA ()
(7) 
where is the point of subtraction in the coordinate space, and is the transverse coordinate of the emitted gluon. In the large limit, the emitted gluon can be viewed as a quarkantiquark pair with a quark leg at transverse coordinate and an antiquark leg at transverse coordinate . The first line on the r.h.s of (7) refers to as the ’running coupling’ contribution and resums all power of corrections to all orders. It has the same structure as the leading order BK equation but with modified kernel due to the running coupling corrections. The modified kernel has two kinds of expressionsBnlo (); KW (), which depends on the scheme choice (see JA () for more discussions on the scheme choice). We adopt the choice derived by Balitsky in Ref.Bnlo (). The kernel of the running coupling BK can be written asBnlo ()
(8) 
We would like to note that the kernel, , in the second line on the r.h.s of (7) is not discussed in this work (see JA () for relevant discussion), since it corresponds to quadratic terms of matrix which are negligible in current study.
In the saturation regime in which the interaction between dipole and hadron is very strong, , and unitarity corrections become important, the quadratic terms in (7) can be neglected in which case one needs only keep the second term in the first line on the r.h.s of (7). The evolution equation including running coupling is given by
(9) 
The expert reader will recognize that Eq.(9) has the same form as the LO BK Eq.(5) with only the kernel modified by running coupling corrections. We use similar approach as the solution to Eq.(5) to analytically solve the rcBK equation.
In the saturation region, the main contribution to the integration on the r.h.s of (9) comes from either or region, see Figure 1. We choose the first regime, . Under this approximation, the kernel becomes
(10)  
with the running coupling at one loop accuracy
(11) 
Substituting the simplified kernel (10) into (9), we can get rcBK equation in the saturation region as:
(12) 
whose solution isWX ()
(13)  
with
(14) 
Comparing to the matrix in the LO case, the matrix in (13) is modified by the running coupling corrections. The quadratic rapidity decrease in exponent at LO BK (see Eq. (6)) is replaced by linear decrease with rapidity in (13), which indicates that the scattering amplitude is slowed down by the running coupling corrections. This result is in agreement with theoretic expectationsBnlo ().
We wish to note that although the running coupling kernel depends on the scheme choice, if one chooses the KovchegovWeigert scheme instead of Balitsky scheme and will find that the analytic solution to the rcBK is independent of the scheme choice. In other words, the running coupling Balitsky and KovchegovWeigert evolution equations give the same analytic solution in the saturation region. These two running coupling evolution equations are equivalent to each other in the saturation region.
Iii The full NLO, DLA, full NLL BK equations and their solutions
In the last section, we discuss the running coupling corrections to the LO BK equation. However, the running coupling (or bubble chain of quark loops) is not the only higher order perturbative corrections to the LO BK equation. The gluon loops also bring a large contributions to the LO BK equation. The full NLO BK equation can be derived by including contributions from the quark and gluon loops as well as from the tree gluon diagrams with quadratic and cubic nonlinearitiesBC ().
iii.1 Full NLO BK equation and its analytic solution
We analytically solve the full NLO BK evolution equation derived by Balitsky and Chirilli in Ref.BC (), which can be written as:
(15) 
where is the renormalization scale, is the first coefficient of the function and is the number of flavors. The full NLO BK equation, Eq.(III.1), shows several remarkable characteristics as compared to the LO BK equation.

The kernel in the first integration receives a correction of order compared to the LO BK kernel. Especially, the terms proportional to come from quark loop contributions (or running coupling corrections). The rest terms contain contributions from gluon loops, which include corrections enhanced by ’double logarithms’. The double logarithms become large in the small dipole size compared to the saturation scale , which can lead to negative value of scattering amplitude LM ().

The second and third integration terms are of order . The double integrations over and refer to two additional partons (except for original quark and antiquark) with transverse coordinate and at the time of scattering.

The second integration does not involve the number of quark flavor, , which indicates that both two additional partons are gluons. This integration term describes the following sequence of evolutions: First the original quarkantiquark dipole () radiates a gluon with transverse coordinate , producing two daughter dipoles () and (); following the dipole () emits another gluon with transverse coordinate , thus creating the dipoles () and (). The cubic term describes the two daughter gluons simultaneous interacting with the target, while the quadratic term expresses the gluon at absent at the time of scattering.

The third integration is proportional to the number of quark flavor, , which indicates that the additional partons are a quark and an antiquark at the time of scattering. Otherwise, the evolutions are the same as the second integration case.
The running coupling terms in (III.1) are scale dependent, we use the scheme developed in Bnlo () to rewrite running coupling part and replace all terms proportional to by the Balitsky running coupling. The kernel in the first integration becomes:
(16) 
We analytically solve Eq.(III.1) in the saturation region where the scattering is strong, thus the scattering amplitude approaches unit (), in other words the matrix closes to zero at this regime. Therefore, the nonlinear terms in (III.1) can be neglected in the saturation region. One finds that the full NLO BK equation reduces to
(17) 
Note that in the saturation region the full NLO BK equation (17) has similar structure as the rcBK equation (9) as well as LO BK equation (5), but with the kernel modified by quark and gluon loops. To solve Eq.(17), we can choose saturation region either or , see Figure 1. In order to compare with running coupling case, we select the first region, . The becomes
(18)  
In the above equation, we neglect the second term in the first line on the r.h.s due to much larger than . The double logarithmic term also can be neglected. One can see that the second term in the second line on the r.h.s of (18) comes from the gluon loop correction.
Inserting the reduced kernel (18) into (17), the full NLO BK equation in the saturation region can be written as
(19) 
and has the following solution
(20)  
with and . Interestingly, the linear rapidity decrease in the exponent, (13), is now replaced by the rapidity raised to the power of 3/2, which indicates that gluon loop corrections compensate part of the decrease made by quark loops in the saturation region. In other words, the evolution speed of the scattering amplitude is slowed down in the full NLO BK compared to LO BK (quadratic rapidity decrease, see Eq.(6)), but it is decreasing not as much as the rcBK case.
iii.2 Full NLL, DLA BK equations and their analytic solutions
The numerical study of the full NLO BK equation found that its solution can turn to negativeLM (), which is not meaningful in physics. This disaster can be traced back to the double logarithmic term which has to be resummed. By an explicit computation of the Feynman diagrams in light cone perturbation theory, the authors in Ref.IM () established an effective method to resum double logarithms to all orders and got a DLA BK equation which solves the problem of negative value of the scattering amplitude. The numerical solutions to the DLA BK equation show that the resummation of large double logarithmic corrections plays a significant role in stabilizing and slowing down the evolution. However, the DLA equation does not include contributions from quark and gluon loop corrections. It is possible to promote the DLA equation to full NLL accuracy by adding the NLO BK corrections derived by Balitsky and Chirilli in BC (). The kernel including double transverse logarithms (DTL) and full NLO corrections can be written as:
(21)  
with the DLA kernelIM ()
(22) 
where is the Bessel function, and . Note that the explicit double logarithmic term as the one in (III.1) is disappeared in (21), since it has already resummed into the kernel .
Let’s analyze the DLA kernel in the saturation region where or , see Figure 1. To keep the consistency with previous two sections, we work in the region at this section. The tends towards zero when , thus , which corroborates the statement that the double logarithm only important in phasespace where the scattering is weakIM (). In the saturation region, the double logarithmic corrections can be neglected. The becomes
(23) 
which is the same as the simplified full NLO kernel in (18). The same kernel implies that the NLL BK equation in saturation region shares the exact the same solution as the full NLO BK equation. Its solution can be written as
(24)  
Let’s discuss the solution to the DLA BK equation derived in IM ()
(25) 
As we know the in the saturation region, the solution of the DLA BK should be the same as the LO BK equation, see (6). We can see that the double transverse logarithmic correction is trivial in the saturation region and can be neglected, which is in agreement with the theoretical expectationsIM (). In the next section, we numerically solve the LO, DLA, rc and NLL BK equations to test the analytic results.
Iv Numerical Solution
To test the outcomes what we get in the above sections for the analytic solutions to the LO, full NLO, and full NLL BK equations. We numerically solve these equations in this section. These BK equations are integrodifferential equations, they can numerically straightforward to solve on a lattice. The impact parameter dependence of scattering amplitude, , is neglected, therefore the amplitude does not depend on angle, . We can solve at discrete values of transverse dipole size. The BK equation can be viewed as a set of differential equations which can be solved by using GNU Scientific Library (GSL). The GSL computes the numerical integrals using adaptive integration routines, interpolates data points using the cubic spline interpolation, and solves the differential equations based on RungeKutta method.
The initial condition is needed to solve the BK equation. In this study, we take the initial condition given by the McLerranVenugopalan model (MV)MV ():
(26) 
Note that as we are interested in the rapidity dependence of the amplitude and the DLA impact on the amplitude in the saturation region, the initial condition is not resumed when solving the NLL BK equation in current work.
For the qualitative properties of BK equation studied in this work, the precise value of the running coupling constant is not important. To specific, we evaluate the running coupling in terms of a popular oneloop expression
(27) 
with . The value of the fudge factor can be obtained by fitting to the HERA data, its value is taken from Ref.JL (). We freeze the coupling to a fixed value for larger dipole size, , to regularize the infrared behavior.
The left hand panel of Fig.2 shows the evolution of the dipole scattering amplitude with fixed coupling for 5 different rapidities. The evolution is faster for fixed coupling than the DLA BK equation in the nonsaturation region, especially for large rapidities. While it seems that the DLA does not contribute to the evolution in the saturation region, as one can see the numerical LO and DLA results overlap each other from the inner zooming in diagram. These numerical outcomes are in agreement with the analytic results obtained in previous sections.
To demonstrate the change of rapidity dependence of scattering amplitude from running coupling to full NLL case, we numerically solve the NLL BK equation which includes corrections from full NLO and double logarithmic resummation. The right hand panel of Fig.2 shows the solutions of rc and full NLL BK equations. The evolution is faster for the running coupling than the NLL in the nonsaturation region where the dipole size is small and the scattering is weak. This outcome implies that the evolution is further slowed down on top of running coupling by the corrections from gluon loops and double logarithmic resummation in the nonsaturation region. However, in the saturation region where the DLA can be neglected, we get opposite result, the evolution is slower for the running coupling than the NLL, one can see that the value of the amplitude in NLL larger than the one in running coupling case for the same rapidity from the inner diagram on the right hand panel of Fig.2. This result can be traced back to the gluon loop contributions in the full NLO case, see the second line in (III.1). The gluon loop corrections take the rapidity dependence of the exponent in matrix from linear to rapidity raised to the power of 3/2 (see Eqs.(13) and (20)), which indicate that the evolution in running coupling case is relative slow compared to the NLL one, but the evolutions for both of them are slower than the LO BK. There is a crossover between the saturation and the nonsaturation regions, where the dipole size is around . Beyond this crossover region, the evolution is dominated by the NLL corrections which slow down the evolution.
Acknowledgements.
WX thanks the Department of Physics of The University of Colorado Boulder for the hospitality during the early stages of this work. This work is supported by National Natural Science Foundation of China under Grant No.11305040, the Education Department of Guizhou Province Talent Fund under Grant No.[2015]5508, the Science and Technology Department of Guizhou Province Fund Nos.[2015]2114.References
 (1) I. Balitsky, Nucl. Phys. B463 (1996) 99;
 (2) J. JalilianMarian, A. Kovner, A. Leonidov, and H. Weigert, Nucl. Phys. B504 (1997) 415;
 (3) J. JalilianMarian, A. Kovner, A. Leonidov, and H. Weigert, Phys. Rev. D59 (1998) 014014;
 (4) E. Iancu, A. Leonidov, and L.D. McLerran, Nucl. Phys. A692 (2001) 583;
 (5) E. Ferreiro, E. Iancu, A. Leonidov, and L.D. McLerran, Nucl. Phys. A703 (2002) 489;
 (6) Yu.V. Kovchegov, Phys. Rev. D60 (1999) 034008; ibid. D61 (1999) 074018.
 (7) J.L. Albacete, N. Armesto, J.G. Milhano, and C.A. Salgado, Phys. Rev. D80 (2009) 034031;
 (8) J.L. Albacete, N. Armesto, J.G. Milhano, P. QuirogaArias, and C.A. Salgado, Eur. Phys. J. C71 (2011) 1705;
 (9) T. Lappi, and H. Mntysaari, Phys. Rev. D88 (2013) 114020;
 (10) G.A. Chirilli, B.W. Xiao, and F. Yuan, Phys. Rev. Lett. 108 (2012) 122301;
 (11) Yu.V. Kovchegov and H. Weigert, Nucl. Phys. A784 (2007) 188;
 (12) I. Balitsky, Phys. Rev. D75 (2007) 014001;
 (13) I. Balitsky, and G.A. Chirilli, Phys. Rev. D77 (2008) 014019;
 (14) T. Lappi, and H. Mntysaari, Phys. Rev. D91 (2015) 074016;
 (15) E. Iancu, J. Madrigal, A. Mueller, G. Soyez, and D. Triantafyllopoulos, Phys. Lett. B744 (2015) 293;
 (16) T. Lappi, and H. Mntysaari, Phys. Rev. D93 (2016) 094004;
 (17) E. Iancu, J. Madrigal, A. Mueller, G. Soyez, and D. Triantafyllopoulos, Phys. Lett. B750 (2015) 643;
 (18) J.L. Albacete, Nucl. Phys. A957 (2017) 71;
 (19) A.H. Rezaeian, Phys. Lett. B718 (2013) 1058;
 (20) A. Stasto, B.W. Xiao, and F. Yuan, Phys. Lett. B716 (2012) 430;
 (21) E. Iancu, K. Itakura, and S. Munier, Phys. Lett. B590 (2004) 199;
 (22) W. Xiang, Phys. Rev. D79 (2009) 014012;
 (23) A.H. Mueller, Nucl. Phys. B415 (1994) 373;
 (24) A.H. Mueller, hepph/0111244;
 (25) E. Levin, and K. Tuchin, Nucl. Phys. B573 (2000) 83;
 (26) J.L. Albacete, and Yu.V. Kovchegov, Phys. Rev. D75 (2007) 125021;
 (27) L.D. McLerran, and R. Venugopalan, Phys. Rev. D49 (1994) 2233;