Solving The High Energy Evolution Equation
Including Running Coupling Corrections
Javier L. Albacete∗ and Yuri V. Kovchegov†
Department of Physics, The Ohio State University
Columbus, OH 43210,USA
April 2007
Abstract
We study the solution of the nonlinear BK evolution equation with the recently calcu-
lated running coupling corrections [1, 2]. Performing a numerical solution we confirm the
earlier result of [3] (obtained by exploring several possible scales for the running coupling)
that the high energy evolution with the running coupling leads to a universal scaling
behavior for the dipole-nucleus scattering amplitude, which is independent of the initial
conditions. It is important to stress that the running coupling corrections calculated
recently significantly change the shape of the scaling function as compared to the fixed
coupling case, in particular leading to a considerable increase in the anomalous dimension
and to a slow-down of the evolution with rapidity. We then concentrate on elucidating
the differences between the two recent calculations of the running coupling corrections.
We explain that the difference is due to an extra contribution to the evolution kernel,
referred to as the subtraction term, which arises when running coupling corrections are
included. These subtraction terms were neglected in both recent calculations. We evalu-
ate numerically the subtraction terms for both calculations, and demonstrate that when
the subtraction terms are added back to the evolution kernels obtained in the two works
the resulting dipole amplitudes agree with each other! We then use the complete running
coupling kernel including the subtraction term to find the numerical solution of the re-
sulting full non-linear evolution equation with the running coupling corrections. Again
the scaling regime is recovered at very large rapidity with the scaling function unaltered
by the subtraction term.
∗e-mail: albacete@mps.ohio-state.edu
†e-mail: yuri@mps.ohio-state.edu
http://arxiv.org/abs/0704.0612v2
1 Introduction
Recently our understanding of the linear Balitsky–Fadin–Kuraev–Lipatov (BFKL) [4, 5] and
non-linear Jalilian-Marian–Iancu–McLerran–Weigert–Leonidov–Kovner (JIMWLK) [6–13] and
Balitsky-Kovchegov (BK) [14–18] small-x evolution equations in the Color Glass Condensate
[6–29] has been improved due to the completion of the calculations determining the scale of the
running coupling in the evolution kernel in [1,2,30,31]. The calculations in [1,2] proceeded by
including αsNf corrections into the evolution kernel and by then completing Nf to the complete
one-loop QCD beta-function by replacing Nf → −6πβ2. Calculation of the αsNf corrections
is particularly easy in the s-channel light-cone perturbation theory formalism [32, 33] used to
derive the BK and JIMWLK equations: there αs Nf corrections are solely due to chains of
quark bubbles placed onto the s-channel gluon lines.
The analytical results of [1,2] are not very concise and could not have been guessed without
an explicit calculation. After finding αsNf corrections, the obtained contributions had to be
divided into the running coupling part, which has a form of a running coupling correction
to the leading-order (LO) JIMWLK or BK kernel, and into the ”subtraction piece”, which
would bring in new structures into the kernel. Such separation had to be done both in [1] and
in [2]. Unfortunately, there appears to be no unique way to perform this separation: it is not
surprising, therefore, that it was done differently in both papers [1, 2]. This resulted in two
different running coupling terms, shown below in Eqs. (35) and (36) along with Eqs. (7) and
(8). Such a discrepancy has led to a misconception in the community that the calculations
of [1] and of [2] disagree at some fundamental level.
Indeed to compare the results of [1] and [2] one has to undo the separation into the running
coupling and subtraction terms: combining both terms one should compare full kernels of the
evolution equation obtained in [1, 2]. There is another more physical reason to perform such
comparison: in principle, there is no small parameter making the subtraction term smaller than
the running coupling term and thus justifying neglecting the former compared to the latter.
Even the labeling of one term as “running coupling” piece is somewhat misleading, since it may
give an impression that the neglected subtraction term has no running coupling corrections in
it. As was shown in [31] both terms actually contribute to the running coupling corrections to
the BFKL equation (if one uses the separation of [2] to define the terms).
In this paper we perform numerical analysis of the BK evolution equation with the αsNf
corrections resummed to all orders and with Nf completed to the QCD beta-function, Nf →
−6πβ2, with β2 given in Eq. (20). We first solve the BK equation keeping the running coupling
term only, with the kernels given by Eqs. (7) and (8). Indeed the solutions we find this way
are different from each other. We then evaluate the subtraction terms for both cases and show
that inclusion of subtraction terms puts the results of [1] and [2] in perfect agreement with each
other! We complete our analysis by solving the BK equation with the full kernel including both
the running coupling and subtraction terms.
This work is structured as follows. Section 2 begins with Sect. 2.1 in which we review
the αsNf corrections to the dipole scattering amplitude evolution equation recently derived
in [1, 2] and the subtraction method employed in both works to separate the running coupling
contributions from the subtraction terms. We discuss the scheme dependence of the running
coupling terms introduced by this separation. We proceed in Sect. 2.2 by deriving the explicit
expressions for the subtracted terms. The calculation is based on the results of [2]. Our
analytical results are summarized in Sect. 2.3, where we give the explicit final expression
for the kernel of the subtraction term in Eq. (39), which, combined with Eq. (38) gives us the
subtraction terms (40) and (41) for the subtractions performed in [1] and in [2] correspondingly.
In Sect. 3 we explain the numerical method we use to solve the evolution equations. We
also list the initial conditions used, along with the definition of the saturation scale employed.
Throughout the paper we will avoid the important question of the Landau pole and the contri-
bution of renormalons to small-x evolution. As we explain in Sect. 3, we will simply “freeze”
the running coupling at a constant value in the infrared. For a detailed study of the renormalon
effects in the non-linear evolution we refer the readers to [30].
Our numerical results are presented in Sect. 4. By solving the evolution equations with the
running coupling term only in Sect. 4.1 we show that the resulting dipole amplitude differs
significantly from the fixed coupling case. We also observe that the amplitude obtained by
solving the equation obtained in [2] is very close to the result of solving the BK evolution
with a postulated parent-dipole running of the coupling constant. Both these amplitudes are
quite different from the solution of the equation derived in [1], as one can see from Fig. 4.
In spite of that, all three evolution equations studied (the ones derived in [1], [2] and the
parent-dipole running coupling model) give approximately identical scaling function for the
dipole amplitude at high rapidity, as demonstrated in Fig. 5 in Sect. 4.2. It is worth noting
that, as can be seen from Fig. 6, the anomalous dimension we extracted from our solution
is γ ≈ 0.85, which is different from the fixed coupling anomalous dimension of γ ≈ 0.64.
The former anomalous dimension also appears to disagree with the predictions of analytical
approximations to the behavior of the dipole amplitude with running coupling from [34–38]. In
Sect. 4.3 we numerically evaluate the subtraction terms for both [1] and [2] and show that their
contributions are important, as shown in Fig. 7. However, subtraction terms decrease with
increasing rapidity, such that at high enough rapidities their relative contribution becomes
small (see Fig. 8). In Fig. 9 we show that inclusion of subtraction terms makes the results of [1]
and [2] agree with each other. Finally, the numerical solution of the full (all orders in αs β2)
evolution equation including both the running coupling and subtraction terms is performed in
Sect. 4.4. The results are shown in Fig. 10. All the main features of the evolution with the
running coupling are preserved in the full solution: the growth of the dipole amplitude and
of the saturation scale with rapidity is slowed down (for the latter see Fig. 11). The scaling
function of Fig. 5 is unaltered by the subtraction term, as shown in Fig. 12.
We summarize and discuss our main conclusions in Sect. 5.
2 Scheme dependence
2.1 Inclusion of running coupling corrections: general concepts
The BK evolution equation for the dipole scattering matrix reads
∂S(x0, x1; Y )
d2z K(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] , (1)
where
K(x0, x1, z) =
r21 r
is the kernel of the evolution. Here transverse two-dimensional vectors x0 and x1 denote the
transverse coordinates of the quark and the anti-quark in the parent dipole, while z is the
position of the gluon produced in one step of evolution [39–42]. We have introduced the
notation r = x0 − x1, r1 = x0 − z, r2 = z − x1 for the sizes of the parent and of the new
(daughter) dipoles created by one step of the evolution. The notation r ≡ |r| for all the 2-
dimensional vectors will be also employed throughout the rest of the paper. Eq. (1) admits
a clear physical interpretation: the original parent dipole, when boosted to higher rapidities,
may emit a new gluon which, in the large-Nc limit, is equivalent to a quark-antiquark pair.
Thus, the original dipole splits into two new dipoles sharing a common transverse coordinate:
the transverse position of the emitted gluon, z. The nonlinear term in the right hand side of
Eq. (1) accounts for either one of the two new dipoles interacting with the target, along with
the possibilities of only one dipole interacting or no interaction at all, while the subtracted
linear term reflects virtual corrections. The kernel of the evolution is just the probability of
one gluon emission calculated at leading logarithmic accuracy in αs ln(1/xB), where xB is the
fraction of momentum carried by the emitted gluon [39–42].
Under the eikonal approximation the dipole scattering matrix off a hadronic target at a
fixed rapidity is given by the average over the hadron field configurations of Wilson lines V
calculated along fixed transverse coordinates (those of the quark and of the antiquark). More
specifically
S(x0, x1; Y ) =
V (x0)V
†(x1)
〉 . (3)
Hence, the integrand of Eq. (1) can be regarded as a three point function in the sense that
the gluon fields of the target are evaluated at three different transverse positions, those of the
original quark and antiquark plus the one of the emitted gluon.
However, the inclusion of higher order corrections to the evolution equation via all order
resummation of αsNf contributions as recently derived in [1, 2] brings in new physical chan-
nels that modify the three point structure of the leading-log equation. The dipole structure
generated under evolution by diagrams like the one depicted in Fig. 1A (for more detailed dis-
cussion of the diagrammatic content of the high order corrections see [2]) is identical to the one
previously discussed for the leading order equation, the only novelty being that the propagator
of the emitted gluon is now dressed with quark loops, modifying the emission probability but
leaving untouched the interaction terms. On the contrary, diagrams like the one in Fig. 1B
in which a quark-antiquark pair (rather than a gluon) is added to the evolved wave function
modify the interaction structure of the evolution equation. The evolution of the parent dipole
scattering matrix driven by these kind of terms is proportional to the scattering matrix of the
two newly created dipoles (the one formed by the original quark and the new antiquark and
vice versa), ∼ S(x0, z1)S(z2, x1). This term depends on four different transverse coordinates,
i.e., it is a four point function and, therefore, its contribution to the evolution equation cannot
be accounted for by a mere modification of the emission kernel of the leading order equation.
To discuss in more detail the modifications introduced by the high order corrections, we
find it useful to rewrite the evolution equation in the following, rather general way
∂S(x0, x1; Y )
= F [S(x0, x1; Y )] (4)
where F is a functional of the dipole scattering matrix which for the original derivation of the
equation is given by the right hand side of Eq. (1). In general it can be decomposed into two
Figure 1: Schematic representation of the diagrams contributing to quark-NLO evolution.
pieces
F [S] = R [S]− S [S] . (5)
The first term, R, which we will call the ’running coupling’ contribution, gathers all the
higher order in αsNf corrections to the evolution that can be recast in a functional form that
looks identical to the leading order one but with a modified kernel, K̃, which includes all the
terms setting the scale for the running coupling:
R [S(x0, x1; Y )] =
d2z K̃(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] . (6)
The second term, S, henceforth referred to as the ’subtraction’ contribution, encodes those
contributions that depart from the three point structure of the leading-log equation. The
explicit derivation and expressions for this term are presented in the next section. The relative
minus sign between the two terms in Eq. (5) has been introduced for latter convenience.
Importantly, the decomposition of F into running coupling and subtraction contributions,
although constrained by unitarity arguments, is not unique. Two different separation schemes
have been proposed in [1, 2]. They are both based on a similar strategy, sketched in Fig. 2,
that can be summarized as follows. The newly created quark-antiquark pair added to the
wave function in the diagrams Fig. 1B is shrunk to a point, called the subtraction point, by
integrating out one of the coordinates in the dipole-qq̄ wave function, rendering the previously
discussed four point nature of these contributions into a three point one. This integrated
three point contribution is added to the running coupling contribution, whereas the original
four point term minus its integrated version are assigned to the subtraction contribution. The
divergence between the two approaches stems from the choice of the subtraction point. In
the subtraction scheme proposed by Balitsky in [1] the subtraction point is chosen to be the
transverse coordinate of either the quark, z2, or the antiquark, z1. The kernel for the running
coupling functional, Eq. (6), obtained in this way is
K̃Bal(r, r1, r2) =
Nc αs(r
r21 r
. (7)
On the other hand, in the subtraction procedure followed in [2] (which we will refer to as
KW) the zero size quark-antiquark pair is fixed at the transverse coordinate of the gluon,
x 0 x 0
Figure 2: Schematic representation of the subtraction procedure.
z = αz1 + (1 − α)z2, where α is the fraction of the gluon’s longitudinal momentum carried by
the quark, yielding the following expression for the kernel of the running coupling contribution:
K̃KW(r, r1, r2) =
− 2 αs(r
1)αs(r
αs(R2)
r1 · r2
r21 r
+ αs(r
, (8)
where
R2(r, r1, r2) = r1 r2
r1·r2
. (9)
As we shall discuss later, the scheme dependence originated by the choice of the subtraction
point is substantial and has an important effect in the solutions of the evolution equation when
only the running contribution is taken into account.
In our numerical study we will also consider the following ad hoc prescription for the kernel
of the running coupling functional in which the scale for the running of the coupling is set to
be the size of the parent dipole
K̃pd(r, r1, r2) =
Nc αs(r
r21 r
. (10)
This prescription is useful as a benchmark used to compare with previous numerical [3] and
analytical works [34, 35, 43] where this ansatz was used.
2.2 Derivation of the subtraction term
We begin by considering the NLO contribution to the kernel of the JIMWLK and BK evolution
equations with the s-channel gluon splitting into a quark-antiquark pair, which then interacts
with the target, as shown on the left hand side of Fig. 3. The contribution of this diagram has
been calculated in [2]. The resulting JIMWLK kernel is [2]
KNLO1 (x0, x1; z1, z2) = 4Nf
(2π)2
(2π)2
(2π)2
(2π)2
e−iq·(z−x0)+iq
′·(z−x
)−i(k−k′)·z
q2q′2
(1− 2α)2q · k k′ · q′ + q · q′ k · k′ − q · k′ k · q′
k2 + q2α(1− α)
k′2 + q′2α(1− α)
2α (1− α) (1− 2α)
k2 + q2α(1− α)
k′2 + q′2α(1− α)
k · q
k′ · q′
4α2 (1− α)2
k2 + q2α(1− α)
k′2 + q′2α(1− α)
 . (11)
The momentum labels in the above equation are explained on the left hand side of Fig. 3.
If k1 and k2 are the transverse momenta of the quark and of the antiquark in the produced
pair as shown in Fig. 3, then the transverse momentum of the gluon is q = k1 + k2. The
other transverse momentum we use is k = k1(1 − α) − k2α, where α is the fraction the of
gluon’s “plus” momentum carried by the quark, α ≡ k1+/(k1+ + k2+). The prime over the
transverse momentum denotes the momentum of the same particle in the complex conjugate
amplitude. For instance q′ is the momentum of the s-channel gluon in the complex conjugate
amplitude. Finally, z1 and z2 denote the transverse coordinates of the quark and the antiquark.
In Eq. (11) we use z12 = z1−z2 (the transverse separation between the quark and the antiquark)
and z = α z1 + (1− α) z2 (the transverse coordinate of the gluon).
Figure 3: A lowest order leading-Nf NLO correction which gives rise to the subtraction term
is shown on the left. The same diagram with the gluon lines “dressed” by chains of fermion
bubbles, as shown on the right, gives the full (resumming all powers of αµNf ) contribution to
the subtraction term. Calculation of the subtraction term is pictured in Fig. 2.
To obtain the BK kernel from Eq. (11) one should sum over all possible emissions of the
gluon off the quark and antiquark lines in the incoming dipole both in the amplitude and in
the complex conjugate amplitude, which is accomplished by
KNLO1 (x0, x1; z1, z2) = CF
m,n=0
(−1)m+n KNLO1 (xm, xn; z1, z2). (12)
Below we will label the JIMWLK kernel by calligraphic letter K and the corresponding BK
kernel by K.
The contribution of the kernel from Eq. (12) to the right hand side of the NLO version of
Eq. (1) is given by the following term
d2z1 d
1 (x0, x1; z1, z2)S(x0, z1, Y )S(z2, x1, Y ) (13)
with αµ the bare coupling.
As shown in Fig. 2, at the NLO level, the subtraction term introduced in Eq. (5), is then
defined by
SNLO[S] = α2µ
d2z1 d
2z2 K
1 (x0, x1; z1, z2)
× [S(x0, w, Y )S(w, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )], (14)
where w is the point of subtraction in the transverse coordinate space. In [1] it was chosen to
be equal to the transverse coordinate of either the quark or the antiquark,
w = z1 or w = z2, (15)
as both choices lead to the same subtraction term SBalNLO[S]:
SBalNLO[S] =
d2z1 d
2z2 K
1 (x0, x1; z1, z2)
× [S(x0, z1, Y )S(z1, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] . (16)
In [2] the subtraction point was chosen to be the transverse coordinate of the gluon z,
w = z = α z1 + (1− α) z2. (17)
This leads to the following subtraction term, which we denote SKWNLO[S]:
SKWNLO[S] =
d2z1 d
1 (x0, x1; z1, z2)
× [S(x0, z, Y )S(z, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] . (18)
Indeed the complete kernel in Eq. (5) is independent of the choice of w. However, since the
subtraction term of Eq. (14) was neglected both in [1] and in [2], different choices of w led to
different expressions for the remaining running coupling part R[S], i.e., to different answers as
far as investigations in [1] and in [2] were concerned. Different choice of w is the main source
of the discrepancy of final answers of [1] and [2], though it does not imply any disagreement in
the full expression (5).
Our goal in this Section is to evaluate KNLO1 (x0, x1; z1, z2) from Eq. (11) including the
running coupling corrections. The s-channel light-cone perturbation theory formalism makes
such inclusion simple [2]: all we have to do is include infinite chains of quark bubbles on the
gluon lines in the amplitude and in the complex conjugate amplitude, as depicted on the right
hand side of Fig. 3. Performing calculations similar to those done in [2] one arrives at
K ❣1 (x0, x1; z1, z2) = 4Nf
(2π)2
(2π)2
(2π)2
(2π)2
e−iq·(z−x0)+iq
′·(z−x
)−i(k−k′)·z
q2q′2
(1− 2α)2q · k k′ · q′ + q · q′ k · k′ − q · k′ k · q′
k2 + q2α(1− α)
k′2 + q′2α(1− α)
2α (1− α) (1− 2α)
k2 + q2α(1− α)
k′2 + q′2α(1− α)
k · q
k′ · q′
4α2 (1− α)2
k2 + q2α(1− α)
k′2 + q′2α(1− α)
1 + αµβ2 ln
q2 e−5/3
1 + αµβ2 ln
q′2 e−5/3
) (19)
where K ❣1 denotes the kernel with the running coupling corrections resummed to all orders.
Just like in [2, 31], here we will use MS renormalization scheme. Inclusion of fermion bubble
chains generated two denominators at the end of Eq. (19), which is its only difference from
Eq. (11). Here
11Nc − 2Nf
. (20)
Now we have to perform the transverse momentum integrals in Eq. (19). First we expand
the denominators at the end of Eq. (19) into a power series and rewrite Eq. (19) as
K ❣1 (x0, x1; z1, z2) = 4Nf
n,m=0
(−αµβ2)n+m
(2π)2
(2π)2
(2π)2
(2π)2
e−iq·(z−x0)+iq
′·(z−x
)−i(k−k′)·z
q2q′2
(1− 2α)2q · k k′ · q′ + q · q′ k · k′ − q · k′ k · q′
k2 + q2α(1− α)
k′2 + q′2α(1− α)
2α (1− α) (1− 2α)
k2 + q2α(1− α)
k′2 + q′2α(1− α)
k · q
k′ · q′
4α2 (1− α)2
k2 + q2α(1− α)
k′2 + q′2α(1− α)
λ=λ′=0
where we have defined µ2 = µ2
e5/3 to make the expressions more compact.
Indeed we can not always expand the denominators of Eq. (19) into a geometric series
employed in Eq. (21), but one has to remember that the summation of bubble chain diagrams
shown on the right side of Fig. 3 gives one the geometric series. Hence the geometric series
come first: later they are absorbed into the denominators shown in Eq. (19), which is an
approximation not valid for all q and q′. Therefore, by keeping the geometric series in Eq. (21)
we are not making any approximations. In general, in what follows we are not going to keep
track of the issues of convergence of perturbation series. The contribution of renormalons to
non-linear small-x evolution was thoroughly investigated in [30] and was found to be significant
at low Q2. We refer the interested reader to [30] for more details on this issue.
Using the following formulas
(2π)2
e−ik·z
k2 + q2
K0(q z) (22)
(2π)2
e−ik·z
k2 + q2
q K1(q z) (23)
we can now perform the k- and k′-integrals in Eq. (21). Integrating over the angles of q and q′
as well yields
K ❣1 (x0, x1; z1, z2) =
(2π)4
n,m=0
(−αµβ2)n+m
dq q dq′ q′
z212 |z − x0| |z − x1|
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
× J1(q |z − x0|)K1(z12 q
α ᾱ) J1(q
′ |z − x1|)K1(z12 q′
α ᾱ) + 2α ᾱ (α− ᾱ)
z12 · (z − x0)
z12 |z − x0|
J1(q |z − x0|)K1(z12 q
α ᾱ) J0(q
′ |z − x1|)K0(z12 q′
α ᾱ)
z12 · (z − x1)
z12 |z − x1|
J0(q |z − x0|)K0(z12 q
α ᾱ) J1(q
′ |z − x1|)K1(z12 q′
α ᾱ)
+ 4α2 ᾱ2 J0(q |z − x0|)K0(z12 q
α ᾱ) J0(q
′ |z − x1|)K0(z12 q′
α ᾱ)
λ=λ′=0
. (24)
We have defined
ᾱ = 1− α (25)
for brevity. Now the integrals over q and q′ can be carried out to give
K ❣1 (x0, x1; z1, z2) =
(2π)4
n,m=0
(−αµβ2)n+m
z212 µ
2 α ᾱ
)λ+λ′
Γ2(1 + λ)
Γ2(1 + λ′)
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
(1 + λ) (1 + λ′)
z812 (α ᾱ)
1 + λ, 2 + λ; 2;−|z − x0|
α ᾱ z212
1 + λ′, 2 + λ′; 2;−|z − x1|
α ᾱ z212
α− ᾱ
z12 · (z − x0)
(1 + λ)F
1 + λ, 2 + λ; 2;−|z − x0|
α ᾱ z212
1 + λ′, 1 + λ′; 1;−|z − x1|
α ᾱ z212
z12 · (z − x1)
1 + λ, 1 + λ; 1;−
|z − x0|2
α ᾱ z212
(1 + λ′)F
1 + λ′, 2 + λ′; 2;−
|z − x1|2
α ᾱ z212
1 + λ, 1 + λ; 1;−
|z − x0|2
α ᾱ z212
1 + λ′, 1 + λ′; 1;−
|z − x1|2
α ᾱ z212
λ=λ′=0
. (26)
Unfortunately further simplification of the expression in Eq. (26) is impossible without approx-
imations. The series resulting from summation over n and m are likely to be divergent due to
renormalons. As we mentioned before, here we neglect the renormalon problem referring the
reader to [30]. Similar to how it was done in [2] we are not going to attempt to resum the series
exactly: instead we will calculate the next-to-leading order terms and assume that with a good
accuracy they give us the scale(s) of the running coupling constant. This procedure is similar
to the well-known prescription due to Brodsky, Lepage and Mackenzie [44].
Using the Taylor-expansions of hypergeometric functions
F (1 + λ, 2 + λ; 2; z) =
− λ 1
1 + ln(1− z) + 1
ln(1− z)
+ o(λ2). (27)
F (1 + λ, 1 + λ; 1; z) =
ln (1− z) + o(λ2) (28)
after some algebra we obtain
K ❣1 (x0, x1; z1, z2) =
[α (z1 − x0)2 + ᾱ (z2 − x0)2] [α (z1 − x1)2 + ᾱ (z2 − x0)2] z412
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
1− αµ β2 ln
R2T (x0)µ
+ o(α2µ)
1− αµ β2 ln
R2T (x1)µ
+ o(α2µ)
+ 2α ᾱ (α− ᾱ) z212
z12 · (z − x0)
1− αµ β2 ln
R2T (x0)µ
+ o(α2µ)
1− αµ β2 ln
R2L(x1)µ
+ o(α2µ)
+ z12 · (z − x1)
1− αµ β2 ln
R2L(x0)µ
+ o(α2µ)
1− αµ β2 ln
R2T (x1)µ
+ o(α2µ)
+4α2 ᾱ2 z412
1− αµ β2 ln
R2L(x0)µ
+ o(α2µ)
1− αµ β2 ln
R2L(x1)µ
+ o(α2µ)
In arriving at Eq. (29) we employed functions RT (x) and RL(x), which have dimensions of
transverse coordinates and are defined by
R2T (x)µ
4 e−2γ−5/3
[α (z1 − x)2 + ᾱ (z2 − x)2]µ2MS
α ᾱ z212
(z − x)2
α (z1 − x)2 + ᾱ (z2 − x)2
α ᾱ z212
R2L(x)µ
4 e−2γ−5/3
[α (z1 − x)2 + ᾱ (z2 − x)2]µ2MS
α (z1 − x)2 + ᾱ (z2 − x)2
α ᾱ z212
The subscripts T and L stand for transverse and longitudinal gluon polarizations which give
rise to the two different functions under the logarithm.
Recombining the series in Eq. (29) into physical running couplings finally yields
α2µK ❣1 (x0, x1; z1, z2) =
[α (z1 − x0)2 + ᾱ (z2 − x0)2] [α (z1 − x1)2 + ᾱ (z2 − x0)2] z412
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
R2T (x0)
R2T (x1)
+ 2α ᾱ (α− ᾱ) z212
z12 · (z − x0) αs
R2T (x0)
R2L(x1)
+z12 · (z − x1)αs
R2L(x0)
R2T (x1)
+ 4α2 ᾱ2 z412 αs
R2L(x0)
R2L(x1)
with the physical running coupling in the MS scheme given by
αs(1/R
1 + αµβ2 ln
R2 µ2
) . (33)
Eq. (32) is the contribution to the JIMWLK evolution kernel of the resummed diagram on
the right hand side of Fig. 3.
2.3 Brief summary of analytical results
Let us briefly summarize our analytical results. The non-linear small-x evolution equation with
the running coupling corrections included reads
∂S(x0, x1; Y )
= R [S]− S [S] . (34)
The first term on the right hand side of Eq. (34) is referred to as the running coupling
contribution. It was calculated independently in [1] and in [2]: the results of those calculations
are given above in Eqs. (7) and (8) correspondingly, which have to be combined with Eq. (6)
to obtain
RBal [S] =
d2z K̃Bal(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] (35)
RKW [S] =
d2z K̃KW(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] . (36)
One notices immediately that RBal [S] calculated in [1] is different from RKW [S] calculated
in [2] due to the difference in the kernels K̃Bal and K̃KW in Eqs. (7) and (8). However that
does not imply disagreement between the calculations of [1] and [2]: after all, it is the full
kernel on the right of Eq. (34), R [S]− S [S], that needs to be compared. To do that one has
to calculate the second term on the right hand side of Eq. (34).
The second term on the right hand side of Eq. (34) is referred to as the subtraction contri-
bution. It is given by
S[S] = α2µ
d2z1 d
2z2K ❣1 (x0, x1; z1, z2) [S(x0, w, Y )S(w, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )]
with the resummed BK kernel
K ❣1 (x0, x1; z1, z2) = CF
m,n=0
(−1)m+n K ❣1 (xm, xn; z1, z2). (38)
The resummed JIMWLK kernel K ❣1 (xm, xn; z1, z2) is given by Eq. (32), along with Eqs. (30)
and (31) defining the scales of the running couplings. In the numerical solution below we will
replace Nf → −6πβ2 in its prefactor, obtaining
α2µK ❣1 (x0, x1; z1, z2) = −
[α (z1 − x0)2 + ᾱ (z2 − x0)2] [α (z1 − x1)2 + ᾱ (z2 − x0)2] z412
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
R2T (x0)
R2T (x1)
+ 2α ᾱ (α− ᾱ) z212
z12 · (z − x0) αs
R2T (x0)
R2L(x1)
+z12 · (z − x1)αs
R2L(x0)
R2T (x1)
+ 4α2 ᾱ2 z412 αs
R2L(x0)
R2L(x1)
This substitution is the same as for all other factors ofNf . The same substitution was performed
in [2] to calculate the running coupling term. In fact, as was shown in [31], the linear part of the
subtraction term (calculated using the prescription of [2]) contributes to the running coupling
corrections to the BFKL equation. Therefore, in that case, the factor of Nf in front of Eq. (32)
is definitely a part of the beta-function. Hence the replacement Nf → −6πβ2 is justified even
in the subtraction term. Once again, in the numerical solution below we will use Eq. (39) along
with Eq. (38) in Eq. (37) to calculate the subtraction term S[S].
Substituting w = z1 (or, equivalently, w = z2) in Eq. (37) would yield the subtraction term
SBal[S] =α2µ
d2z1 d
2z2K ❣1 (x0, x1; z1, z2)
× [S(x0, z1, Y )S(z1, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] (40)
which has to be subtracted from RBal [S] calculated in [1] and given by Eq. (35) to obtain the
complete evolution equation resumming all orders of αsNf in the kernel.
Substituting w = z = α z1 + (1− α) z2 in Eq. (37) yields
SKW[S] =α2µ
d2z1 d
2z2K ❣1 (x0, x1; z1, z2)
× [S(x0, z, Y )S(z, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] (41)
which has to be subtracted from RKW [S] calculated in [2] and given in Eq. (36) again to obtain
the complete evolution equation resumming all orders of αsNf in the kernel. We checked
explicitly by performing analytic calculations that the two evolution equations obtained this
way agree at the NLO and NNLO. Below we will check the agreement of the two calculations
to all orders by performing a numerical analysis of the solutions of these equations.
The above discussion demonstrates that the separation of the evolution kernel into the
running coupling and subtraction pieces, as done in Eq. (34), is somewhat artificial, and has
no small parameter justifying one or another separation prescription. Therefore, the small-x
evolution equation including all running coupling (or, more precisely, αsNf ) corrections should
combine both terms in Eq. (34). Below we will solve such evolution equation numerically to
obtain the full small-x evolution with the running coupling.
3 Numerical setup and initial conditions
In our numerical study we consider the translational invariant approximation in which the
scattering matrix is independent of the impact parameter of the collision, i.e., S = S(r, Y ).
To solve the integro-differential equations, corresponding to the BK equation with running
coupling we employ a second-order Runge-Kutta method with a step size in rapidity ∆Y = 0.1.
We discretize the variable |r| into 800 points equally separated in logarithmic space between
rmin = 10
−8 and rmax = 50. Throughout this paper, the units of r will be GeV
−1, and those
of Qs will be GeV. All the integrals have been performed using improved adaptative Gaussian
quadrature methods. The accuracy of this numerical method has been checked in [3] to be
better than a 4% in all the r range.
We consider three different initial conditions for the dipole scattering amplitude, N(r, Y ) =
1− S(r, Y ). The first one is taken from the McLerran-Venugopalan (MV) model [22, 23]:
NMV (r, Y = 0) = 1− exp
. (42)
where a constant term has been added to the argument of the logarithm in the exponent in
order to regularize it for large values of r. The other two initial conditions are given by
NAN (r) = 1− exp
−(r Q
, (43)
with γ = 0.6 and γ = 0.8. These two last initial conditions will be referred hereinafter as AN06
and AN08 respectively. The interest in this ansatz, reminiscent of the Golec-Biernat–Wusthoff
model [45], is that the small-r behavior NAN ∝ r2γ corresponds to an anomalous dimension
1− γ of the unintegrated gluon distribution at large transverse momentum. (AN labels initial
conditions with anomalous dimension.) Our choices γ = 0.6 and γ = 0.8 can be motivated
a posteriori by the observation that the anomalous dimension of the evolved BK solution for
running coupling lies in between those two values and the one for the MV initial condition,
γ ≈ 1 (see Section 4.2). Thus, the choice of distinct initial conditions allows us to better
track the onset of the expected asymptotic universal behavior that is eventually reached at
high energies and to study the influence of the pre-asymptotic, non-universal corrections to the
solutions of the evolution equations. To completely determine our initial conditions, we set
Q′s = 1 GeV at Y=0 in Eqs. (42) and (43) and put Λ = 0.2 GeV. Although Q
s is normally
identified with the saturation scale, our definition of the saturation scale through the rest of
the paper will be purely pragmatical and given by the condition
N(r = 1/Qs(Y ), Y ) = κ, (44)
with κ = 0.5. We have checked that this choice of κ, albeit arbitrary, does not affect any of the
major conclusions to be drawn in the rest of the paper.
Finally, in order to avoid the Landau pole and to regularize the running coupling at large
transverse sizes we stick to the following procedure: for small transverse distances r < rfr, with
rfr defined by αs(1/r
fr) = 0.5, the running coupling is given by the one loop expression
αs(1/r
β2 ln
r2 Λ2
) (45)
with Nf = 3 and Λ = 0.2 GeV, whereas for larger sizes, r > rfr, we “freeze” the coupling at a
fixed value αs = 0.5. A detailed study of the role of Landau pole in non-linear small-x evolution
is given in [30].
4 Results
In this Section, we discuss our numerical results and how they compare to previous numerical
work and analytical estimates.
4.1 Running coupling
Fig. 4 shows the solutions of the evolution equation when only the running coupling contribution
is taken into account, i.e., neglecting the subtraction term in Eq. (34), for different initial
conditions and for the three schemes considered in this work: Balitsky’s, given by Eqs. (7) and
(35), KW, given by Eqs. (8) and (36), and the the ad hoc parent dipole implementation of the
running coupling, shown in Eq. (10).
MV init. cond.     Y=0,5,15,30KW
Balitsky
parent dipole
init. cond.
AN08 init. cond.     Y=0,5,15,30
AN06 init. cond.     Y=0,5,15,30
r )-1(GeV
Figure 4: Solutions of the BK equation at rapidities Y=0, 5, 15 and 30 (curves are labeled from
right to left) for the three running coupling schemes considered in this work: KW (solid line),
Balitsky (dashed line) and parent dipole (dashed-dotted lines). The initial conditions are MV
(top), AN08 (middle) and AN06 (bottom).
As previously observed in [3, 46], the most relevant effect of including running coupling
corrections in the evolution equation is a considerable reduction in the speed of the evolution
with respect to the fixed coupling case. This is a common feature of the different running
coupling schemes studied here and of other phenomenological ones considered in the literature
(a detailed comparison between the solutions for fixed coupling evolution and for parent dipole
running coupling can be found e.g. in [3]). This is not a surprising result, since a generic effect
of the running of the coupling is to suppress the emission of small transverse size dipoles, which
is the leading mechanism driving the evolution.
However, despite this common feature of the running coupling solutions, significant dif-
ferences are found between the solutions obtained under different schemes as we infer from
Fig. 4. In particular, the evolution is much faster with the KW prescription than with that of
Balitsky. Equivalently, the KW prescription yields a stronger growth of the saturation scale
with rapidity/energy than Balitsky’s. Moreover, the solutions obtained when the parent dipole
prescription is used lay much closer to those obtained within the KW scheme than to the ones
obtained when Balitsky’s scheme is applied, contrary to what was suggested in [1]. As argued
before, the differences observed in the solutions obtained using the two subtraction schemes
are entirely due to neglecting the subtraction contribution and reflect the arbitrariness of the
separation procedure.
4.2 Geometric scaling
It has been found in previous analytical [34, 36, 47] and numerical studies on the solutions
of the BK equation at leading order [3, 48–50] and for different heuristic implementations of
next-to-leading order corrections [3,46], including the parent dipole prescription for the running
coupling also considered in this work, that the solutions of the evolution equation at high enough
rapidities are no longer a function of two separate variables r and Y , but rather they depend
on a single scaling variable, τ = r Qs(Y ). This feature of the evolution, commonly referred to
as geometric scaling, is an exact property of the solutions for fixed coupling evolution due to
the conformal invariance of the leading-log kernel, and has become one of the key connections
between the saturation based formalisms and the phenomenology of heavy ion collisions and
deep inelastic scattering experiments [51–57].
It can be seen from Fig. 5 that the solutions of the BK equation with the running coupling
terms discussed in the previous section also exhibit the property of scaling, in agreement with
the analytical study carried out in [38], shown by the fact that the rescaled high rapidity
solutions lay on a single curve which is independent of both the running coupling scheme and
of the initial condition. The scaling behavior of the solution is observed in the whole τ range
studied in this work, including the saturation region, τ > 1. The tiny deviations from a pure
scaling behavior observed in Fig. 5 may be attributed to the fact that the full asymptotic
behavior is reached at even larger rapidities (Y & 80, [3]) than those achieved by the numerical
solution performed in this work.
Remarkably, the scaling function for both KW and Balitsky’s scheme coincides with the one
obtained with the parent dipole prescription, up to the above mentioned scaling violations. It
has been observed in [3, 46, 48] that the scaling function differs significantly in the fixed and
running coupling cases. Following that work, and to make a more quantitative study of the
scaling property, we fitted our solutions to the functional form [34]
f(τ) = a τ 2 γ
ln τ 2 + b
, (46)
with a, b and γ free parameters, within a fixed window below the saturation region, τ ∈
[10−5, 0.1]. Noticeably, at large enough rapidities the whole fitting window lays within the
geometric scaling window proposed in [36]: (Λ/Qs(Y )) < τ < 1, where Λ is some initial scale.
-210 -110 1
MV init. cond.   Y=0,40
Balitsky
parent dipole
fixed coupling
init. cond.
τ -110 1 10
AN06 init. cond.   Y=0,40
Figure 5: Solutions of the BK equation at rapidities Y=0 and 40 for KW (solid line), Balitsky
(dashed line) and parent dipole (dashed-dotted lines) schemes plotted versus the scaling variable
τ = rQs(Y ). The asymptotic solution obtained with fixed coupling αs = 0.2 at Y = 40 in [3]
is shown (black dashed-dotted line) for comparison. The initial conditions are MV (left) and
AN06 (right).
The value of γ extracted from the fits at rapidity Y = 40 lays in between γ ∼ 0.8 and γ ∼ 0.9.
This conclusion holds for the three initial conditions used here: the anomalous dimension
seems to converge to some intermediate value, in agreement with the value found in [3], for
asymptotic running coupling solutions (γ ∼ 0.85 at Y = 70). This result for anomalous
dimension is very far away from the value obtained in [3] for fixed coupling solutions (γ ∼ 0.64
at Y = 70) and from the predicted anomalous dimension for both running and fixed coupling
solutions from analytical studies of the equation based on saddle point techniques [34–38],
γc = χ(γc)/χ
′(γc) = 0.6275, where χ is the leading-log BFKL kernel.
It might be argued that the numerical value of the anomalous dimension extracted from our
fits is conditioned by the choice of the fitting function and by the fitting interval. Actually, it
was shown in [3] that the solutions of the evolution could be well fitted by other functional forms,
including the double-leading-log solution of BFKL, within a similar fitting region to the one
considered in this work. On the other hand, several phenomenological parameterizations of the
solution of the evolution have been proposed in [54–57] and have successfully confronted HERA
10 -410
10 -210 -110 1
 2*0.6τ ~ 
Fixed coupling
 2*0.85τ ~ 
Running coupling
Figure 6: Asymptotic solutions (Y=40) of the evolution equation for running coupling (solid
line) and fixed coupling with αs = 0.2 (dashed line). A fit to a power-law function aτ
2γ in the
region τ ∈ [10−6, 10−2] yields γ ≈ 0.85 for the running coupling solution and γ ≈ 0.6 for the
fixed coupling one.
and RHIC experimental data. There, the dipole scattering amplitude at arbitrary rapidity is
assumed to be given by a functional form analogous to our ansatz for the initial condition
Eq. (43), but allowing for geometric scaling violations by replacing γ → γ(r, Y ). The value
of the anomalous dimension at r = 1/Qs and/or for Y → ∞ is fixed to be the BFKL saddle
point, γc ∼ 0.63 (the saddle point value considered in [57] is slightly different, γ ∼ 0.53), while
the value γ = 1 is recovered in the limit r → ∞ at any finite rapidity. The success of these
phenomenological works supports the claim that the anomalous dimension of the solution is
given by the BFKL saddle point, in agreement with the above mentioned analytical predictions.
However, the relevant values of momenta probed at current phenomenological applications are
very distinct from the fitting region considered here. For example, the inclusive structure
function measured in HERA is fitted in [54,56] within the region 0.045 GeV2 < Q2 < 45 GeV2,
whereas charged hadron pt spectra in dAu collisions is well reproduced by [55–57] in the region
1 GeV < pt < 4.5 GeV. Note that, for both sets of data, the measured regions overlap with
the deeply saturated domain of the solution. On the contrary, our fitting region 10−5 < τ <
1 corresponds to values of momenta ∼ 10Qs(Y ) < pt < 105Qs(Y ) (always well above the
saturation scale), with Qs(Y = 40) ∼ 500÷1000 GeV for the different running coupling schemes
considered and, therefore, has no overlap with the kinematic regions measured experimentally,
since we scrutinize a momentum region strongly shifted to the ultraviolet compared to currently
available data. Moreover, it should be noticed that the rapidity interval covered by both sets
of experimental data is ∆Y < 4 in both cases, while we study the solutions of the evolution
(Y=0)]MVSub[N
-410 -310 -210 -110 1
(Y=30)]MVSub[N
(Y=0)]
Sub[N
10 -210 -110 1
(Y=30)]
Sub[N
(Y=0)]
Sub[N
10 -210 -110 1
(Y=30)]
Sub[N
Figure 7: Subtraction contribution calculated in the KW scheme (triangles) and in Balitsky’s
(stars). The trial functions correspond to the solutions of the evolution under Balitsky run-
ning coupling scheme at rapidities Y = 0, 30 for MV (left), AN08 (center) and AN06 initial
conditions.
at asymptotic rapidities, Y ∼ 40. We have checked that shifting our fitting region to larger
values of τ (smaller momentum) would bring the value of γ extracted from our fits closer to the
saddle point BFKL one, since the transition from the ultraviolet region to the deeply saturated
domain of the scaling solution is realized by a locally less steeper function (see Figs. (5) and
(12)). Therefore, there is no contradiction at all between the success of the phenomenological
parameterizations of the solutions and the results reported here.
With the above clarifications we reach the following conclusion: the asymptotic scaling
solutions corresponding to fixed and running coupling evolution are intrinsically different in
the whole r-range. This is emphasized in Fig. 6, where we represent the scaling solutions in
a log scale for τ < 1. It is clear that the tail of the distribution falls off with decreasing τ
much steeper for the running coupling solution than for the fixed coupling one. A fit to a
pure power-law function, f = a τ 2γ , in the region τ ∈ [10−6, 10−2] yields γ ∼ 0.85 for the
running coupling and γ ∼ 0.61 for the fixed coupling solution. The differences between fixed
and running coupling solutions at τ > 1 are evident from Fig. 5. This is a puzzling result that
remains to be understood from purely analytical methods.
(Y=0)]MV[N
10 -210 -110 1
(Y=30)]MV[N
(Y=0)]
10 -210 -110 1
(Y=30)]AN08[N
(Y=0)]
10 -210 -110 1
(Y=30)]AN06[N
Figure 8: Ratio of the subtraction over the running terms, D(r, Y ) = S[N(r, Y )]/R[N(r, Y )],
calculated in both KW (triangles) and Balitsky (stars) schemes for MV (left), AN08 (middle)
and AN06 (right) initial conditions at rapidities Y=0 (top) and Y=30 (bottom).
4.3 Subtraction Term
Before attempting to solve the complete evolution equation, and in order to gain insight in the
nature and structure of the subtraction contribution, we first evaluate the subtraction functional
for both Balitsky, Eq. (40), and KW, Eq. (41), schemes using a set of trial functions for S which
we choose to consist of the solutions of the evolution equation with the running coupling in
Balitsky’s scheme at different rapidities and of the three initial conditions considered above in
this work.
Two main remarks can be made about our results, shown in Fig. 7:
i) For all the trial functions considered in this work, the subtraction contribution is much
larger in the KW scheme than in Balitsky’s. A plausible explanation for this is that Balitsky’s
subtraction contribution, Eq. (40), when expanded in terms of dipole scattering amplitudes,
N = 1 − S, reduces to a sum of non-linear terms, since all the linear terms in the expansion
cancel each other due to the z1 ↔ z2 symmetry of the kernel, whereas in the KW case no such
cancellation happens and the subtraction contribution, Eq. (41), also includes linear terms,
which are dominant over the non-linear ones in the non-saturated domain where N ≪ 1.
ii) The subtraction contribution S has the same sign as the running coupling contribution
R in the whole τ range which, together with the relative minus sign assigned to the subtraction
term in Eq. (34), implies that the proper inclusion of the subtraction term reduces the value
of the functional that governs the evolution, F . In other words: the subtraction contribution
tends to systematically slow down the evolution, as we shall explicitly confirm in the next
subsection.
To better quantify the size of the subtraction contribution, we plot the ratio D(r, Y ) ≡
S[N(r, Y )]/R[N(r, Y )] in Fig. 8. At Y = 0, the relative weight of the subtraction contribution
with respect to the running one within the KW scheme and for a MV initial condition goes
from a D ∼ 0.4 at small τ to D ∼ 1 at τ ∼ 1. The same ratio for the Balitsky scheme takes
significantly smaller values: it goes from D ∼ 0.1 at small τ to D ∼ 0.4 for τ ∼ 1. As the
evolved solutions get closer to the scaling function, i.e. for larger rapidities, the r dependence of
the ratio becomes flatter and its overall normalization goes down to an approximately constant
value D ∼ 0.15 for the KW scheme and D ∼ 0.025 for that of Balitsky. This behavior remains
unaltered when going from rapidity Y = 20 to Y = 30, which suggests that the ratio may
saturate to a fixed value in the asymptotic region.
-210 -110 1
KWRun
BalRun
KW(Run-Sub)
(Run-Sub)
(Y=0)MVN
τ -210 -110 1
(Y=30)MVN
Figure 9: Total kernel F = R−S calculated under Balitsky’s scheme, Eqs (7) and (40), (solid
line) and under the KW scheme, Eqs (8) and (41), (dashed line). The overlap of the two lines
shows the agreement between the two calculations. Triangles stand for the running coupling
term calculated in the KW approach, RKW, while stars stand for the running coupling term
under Balitsky’s scheme, RBal. The trial functions N(r, Y ) correspond to the solution of the
evolution with only running coupling under Balitsky’s scheme at Y=0 (left) and Y=30 (right)
for a MV initial condition.
-110 1
1.2 MV i.c.
Y=0,3,10
KWRun
BalRun
F=Run-Sub
init. cond.
r )-1(GeV
-110 1
scaling function i.c.
Y=0,3,10
r )-1(GeV
Figure 10: Solutions of the complete (all orders in αs β2) evolution equation given in Eq. (34)
(solid lines), and of the equation with Balitsky’s (dashed lines) and KW’s (dashed-dotted)
running coupling schemes at rapidities Y = 0, 5 and 10. Left plot uses MV initial condition.
The right plot employs the initial condition given by the dipole amplitude at rapidity Y = 35
evolved using Balitsky’s running coupling scheme and with r-dependence rescaled down such
that Qs = Q
s = 1 GeV.
Finally, we have checked that combining the subtraction and running coupling contributions
for both schemes adds up to the same result. This is shown in Fig. 9, where we plot the value
of the total functional F = R − S calculated under the KW scheme (Eqs. (8) and (36) for
the running coupling term, R, and Eq. (41) for the subtraction term, S) and under Balitsky’s
scheme (Eqs. (7) and (35) for the running coupling term and Eq. (40) for the subtraction term).
The two results coincide within the estimation of the numerical accuracy previously discussed.
The agreement between the two results is better in the small-τ region, where the two curves lay
almost on top of each other. In the saturation region, τ & 1, the agreement is slightly worse,
although the differences between the values of F calculated in both schemes is still much less
than the differences between the running coupling terms themselves. This slight remaining
disagreement between the Balitsky’s and KW prescriptions may also be due to inaccuracies in
a Fourier transform of a geometric series performed in arriving at Eq. (39). This result serves
as a cross-check of our numerical method and as an additional confirmation of the agreement
of the independent calculations derived in [1, 2].
0 1 2 3 4 5 6 7 8 9
MV i.c.KWRun
BalRun
F=Run-Sub
)2(GeV
1 2 3 4 5 6 7 8 9
scaling function i.c.
Figure 11: Saturation scale corresponding to the solutions plotted in Fig. 10.
4.4 Complete running coupling BK equation
In this section we calculate the solutions of the complete evolution equation, Eq. (34), including
both the running and subtraction terms obtained by the all-orders αs Nf resummation and by
the Nf → −6πβ2 replacement. Since the numerical evaluation of the subtraction contribution
at each point of the grid and each step of the evolution would require an exceedingly large
amount of CPU time consumption, the strategy followed to include it in the evolution equation
consists of calculating such contribution only in a small set of grid points at each step of the
evolution, which we fixed at n = 16, between the points r1 and r2, which are determined at each
step of the evolution by the conditions N(Y, r1) = 10
−9, and N(Y, r2) = 0.99, and then using
power-law interpolation and extrapolation to the other points of the grid. Both the running
and subtracted terms are calculated within Balitsky scheme. This procedure is motivated by
the fact that, as discussed in the previous section, the subtraction contribution can be regarded
as a small perturbation with respect to the running coupling term within Balitsky’s scheme and
by the fact that it is a rather smooth function that can be well fitted by a power-law function in
most of the r-range. The accuracy of this procedure has been checked by doubling the number
of points at which the subtraction contribution is calculated at each step of the evolution, i.e.
by setting n=32. At Y=2, the differences between the solutions obtained with the two above
mentioned choices for n were less than a 8% in the tail of the solution, r < r1, and less than a
3% for r > r1.
The results of the evolution calculated in this way and using MV and rescaled asymptotic
running coupling solution (Y=35) as initial conditions are plotted in Fig. 10. They confirm
-210 -110 1 10
Y=10F=Run-Sub
(Y=35)  MV BALi.c.=N
Figure 12: Rescaled solutions given by the complete αs β2-evolution equation (solid line) and for
KW (dashed-dotted line) and Balitsky’s (dashed line) running coupling schemes at Y = 10. The
initial condition corresponds to the dipole amplitude at rapidity Y = 35 evolved using Balitsky’s
running coupling scheme and with r-dependence rescaled down such that Qs = Q
s = 1 GeV.
the expectations raised in the previous Subsection: the inclusion of the subtraction terms
considerably slows down the evolution with respect to the sole consideration of the running
coupling contributions. Moreover, the reduction in the speed of the wave front is much larger
for the KW scheme than for that of Balitsky one for both initial conditions. However, the closer
the initial condition is to the asymptotic running coupling scaling function, the smaller are the
effects of the subtraction contribution. These features can be better quantified by inspecting
the rapidity dependence of the saturation scale generated by the evolution, plotted in Fig. 11.
At rapidity Y = 10 the ratio of the saturation scale Qs yielded by the KW scheme to Qs given
by the complete αs β2-evolution equation is a factor of ∼ 2.5 for the MV initial condition and
a factor of ∼ 2.1 for the asymptotic running coupling initial condition. At the same rapidity,
the ratio of the saturation scale obtained under Balitsky’s scheme to Qs corresponding to the
complete αs β2-evolution is ∼ 1.25 for the MV initial condition and ∼ 1.15 for the scaling
function initial condition. Thus, in spite of the smallness of the ratio of the subtraction terms
to the running coupling contributions at high rapidity, which is ∼ 0.025 for Balitsky’s and
∼ 0.15 for KW scheme at Y = 30 (see bottom plots in Fig. 8), the proper inclusion of the
subtraction term results in fairly sizable effects in the solutions of the evolution equation.
Finally, we notice that the scaling behavior of the solution is not affected by the subtraction
term. This is seen in figure Fig. 12, where we evolve starting from an initial condition already
close to the running coupling scaling function and plot the solutions of the evolution equation
obtained with just running coupling terms (see Section 4.1) and the solution of the complete
αs β2-evolution at rapidity Y = 10. It is clear that, within the numerical accuracy, no departure
from the scaling behavior is observed. Therefore the main effect of a proper consideration of the
subtraction term is the one of reducing the speed of the evolution. It does not violate or modify
the geometric scaling property of the solutions established in Section 4.2. In our understanding
geometric scaling appears to persist when the running coupling effects are included because,
at high enough rapidity Qs(Y ) ≫ Λ, such that the new (from the LO standpoint) momentum
scale Λ introduced by the running coupling can be safely neglected. Hence the dynamics is
again characterized by a single momentum scale Qs(Y ). At the same time running coupling
does modify the evolution kernel, leading to a different shape of the scaling function.
5 Conclusions
In this paper we have taken into account all corrections to the kernels of the non-linear JIMWLK
and BK evolution equations containing powers of αs Nf . We reiterated the fact that the sep-
aration of the resulting kernel resumming all powers of αs Nf into the running coupling and
subtraction parts, as done in the previous calculations of [1, 2], is not justified parametrically.
We have then performed numerical analysis with the following conclusions.
• First we solved the evolution equations derived in [1] and [2] keeping only the running
coupling part or the evolution kernel and neglecting the subtraction term. Comparing to
the results for fixed coupling obtained in [3] we confirmed the conclusion reached in [3]
that the growth with rapidity is substantially reduced when running coupling corrections
are included. The results for three different initial conditions are shown in Fig. 4. We
observe that the solution of the equation derived in [2] differs significantly from that
derived in [1], but agrees (with good numerical accuracy) with the solution of the BK
evolution equation with the coupling running at the parent dipole size. (The latter is just
a model of the running coupling not resulting from any calculations, which we plot for
illustrative purposes.) We also observe that at sufficiently high rapidity both equations
from [1] and from [2] give us the same scaling function for the dipole amplitude N(r, Y )
as a function of r Qs(Y ), which is also in agreement with the scaling function given by
the parent dipole running, as shown in Fig. 5. The fact that the scaling is preserved when
the running coupling corrections are included was previously established in [3], though
for models of running coupling only. The shape of the scaling function is very different
from that obtained from the fixed coupling evolution equations. In particular, we found
that for dipole sizes below 0.1/Qs the anomalous dimension of the scaling function in the
running coupling case becomes γ ≈ 0.85 (see Fig. 6). This is different from the result of
several analytical estimates [34–38], which expect the anomalous dimension not to change
when running coupling corrections are included and to remain at its fixed coupling value
of γ ≈ 0.63.
• We have then evaluated the subtraction term for both calculations performed in [1] and
[2]. We demonstrated that subtracting the subtraction terms from the running coupling
terms makes the full answer agree for both calculations of [1] and [2], as shown in Fig. 9
for the right hand side of the evolution equation. It turns out that the subtraction
term SBal[S], which has to be subtracted from the result of [1], is systematically smaller
than SKW[S], to be subtracted from the result of [2], over the whole rapidity range
studied here. This implies that the result of [1] should have a smaller correction than
the result of [2] and is thus closer to the full answer. The subtraction terms SBal[S]
and SKW[S] are plotted in Fig. 7 as functions of the dipole size r for different values of
rapidity. Their relative contributions to the evolution kernel are shown in Fig. 8, where
we plotted the subtraction functional divided by the running coupling functional. From
those figures we conclude that both the magnitude of these extra terms and their relative
contribution to the evolution kernel decrease with increasing rapidity. Hence, while at
”moderate” rapidities (the ones closer to realistic experimental values) the subtraction
term is important for both calculations [1, 2], it becomes increasingly less important at
asymptotically large rapidities. The physics is easy to understand: the subtraction terms
are o(α2s), while the running coupling part of the kernel is o(αs). Hence, if we suppose
that the effective value of the coupling is given by its magnitude at the saturation scale
Qs(Y ), then, as rapidity increases, the coupling would decrease, making the subtraction
term much smaller than the running coupling term. Indeed, while at asymptotically high
rapidities the assumption of [1,2] that the subtraction term could be neglected is justified,
making the results of [1] and [2] agree with each other, for rapidities relevant to modern
days experiments the subtraction term is numerically important.
• With the last conclusion in mind we continued by numerically solving the full evolution
equation resumming all powers of αs Nf in the evolution kernel, which now would combine
both the running coupling and the subtraction terms. The five-dimensional integral in
the subtraction term (37) made obtaining this solution rather difficult. The outcome of
the calculation is shown in Fig. 10. All the main conclusions stated above were again
confirmed by the solution of the full equation. At asymptotically high rapidity scaling
regime is recovered, as can be seen from Fig. 12. As the subtraction term is less important
in that regime, the scaling function appears to be the same as in the case of having only
the running coupling term in the kernel. The anomalous dimension again turns out
to be γ ≈ 0.85, in disagreement with the analytical expectations of [34–38]. However,
the scaling of the saturation scale with rapidity appears to be in agreement with the
expectations of analytical work of [34, 35, 38], as shown in Fig. 11.
We conclude by observing that the knowledge of the non-linear small-x evolution equation
with all the running coupling corrections included brings us to an unprecedented level of pre-
cision allowing for a much more detailed comparison with experiments than was ever possible
before.
Acknowledgments
We would like to thank Heribert Weigert for many informative and helpful discussions at
the beginning of this work. A portion of the performed work was motivated by stimulating
discussions with Robi Peschanski, which we gratefully acknowledge.
This research is sponsored in part by the U.S. Department of Energy under Grant No. DE-
FG02-05ER41377. This work was supported in part by an allocation of computing time from
the Ohio Supercomputer Center.
References
[1] I. I. Balitsky, Quark Contribution to the Small-x Evolution of Color Dipole, Phys. Rev. D
75 (2007) 014001, [hep-ph/0609105].
[2] Y. Kovchegov and H. Weigert, Triumvirate of Running Couplings in Small-x Evolution,
Nucl. Phys. A 784 (2007) 188–226, [hep-ph/0609090].
[3] J. L. Albacete, N. Armesto, J. G. Milhano, C. A. Salgado, and U. A. Wiedemann,
Numerical analysis of the Balitsky-Kovchegov equation with running coupling:
Dependence of the saturation scale on nuclear size and rapidity, Phys. Rev. D71 (2005)
014003, [hep-ph/0408216].
[4] E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, The Pomeranchuk singularity in
non-Abelian gauge theories, Sov. Phys. JETP 45 (1977) 199–204.
[5] Y. Y. Balitsky and L. N. Lipatov Sov. J. Nucl. Phys. 28 (1978) 822.
[6] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, The BFKL equation from
the Wilson renormalization group, Nucl. Phys. B504 (1997) 415–431, [hep-ph/9701284].
[7] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, The Wilson renormalization
group for low x physics: Towards the high density regime, Phys. Rev. D59 (1998) 014014,
[hep-ph/9706377].
[8] J. Jalilian-Marian, A. Kovner, and H. Weigert, The Wilson renormalization group for low
x physics: Gluon evolution at finite parton density, Phys. Rev. D59 (1998) 014015,
[hep-ph/9709432].
[9] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Unitarization of gluon
distribution in the doubly logarithmic regime at high density, Phys. Rev. D59 (1999)
034007, [hep-ph/9807462].
[10] A. Kovner, J. G. Milhano, and H. Weigert, Relating different approaches to nonlinear
QCD evolution at finite gluon density, Phys. Rev. D62 (2000) 114005, [hep-ph/0004014].
[11] H. Weigert, Unitarity at small Bjorken x, Nucl. Phys. A703 (2002) 823–860,
[hep-ph/0004044].
[12] E. Iancu, A. Leonidov, and L. D. McLerran, Nonlinear gluon evolution in the color glass
condensate. I, Nucl. Phys. A692 (2001) 583–645, [hep-ph/0011241].
[13] E. Ferreiro, E. Iancu, A. Leonidov, and L. McLerran, Nonlinear gluon evolution in the
color glass condensate. II, Nucl. Phys. A703 (2002) 489–538, [hep-ph/0109115].
[14] I. Balitsky, Operator expansion for high-energy scattering, Nucl. Phys. B463 (1996)
99–160, [hep-ph/9509348].
[15] I. Balitsky, Operator expansion for diffractive high-energy scattering, hep-ph/9706411.
http://arxiv.org/abs/hep-ph/0609105
http://xxx.lanl.gov/abs/hep-ph/0609105
http://arxiv.org/abs/hep-ph/0609090
http://xxx.lanl.gov/abs/hep-ph/0609090
http://arxiv.org/abs/hep-ph/0408216
http://xxx.lanl.gov/abs/hep-ph/0408216
http://arxiv.org/abs/hep-ph/9701284
http://xxx.lanl.gov/abs/hep-ph/9701284
http://arxiv.org/abs/hep-ph/9706377
http://xxx.lanl.gov/abs/hep-ph/9706377
http://arxiv.org/abs/hep-ph/9709432
http://xxx.lanl.gov/abs/hep-ph/9709432
http://arxiv.org/abs/hep-ph/9807462
http://xxx.lanl.gov/abs/hep-ph/9807462
http://arxiv.org/abs/hep-ph/0004014
http://xxx.lanl.gov/abs/hep-ph/0004014
http://arxiv.org/abs/hep-ph/0004044
http://xxx.lanl.gov/abs/hep-ph/0004044
http://arxiv.org/abs/hep-ph/0011241
http://xxx.lanl.gov/abs/hep-ph/0011241
http://arxiv.org/abs/hep-ph/0109115
http://xxx.lanl.gov/abs/hep-ph/0109115
http://arxiv.org/abs/hep-ph/9509348
http://xxx.lanl.gov/abs/hep-ph/9509348
http://arxiv.org/abs/hep-ph/9706411
http://xxx.lanl.gov/abs/hep-ph/9706411
[16] I. Balitsky, Factorization and high-energy effective action, Phys. Rev. D60 (1999)
014020, [hep-ph/9812311].
[17] Y. V. Kovchegov, Small-x F2 structure function of a nucleus including multiple pomeron
exchanges, Phys. Rev. D60 (1999) 034008, [hep-ph/9901281].
[18] Y. V. Kovchegov, Unitarization of the BFKL pomeron on a nucleus, Phys. Rev. D61
(2000) 074018, [hep-ph/9905214].
[19] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Singlet structure function at small x:
Unitarization of gluon ladders, Nucl. Phys. B188 (1981) 555–576.
[20] A. H. Mueller and J.-w. Qiu, Gluon recombination and shadowing at small values of x,
Nucl. Phys. B268 (1986) 427.
[21] L. D. McLerran and R. Venugopalan, Green’s functions in the color field of a large
nucleus, Phys. Rev. D50 (1994) 2225–2233, [hep-ph/9402335].
[22] L. D. McLerran and R. Venugopalan, Gluon distribution functions for very large nuclei
at small transverse momentum, Phys. Rev. D49 (1994) 3352–3355, [hep-ph/9311205].
[23] L. D. McLerran and R. Venugopalan, Computing quark and gluon distribution functions
for very large nuclei, Phys. Rev. D49 (1994) 2233–2241, [hep-ph/9309289].
[24] Y. V. Kovchegov, Non-Abelian Weizsaecker-Williams field and a two- dimensional
effective color charge density for a very large nucleus, Phys. Rev. D54 (1996) 5463–5469,
[hep-ph/9605446].
[25] Y. V. Kovchegov, Quantum structure of the non-Abelian Weizsaecker-Williams field for a
very large nucleus, Phys. Rev. D55 (1997) 5445–5455, [hep-ph/9701229].
[26] J. Jalilian-Marian, A. Kovner, L. D. McLerran, and H. Weigert, The intrinsic glue
distribution at very small x, Phys. Rev. D55 (1997) 5414–5428, [hep-ph/9606337].
[27] E. Iancu and R. Venugopalan, The color glass condensate and high energy scattering in
QCD, hep-ph/0303204.
[28] H. Weigert, Evolution at small xbj: The Color Glass Condensate, Prog. Part. Nucl. Phys.
55 (2005) 461–565, [hep-ph/0501087].
[29] J. Jalilian-Marian and Y. V. Kovchegov, Saturation physics and deuteron gold collisions
at RHIC, Prog. Part. Nucl. Phys. 56 (2006) 104–231, [hep-ph/0505052].
[30] E. Gardi, J. Kuokkanen, K. Rummukainen, and H. Weigert, Running coupling and power
corrections in nonlinear evolution at the high-energy limit, Nucl. Phys. A784 (2007)
282–340, [hep-ph/0609087].
[31] Y. V. Kovchegov and H. Weigert, Quark loop contribution to BFKL evolution: Running
coupling and leading-N(f) NLO intercept, Nucl. Phys. A789, 260 (2007),
[hep-ph/0612071].
http://arxiv.org/abs/hep-ph/9812311
http://xxx.lanl.gov/abs/hep-ph/9812311
http://arxiv.org/abs/hep-ph/9901281
http://xxx.lanl.gov/abs/hep-ph/9901281
http://arxiv.org/abs/hep-ph/9905214
http://xxx.lanl.gov/abs/hep-ph/9905214
http://arxiv.org/abs/hep-ph/9402335
http://xxx.lanl.gov/abs/hep-ph/9402335
http://arxiv.org/abs/hep-ph/9311205
http://xxx.lanl.gov/abs/hep-ph/9311205
http://arxiv.org/abs/hep-ph/9309289
http://xxx.lanl.gov/abs/hep-ph/9309289
http://arxiv.org/abs/hep-ph/9605446
http://xxx.lanl.gov/abs/hep-ph/9605446
http://arxiv.org/abs/hep-ph/9701229
http://xxx.lanl.gov/abs/hep-ph/9701229
http://arxiv.org/abs/hep-ph/9606337
http://xxx.lanl.gov/abs/hep-ph/9606337
http://arxiv.org/abs/hep-ph/0303204
http://xxx.lanl.gov/abs/hep-ph/0303204
http://arxiv.org/abs/hep-ph/0501087
http://xxx.lanl.gov/abs/hep-ph/0501087
http://arxiv.org/abs/hep-ph/0505052
http://xxx.lanl.gov/abs/hep-ph/0505052
http://arxiv.org/abs/hep-ph/0609087
http://xxx.lanl.gov/abs/hep-ph/0609087
http://arxiv.org/abs/hep-ph/0612071
http://xxx.lanl.gov/abs/hep-ph/0612071
[32] G. P. Lepage and S. J. Brodsky, Exclusive processes in perturbative quantum
chromodynamics, Phys. Rev. D22 (1980) 2157.
[33] S. J. Brodsky, H.-C. Pauli, and S. S. Pinsky, Quantum chromodynamics and other field
theories on the light cone, Phys. Rept. 301 (1998) 299–486, [hep-ph/9705477].
[34] A. H. Mueller and D. N. Triantafyllopoulos, The energy dependence of the saturation
momentum, Nucl. Phys. B640 (2002) 331–350, [hep-ph/0205167].
[35] D. N. Triantafyllopoulos, The energy dependence of the saturation momentum from RG
improved BFKL evolution, Nucl. Phys. B648 (2003) 293–316, [hep-ph/0209121].
[36] E. Iancu, K. Itakura, and L. McLerran, Geometric scaling above the saturation scale,
Nucl. Phys. A708 (2002) 327–352, [hep-ph/0203137].
[37] S. Munier and R. Peschanski, Universality and tree structure of high energy QCD, Phys.
Rev. D70 (2004) 077503, [hep-ph/0401215].
[38] G. Beuf and R. Peschanski, Universality of QCD traveling-waves with running coupling,
hep-ph/0702131.
[39] A. H. Mueller, Soft gluons in the infinite momentum wave function and the BFKL
pomeron, Nucl. Phys. B415 (1994) 373–385.
[40] A. H. Mueller and B. Patel, Single and double BFKL pomeron exchange and a dipole
picture of high-energy hard processes, Nucl. Phys. B425 (1994) 471–488,
[hep-ph/9403256].
[41] A. H. Mueller, Unitarity and the BFKL pomeron, Nucl. Phys. B437 (1995) 107–126,
[hep-ph/9408245].
[42] Z. Chen and A. H. Mueller, The dipole picture of high-energy scattering, the BFKL
equation and many gluon compound states, Nucl. Phys. B451 (1995) 579–604.
[43] S. Munier and R. Peschanski, Traveling wave fronts and the transition to saturation,
Phys. Rev. D69 (2004) 034008, [hep-ph/0310357].
[44] S. J. Brodsky, G. P. Lepage, and P. B. Mackenzie, On the elimination of scale
ambiguities in perturbative quantum chromodynamics, Phys. Rev. D28 (1983) 228.
[45] K. Golec-Biernat and M. Wüsthoff, Saturation effects in deep inelastic scattering at low
Q2 and its implications on diffraction, Phys. Rev. D59 (1998) 014017, [hep-ph/9807513].
[46] M. A. Braun, Pomeron fan diagrams with an infrared cutoff and running coupling, Phys.
Lett. B 576 (2003) 115, [hep-ph/0308320].
[47] S. Munier and R. Peschanski, Geometric scaling as traveling waves, Phys. Rev. Lett. 91
(2003) 232001, [hep-ph/0309177].
http://arxiv.org/abs/hep-ph/9705477
http://xxx.lanl.gov/abs/hep-ph/9705477
http://arxiv.org/abs/hep-ph/0205167
http://xxx.lanl.gov/abs/hep-ph/0205167
http://arxiv.org/abs/hep-ph/0209121
http://xxx.lanl.gov/abs/hep-ph/0209121
http://arxiv.org/abs/hep-ph/0203137
http://xxx.lanl.gov/abs/hep-ph/0203137
http://arxiv.org/abs/hep-ph/0401215
http://xxx.lanl.gov/abs/hep-ph/0401215
http://arxiv.org/abs/hep-ph/0702131
http://xxx.lanl.gov/abs/hep-ph/0702131
http://arxiv.org/abs/hep-ph/9403256
http://xxx.lanl.gov/abs/hep-ph/9403256
http://arxiv.org/abs/hep-ph/9408245
http://xxx.lanl.gov/abs/hep-ph/9408245
http://arxiv.org/abs/hep-ph/0310357
http://xxx.lanl.gov/abs/hep-ph/0310357
http://arxiv.org/abs/hep-ph/9807513
http://xxx.lanl.gov/abs/hep-ph/9807513
http://arxiv.org/abs/hep-ph/0308320
http://xxx.lanl.gov/abs/hep-ph/0308320
http://arxiv.org/abs/hep-ph/0309177
http://xxx.lanl.gov/abs/hep-ph/0309177
[48] J. L. Albacete, N. Armesto, A. Kovner, C. A. Salgado, and U. A. Wiedemann, Energy
dependence of the Cronin effect from non-linear QCD evolution, Phys. Rev. Lett. 92
(2004) 082001, [hep-ph/0307179].
[49] M. Lublinsky, Scaling phenomena from non-linear evolution in high energy DIS, Eur.
Phys. J. C21 (2001) 513–519, [hep-ph/0106112].
[50] N. Armesto and M. A. Braun, Parton densities and dipole cross-sections at small x in
large nuclei, Eur. Phys. J. C20 (2001) 517–522, [hep-ph/0104038].
[51] A. M. Stasto, K. Golec-Biernat, and J. Kwiecinski, Geometric scaling for the total γ∗p
cross-section in the low x region, Phys. Rev. Lett. 86 (2001) 596–599, [hep-ph/0007192].
[52] N. Armesto, C. A. Salgado, and U. A. Wiedemann, Relating high-energy lepton hadron,
proton nucleus and nucleus nucleus collisions through geometric scaling,
hep-ph/0407018.
[53] J. L. Albacete, N. Armesto, J. G. Milhano, C. A. Salgado, and U. A. Wiedemann,
Nuclear size and rapidity dependence of the saturation scale from QCD evolution and
experimental data, Eur. Phys. J. C43 (2005) 353–360, [hep-ph/0502167].
[54] E. Iancu, K. Itakura, and S. Munier, Saturation and BFKL dynamics in the HERA data
at small x, Phys. Lett. B590 (2004) 199–208, [hep-ph/0310338].
[55] A. Dumitru, A. Hayashigaki, and J. Jalilian-Marian, Geometric scaling violations in the
central rapidity region of d + Au collisions at RHIC, Nucl. Phys. A770 (2006) 57–70,
[hep-ph/0512129].
[56] V. P. Goncalves, M. S. Kugeratski, M. V. T. Machado, and F. S. Navarra, Saturation
physics at HERA and RHIC: An unified description, Phys. Lett. B643 (2006) 273–278,
[hep-ph/0608063].
[57] D. Kharzeev, Y. V. Kovchegov, and K. Tuchin, Nuclear modification factor in d + Au
collisions: Onset of suppression in the color glass condensate, Phys. Lett. B599 (2004)
23–31, [hep-ph/0405045].
http://arxiv.org/abs/hep-ph/0307179
http://xxx.lanl.gov/abs/hep-ph/0307179
http://arxiv.org/abs/hep-ph/0106112
http://xxx.lanl.gov/abs/hep-ph/0106112
http://arxiv.org/abs/hep-ph/0104038
http://xxx.lanl.gov/abs/hep-ph/0104038
http://arxiv.org/abs/hep-ph/0007192
http://xxx.lanl.gov/abs/hep-ph/0007192
http://arxiv.org/abs/hep-ph/0407018
http://xxx.lanl.gov/abs/hep-ph/0407018
http://arxiv.org/abs/hep-ph/0502167
http://xxx.lanl.gov/abs/hep-ph/0502167
http://arxiv.org/abs/hep-ph/0310338
http://xxx.lanl.gov/abs/hep-ph/0310338
http://arxiv.org/abs/hep-ph/0512129
http://xxx.lanl.gov/abs/hep-ph/0512129
http://arxiv.org/abs/hep-ph/0608063
http://xxx.lanl.gov/abs/hep-ph/0608063
http://arxiv.org/abs/hep-ph/0405045
http://xxx.lanl.gov/abs/hep-ph/0405045
	Introduction
	Scheme dependence
	Inclusion of running coupling corrections: general concepts
	Derivation of the subtraction term
	Brief summary of analytical results
	Numerical setup and initial conditions
	Results
	Running coupling
	Geometric scaling
	Subtraction Term
	Complete running coupling BK equation
	Conclusions

  We study the solution of the nonlinear BK evolution equation with the
recently calculated running coupling corrections [hep-ph/0609105,
hep-ph/0609090]. Performing a numerical solution we confirm the earlier result
of [hep-ph/0408216] that the high energy evolution with the running coupling
leads to a universal scaling behavior for the dipole scattering amplitude. The
running coupling corrections calculated recently significantly change the shape
of the scaling function as compared to the fixed coupling case leading to a
considerable increase in the anomalous dimension and to a slow-down of the
evolution with rapidity. The difference between the two recent calculations is
due to an extra contribution to the evolution kernel, referred to as the
subtraction term, which arises when running coupling corrections are included.
These subtraction terms were neglected in both recent calculations. We evaluate
numerically the subtraction terms for both calculations, and demonstrate that
when the subtraction terms are added back to the evolution kernels obtained in
the two works the resulting dipole amplitudes agree with each other! We then
use the complete running coupling kernel including the subtraction term to find
the numerical solution of the resulting full non-linear evolution equation with
the running coupling corrections. Again the scaling regime is recovered at very
large rapidity.

Introduction
Recently our understanding of the linear Balitsky–Fadin–Kuraev–Lipatov (BFKL) [4, 5] and
non-linear Jalilian-Marian–Iancu–McLerran–Weigert–Leonidov–Kovner (JIMWLK) [6–13] and
Balitsky-Kovchegov (BK) [14–18] small-x evolution equations in the Color Glass Condensate
[6–29] has been improved due to the completion of the calculations determining the scale of the
running coupling in the evolution kernel in [1,2,30,31]. The calculations in [1,2] proceeded by
including αsNf corrections into the evolution kernel and by then completing Nf to the complete
one-loop QCD beta-function by replacing Nf → −6πβ2. Calculation of the αsNf corrections
is particularly easy in the s-channel light-cone perturbation theory formalism [32, 33] used to
derive the BK and JIMWLK equations: there αs Nf corrections are solely due to chains of
quark bubbles placed onto the s-channel gluon lines.
The analytical results of [1,2] are not very concise and could not have been guessed without
an explicit calculation. After finding αsNf corrections, the obtained contributions had to be
divided into the running coupling part, which has a form of a running coupling correction
to the leading-order (LO) JIMWLK or BK kernel, and into the ”subtraction piece”, which
would bring in new structures into the kernel. Such separation had to be done both in [1] and
in [2]. Unfortunately, there appears to be no unique way to perform this separation: it is not
surprising, therefore, that it was done differently in both papers [1, 2]. This resulted in two
different running coupling terms, shown below in Eqs. (35) and (36) along with Eqs. (7) and
(8). Such a discrepancy has led to a misconception in the community that the calculations
of [1] and of [2] disagree at some fundamental level.
Indeed to compare the results of [1] and [2] one has to undo the separation into the running
coupling and subtraction terms: combining both terms one should compare full kernels of the
evolution equation obtained in [1, 2]. There is another more physical reason to perform such
comparison: in principle, there is no small parameter making the subtraction term smaller than
the running coupling term and thus justifying neglecting the former compared to the latter.
Even the labeling of one term as “running coupling” piece is somewhat misleading, since it may
give an impression that the neglected subtraction term has no running coupling corrections in
it. As was shown in [31] both terms actually contribute to the running coupling corrections to
the BFKL equation (if one uses the separation of [2] to define the terms).
In this paper we perform numerical analysis of the BK evolution equation with the αsNf
corrections resummed to all orders and with Nf completed to the QCD beta-function, Nf →
−6πβ2, with β2 given in Eq. (20). We first solve the BK equation keeping the running coupling
term only, with the kernels given by Eqs. (7) and (8). Indeed the solutions we find this way
are different from each other. We then evaluate the subtraction terms for both cases and show
that inclusion of subtraction terms puts the results of [1] and [2] in perfect agreement with each
other! We complete our analysis by solving the BK equation with the full kernel including both
the running coupling and subtraction terms.
This work is structured as follows. Section 2 begins with Sect. 2.1 in which we review
the αsNf corrections to the dipole scattering amplitude evolution equation recently derived
in [1, 2] and the subtraction method employed in both works to separate the running coupling
contributions from the subtraction terms. We discuss the scheme dependence of the running
coupling terms introduced by this separation. We proceed in Sect. 2.2 by deriving the explicit
expressions for the subtracted terms. The calculation is based on the results of [2]. Our
analytical results are summarized in Sect. 2.3, where we give the explicit final expression
for the kernel of the subtraction term in Eq. (39), which, combined with Eq. (38) gives us the
subtraction terms (40) and (41) for the subtractions performed in [1] and in [2] correspondingly.
In Sect. 3 we explain the numerical method we use to solve the evolution equations. We
also list the initial conditions used, along with the definition of the saturation scale employed.
Throughout the paper we will avoid the important question of the Landau pole and the contri-
bution of renormalons to small-x evolution. As we explain in Sect. 3, we will simply “freeze”
the running coupling at a constant value in the infrared. For a detailed study of the renormalon
effects in the non-linear evolution we refer the readers to [30].
Our numerical results are presented in Sect. 4. By solving the evolution equations with the
running coupling term only in Sect. 4.1 we show that the resulting dipole amplitude differs
significantly from the fixed coupling case. We also observe that the amplitude obtained by
solving the equation obtained in [2] is very close to the result of solving the BK evolution
with a postulated parent-dipole running of the coupling constant. Both these amplitudes are
quite different from the solution of the equation derived in [1], as one can see from Fig. 4.
In spite of that, all three evolution equations studied (the ones derived in [1], [2] and the
parent-dipole running coupling model) give approximately identical scaling function for the
dipole amplitude at high rapidity, as demonstrated in Fig. 5 in Sect. 4.2. It is worth noting
that, as can be seen from Fig. 6, the anomalous dimension we extracted from our solution
is γ ≈ 0.85, which is different from the fixed coupling anomalous dimension of γ ≈ 0.64.
The former anomalous dimension also appears to disagree with the predictions of analytical
approximations to the behavior of the dipole amplitude with running coupling from [34–38]. In
Sect. 4.3 we numerically evaluate the subtraction terms for both [1] and [2] and show that their
contributions are important, as shown in Fig. 7. However, subtraction terms decrease with
increasing rapidity, such that at high enough rapidities their relative contribution becomes
small (see Fig. 8). In Fig. 9 we show that inclusion of subtraction terms makes the results of [1]
and [2] agree with each other. Finally, the numerical solution of the full (all orders in αs β2)
evolution equation including both the running coupling and subtraction terms is performed in
Sect. 4.4. The results are shown in Fig. 10. All the main features of the evolution with the
running coupling are preserved in the full solution: the growth of the dipole amplitude and
of the saturation scale with rapidity is slowed down (for the latter see Fig. 11). The scaling
function of Fig. 5 is unaltered by the subtraction term, as shown in Fig. 12.
We summarize and discuss our main conclusions in Sect. 5.
2 Scheme dependence
2.1 Inclusion of running coupling corrections: general concepts
The BK evolution equation for the dipole scattering matrix reads
∂S(x0, x1; Y )
d2z K(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] , (1)
where
K(x0, x1, z) =
r21 r
is the kernel of the evolution. Here transverse two-dimensional vectors x0 and x1 denote the
transverse coordinates of the quark and the anti-quark in the parent dipole, while z is the
position of the gluon produced in one step of evolution [39–42]. We have introduced the
notation r = x0 − x1, r1 = x0 − z, r2 = z − x1 for the sizes of the parent and of the new
(daughter) dipoles created by one step of the evolution. The notation r ≡ |r| for all the 2-
dimensional vectors will be also employed throughout the rest of the paper. Eq. (1) admits
a clear physical interpretation: the original parent dipole, when boosted to higher rapidities,
may emit a new gluon which, in the large-Nc limit, is equivalent to a quark-antiquark pair.
Thus, the original dipole splits into two new dipoles sharing a common transverse coordinate:
the transverse position of the emitted gluon, z. The nonlinear term in the right hand side of
Eq. (1) accounts for either one of the two new dipoles interacting with the target, along with
the possibilities of only one dipole interacting or no interaction at all, while the subtracted
linear term reflects virtual corrections. The kernel of the evolution is just the probability of
one gluon emission calculated at leading logarithmic accuracy in αs ln(1/xB), where xB is the
fraction of momentum carried by the emitted gluon [39–42].
Under the eikonal approximation the dipole scattering matrix off a hadronic target at a
fixed rapidity is given by the average over the hadron field configurations of Wilson lines V
calculated along fixed transverse coordinates (those of the quark and of the antiquark). More
specifically
S(x0, x1; Y ) =
V (x0)V
†(x1)
〉 . (3)
Hence, the integrand of Eq. (1) can be regarded as a three point function in the sense that
the gluon fields of the target are evaluated at three different transverse positions, those of the
original quark and antiquark plus the one of the emitted gluon.
However, the inclusion of higher order corrections to the evolution equation via all order
resummation of αsNf contributions as recently derived in [1, 2] brings in new physical chan-
nels that modify the three point structure of the leading-log equation. The dipole structure
generated under evolution by diagrams like the one depicted in Fig. 1A (for more detailed dis-
cussion of the diagrammatic content of the high order corrections see [2]) is identical to the one
previously discussed for the leading order equation, the only novelty being that the propagator
of the emitted gluon is now dressed with quark loops, modifying the emission probability but
leaving untouched the interaction terms. On the contrary, diagrams like the one in Fig. 1B
in which a quark-antiquark pair (rather than a gluon) is added to the evolved wave function
modify the interaction structure of the evolution equation. The evolution of the parent dipole
scattering matrix driven by these kind of terms is proportional to the scattering matrix of the
two newly created dipoles (the one formed by the original quark and the new antiquark and
vice versa), ∼ S(x0, z1)S(z2, x1). This term depends on four different transverse coordinates,
i.e., it is a four point function and, therefore, its contribution to the evolution equation cannot
be accounted for by a mere modification of the emission kernel of the leading order equation.
To discuss in more detail the modifications introduced by the high order corrections, we
find it useful to rewrite the evolution equation in the following, rather general way
∂S(x0, x1; Y )
= F [S(x0, x1; Y )] (4)
where F is a functional of the dipole scattering matrix which for the original derivation of the
equation is given by the right hand side of Eq. (1). In general it can be decomposed into two
Figure 1: Schematic representation of the diagrams contributing to quark-NLO evolution.
pieces
F [S] = R [S]− S [S] . (5)
The first term, R, which we will call the ’running coupling’ contribution, gathers all the
higher order in αsNf corrections to the evolution that can be recast in a functional form that
looks identical to the leading order one but with a modified kernel, K̃, which includes all the
terms setting the scale for the running coupling:
R [S(x0, x1; Y )] =
d2z K̃(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] . (6)
The second term, S, henceforth referred to as the ’subtraction’ contribution, encodes those
contributions that depart from the three point structure of the leading-log equation. The
explicit derivation and expressions for this term are presented in the next section. The relative
minus sign between the two terms in Eq. (5) has been introduced for latter convenience.
Importantly, the decomposition of F into running coupling and subtraction contributions,
although constrained by unitarity arguments, is not unique. Two different separation schemes
have been proposed in [1, 2]. They are both based on a similar strategy, sketched in Fig. 2,
that can be summarized as follows. The newly created quark-antiquark pair added to the
wave function in the diagrams Fig. 1B is shrunk to a point, called the subtraction point, by
integrating out one of the coordinates in the dipole-qq̄ wave function, rendering the previously
discussed four point nature of these contributions into a three point one. This integrated
three point contribution is added to the running coupling contribution, whereas the original
four point term minus its integrated version are assigned to the subtraction contribution. The
divergence between the two approaches stems from the choice of the subtraction point. In
the subtraction scheme proposed by Balitsky in [1] the subtraction point is chosen to be the
transverse coordinate of either the quark, z2, or the antiquark, z1. The kernel for the running
coupling functional, Eq. (6), obtained in this way is
K̃Bal(r, r1, r2) =
Nc αs(r
r21 r
. (7)
On the other hand, in the subtraction procedure followed in [2] (which we will refer to as
KW) the zero size quark-antiquark pair is fixed at the transverse coordinate of the gluon,
x 0 x 0
Figure 2: Schematic representation of the subtraction procedure.
z = αz1 + (1 − α)z2, where α is the fraction of the gluon’s longitudinal momentum carried by
the quark, yielding the following expression for the kernel of the running coupling contribution:
K̃KW(r, r1, r2) =
− 2 αs(r
1)αs(r
αs(R2)
r1 · r2
r21 r
+ αs(r
, (8)
where
R2(r, r1, r2) = r1 r2
r1·r2
. (9)
As we shall discuss later, the scheme dependence originated by the choice of the subtraction
point is substantial and has an important effect in the solutions of the evolution equation when
only the running contribution is taken into account.
In our numerical study we will also consider the following ad hoc prescription for the kernel
of the running coupling functional in which the scale for the running of the coupling is set to
be the size of the parent dipole
K̃pd(r, r1, r2) =
Nc αs(r
r21 r
. (10)
This prescription is useful as a benchmark used to compare with previous numerical [3] and
analytical works [34, 35, 43] where this ansatz was used.
2.2 Derivation of the subtraction term
We begin by considering the NLO contribution to the kernel of the JIMWLK and BK evolution
equations with the s-channel gluon splitting into a quark-antiquark pair, which then interacts
with the target, as shown on the left hand side of Fig. 3. The contribution of this diagram has
been calculated in [2]. The resulting JIMWLK kernel is [2]
KNLO1 (x0, x1; z1, z2) = 4Nf
(2π)2
(2π)2
(2π)2
(2π)2
e−iq·(z−x0)+iq
′·(z−x
)−i(k−k′)·z
q2q′2
(1− 2α)2q · k k′ · q′ + q · q′ k · k′ − q · k′ k · q′
k2 + q2α(1− α)
k′2 + q′2α(1− α)
2α (1− α) (1− 2α)
k2 + q2α(1− α)
k′2 + q′2α(1− α)
k · q
k′ · q′
4α2 (1− α)2
k2 + q2α(1− α)
k′2 + q′2α(1− α)
 . (11)
The momentum labels in the above equation are explained on the left hand side of Fig. 3.
If k1 and k2 are the transverse momenta of the quark and of the antiquark in the produced
pair as shown in Fig. 3, then the transverse momentum of the gluon is q = k1 + k2. The
other transverse momentum we use is k = k1(1 − α) − k2α, where α is the fraction the of
gluon’s “plus” momentum carried by the quark, α ≡ k1+/(k1+ + k2+). The prime over the
transverse momentum denotes the momentum of the same particle in the complex conjugate
amplitude. For instance q′ is the momentum of the s-channel gluon in the complex conjugate
amplitude. Finally, z1 and z2 denote the transverse coordinates of the quark and the antiquark.
In Eq. (11) we use z12 = z1−z2 (the transverse separation between the quark and the antiquark)
and z = α z1 + (1− α) z2 (the transverse coordinate of the gluon).
Figure 3: A lowest order leading-Nf NLO correction which gives rise to the subtraction term
is shown on the left. The same diagram with the gluon lines “dressed” by chains of fermion
bubbles, as shown on the right, gives the full (resumming all powers of αµNf ) contribution to
the subtraction term. Calculation of the subtraction term is pictured in Fig. 2.
To obtain the BK kernel from Eq. (11) one should sum over all possible emissions of the
gluon off the quark and antiquark lines in the incoming dipole both in the amplitude and in
the complex conjugate amplitude, which is accomplished by
KNLO1 (x0, x1; z1, z2) = CF
m,n=0
(−1)m+n KNLO1 (xm, xn; z1, z2). (12)
Below we will label the JIMWLK kernel by calligraphic letter K and the corresponding BK
kernel by K.
The contribution of the kernel from Eq. (12) to the right hand side of the NLO version of
Eq. (1) is given by the following term
d2z1 d
1 (x0, x1; z1, z2)S(x0, z1, Y )S(z2, x1, Y ) (13)
with αµ the bare coupling.
As shown in Fig. 2, at the NLO level, the subtraction term introduced in Eq. (5), is then
defined by
SNLO[S] = α2µ
d2z1 d
2z2 K
1 (x0, x1; z1, z2)
× [S(x0, w, Y )S(w, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )], (14)
where w is the point of subtraction in the transverse coordinate space. In [1] it was chosen to
be equal to the transverse coordinate of either the quark or the antiquark,
w = z1 or w = z2, (15)
as both choices lead to the same subtraction term SBalNLO[S]:
SBalNLO[S] =
d2z1 d
2z2 K
1 (x0, x1; z1, z2)
× [S(x0, z1, Y )S(z1, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] . (16)
In [2] the subtraction point was chosen to be the transverse coordinate of the gluon z,
w = z = α z1 + (1− α) z2. (17)
This leads to the following subtraction term, which we denote SKWNLO[S]:
SKWNLO[S] =
d2z1 d
1 (x0, x1; z1, z2)
× [S(x0, z, Y )S(z, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] . (18)
Indeed the complete kernel in Eq. (5) is independent of the choice of w. However, since the
subtraction term of Eq. (14) was neglected both in [1] and in [2], different choices of w led to
different expressions for the remaining running coupling part R[S], i.e., to different answers as
far as investigations in [1] and in [2] were concerned. Different choice of w is the main source
of the discrepancy of final answers of [1] and [2], though it does not imply any disagreement in
the full expression (5).
Our goal in this Section is to evaluate KNLO1 (x0, x1; z1, z2) from Eq. (11) including the
running coupling corrections. The s-channel light-cone perturbation theory formalism makes
such inclusion simple [2]: all we have to do is include infinite chains of quark bubbles on the
gluon lines in the amplitude and in the complex conjugate amplitude, as depicted on the right
hand side of Fig. 3. Performing calculations similar to those done in [2] one arrives at
K ❣1 (x0, x1; z1, z2) = 4Nf
(2π)2
(2π)2
(2π)2
(2π)2
e−iq·(z−x0)+iq
′·(z−x
)−i(k−k′)·z
q2q′2
(1− 2α)2q · k k′ · q′ + q · q′ k · k′ − q · k′ k · q′
k2 + q2α(1− α)
k′2 + q′2α(1− α)
2α (1− α) (1− 2α)
k2 + q2α(1− α)
k′2 + q′2α(1− α)
k · q
k′ · q′
4α2 (1− α)2
k2 + q2α(1− α)
k′2 + q′2α(1− α)
1 + αµβ2 ln
q2 e−5/3
1 + αµβ2 ln
q′2 e−5/3
) (19)
where K ❣1 denotes the kernel with the running coupling corrections resummed to all orders.
Just like in [2, 31], here we will use MS renormalization scheme. Inclusion of fermion bubble
chains generated two denominators at the end of Eq. (19), which is its only difference from
Eq. (11). Here
11Nc − 2Nf
. (20)
Now we have to perform the transverse momentum integrals in Eq. (19). First we expand
the denominators at the end of Eq. (19) into a power series and rewrite Eq. (19) as
K ❣1 (x0, x1; z1, z2) = 4Nf
n,m=0
(−αµβ2)n+m
(2π)2
(2π)2
(2π)2
(2π)2
e−iq·(z−x0)+iq
′·(z−x
)−i(k−k′)·z
q2q′2
(1− 2α)2q · k k′ · q′ + q · q′ k · k′ − q · k′ k · q′
k2 + q2α(1− α)
k′2 + q′2α(1− α)
2α (1− α) (1− 2α)
k2 + q2α(1− α)
k′2 + q′2α(1− α)
k · q
k′ · q′
4α2 (1− α)2
k2 + q2α(1− α)
k′2 + q′2α(1− α)
λ=λ′=0
where we have defined µ2 = µ2
e5/3 to make the expressions more compact.
Indeed we can not always expand the denominators of Eq. (19) into a geometric series
employed in Eq. (21), but one has to remember that the summation of bubble chain diagrams
shown on the right side of Fig. 3 gives one the geometric series. Hence the geometric series
come first: later they are absorbed into the denominators shown in Eq. (19), which is an
approximation not valid for all q and q′. Therefore, by keeping the geometric series in Eq. (21)
we are not making any approximations. In general, in what follows we are not going to keep
track of the issues of convergence of perturbation series. The contribution of renormalons to
non-linear small-x evolution was thoroughly investigated in [30] and was found to be significant
at low Q2. We refer the interested reader to [30] for more details on this issue.
Using the following formulas
(2π)2
e−ik·z
k2 + q2
K0(q z) (22)
(2π)2
e−ik·z
k2 + q2
q K1(q z) (23)
we can now perform the k- and k′-integrals in Eq. (21). Integrating over the angles of q and q′
as well yields
K ❣1 (x0, x1; z1, z2) =
(2π)4
n,m=0
(−αµβ2)n+m
dq q dq′ q′
z212 |z − x0| |z − x1|
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
× J1(q |z − x0|)K1(z12 q
α ᾱ) J1(q
′ |z − x1|)K1(z12 q′
α ᾱ) + 2α ᾱ (α− ᾱ)
z12 · (z − x0)
z12 |z − x0|
J1(q |z − x0|)K1(z12 q
α ᾱ) J0(q
′ |z − x1|)K0(z12 q′
α ᾱ)
z12 · (z − x1)
z12 |z − x1|
J0(q |z − x0|)K0(z12 q
α ᾱ) J1(q
′ |z − x1|)K1(z12 q′
α ᾱ)
+ 4α2 ᾱ2 J0(q |z − x0|)K0(z12 q
α ᾱ) J0(q
′ |z − x1|)K0(z12 q′
α ᾱ)
λ=λ′=0
. (24)
We have defined
ᾱ = 1− α (25)
for brevity. Now the integrals over q and q′ can be carried out to give
K ❣1 (x0, x1; z1, z2) =
(2π)4
n,m=0
(−αµβ2)n+m
z212 µ
2 α ᾱ
)λ+λ′
Γ2(1 + λ)
Γ2(1 + λ′)
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
(1 + λ) (1 + λ′)
z812 (α ᾱ)
1 + λ, 2 + λ; 2;−|z − x0|
α ᾱ z212
1 + λ′, 2 + λ′; 2;−|z − x1|
α ᾱ z212
α− ᾱ
z12 · (z − x0)
(1 + λ)F
1 + λ, 2 + λ; 2;−|z − x0|
α ᾱ z212
1 + λ′, 1 + λ′; 1;−|z − x1|
α ᾱ z212
z12 · (z − x1)
1 + λ, 1 + λ; 1;−
|z − x0|2
α ᾱ z212
(1 + λ′)F
1 + λ′, 2 + λ′; 2;−
|z − x1|2
α ᾱ z212
1 + λ, 1 + λ; 1;−
|z − x0|2
α ᾱ z212
1 + λ′, 1 + λ′; 1;−
|z − x1|2
α ᾱ z212
λ=λ′=0
. (26)
Unfortunately further simplification of the expression in Eq. (26) is impossible without approx-
imations. The series resulting from summation over n and m are likely to be divergent due to
renormalons. As we mentioned before, here we neglect the renormalon problem referring the
reader to [30]. Similar to how it was done in [2] we are not going to attempt to resum the series
exactly: instead we will calculate the next-to-leading order terms and assume that with a good
accuracy they give us the scale(s) of the running coupling constant. This procedure is similar
to the well-known prescription due to Brodsky, Lepage and Mackenzie [44].
Using the Taylor-expansions of hypergeometric functions
F (1 + λ, 2 + λ; 2; z) =
− λ 1
1 + ln(1− z) + 1
ln(1− z)
+ o(λ2). (27)
F (1 + λ, 1 + λ; 1; z) =
ln (1− z) + o(λ2) (28)
after some algebra we obtain
K ❣1 (x0, x1; z1, z2) =
[α (z1 − x0)2 + ᾱ (z2 − x0)2] [α (z1 − x1)2 + ᾱ (z2 − x0)2] z412
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
1− αµ β2 ln
R2T (x0)µ
+ o(α2µ)
1− αµ β2 ln
R2T (x1)µ
+ o(α2µ)
+ 2α ᾱ (α− ᾱ) z212
z12 · (z − x0)
1− αµ β2 ln
R2T (x0)µ
+ o(α2µ)
1− αµ β2 ln
R2L(x1)µ
+ o(α2µ)
+ z12 · (z − x1)
1− αµ β2 ln
R2L(x0)µ
+ o(α2µ)
1− αµ β2 ln
R2T (x1)µ
+ o(α2µ)
+4α2 ᾱ2 z412
1− αµ β2 ln
R2L(x0)µ
+ o(α2µ)
1− αµ β2 ln
R2L(x1)µ
+ o(α2µ)
In arriving at Eq. (29) we employed functions RT (x) and RL(x), which have dimensions of
transverse coordinates and are defined by
R2T (x)µ
4 e−2γ−5/3
[α (z1 − x)2 + ᾱ (z2 − x)2]µ2MS
α ᾱ z212
(z − x)2
α (z1 − x)2 + ᾱ (z2 − x)2
α ᾱ z212
R2L(x)µ
4 e−2γ−5/3
[α (z1 − x)2 + ᾱ (z2 − x)2]µ2MS
α (z1 − x)2 + ᾱ (z2 − x)2
α ᾱ z212
The subscripts T and L stand for transverse and longitudinal gluon polarizations which give
rise to the two different functions under the logarithm.
Recombining the series in Eq. (29) into physical running couplings finally yields
α2µK ❣1 (x0, x1; z1, z2) =
[α (z1 − x0)2 + ᾱ (z2 − x0)2] [α (z1 − x1)2 + ᾱ (z2 − x0)2] z412
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
R2T (x0)
R2T (x1)
+ 2α ᾱ (α− ᾱ) z212
z12 · (z − x0) αs
R2T (x0)
R2L(x1)
+z12 · (z − x1)αs
R2L(x0)
R2T (x1)
+ 4α2 ᾱ2 z412 αs
R2L(x0)
R2L(x1)
with the physical running coupling in the MS scheme given by
αs(1/R
1 + αµβ2 ln
R2 µ2
) . (33)
Eq. (32) is the contribution to the JIMWLK evolution kernel of the resummed diagram on
the right hand side of Fig. 3.
2.3 Brief summary of analytical results
Let us briefly summarize our analytical results. The non-linear small-x evolution equation with
the running coupling corrections included reads
∂S(x0, x1; Y )
= R [S]− S [S] . (34)
The first term on the right hand side of Eq. (34) is referred to as the running coupling
contribution. It was calculated independently in [1] and in [2]: the results of those calculations
are given above in Eqs. (7) and (8) correspondingly, which have to be combined with Eq. (6)
to obtain
RBal [S] =
d2z K̃Bal(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] (35)
RKW [S] =
d2z K̃KW(x0, x1, z) [S(x0, z; Y )S(z, x1; Y )− S(x0, x1; Y )] . (36)
One notices immediately that RBal [S] calculated in [1] is different from RKW [S] calculated
in [2] due to the difference in the kernels K̃Bal and K̃KW in Eqs. (7) and (8). However that
does not imply disagreement between the calculations of [1] and [2]: after all, it is the full
kernel on the right of Eq. (34), R [S]− S [S], that needs to be compared. To do that one has
to calculate the second term on the right hand side of Eq. (34).
The second term on the right hand side of Eq. (34) is referred to as the subtraction contri-
bution. It is given by
S[S] = α2µ
d2z1 d
2z2K ❣1 (x0, x1; z1, z2) [S(x0, w, Y )S(w, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )]
with the resummed BK kernel
K ❣1 (x0, x1; z1, z2) = CF
m,n=0
(−1)m+n K ❣1 (xm, xn; z1, z2). (38)
The resummed JIMWLK kernel K ❣1 (xm, xn; z1, z2) is given by Eq. (32), along with Eqs. (30)
and (31) defining the scales of the running couplings. In the numerical solution below we will
replace Nf → −6πβ2 in its prefactor, obtaining
α2µK ❣1 (x0, x1; z1, z2) = −
[α (z1 − x0)2 + ᾱ (z2 − x0)2] [α (z1 − x1)2 + ᾱ (z2 − x0)2] z412
− 4α ᾱ z12 · (z − x0) z12 · (z − x1) + z212 (z − x0) · (z − x1)
R2T (x0)
R2T (x1)
+ 2α ᾱ (α− ᾱ) z212
z12 · (z − x0) αs
R2T (x0)
R2L(x1)
+z12 · (z − x1)αs
R2L(x0)
R2T (x1)
+ 4α2 ᾱ2 z412 αs
R2L(x0)
R2L(x1)
This substitution is the same as for all other factors ofNf . The same substitution was performed
in [2] to calculate the running coupling term. In fact, as was shown in [31], the linear part of the
subtraction term (calculated using the prescription of [2]) contributes to the running coupling
corrections to the BFKL equation. Therefore, in that case, the factor of Nf in front of Eq. (32)
is definitely a part of the beta-function. Hence the replacement Nf → −6πβ2 is justified even
in the subtraction term. Once again, in the numerical solution below we will use Eq. (39) along
with Eq. (38) in Eq. (37) to calculate the subtraction term S[S].
Substituting w = z1 (or, equivalently, w = z2) in Eq. (37) would yield the subtraction term
SBal[S] =α2µ
d2z1 d
2z2K ❣1 (x0, x1; z1, z2)
× [S(x0, z1, Y )S(z1, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] (40)
which has to be subtracted from RBal [S] calculated in [1] and given by Eq. (35) to obtain the
complete evolution equation resumming all orders of αsNf in the kernel.
Substituting w = z = α z1 + (1− α) z2 in Eq. (37) yields
SKW[S] =α2µ
d2z1 d
2z2K ❣1 (x0, x1; z1, z2)
× [S(x0, z, Y )S(z, x1, Y )− S(x0, z1, Y )S(z2, x1, Y )] (41)
which has to be subtracted from RKW [S] calculated in [2] and given in Eq. (36) again to obtain
the complete evolution equation resumming all orders of αsNf in the kernel. We checked
explicitly by performing analytic calculations that the two evolution equations obtained this
way agree at the NLO and NNLO. Below we will check the agreement of the two calculations
to all orders by performing a numerical analysis of the solutions of these equations.
The above discussion demonstrates that the separation of the evolution kernel into the
running coupling and subtraction pieces, as done in Eq. (34), is somewhat artificial, and has
no small parameter justifying one or another separation prescription. Therefore, the small-x
evolution equation including all running coupling (or, more precisely, αsNf ) corrections should
combine both terms in Eq. (34). Below we will solve such evolution equation numerically to
obtain the full small-x evolution with the running coupling.
3 Numerical setup and initial conditions
In our numerical study we consider the translational invariant approximation in which the
scattering matrix is independent of the impact parameter of the collision, i.e., S = S(r, Y ).
To solve the integro-differential equations, corresponding to the BK equation with running
coupling we employ a second-order Runge-Kutta method with a step size in rapidity ∆Y = 0.1.
We discretize the variable |r| into 800 points equally separated in logarithmic space between
rmin = 10
−8 and rmax = 50. Throughout this paper, the units of r will be GeV
−1, and those
of Qs will be GeV. All the integrals have been performed using improved adaptative Gaussian
quadrature methods. The accuracy of this numerical method has been checked in [3] to be
better than a 4% in all the r range.
We consider three different initial conditions for the dipole scattering amplitude, N(r, Y ) =
1− S(r, Y ). The first one is taken from the McLerran-Venugopalan (MV) model [22, 23]:
NMV (r, Y = 0) = 1− exp
. (42)
where a constant term has been added to the argument of the logarithm in the exponent in
order to regularize it for large values of r. The other two initial conditions are given by
NAN (r) = 1− exp
−(r Q
, (43)
with γ = 0.6 and γ = 0.8. These two last initial conditions will be referred hereinafter as AN06
and AN08 respectively. The interest in this ansatz, reminiscent of the Golec-Biernat–Wusthoff
model [45], is that the small-r behavior NAN ∝ r2γ corresponds to an anomalous dimension
1− γ of the unintegrated gluon distribution at large transverse momentum. (AN labels initial
conditions with anomalous dimension.) Our choices γ = 0.6 and γ = 0.8 can be motivated
a posteriori by the observation that the anomalous dimension of the evolved BK solution for
running coupling lies in between those two values and the one for the MV initial condition,
γ ≈ 1 (see Section 4.2). Thus, the choice of distinct initial conditions allows us to better
track the onset of the expected asymptotic universal behavior that is eventually reached at
high energies and to study the influence of the pre-asymptotic, non-universal corrections to the
solutions of the evolution equations. To completely determine our initial conditions, we set
Q′s = 1 GeV at Y=0 in Eqs. (42) and (43) and put Λ = 0.2 GeV. Although Q
s is normally
identified with the saturation scale, our definition of the saturation scale through the rest of
the paper will be purely pragmatical and given by the condition
N(r = 1/Qs(Y ), Y ) = κ, (44)
with κ = 0.5. We have checked that this choice of κ, albeit arbitrary, does not affect any of the
major conclusions to be drawn in the rest of the paper.
Finally, in order to avoid the Landau pole and to regularize the running coupling at large
transverse sizes we stick to the following procedure: for small transverse distances r < rfr, with
rfr defined by αs(1/r
fr) = 0.5, the running coupling is given by the one loop expression
αs(1/r
β2 ln
r2 Λ2
) (45)
with Nf = 3 and Λ = 0.2 GeV, whereas for larger sizes, r > rfr, we “freeze” the coupling at a
fixed value αs = 0.5. A detailed study of the role of Landau pole in non-linear small-x evolution
is given in [30].
4 Results
In this Section, we discuss our numerical results and how they compare to previous numerical
work and analytical estimates.
4.1 Running coupling
Fig. 4 shows the solutions of the evolution equation when only the running coupling contribution
is taken into account, i.e., neglecting the subtraction term in Eq. (34), for different initial
conditions and for the three schemes considered in this work: Balitsky’s, given by Eqs. (7) and
(35), KW, given by Eqs. (8) and (36), and the the ad hoc parent dipole implementation of the
running coupling, shown in Eq. (10).
MV init. cond.     Y=0,5,15,30KW
Balitsky
parent dipole
init. cond.
AN08 init. cond.     Y=0,5,15,30
AN06 init. cond.     Y=0,5,15,30
r )-1(GeV
Figure 4: Solutions of the BK equation at rapidities Y=0, 5, 15 and 30 (curves are labeled from
right to left) for the three running coupling schemes considered in this work: KW (solid line),
Balitsky (dashed line) and parent dipole (dashed-dotted lines). The initial conditions are MV
(top), AN08 (middle) and AN06 (bottom).
As previously observed in [3, 46], the most relevant effect of including running coupling
corrections in the evolution equation is a considerable reduction in the speed of the evolution
with respect to the fixed coupling case. This is a common feature of the different running
coupling schemes studied here and of other phenomenological ones considered in the literature
(a detailed comparison between the solutions for fixed coupling evolution and for parent dipole
running coupling can be found e.g. in [3]). This is not a surprising result, since a generic effect
of the running of the coupling is to suppress the emission of small transverse size dipoles, which
is the leading mechanism driving the evolution.
However, despite this common feature of the running coupling solutions, significant dif-
ferences are found between the solutions obtained under different schemes as we infer from
Fig. 4. In particular, the evolution is much faster with the KW prescription than with that of
Balitsky. Equivalently, the KW prescription yields a stronger growth of the saturation scale
with rapidity/energy than Balitsky’s. Moreover, the solutions obtained when the parent dipole
prescription is used lay much closer to those obtained within the KW scheme than to the ones
obtained when Balitsky’s scheme is applied, contrary to what was suggested in [1]. As argued
before, the differences observed in the solutions obtained using the two subtraction schemes
are entirely due to neglecting the subtraction contribution and reflect the arbitrariness of the
separation procedure.
4.2 Geometric scaling
It has been found in previous analytical [34, 36, 47] and numerical studies on the solutions
of the BK equation at leading order [3, 48–50] and for different heuristic implementations of
next-to-leading order corrections [3,46], including the parent dipole prescription for the running
coupling also considered in this work, that the solutions of the evolution equation at high enough
rapidities are no longer a function of two separate variables r and Y , but rather they depend
on a single scaling variable, τ = r Qs(Y ). This feature of the evolution, commonly referred to
as geometric scaling, is an exact property of the solutions for fixed coupling evolution due to
the conformal invariance of the leading-log kernel, and has become one of the key connections
between the saturation based formalisms and the phenomenology of heavy ion collisions and
deep inelastic scattering experiments [51–57].
It can be seen from Fig. 5 that the solutions of the BK equation with the running coupling
terms discussed in the previous section also exhibit the property of scaling, in agreement with
the analytical study carried out in [38], shown by the fact that the rescaled high rapidity
solutions lay on a single curve which is independent of both the running coupling scheme and
of the initial condition. The scaling behavior of the solution is observed in the whole τ range
studied in this work, including the saturation region, τ > 1. The tiny deviations from a pure
scaling behavior observed in Fig. 5 may be attributed to the fact that the full asymptotic
behavior is reached at even larger rapidities (Y & 80, [3]) than those achieved by the numerical
solution performed in this work.
Remarkably, the scaling function for both KW and Balitsky’s scheme coincides with the one
obtained with the parent dipole prescription, up to the above mentioned scaling violations. It
has been observed in [3, 46, 48] that the scaling function differs significantly in the fixed and
running coupling cases. Following that work, and to make a more quantitative study of the
scaling property, we fitted our solutions to the functional form [34]
f(τ) = a τ 2 γ
ln τ 2 + b
, (46)
with a, b and γ free parameters, within a fixed window below the saturation region, τ ∈
[10−5, 0.1]. Noticeably, at large enough rapidities the whole fitting window lays within the
geometric scaling window proposed in [36]: (Λ/Qs(Y )) < τ < 1, where Λ is some initial scale.
-210 -110 1
MV init. cond.   Y=0,40
Balitsky
parent dipole
fixed coupling
init. cond.
τ -110 1 10
AN06 init. cond.   Y=0,40
Figure 5: Solutions of the BK equation at rapidities Y=0 and 40 for KW (solid line), Balitsky
(dashed line) and parent dipole (dashed-dotted lines) schemes plotted versus the scaling variable
τ = rQs(Y ). The asymptotic solution obtained with fixed coupling αs = 0.2 at Y = 40 in [3]
is shown (black dashed-dotted line) for comparison. The initial conditions are MV (left) and
AN06 (right).
The value of γ extracted from the fits at rapidity Y = 40 lays in between γ ∼ 0.8 and γ ∼ 0.9.
This conclusion holds for the three initial conditions used here: the anomalous dimension
seems to converge to some intermediate value, in agreement with the value found in [3], for
asymptotic running coupling solutions (γ ∼ 0.85 at Y = 70). This result for anomalous
dimension is very far away from the value obtained in [3] for fixed coupling solutions (γ ∼ 0.64
at Y = 70) and from the predicted anomalous dimension for both running and fixed coupling
solutions from analytical studies of the equation based on saddle point techniques [34–38],
γc = χ(γc)/χ
′(γc) = 0.6275, where χ is the leading-log BFKL kernel.
It might be argued that the numerical value of the anomalous dimension extracted from our
fits is conditioned by the choice of the fitting function and by the fitting interval. Actually, it
was shown in [3] that the solutions of the evolution could be well fitted by other functional forms,
including the double-leading-log solution of BFKL, within a similar fitting region to the one
considered in this work. On the other hand, several phenomenological parameterizations of the
solution of the evolution have been proposed in [54–57] and have successfully confronted HERA
10 -410
10 -210 -110 1
 2*0.6τ ~ 
Fixed coupling
 2*0.85τ ~ 
Running coupling
Figure 6: Asymptotic solutions (Y=40) of the evolution equation for running coupling (solid
line) and fixed coupling with αs = 0.2 (dashed line). A fit to a power-law function aτ
2γ in the
region τ ∈ [10−6, 10−2] yields γ ≈ 0.85 for the running coupling solution and γ ≈ 0.6 for the
fixed coupling one.
and RHIC experimental data. There, the dipole scattering amplitude at arbitrary rapidity is
assumed to be given by a functional form analogous to our ansatz for the initial condition
Eq. (43), but allowing for geometric scaling violations by replacing γ → γ(r, Y ). The value
of the anomalous dimension at r = 1/Qs and/or for Y → ∞ is fixed to be the BFKL saddle
point, γc ∼ 0.63 (the saddle point value considered in [57] is slightly different, γ ∼ 0.53), while
the value γ = 1 is recovered in the limit r → ∞ at any finite rapidity. The success of these
phenomenological works supports the claim that the anomalous dimension of the solution is
given by the BFKL saddle point, in agreement with the above mentioned analytical predictions.
However, the relevant values of momenta probed at current phenomenological applications are
very distinct from the fitting region considered here. For example, the inclusive structure
function measured in HERA is fitted in [54,56] within the region 0.045 GeV2 < Q2 < 45 GeV2,
whereas charged hadron pt spectra in dAu collisions is well reproduced by [55–57] in the region
1 GeV < pt < 4.5 GeV. Note that, for both sets of data, the measured regions overlap with
the deeply saturated domain of the solution. On the contrary, our fitting region 10−5 < τ <
1 corresponds to values of momenta ∼ 10Qs(Y ) < pt < 105Qs(Y ) (always well above the
saturation scale), with Qs(Y = 40) ∼ 500÷1000 GeV for the different running coupling schemes
considered and, therefore, has no overlap with the kinematic regions measured experimentally,
since we scrutinize a momentum region strongly shifted to the ultraviolet compared to currently
available data. Moreover, it should be noticed that the rapidity interval covered by both sets
of experimental data is ∆Y < 4 in both cases, while we study the solutions of the evolution
(Y=0)]MVSub[N
-410 -310 -210 -110 1
(Y=30)]MVSub[N
(Y=0)]
Sub[N
10 -210 -110 1
(Y=30)]
Sub[N
(Y=0)]
Sub[N
10 -210 -110 1
(Y=30)]
Sub[N
Figure 7: Subtraction contribution calculated in the KW scheme (triangles) and in Balitsky’s
(stars). The trial functions correspond to the solutions of the evolution under Balitsky run-
ning coupling scheme at rapidities Y = 0, 30 for MV (left), AN08 (center) and AN06 initial
conditions.
at asymptotic rapidities, Y ∼ 40. We have checked that shifting our fitting region to larger
values of τ (smaller momentum) would bring the value of γ extracted from our fits closer to the
saddle point BFKL one, since the transition from the ultraviolet region to the deeply saturated
domain of the scaling solution is realized by a locally less steeper function (see Figs. (5) and
(12)). Therefore, there is no contradiction at all between the success of the phenomenological
parameterizations of the solutions and the results reported here.
With the above clarifications we reach the following conclusion: the asymptotic scaling
solutions corresponding to fixed and running coupling evolution are intrinsically different in
the whole r-range. This is emphasized in Fig. 6, where we represent the scaling solutions in
a log scale for τ < 1. It is clear that the tail of the distribution falls off with decreasing τ
much steeper for the running coupling solution than for the fixed coupling one. A fit to a
pure power-law function, f = a τ 2γ , in the region τ ∈ [10−6, 10−2] yields γ ∼ 0.85 for the
running coupling and γ ∼ 0.61 for the fixed coupling solution. The differences between fixed
and running coupling solutions at τ > 1 are evident from Fig. 5. This is a puzzling result that
remains to be understood from purely analytical methods.
(Y=0)]MV[N
10 -210 -110 1
(Y=30)]MV[N
(Y=0)]
10 -210 -110 1
(Y=30)]AN08[N
(Y=0)]
10 -210 -110 1
(Y=30)]AN06[N
Figure 8: Ratio of the subtraction over the running terms, D(r, Y ) = S[N(r, Y )]/R[N(r, Y )],
calculated in both KW (triangles) and Balitsky (stars) schemes for MV (left), AN08 (middle)
and AN06 (right) initial conditions at rapidities Y=0 (top) and Y=30 (bottom).
4.3 Subtraction Term
Before attempting to solve the complete evolution equation, and in order to gain insight in the
nature and structure of the subtraction contribution, we first evaluate the subtraction functional
for both Balitsky, Eq. (40), and KW, Eq. (41), schemes using a set of trial functions for S which
we choose to consist of the solutions of the evolution equation with the running coupling in
Balitsky’s scheme at different rapidities and of the three initial conditions considered above in
this work.
Two main remarks can be made about our results, shown in Fig. 7:
i) For all the trial functions considered in this work, the subtraction contribution is much
larger in the KW scheme than in Balitsky’s. A plausible explanation for this is that Balitsky’s
subtraction contribution, Eq. (40), when expanded in terms of dipole scattering amplitudes,
N = 1 − S, reduces to a sum of non-linear terms, since all the linear terms in the expansion
cancel each other due to the z1 ↔ z2 symmetry of the kernel, whereas in the KW case no such
cancellation happens and the subtraction contribution, Eq. (41), also includes linear terms,
which are dominant over the non-linear ones in the non-saturated domain where N ≪ 1.
ii) The subtraction contribution S has the same sign as the running coupling contribution
R in the whole τ range which, together with the relative minus sign assigned to the subtraction
term in Eq. (34), implies that the proper inclusion of the subtraction term reduces the value
of the functional that governs the evolution, F . In other words: the subtraction contribution
tends to systematically slow down the evolution, as we shall explicitly confirm in the next
subsection.
To better quantify the size of the subtraction contribution, we plot the ratio D(r, Y ) ≡
S[N(r, Y )]/R[N(r, Y )] in Fig. 8. At Y = 0, the relative weight of the subtraction contribution
with respect to the running one within the KW scheme and for a MV initial condition goes
from a D ∼ 0.4 at small τ to D ∼ 1 at τ ∼ 1. The same ratio for the Balitsky scheme takes
significantly smaller values: it goes from D ∼ 0.1 at small τ to D ∼ 0.4 for τ ∼ 1. As the
evolved solutions get closer to the scaling function, i.e. for larger rapidities, the r dependence of
the ratio becomes flatter and its overall normalization goes down to an approximately constant
value D ∼ 0.15 for the KW scheme and D ∼ 0.025 for that of Balitsky. This behavior remains
unaltered when going from rapidity Y = 20 to Y = 30, which suggests that the ratio may
saturate to a fixed value in the asymptotic region.
-210 -110 1
KWRun
BalRun
KW(Run-Sub)
(Run-Sub)
(Y=0)MVN
τ -210 -110 1
(Y=30)MVN
Figure 9: Total kernel F = R−S calculated under Balitsky’s scheme, Eqs (7) and (40), (solid
line) and under the KW scheme, Eqs (8) and (41), (dashed line). The overlap of the two lines
shows the agreement between the two calculations. Triangles stand for the running coupling
term calculated in the KW approach, RKW, while stars stand for the running coupling term
under Balitsky’s scheme, RBal. The trial functions N(r, Y ) correspond to the solution of the
evolution with only running coupling under Balitsky’s scheme at Y=0 (left) and Y=30 (right)
for a MV initial condition.
-110 1
1.2 MV i.c.
Y=0,3,10
KWRun
BalRun
F=Run-Sub
init. cond.
r )-1(GeV
-110 1
scaling function i.c.
Y=0,3,10
r )-1(GeV
Figure 10: Solutions of the complete (all orders in αs β2) evolution equation given in Eq. (34)
(solid lines), and of the equation with Balitsky’s (dashed lines) and KW’s (dashed-dotted)
running coupling schemes at rapidities Y = 0, 5 and 10. Left plot uses MV initial condition.
The right plot employs the initial condition given by the dipole amplitude at rapidity Y = 35
evolved using Balitsky’s running coupling scheme and with r-dependence rescaled down such
that Qs = Q
s = 1 GeV.
Finally, we have checked that combining the subtraction and running coupling contributions
for both schemes adds up to the same result. This is shown in Fig. 9, where we plot the value
of the total functional F = R − S calculated under the KW scheme (Eqs. (8) and (36) for
the running coupling term, R, and Eq. (41) for the subtraction term, S) and under Balitsky’s
scheme (Eqs. (7) and (35) for the running coupling term and Eq. (40) for the subtraction term).
The two results coincide within the estimation of the numerical accuracy previously discussed.
The agreement between the two results is better in the small-τ region, where the two curves lay
almost on top of each other. In the saturation region, τ & 1, the agreement is slightly worse,
although the differences between the values of F calculated in both schemes is still much less
than the differences between the running coupling terms themselves. This slight remaining
disagreement between the Balitsky’s and KW prescriptions may also be due to inaccuracies in
a Fourier transform of a geometric series performed in arriving at Eq. (39). This result serves
as a cross-check of our numerical method and as an additional confirmation of the agreement
of the independent calculations derived in [1, 2].
0 1 2 3 4 5 6 7 8 9
MV i.c.KWRun
BalRun
F=Run-Sub
)2(GeV
1 2 3 4 5 6 7 8 9
scaling function i.c.
Figure 11: Saturation scale corresponding to the solutions plotted in Fig. 10.
4.4 Complete running coupling BK equation
In this section we calculate the solutions of the complete evolution equation, Eq. (34), including
both the running and subtraction terms obtained by the all-orders αs Nf resummation and by
the Nf → −6πβ2 replacement. Since the numerical evaluation of the subtraction contribution
at each point of the grid and each step of the evolution would require an exceedingly large
amount of CPU time consumption, the strategy followed to include it in the evolution equation
consists of calculating such contribution only in a small set of grid points at each step of the
evolution, which we fixed at n = 16, between the points r1 and r2, which are determined at each
step of the evolution by the conditions N(Y, r1) = 10
−9, and N(Y, r2) = 0.99, and then using
power-law interpolation and extrapolation to the other points of the grid. Both the running
and subtracted terms are calculated within Balitsky scheme. This procedure is motivated by
the fact that, as discussed in the previous section, the subtraction contribution can be regarded
as a small perturbation with respect to the running coupling term within Balitsky’s scheme and
by the fact that it is a rather smooth function that can be well fitted by a power-law function in
most of the r-range. The accuracy of this procedure has been checked by doubling the number
of points at which the subtraction contribution is calculated at each step of the evolution, i.e.
by setting n=32. At Y=2, the differences between the solutions obtained with the two above
mentioned choices for n were less than a 8% in the tail of the solution, r < r1, and less than a
3% for r > r1.
The results of the evolution calculated in this way and using MV and rescaled asymptotic
running coupling solution (Y=35) as initial conditions are plotted in Fig. 10. They confirm
-210 -110 1 10
Y=10F=Run-Sub
(Y=35)  MV BALi.c.=N
Figure 12: Rescaled solutions given by the complete αs β2-evolution equation (solid line) and for
KW (dashed-dotted line) and Balitsky’s (dashed line) running coupling schemes at Y = 10. The
initial condition corresponds to the dipole amplitude at rapidity Y = 35 evolved using Balitsky’s
running coupling scheme and with r-dependence rescaled down such that Qs = Q
s = 1 GeV.
the expectations raised in the previous Subsection: the inclusion of the subtraction terms
considerably slows down the evolution with respect to the sole consideration of the running
coupling contributions. Moreover, the reduction in the speed of the wave front is much larger
for the KW scheme than for that of Balitsky one for both initial conditions. However, the closer
the initial condition is to the asymptotic running coupling scaling function, the smaller are the
effects of the subtraction contribution. These features can be better quantified by inspecting
the rapidity dependence of the saturation scale generated by the evolution, plotted in Fig. 11.
At rapidity Y = 10 the ratio of the saturation scale Qs yielded by the KW scheme to Qs given
by the complete αs β2-evolution equation is a factor of ∼ 2.5 for the MV initial condition and
a factor of ∼ 2.1 for the asymptotic running coupling initial condition. At the same rapidity,
the ratio of the saturation scale obtained under Balitsky’s scheme to Qs corresponding to the
complete αs β2-evolution is ∼ 1.25 for the MV initial condition and ∼ 1.15 for the scaling
function initial condition. Thus, in spite of the smallness of the ratio of the subtraction terms
to the running coupling contributions at high rapidity, which is ∼ 0.025 for Balitsky’s and
∼ 0.15 for KW scheme at Y = 30 (see bottom plots in Fig. 8), the proper inclusion of the
subtraction term results in fairly sizable effects in the solutions of the evolution equation.
Finally, we notice that the scaling behavior of the solution is not affected by the subtraction
term. This is seen in figure Fig. 12, where we evolve starting from an initial condition already
close to the running coupling scaling function and plot the solutions of the evolution equation
obtained with just running coupling terms (see Section 4.1) and the solution of the complete
αs β2-evolution at rapidity Y = 10. It is clear that, within the numerical accuracy, no departure
from the scaling behavior is observed. Therefore the main effect of a proper consideration of the
subtraction term is the one of reducing the speed of the evolution. It does not violate or modify
the geometric scaling property of the solutions established in Section 4.2. In our understanding
geometric scaling appears to persist when the running coupling effects are included because,
at high enough rapidity Qs(Y ) ≫ Λ, such that the new (from the LO standpoint) momentum
scale Λ introduced by the running coupling can be safely neglected. Hence the dynamics is
again characterized by a single momentum scale Qs(Y ). At the same time running coupling
does modify the evolution kernel, leading to a different shape of the scaling function.
5 Conclusions
In this paper we have taken into account all corrections to the kernels of the non-linear JIMWLK
and BK evolution equations containing powers of αs Nf . We reiterated the fact that the sep-
aration of the resulting kernel resumming all powers of αs Nf into the running coupling and
subtraction parts, as done in the previous calculations of [1, 2], is not justified parametrically.
We have then performed numerical analysis with the following conclusions.
• First we solved the evolution equations derived in [1] and [2] keeping only the running
coupling part or the evolution kernel and neglecting the subtraction term. Comparing to
the results for fixed coupling obtained in [3] we confirmed the conclusion reached in [3]
that the growth with rapidity is substantially reduced when running coupling corrections
are included. The results for three different initial conditions are shown in Fig. 4. We
observe that the solution of the equation derived in [2] differs significantly from that
derived in [1], but agrees (with good numerical accuracy) with the solution of the BK
evolution equation with the coupling running at the parent dipole size. (The latter is just
a model of the running coupling not resulting from any calculations, which we plot for
illustrative purposes.) We also observe that at sufficiently high rapidity both equations
from [1] and from [2] give us the same scaling function for the dipole amplitude N(r, Y )
as a function of r Qs(Y ), which is also in agreement with the scaling function given by
the parent dipole running, as shown in Fig. 5. The fact that the scaling is preserved when
the running coupling corrections are included was previously established in [3], though
for models of running coupling only. The shape of the scaling function is very different
from that obtained from the fixed coupling evolution equations. In particular, we found
that for dipole sizes below 0.1/Qs the anomalous dimension of the scaling function in the
running coupling case becomes γ ≈ 0.85 (see Fig. 6). This is different from the result of
several analytical estimates [34–38], which expect the anomalous dimension not to change
when running coupling corrections are included and to remain at its fixed coupling value
of γ ≈ 0.63.
• We have then evaluated the subtraction term for both calculations performed in [1] and
[2]. We demonstrated that subtracting the subtraction terms from the running coupling
terms makes the full answer agree for both calculations of [1] and [2], as shown in Fig. 9
for the right hand side of the evolution equation. It turns out that the subtraction
term SBal[S], which has to be subtracted from the result of [1], is systematically smaller
than SKW[S], to be subtracted from the result of [2], over the whole rapidity range
studied here. This implies that the result of [1] should have a smaller correction than
the result of [2] and is thus closer to the full answer. The subtraction terms SBal[S]
and SKW[S] are plotted in Fig. 7 as functions of the dipole size r for different values of
rapidity. Their relative contributions to the evolution kernel are shown in Fig. 8, where
we plotted the subtraction functional divided by the running coupling functional. From
those figures we conclude that both the magnitude of these extra terms and their relative
contribution to the evolution kernel decrease with increasing rapidity. Hence, while at
”moderate” rapidities (the ones closer to realistic experimental values) the subtraction
term is important for both calculations [1, 2], it becomes increasingly less important at
asymptotically large rapidities. The physics is easy to understand: the subtraction terms
are o(α2s), while the running coupling part of the kernel is o(αs). Hence, if we suppose
that the effective value of the coupling is given by its magnitude at the saturation scale
Qs(Y ), then, as rapidity increases, the coupling would decrease, making the subtraction
term much smaller than the running coupling term. Indeed, while at asymptotically high
rapidities the assumption of [1,2] that the subtraction term could be neglected is justified,
making the results of [1] and [2] agree with each other, for rapidities relevant to modern
days experiments the subtraction term is numerically important.
• With the last conclusion in mind we continued by numerically solving the full evolution
equation resumming all powers of αs Nf in the evolution kernel, which now would combine
both the running coupling and the subtraction terms. The five-dimensional integral in
the subtraction term (37) made obtaining this solution rather difficult. The outcome of
the calculation is shown in Fig. 10. All the main conclusions stated above were again
confirmed by the solution of the full equation. At asymptotically high rapidity scaling
regime is recovered, as can be seen from Fig. 12. As the subtraction term is less important
in that regime, the scaling function appears to be the same as in the case of having only
the running coupling term in the kernel. The anomalous dimension again turns out
to be γ ≈ 0.85, in disagreement with the analytical expectations of [34–38]. However,
the scaling of the saturation scale with rapidity appears to be in agreement with the
expectations of analytical work of [34, 35, 38], as shown in Fig. 11.
We conclude by observing that the knowledge of the non-linear small-x evolution equation
with all the running coupling corrections included brings us to an unprecedented level of pre-
cision allowing for a much more detailed comparison with experiments than was ever possible
before.
Acknowledgments
We would like to thank Heribert Weigert for many informative and helpful discussions at
the beginning of this work. A portion of the performed work was motivated by stimulating
discussions with Robi Peschanski, which we gratefully acknowledge.
This research is sponsored in part by the U.S. Department of Energy under Grant No. DE-
FG02-05ER41377. This work was supported in part by an allocation of computing time from
the Ohio Supercomputer Center.
References
[1] I. I. Balitsky, Quark Contribution to the Small-x Evolution of Color Dipole, Phys. Rev. D
75 (2007) 014001, [hep-ph/0609105].
[2] Y. Kovchegov and H. Weigert, Triumvirate of Running Couplings in Small-x Evolution,
Nucl. Phys. A 784 (2007) 188–226, [hep-ph/0609090].
[3] J. L. Albacete, N. Armesto, J. G. Milhano, C. A. Salgado, and U. A. Wiedemann,
Numerical analysis of the Balitsky-Kovchegov equation with running coupling:
Dependence of the saturation scale on nuclear size and rapidity, Phys. Rev. D71 (2005)
014003, [hep-ph/0408216].
[4] E. A. Kuraev, L. N. Lipatov, and V. S. Fadin, The Pomeranchuk singularity in
non-Abelian gauge theories, Sov. Phys. JETP 45 (1977) 199–204.
[5] Y. Y. Balitsky and L. N. Lipatov Sov. J. Nucl. Phys. 28 (1978) 822.
[6] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, The BFKL equation from
the Wilson renormalization group, Nucl. Phys. B504 (1997) 415–431, [hep-ph/9701284].
[7] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, The Wilson renormalization
group for low x physics: Towards the high density regime, Phys. Rev. D59 (1998) 014014,
[hep-ph/9706377].
[8] J. Jalilian-Marian, A. Kovner, and H. Weigert, The Wilson renormalization group for low
x physics: Gluon evolution at finite parton density, Phys. Rev. D59 (1998) 014015,
[hep-ph/9709432].
[9] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Unitarization of gluon
distribution in the doubly logarithmic regime at high density, Phys. Rev. D59 (1999)
034007, [hep-ph/9807462].
[10] A. Kovner, J. G. Milhano, and H. Weigert, Relating different approaches to nonlinear
QCD evolution at finite gluon density, Phys. Rev. D62 (2000) 114005, [hep-ph/0004014].
[11] H. Weigert, Unitarity at small Bjorken x, Nucl. Phys. A703 (2002) 823–860,
[hep-ph/0004044].
[12] E. Iancu, A. Leonidov, and L. D. McLerran, Nonlinear gluon evolution in the color glass
condensate. I, Nucl. Phys. A692 (2001) 583–645, [hep-ph/0011241].
[13] E. Ferreiro, E. Iancu, A. Leonidov, and L. McLerran, Nonlinear gluon evolution in the
color glass condensate. II, Nucl. Phys. A703 (2002) 489–538, [hep-ph/0109115].
[14] I. Balitsky, Operator expansion for high-energy scattering, Nucl. Phys. B463 (1996)
99–160, [hep-ph/9509348].
[15] I. Balitsky, Operator expansion for diffractive high-energy scattering, hep-ph/9706411.
http://arxiv.org/abs/hep-ph/0609105
http://xxx.lanl.gov/abs/hep-ph/0609105
http://arxiv.org/abs/hep-ph/0609090
http://xxx.lanl.gov/abs/hep-ph/0609090
http://arxiv.org/abs/hep-ph/0408216
http://xxx.lanl.gov/abs/hep-ph/0408216
http://arxiv.org/abs/hep-ph/9701284
http://xxx.lanl.gov/abs/hep-ph/9701284
http://arxiv.org/abs/hep-ph/9706377
http://xxx.lanl.gov/abs/hep-ph/9706377
http://arxiv.org/abs/hep-ph/9709432
http://xxx.lanl.gov/abs/hep-ph/9709432
http://arxiv.org/abs/hep-ph/9807462
http://xxx.lanl.gov/abs/hep-ph/9807462
http://arxiv.org/abs/hep-ph/0004014
http://xxx.lanl.gov/abs/hep-ph/0004014
http://arxiv.org/abs/hep-ph/0004044
http://xxx.lanl.gov/abs/hep-ph/0004044
http://arxiv.org/abs/hep-ph/0011241
http://xxx.lanl.gov/abs/hep-ph/0011241
http://arxiv.org/abs/hep-ph/0109115
http://xxx.lanl.gov/abs/hep-ph/0109115
http://arxiv.org/abs/hep-ph/9509348
http://xxx.lanl.gov/abs/hep-ph/9509348
http://arxiv.org/abs/hep-ph/9706411
http://xxx.lanl.gov/abs/hep-ph/9706411
[16] I. Balitsky, Factorization and high-energy effective action, Phys. Rev. D60 (1999)
014020, [hep-ph/9812311].
[17] Y. V. Kovchegov, Small-x F2 structure function of a nucleus including multiple pomeron
exchanges, Phys. Rev. D60 (1999) 034008, [hep-ph/9901281].
[18] Y. V. Kovchegov, Unitarization of the BFKL pomeron on a nucleus, Phys. Rev. D61
(2000) 074018, [hep-ph/9905214].
[19] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Singlet structure function at small x:
Unitarization of gluon ladders, Nucl. Phys. B188 (1981) 555–576.
[20] A. H. Mueller and J.-w. Qiu, Gluon recombination and shadowing at small values of x,
Nucl. Phys. B268 (1986) 427.
[21] L. D. McLerran and R. Venugopalan, Green’s functions in the color field of a large
nucleus, Phys. Rev. D50 (1994) 2225–2233, [hep-ph/9402335].
[22] L. D. McLerran and R. Venugopalan, Gluon distribution functions for very large nuclei
at small transverse momentum, Phys. Rev. D49 (1994) 3352–3355, [hep-ph/9311205].
[23] L. D. McLerran and R. Venugopalan, Computing quark and gluon distribution functions
for very large nuclei, Phys. Rev. D49 (1994) 2233–2241, [hep-ph/9309289].
[24] Y. V. Kovchegov, Non-Abelian Weizsaecker-Williams field and a two- dimensional
effective color charge density for a very large nucleus, Phys. Rev. D54 (1996) 5463–5469,
[hep-ph/9605446].
[25] Y. V. Kovchegov, Quantum structure of the non-Abelian Weizsaecker-Williams field for a
very large nucleus, Phys. Rev. D55 (1997) 5445–5455, [hep-ph/9701229].
[26] J. Jalilian-Marian, A. Kovner, L. D. McLerran, and H. Weigert, The intrinsic glue
distribution at very small x, Phys. Rev. D55 (1997) 5414–5428, [hep-ph/9606337].
[27] E. Iancu and R. Venugopalan, The color glass condensate and high energy scattering in
QCD, hep-ph/0303204.
[28] H. Weigert, Evolution at small xbj: The Color Glass Condensate, Prog. Part. Nucl. Phys.
55 (2005) 461–565, [hep-ph/0501087].
[29] J. Jalilian-Marian and Y. V. Kovchegov, Saturation physics and deuteron gold collisions
at RHIC, Prog. Part. Nucl. Phys. 56 (2006) 104–231, [hep-ph/0505052].
[30] E. Gardi, J. Kuokkanen, K. Rummukainen, and H. Weigert, Running coupling and power
corrections in nonlinear evolution at the high-energy limit, Nucl. Phys. A784 (2007)
282–340, [hep-ph/0609087].
[31] Y. V. Kovchegov and H. Weigert, Quark loop contribution to BFKL evolution: Running
coupling and leading-N(f) NLO intercept, Nucl. Phys. A789, 260 (2007),
[hep-ph/0612071].
http://arxiv.org/abs/hep-ph/9812311
http://xxx.lanl.gov/abs/hep-ph/9812311
http://arxiv.org/abs/hep-ph/9901281
http://xxx.lanl.gov/abs/hep-ph/9901281
http://arxiv.org/abs/hep-ph/9905214
http://xxx.lanl.gov/abs/hep-ph/9905214
http://arxiv.org/abs/hep-ph/9402335
http://xxx.lanl.gov/abs/hep-ph/9402335
http://arxiv.org/abs/hep-ph/9311205
http://xxx.lanl.gov/abs/hep-ph/9311205
http://arxiv.org/abs/hep-ph/9309289
http://xxx.lanl.gov/abs/hep-ph/9309289
http://arxiv.org/abs/hep-ph/9605446
http://xxx.lanl.gov/abs/hep-ph/9605446
http://arxiv.org/abs/hep-ph/9701229
http://xxx.lanl.gov/abs/hep-ph/9701229
http://arxiv.org/abs/hep-ph/9606337
http://xxx.lanl.gov/abs/hep-ph/9606337
http://arxiv.org/abs/hep-ph/0303204
http://xxx.lanl.gov/abs/hep-ph/0303204
http://arxiv.org/abs/hep-ph/0501087
http://xxx.lanl.gov/abs/hep-ph/0501087
http://arxiv.org/abs/hep-ph/0505052
http://xxx.lanl.gov/abs/hep-ph/0505052
http://arxiv.org/abs/hep-ph/0609087
http://xxx.lanl.gov/abs/hep-ph/0609087
http://arxiv.org/abs/hep-ph/0612071
http://xxx.lanl.gov/abs/hep-ph/0612071
[32] G. P. Lepage and S. J. Brodsky, Exclusive processes in perturbative quantum
chromodynamics, Phys. Rev. D22 (1980) 2157.
[33] S. J. Brodsky, H.-C. Pauli, and S. S. Pinsky, Quantum chromodynamics and other field
theories on the light cone, Phys. Rept. 301 (1998) 299–486, [hep-ph/9705477].
[34] A. H. Mueller and D. N. Triantafyllopoulos, The energy dependence of the saturation
momentum, Nucl. Phys. B640 (2002) 331–350, [hep-ph/0205167].
[35] D. N. Triantafyllopoulos, The energy dependence of the saturation momentum from RG
improved BFKL evolution, Nucl. Phys. B648 (2003) 293–316, [hep-ph/0209121].
[36] E. Iancu, K. Itakura, and L. McLerran, Geometric scaling above the saturation scale,
Nucl. Phys. A708 (2002) 327–352, [hep-ph/0203137].
[37] S. Munier and R. Peschanski, Universality and tree structure of high energy QCD, Phys.
Rev. D70 (2004) 077503, [hep-ph/0401215].
[38] G. Beuf and R. Peschanski, Universality of QCD traveling-waves with running coupling,
hep-ph/0702131.
[39] A. H. Mueller, Soft gluons in the infinite momentum wave function and the BFKL
pomeron, Nucl. Phys. B415 (1994) 373–385.
[40] A. H. Mueller and B. Patel, Single and double BFKL pomeron exchange and a dipole
picture of high-energy hard processes, Nucl. Phys. B425 (1994) 471–488,
[hep-ph/9403256].
[41] A. H. Mueller, Unitarity and the BFKL pomeron, Nucl. Phys. B437 (1995) 107–126,
[hep-ph/9408245].
[42] Z. Chen and A. H. Mueller, The dipole picture of high-energy scattering, the BFKL
equation and many gluon compound states, Nucl. Phys. B451 (1995) 579–604.
[43] S. Munier and R. Peschanski, Traveling wave fronts and the transition to saturation,
Phys. Rev. D69 (2004) 034008, [hep-ph/0310357].
[44] S. J. Brodsky, G. P. Lepage, and P. B. Mackenzie, On the elimination of scale
ambiguities in perturbative quantum chromodynamics, Phys. Rev. D28 (1983) 228.
[45] K. Golec-Biernat and M. Wüsthoff, Saturation effects in deep inelastic scattering at low
Q2 and its implications on diffraction, Phys. Rev. D59 (1998) 014017, [hep-ph/9807513].
[46] M. A. Braun, Pomeron fan diagrams with an infrared cutoff and running coupling, Phys.
Lett. B 576 (2003) 115, [hep-ph/0308320].
[47] S. Munier and R. Peschanski, Geometric scaling as traveling waves, Phys. Rev. Lett. 91
(2003) 232001, [hep-ph/0309177].
http://arxiv.org/abs/hep-ph/9705477
http://xxx.lanl.gov/abs/hep-ph/9705477
http://arxiv.org/abs/hep-ph/0205167
http://xxx.lanl.gov/abs/hep-ph/0205167
http://arxiv.org/abs/hep-ph/0209121
http://xxx.lanl.gov/abs/hep-ph/0209121
http://arxiv.org/abs/hep-ph/0203137
http://xxx.lanl.gov/abs/hep-ph/0203137
http://arxiv.org/abs/hep-ph/0401215
http://xxx.lanl.gov/abs/hep-ph/0401215
http://arxiv.org/abs/hep-ph/0702131
http://xxx.lanl.gov/abs/hep-ph/0702131
http://arxiv.org/abs/hep-ph/9403256
http://xxx.lanl.gov/abs/hep-ph/9403256
http://arxiv.org/abs/hep-ph/9408245
http://xxx.lanl.gov/abs/hep-ph/9408245
http://arxiv.org/abs/hep-ph/0310357
http://xxx.lanl.gov/abs/hep-ph/0310357
http://arxiv.org/abs/hep-ph/9807513
http://xxx.lanl.gov/abs/hep-ph/9807513
http://arxiv.org/abs/hep-ph/0308320
http://xxx.lanl.gov/abs/hep-ph/0308320
http://arxiv.org/abs/hep-ph/0309177
http://xxx.lanl.gov/abs/hep-ph/0309177
[48] J. L. Albacete, N. Armesto, A. Kovner, C. A. Salgado, and U. A. Wiedemann, Energy
dependence of the Cronin effect from non-linear QCD evolution, Phys. Rev. Lett. 92
(2004) 082001, [hep-ph/0307179].
[49] M. Lublinsky, Scaling phenomena from non-linear evolution in high energy DIS, Eur.
Phys. J. C21 (2001) 513–519, [hep-ph/0106112].
[50] N. Armesto and M. A. Braun, Parton densities and dipole cross-sections at small x in
large nuclei, Eur. Phys. J. C20 (2001) 517–522, [hep-ph/0104038].
[51] A. M. Stasto, K. Golec-Biernat, and J. Kwiecinski, Geometric scaling for the total γ∗p
cross-section in the low x region, Phys. Rev. Lett. 86 (2001) 596–599, [hep-ph/0007192].
[52] N. Armesto, C. A. Salgado, and U. A. Wiedemann, Relating high-energy lepton hadron,
proton nucleus and nucleus nucleus collisions through geometric scaling,
hep-ph/0407018.
[53] J. L. Albacete, N. Armesto, J. G. Milhano, C. A. Salgado, and U. A. Wiedemann,
Nuclear size and rapidity dependence of the saturation scale from QCD evolution and
experimental data, Eur. Phys. J. C43 (2005) 353–360, [hep-ph/0502167].
[54] E. Iancu, K. Itakura, and S. Munier, Saturation and BFKL dynamics in the HERA data
at small x, Phys. Lett. B590 (2004) 199–208, [hep-ph/0310338].
[55] A. Dumitru, A. Hayashigaki, and J. Jalilian-Marian, Geometric scaling violations in the
central rapidity region of d + Au collisions at RHIC, Nucl. Phys. A770 (2006) 57–70,
[hep-ph/0512129].
[56] V. P. Goncalves, M. S. Kugeratski, M. V. T. Machado, and F. S. Navarra, Saturation
physics at HERA and RHIC: An unified description, Phys. Lett. B643 (2006) 273–278,
[hep-ph/0608063].
[57] D. Kharzeev, Y. V. Kovchegov, and K. Tuchin, Nuclear modification factor in d + Au
collisions: Onset of suppression in the color glass condensate, Phys. Lett. B599 (2004)
23–31, [hep-ph/0405045].
http://arxiv.org/abs/hep-ph/0307179
http://xxx.lanl.gov/abs/hep-ph/0307179
http://arxiv.org/abs/hep-ph/0106112
http://xxx.lanl.gov/abs/hep-ph/0106112
http://arxiv.org/abs/hep-ph/0104038
http://xxx.lanl.gov/abs/hep-ph/0104038
http://arxiv.org/abs/hep-ph/0007192
http://xxx.lanl.gov/abs/hep-ph/0007192
http://arxiv.org/abs/hep-ph/0407018
http://xxx.lanl.gov/abs/hep-ph/0407018
http://arxiv.org/abs/hep-ph/0502167
http://xxx.lanl.gov/abs/hep-ph/0502167
http://arxiv.org/abs/hep-ph/0310338
http://xxx.lanl.gov/abs/hep-ph/0310338
http://arxiv.org/abs/hep-ph/0512129
http://xxx.lanl.gov/abs/hep-ph/0512129
http://arxiv.org/abs/hep-ph/0608063
http://xxx.lanl.gov/abs/hep-ph/0608063
http://arxiv.org/abs/hep-ph/0405045
http://xxx.lanl.gov/abs/hep-ph/0405045
	Introduction
	Scheme dependence
	Inclusion of running coupling corrections: general concepts
	Derivation of the subtraction term
	Brief summary of analytical results
	Numerical setup and initial conditions
	Results
	Running coupling
	Geometric scaling
	Subtraction Term
	Complete running coupling BK equation
	Conclusions
