Solving The High Energy Evolution Equation Including Running Coupling Corrections

TL;DR

This paper numerically solves the nonlinear BK evolution equation with recent running coupling corrections, revealing a universal scaling behavior and emphasizing the importance of subtraction terms for accurate results.

Contribution

It introduces a comprehensive numerical solution of the BK equation including the subtraction terms, clarifying their impact on the evolution and scaling behavior.

Findings

Running coupling corrections significantly alter the shape of the scaling function.

Subtraction terms are crucial for consistent evolution kernel calculations.

Universal scaling behavior persists at large rapidity with full corrections.

Abstract

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.

Figures

Equations

$$\frac{\partial S({\underline{x}}_0,{\underline{x}}_1;Y)}{\partial Y}=\int d^{2}zK({\underline{x}}_0,{\underline{x}}_1,\underline{z})\left[S({\underline{x}}_0,\underline{z};Y)S(\underline{z},{\underline{x}}_1;Y)-S({\underline{x}}_0,{\underline{x}}_1;Y)\right],$$

$$K({\underline{x}}_0,{\underline{x}}_1,\underline{z})=\frac{{\alpha }_s{N}_c}{2{\pi }^2}\frac{r^2}{r_1^2{r}_2^2}$$

$$S({\underline{x}}_0,{\underline{x}}_1;Y)=\frac{1}{N_c}⟨\text{mbox}tr\left\{V({\underline{x}}_0)V^{†}({\underline{x}}_1)\right\}⟩.$$

$$\frac{\partial S({\underline{x}}_0,{\underline{x}}_1;Y)}{\partial Y}=\mathcal{F}\left[S({\underline{x}}_0,{\underline{x}}_1;Y)\right]$$

$$\mathcal{F}\left[S\right]=\mathcal{R}\left[S\right]-\mathcal{S}\left[S\right].$$

$$\mathcal{R}\left[S({\underline{x}}_0,{\underline{x}}_1;Y)\right]=\int d^{2}z\tilde{K}({\underline{x}}_0,{\underline{x}}_1,\underline{z})\left[S({\underline{x}}_0,\underline{z};Y)S(\underline{z},{\underline{x}}_1;Y)-S({\underline{x}}_0,{\underline{x}}_1;Y)\right].$$

$${\tilde{K}}^{\text{Bal}}(\underline{r},{\underline{r}}_1,{\underline{r}}_2)=\frac{N_c{\alpha }_s(r^2)}{2{\pi }^2}\left[\frac{r^2}{r_1^2{r}_2^2}+\frac{1}{r_1^2}\left(\frac{{\alpha }_s(r_1^2)}{{\alpha }_s(r_2^2)}-1\right)+\frac{1}{r_2^2}\left(\frac{{\alpha }_s(r_2^2)}{{\alpha }_s(r_1^2)}-1\right)\right].$$

$${\tilde{K}}^{\text{KW}}(\underline{r},{\underline{r}}_1,{\underline{r}}_2)=\frac{N_c}{2{\pi }^2}\left[{\alpha }_s(r_1^2)\frac{1}{r_1^2}-2\frac{{\alpha }_s(r_1^2){\alpha }_s(r_2^2)}{{\alpha }_s(R^2)}\frac{{\underline{r}}_1\cdot {\underline{r}}_2}{r_1^2{r}_2^2}+{\alpha }_s(r_2^2)\frac{1}{r_2^2}\right],$$

$$R^2(\underline{r},{\underline{r}}_1,{\underline{r}}_2)=r_1{r}_2{\left(\frac{r_2}{r_1}\right)}^{\frac{r_1^2+r_2^2}{r_1^2-r_2^2}-2\frac{r_1^2{r}_2^2}{{\underline{r}}_1\cdot {\underline{r}}_2}\frac{1}{r_1^2-r_2^2}}.$$

$${\tilde{K}}^{pd}(\underline{r},{\underline{r}}_1,{\underline{r}}_2)=\frac{N_c{\alpha }_s(r^2)}{2{\pi }^2}\frac{r^2}{r_1^2{r}_2^2}.$$

$${\mathcal{K}}_1^{\text{NLO}}({\underline{x}}_0,{\underline{x}}_1;{\underline{z}}_1,{\underline{z}}_2)=$$

$$\displaystyle\times\left[\frac{1}{{\underline{q}}^{2}{\underline{q}}^{\prime 2}}\frac{(1-2\alpha)^{2}{\underline{q}}\cdot{\underline{k}}\ {\underline{k}}^{\prime}\cdot{\underline{q}}^{\prime}+{\underline{q}}\cdot{\underline{q}}^{\prime}\ {\underline{k}}\cdot{\underline{k}}^{\prime}-{\underline{q}}\cdot{\underline{k}}^{\prime}\ {\underline{k}}\cdot{\underline{q}}^{\prime}}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\right.$$

$$\displaystyle+\frac{2\,\alpha\,(1-\alpha)\,(1-2\alpha)}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\,\left(\frac{{\underline{k}}\cdot{\underline{q}}}{{\underline{q}}^{2}}+\frac{{\underline{k}}^{\prime}\cdot{\underline{q}}^{\prime}}{{\underline{q}}^{\prime 2}}\right)$$

$$\displaystyle\left.+\,\frac{4\,\alpha^{2}\,(1-\alpha)^{2}}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\right].$$

$$K_1^{\text{NLO}}({\underline{x}}_0,{\underline{x}}_1;{\underline{z}}_1,{\underline{z}}_2)=C_F\sum _{m,n=0}^1(-1{)}^{m+n}{\mathcal{K}}_1^{\text{NLO}}$$

$${\alpha }_{\mu }^2\int d^2{z}_1{d}^2{z}_2{K}_1^{\text{NLO}}({\underline{x}}_0,{\underline{x}}_1;{\underline{z}}_1,{\underline{z}}_2)S({\underline{x}}_0,{\underline{z}}_1,Y)S({\underline{z}}_2,{\underline{x}}_1,Y)$$

$${\mathcal{S}}_{\text{NLO}}[S]={\alpha }_{\mu }^2$$

$$\times$$

$$\underline{w}={\underline{z}}_1\text{mbox}or\underline{w}={\underline{z}}_2,$$

$${\mathcal{S}}_{\text{NLO}}^{\text{Bal}}[S]=$$

$$\times$$

$$\underline{w}=\underline{z}=\alpha {\underline{z}}_1+(1-\alpha ){\underline{z}}_2.$$

$${\mathcal{S}}_{\text{NLO}}^{\text{KW}}[S]=$$

$$\times$$

$$\displaystyle{\cal K}_{\begin{picture}(10.0,8.0)\put(5.0,5.0){\circle{8.0}} \put(2.9,2.5){{\scriptsize 1}} \end{picture}}({\underline{x}}_{0},{\underline{x}}_{1};{\underline{z}}_{1},{\underline{z}}_{2})\,=$$

$$\displaystyle\times\left[\frac{1}{{\underline{q}}^{2}{\underline{q}}^{\prime 2}}\frac{(1-2\alpha)^{2}{\underline{q}}\cdot{\underline{k}}\ {\underline{k}}^{\prime}\cdot{\underline{q}}^{\prime}+{\underline{q}}\cdot{\underline{q}}^{\prime}\ {\underline{k}}\cdot{\underline{k}}^{\prime}-{\underline{q}}\cdot{\underline{k}}^{\prime}\ {\underline{k}}\cdot{\underline{q}}^{\prime}}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\right.$$

$$\displaystyle+\frac{2\,\alpha\,(1-\alpha)\,(1-2\alpha)}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\,\left(\frac{{\underline{k}}\cdot{\underline{q}}}{{\underline{q}}^{2}}+\frac{{\underline{k}}^{\prime}\cdot{\underline{q}}^{\prime}}{{\underline{q}}^{\prime 2}}\right)$$

$$\displaystyle\left.+\,\frac{4\,\alpha^{2}\,(1-\alpha)^{2}}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\right]$$

$$\times \frac{1}{\left(1+{\alpha }_{\mu }{\beta }_2\ln \frac{{\underline{q}}^2{e}^{-5/3}}{{\mu }_{\overline{\text{MS}}}^2}\right)\left(1+{\alpha }_{\mu }{\beta }_2\ln \frac{{\underline{q}}^{\prime 2}{e}^{-5/3}}{{\mu }_{\overline{\text{MS}}}^2}\right)}$$

$${\beta }_2=\frac{11N_c-2N_f}{12\pi }.$$

$$\displaystyle{\cal K}_{\begin{picture}(10.0,8.0)\put(5.0,5.0){\circle{8.0}} \put(2.9,2.5){{\scriptsize 1}} \end{picture}}({\underline{x}}_{0},{\underline{x}}_{1};{\underline{z}}_{1},{\underline{z}}_{2})\,=$$

$$e^{-i\underline{q}\cdot (\underline{z}-{\underline{x}}_0)+i{\underline{q}}^{\prime }\cdot (\underline{z}-{\underline{x}}_1)-i(\underline{k}-{\underline{k}}^{\prime })\cdot {\underline{z}}_{12}}{\left(\frac{q^2}{{\mu }^2}\right)}^{\lambda }{\left(\frac{q^{\prime 2}}{{\mu }^2}\right)}^{{\lambda }^{\prime }}$$

$$\displaystyle\times\left[\frac{1}{{\underline{q}}^{2}{\underline{q}}^{\prime 2}}\frac{(1-2\alpha)^{2}{\underline{q}}\cdot{\underline{k}}\ {\underline{k}}^{\prime}\cdot{\underline{q}}^{\prime}+{\underline{q}}\cdot{\underline{q}}^{\prime}\ {\underline{k}}\cdot{\underline{k}}^{\prime}-{\underline{q}}\cdot{\underline{k}}^{\prime}\ {\underline{k}}\cdot{\underline{q}}^{\prime}}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\right.$$

$$\displaystyle+\frac{2\,\alpha\,(1-\alpha)\,(1-2\alpha)}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\,\left(\frac{{\underline{k}}\cdot{\underline{q}}}{{\underline{q}}^{2}}+\frac{{\underline{k}}^{\prime}\cdot{\underline{q}}^{\prime}}{{\underline{q}}^{\prime 2}}\right)$$

$$\displaystyle+\left.\left.\left.\frac{4\,\alpha^{2}\,(1-\alpha)^{2}}{\Big{[}{\underline{k}}^{2}+{\underline{q}}^{2}\alpha(1-\alpha)\Big{]}\Big{[}{\underline{k}}^{\prime 2}+{\underline{q}}^{\prime 2}\alpha(1-\alpha)\Big{]}}\right]\right\}\right|_{\lambda=\lambda^{\prime}=0}$$

$$\int \frac{d^{2}k}{(2\pi {)}^2}\frac{e^{-i\underline{k}\cdot \underline{z}}}{k^2+q^2}=\frac{1}{2\pi }{K}_0(qz)$$

$$\int \frac{d^{2}k}{(2\pi {)}^2}{e}^{-i\underline{k}\cdot \underline{z}}\frac{\underline{k}}{k^2+q^2}=\frac{-i}{2\pi }\frac{\underline{z}}{z}q{K}_1(qz)$$

Full text

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, 7, 8, 9, 10, 11, 12, 13] and Balitsky-Kovchegov (BK) [14, 15, 16, 17, 18] small-$x$ evolution equations in the Color Glass Condensate [19, 20, 21, 22, 23, 24, 25, 26, 6, 7, 8, 9, 10, 11, 12, 13, 17, 18, 14, 15, 16, 27, 28, 29] has been improved due to the completion of the calculations determining the scale of the running coupling in the evolution kernel in [1, 30, 2, 31]. The calculations in [1, 2] proceeded by including $\alpha_{s}\,N_{f}$ corrections into the evolution kernel and by then completing $N_{f}$ to the complete one-loop QCD beta-function by replacing $N_{f}\rightarrow-6\pi\beta_{2}$. Calculation of the $\alpha_{s}\,N_{f}$ corrections is particularly easy in the $s$-channel light-cone perturbation theory formalism [32, 33] used to derive the BK and JIMWLK equations: there $\alpha_{s}\,N_{f}$ 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 $\alpha_{s}\,N_{f}$ 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 $\alpha_{s}\,N_{f}$ corrections resummed to all orders and with $N_{f}$ completed to the QCD beta-function, $N_{f}\rightarrow-6\pi\beta_{2}$, with $\beta_{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 $\alpha_{s}\,N_{f}$ 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. (2.3), which, combined with Eq. (38) gives us the subtraction terms (2.3) and (2.3) 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 contribution 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 $\gamma\approx 0.85$, which is different from the fixed coupling anomalous dimension of $\gamma\approx 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, 35, 36, 37, 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 $\alpha_{s}\,\beta_{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

[TABLE]

where

[TABLE]

is the kernel of the evolution. Here transverse two-dimensional vectors ${\underline{x}}_{0}$ and ${\underline{x}}_{1}$ denote the transverse coordinates of the quark and the anti-quark in the parent dipole, while $\underline{z}$ is the position of the gluon produced in one step of evolution [39, 40, 41, 42]. We have introduced the notation $\underline{r}=\underline{x}_{0}-\underline{x}_{1}$, $\underline{r}_{1}=\underline{x}_{0}-\underline{z}$, $\underline{r}_{2}=\underline{z}-\underline{x}_{1}$ for the sizes of the parent and of the new (daughter) dipoles created by one step of the evolution. The notation $r\equiv|\underline{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-$N_{c}$ 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, $\underline{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 $\alpha_{s}\ln(1/x_{B})$, where $x_{B}$ is the fraction of momentum carried by the emitted gluon [39, 40, 41, 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

[TABLE]

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 $\alpha_{s}\,N_{f}$ contributions as recently derived in [1, 2] brings in new physical channels 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 discussion 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), $\sim S(\underline{x}_{0},\underline{z}_{1})S(\underline{z}_{2},\underline{x}_{1})$. 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

[TABLE]

where $\mathcal{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 pieces

[TABLE]

The first term, $\mathcal{R}$, which we will call the ’running coupling’ contribution, gathers all the higher order in $\alpha_{s}\,N_{f}$ 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, $\tilde{K}$, which includes all the terms setting the scale for the running coupling:

[TABLE]

The second term, $\mathcal{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 $\mathcal{F}$ into running coupling and subtraction contributions, although constrained by unitarity arguments, is not unique. Two different separation schemes have been proposed in [2, 1]. 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-$q\bar{q}$ 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, ${\underline{z}}_{2}$, or the antiquark, ${\underline{z}}_{1}$. The kernel for the running coupling functional, Eq. (6), obtained in this way is

[TABLE]

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, $z=\alpha z_{1}+(1-\alpha)z_{2}$, where $\alpha$ 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:

[TABLE]

where

[TABLE]

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

[TABLE]

This prescription is useful as a benchmark used to compare with previous numerical [3] and analytical works [43, 34, 35] 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]

[TABLE]

The momentum labels in the above equation are explained on the left hand side of Fig. 3. If ${\underline{k}}_{1}$ and ${\underline{k}}_{2}$ 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 ${\underline{q}}={\underline{k}}_{1}+{\underline{k}}_{2}$. The other transverse momentum we use is ${\underline{k}}={\underline{k}}_{1}(1-\alpha)-{\underline{k}}_{2}\alpha$, where $\alpha$ is the fraction the of gluon’s “plus” momentum carried by the quark, $\alpha\equiv k_{1+}/(k_{1+}+k_{2+})$. The prime over the transverse momentum denotes the momentum of the same particle in the complex conjugate amplitude. For instance ${\underline{q}}^{\prime}$ is the momentum of the $s$-channel gluon in the complex conjugate amplitude. Finally, ${\underline{z}}_{1}$ and ${\underline{z}}_{2}$ denote the transverse coordinates of the quark and the antiquark. In Eq. (2.2) we use ${\underline{z}}_{12}={\underline{z}}_{1}-{\underline{z}}_{2}$ (the transverse separation between the quark and the antiquark) and ${\underline{z}}=\alpha\,{\underline{z}}_{1}+(1-\alpha)\,{\underline{z}}_{2}$ (the transverse coordinate of the gluon).

To obtain the BK kernel from Eq. (2.2) 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

[TABLE]

Below we will label the JIMWLK kernel by calligraphic letter ${\cal 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

[TABLE]

with $\alpha_{\mu}$ the bare coupling.

As shown in Fig. 2, at the NLO level, the subtraction term introduced in Eq. (5), is then defined by

[TABLE]

where $\underline{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,

[TABLE]

as both choices lead to the same subtraction term ${\mathcal{S}}^{\text{Bal}}_{\text{NLO}}[S]$:

[TABLE]

In [2] the subtraction point was chosen to be the transverse coordinate of the gluon $\underline{z}$,

[TABLE]

This leads to the following subtraction term, which we denote ${\mathcal{S}}^{\text{KW}}_{\text{NLO}}[S]$:

[TABLE]

Indeed the complete kernel in Eq. (5) is independent of the choice of $\underline{w}$. However, since the subtraction term of Eq. (2.2) was neglected both in [1] and in [2], different choices of $\underline{w}$ led to different expressions for the remaining running coupling part ${\cal R}[S]$, i.e., to different answers as far as investigations in [1] and in [2] were concerned. Different choice of $\underline{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 $K_{1}^{\text{NLO}}(\underline{x}_{0},\underline{x}_{1};\underline{z}_{1},\underline{z}_{2})$ from Eq. (2.2) 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

[TABLE]

where {\cal K}_{\begin{picture}(10.0,8.0)\put(5.0,5.0){\circle{8.0}} \put(2.9,2.5){{\scriptsize 1}} \end{picture}} denotes the kernel with the running coupling corrections resummed to all orders. Just like in [2, 31], here we will use $\overline{\text{MS}}$ renormalization scheme. Inclusion of fermion bubble chains generated two denominators at the end of Eq. (2.2), which is its only difference from Eq. (2.2). Here

[TABLE]

Now we have to perform the transverse momentum integrals in Eq. (2.2). First we expand the denominators at the end of Eq. (2.2) into a power series and rewrite Eq. (2.2) as

[TABLE]

where we have defined $\mu^{2}=\mu_{\overline{{\text{MS}}}}^{2}\,e^{5/3}$ to make the expressions more compact.

Indeed we can not always expand the denominators of Eq. (2.2) into a geometric series employed in Eq. (2.2), 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. (2.2), which is an approximation not valid for all $q$ and $q^{\prime}$. Therefore, by keeping the geometric series in Eq. (2.2) 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 $Q^{2}$. We refer the interested reader to [30] for more details on this issue.

Using the following formulas

[TABLE]

and

[TABLE]

we can now perform the $k$- and $k^{\prime}$-integrals in Eq. (2.2). Integrating over the angles of $\underline{q}$ and ${\underline{q}}^{\prime}$ as well yields

[TABLE]

We have defined

[TABLE]

for brevity. Now the integrals over $q$ and $q^{\prime}$ can be carried out to give

[TABLE]

Unfortunately further simplification of the expression in Eq. (2.2) is impossible without approximations. 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

[TABLE]

and

[TABLE]

after some algebra we obtain

[TABLE]

In arriving at Eq. (2.2) we employed functions $R_{T}({\underline{x}})$ and $R_{L}({\underline{x}})$, which have dimensions of transverse coordinates and are defined by

[TABLE]

and

[TABLE]

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. (2.2) into physical running couplings finally yields

[TABLE]

with the physical running coupling in the $\overline{{\text{MS}}}$ scheme given by

[TABLE]

Eq. (2.2) 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

[TABLE]

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

[TABLE]

and

[TABLE]

One notices immediately that $\mathcal{R}^{\text{Bal}}\left[S\right]$ calculated in [1] is different from $\mathcal{R}^{\text{KW}}\left[S\right]$ calculated in [2] due to the difference in the kernels $\tilde{K}^{\text{Bal}}$ and $\tilde{K}^{\text{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), $\mathcal{R}\left[S\right]-\mathcal{S}\left[S\right]$, 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 contribution. It is given by

[TABLE]

with the resummed BK kernel

[TABLE]

The resummed JIMWLK kernel {\cal K}_{\begin{picture}(10.0,8.0)\put(5.0,5.0){\circle{8.0}} \put(2.9,2.5){{\scriptsize 1}} \end{picture}}({\underline{x}}_{m},{\underline{x}}_{n};{\underline{z}}_{1},{\underline{z}}_{2}) is given by Eq. (2.2), along with Eqs. (2.2) and (31) defining the scales of the running couplings. In the numerical solution below we will replace $N_{f}\rightarrow-6\pi\beta_{2}$ in its prefactor, obtaining

[TABLE]

This substitution is the same as for all other factors of $N_{f}$. 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 $N_{f}$ in front of Eq. (2.2) is definitely a part of the beta-function. Hence the replacement $N_{f}\rightarrow-6\pi\beta_{2}$ is justified even in the subtraction term. Once again, in the numerical solution below we will use Eq. (2.3) along with Eq. (38) in Eq. (37) to calculate the subtraction term ${\mathcal{S}}[S]$.

Substituting ${\underline{w}}={\underline{z}}_{1}$ (or, equivalently, ${\underline{w}}={\underline{z}}_{2}$) in Eq. (37) would yield the subtraction term

[TABLE]

which has to be subtracted from $\mathcal{R}^{\text{Bal}}\left[S\right]$ calculated in [1] and given by Eq. (35) to obtain the complete evolution equation resumming all orders of $\alpha_{s}\,N_{f}$ in the kernel.

Substituting ${\underline{w}}={\underline{z}}=\alpha\,{\underline{z}}_{1}+(1-\alpha)\,{\underline{z}}_{2}$ in Eq. (37) yields

[TABLE]

which has to be subtracted from $\mathcal{R}^{\text{KW}}\left[S\right]$ calculated in [2] and given in Eq. (36) again to obtain the complete evolution equation resumming all orders of $\alpha_{s}\,N_{f}$ 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, $\alpha_{s}\,N_{f}$) 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 $\Delta Y=0.1$. We discretize the variable $|\underline{r}|$ into 800 points equally separated in logarithmic space between $r_{min}=10^{-8}$ and $r_{max}=50$. Throughout this paper, the units of $r$ will be GeV*-1*, and those of $Q_{s}$ 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 [23, 22]:

[TABLE]

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

[TABLE]

with $\gamma=0.6$ and $\gamma=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 $N^{AN}\propto r^{2\gamma}$ corresponds to an anomalous dimension $1-\gamma$ of the unintegrated gluon distribution at large transverse momentum. (AN labels initial conditions with anomalous dimension.) Our choices $\gamma=0.6$ and $\gamma=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, $\gamma\approx 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}^{\prime}=1$ GeV at Y=0 in Eqs. (42) and (43) and put $\Lambda=0.2$ GeV. Although $Q_{s}^{\prime}$ 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

[TABLE]

with $\kappa=0.5$. We have checked that this choice of $\kappa$, 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<r_{fr}$, with $r_{fr}$ defined by $\alpha_{s}(1/r_{fr}^{2})=0.5$, the running coupling is given by the one loop expression

[TABLE]

with $N_{f}=3$ and $\Lambda=0.2$ GeV, whereas for larger sizes, $r>r_{fr}$, we “freeze” the coupling at a fixed value $\alpha_{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).

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 differences 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 [36, 34, 47] and numerical studies on the solutions of the BK equation at leading order [3, 48, 49, 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, $\tau=r\,Q_{s}(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, 52, 53, 54, 55, 56, 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 $\tau$ range studied in this work, including the saturation region, $\tau>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\mathchar 13350\relax 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, 48, 46] 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]

[TABLE]

with $a$, $b$ and $\gamma$ free parameters, within a fixed window below the saturation region, $\tau\in[10^{-5},0.1]$. Noticeably, at large enough rapidities the whole fitting window lays within the geometric scaling window proposed in [36]: $(\Lambda/Q_{s}(Y))<\tau<1$, where $\Lambda$ is some initial scale. The value of $\gamma$ extracted from the fits at rapidity $Y=40$ lays in between $\gamma\sim 0.8$ and $\gamma\sim 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 ($\gamma\sim 0.85$ at $Y=70$). This result for anomalous dimension is very far away from the value obtained in [3] for fixed coupling solutions ($\gamma\sim 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, 35, 36, 37, 38], $\gamma_{c}=\chi(\gamma_{c})/\chi^{\prime}(\gamma_{c})=0.6275$, where $\chi$ 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, 55, 56, 57] and have successfully confronted HERA 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 $\gamma\rightarrow\gamma(r,Y)$. The value of the anomalous dimension at $r=1/Q_{s}$ and/or for $Y\rightarrow\infty$ is fixed to be the BFKL saddle point, $\gamma_{c}\sim 0.63$ (the saddle point value considered in [57] is slightly different, $\gamma\sim 0.53$), while the value $\gamma=1$ is recovered in the limit $r\rightarrow\infty$ 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 $<Q^{2}<45$ GeV2, whereas charged hadron $p_{t}$ spectra in dAu collisions is well reproduced by [55, 57, 56] in the region $1$ GeV $<p_{t}<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}<\tau<1$ corresponds to values of momenta $\sim 10\,Q_{s}(Y)<p_{t}<10^{5}\,Q_{s}(Y)$ (always well above the saturation scale), with $Q_{s}(Y=40)\sim 500\div 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 $\Delta Y<4$ in both cases, while we study the solutions of the evolution at asymptotic rapidities, $Y\sim 40$. We have checked that shifting our fitting region to larger values of $\tau$ (smaller momentum) would bring the value of $\gamma$ 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 $\tau<1$. It is clear that the tail of the distribution falls off with decreasing $\tau$ much steeper for the running coupling solution than for the fixed coupling one. A fit to a pure power-law function, $f=a\,\tau^{2\gamma}$, in the region $\tau\in[10^{-6},10^{-2}]$ yields $\gamma\sim 0.85$ for the running coupling and $\gamma\sim 0.61$ for the fixed coupling solution. The differences between fixed and running coupling solutions at $\tau>1$ are evident from Fig. 5. This is a puzzling result that remains to be understood from purely analytical methods.

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. (2.3), and KW, Eq. (2.3), 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. (2.3), 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 ${\underline{z}}_{1}\leftrightarrow{\underline{z}}_{2}$ symmetry of the kernel, whereas in the KW case no such cancellation happens and the subtraction contribution, Eq. (2.3), also includes linear terms, which are dominant over the non-linear ones in the non-saturated domain where $N\ll 1$.

ii) The subtraction contribution $\mathcal{S}$ has the same sign as the running coupling contribution $\mathcal{R}$ in the whole $\tau$ 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, $\mathcal{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 $\mathcal{D}(r,Y)\equiv\mathcal{S}[N(r,Y)]/\mathcal{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 $\mathcal{D}\sim 0.4$ at small $\tau$ to $\mathcal{D}\sim 1$ at $\tau\sim 1$. The same ratio for the Balitsky scheme takes significantly smaller values: it goes from $\mathcal{D}\sim 0.1$ at small $\tau$ to $\mathcal{D}\sim 0.4$ for $\tau\sim 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 $\mathcal{D}\sim 0.15$ for the KW scheme and $\mathcal{D}\sim 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.

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 $\mathcal{F}=\mathcal{R}-\mathcal{S}$ calculated under the KW scheme (Eqs. (8) and (36) for the running coupling term, $\mathcal{R}$, and Eq. (2.3) for the subtraction term, $\mathcal{S}$) and under Balitsky’s scheme (Eqs. (7) and (35) for the running coupling term and Eq. (2.3) 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-$\tau$ region, where the two curves lay almost on top of each other. In the saturation region, $\tau\mathchar 13350\relax 1$, the agreement is slightly worse, although the differences between the values of $\mathcal{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. (2.3). 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].

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 $\alpha_{s}\,N_{f}$ resummation and by the $N_{f}\rightarrow-6\pi\beta_{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 $r_{1}$ and $r_{2}$, which are determined at each step of the evolution by the conditions $N(Y,r_{1})=10^{-9}$, and $N(Y,r_{2})=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<r_{1}$, and less than a $3\%$ for $r>r_{1}$.

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 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 $Q_{s}$ yielded by the KW scheme to $Q_{s}$ given by the complete $\alpha_{s}\,\beta_{2}$-evolution equation is a factor of $\sim 2.5$ for the MV initial condition and a factor of $\sim 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 $Q_{s}$ corresponding to the complete $\alpha_{s}\,\beta_{2}$-evolution is $\sim 1.25$ for the MV initial condition and $\sim 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 $\sim 0.025$ for Balitsky’s and $\sim 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 $\alpha_{s}\,\beta_{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 $Q_{s}(Y)\gg\Lambda$, such that the new (from the LO standpoint) momentum scale $\Lambda$ introduced by the running coupling can be safely neglected. Hence the dynamics is again characterized by a single momentum scale $Q_{s}(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 $\alpha_{s}\,N_{f}$. We reiterated the fact that the separation of the resulting kernel resumming all powers of $\alpha_{s}\,N_{f}$ 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\,Q_{s}(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/Q_{s}$ the anomalous dimension of the scaling function in the running coupling case becomes $\gamma\approx 0.85$ (see Fig. 6). This is different from the result of several analytical estimates [34, 35, 36, 37, 38], which expect the anomalous dimension not to change when running coupling corrections are included and to remain at its fixed coupling value of $\gamma\approx 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 ${\mathcal{S}}^{\text{Bal}}[S]$, which has to be subtracted from the result of [1], is systematically smaller than ${\mathcal{S}}^{\text{KW}}[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 ${\mathcal{S}}^{\text{Bal}}[S]$ and ${\mathcal{S}}^{\text{KW}}[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(\alpha_{s}^{2})$, while the running coupling part of the kernel is $o(\alpha_{s})$. Hence, if we suppose that the effective value of the coupling is given by its magnitude at the saturation scale $Q_{s}(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 $\alpha_{s}\,N_{f}$ 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 $\gamma\approx 0.85$, in disagreement with the analytical expectations of [34, 35, 36, 37, 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 precision 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.