Resumming double logarithms in the QCD evolution of color dipoles

The higher-order perturbative corrections, beyond leading logarithmic accuracy, to the BFKL evolution in QCD at high energy are well known to suffer from a severe lack-of-convergence problem, due to radiative corrections enhanced by double collinear logarithms. Via an explicit calculation of Feynman graphs in light cone (time-ordered) perturbation theory, we show that the corrections enhanced by double logarithms (either energy-collinear, or double collinear) are associated with soft gluon emissions which are strictly ordered in lifetime. These corrections can be resummed to all orders by solving an evolution equation which is non-local in rapidity. This equation can be equivalently rewritten in local form, but with modified kernel and initial conditions, which resum double collinear logs to all orders. We extend this resummation to the next-to-leading order BFKL and BK equations. The first numerical studies of the collinearly-improved BK equation demonstrate the essential role of the resummation in both stabilizing and slowing down the evolution.


Introduction
It is by now well established that the Balitsky-JIMWLK hierarchy 1 [1][2][3][4][5][6][7] and its mean field approximation known as the Balitsky-Kovchegov (BK) equation [8] govern the high-energy evolution of scattering amplitudes in presence of non-linear effects (multiple scattering and gluon saturation) responsible for unitarization. Some of the most remarkable recent developments in that context refer to the first calculations of the next-to-leading order (NLO) corrections [9][10][11] to the B-JIMWLK and BK equations. These new developments parallel and extend previous efforts, towards the end of nineties, which established the NLO version [12][13][14][15][16][17] of the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [18][19][20] the linearized version of the BK equation which applies so long as the scattering is weak. Although the BFKL and B-JIMWLK equations are based on a common evolution mechanism, they differ in the way how they treat the scattering problem: the BFKL equation deals only with single scattering, as appropriate for a dilute target, whereas the B-JIMWLK hierarchy includes the interplay between evolution and multiple scatterings. The former is usually written in transverse momentum space, as an equation for the unintegrated gluon distribution, while the latter is formulated in terms of transverse coordinates (better suited for implementing the eikonal approximation) and keeps trace of the multiple scattering of the individual partons in the projectile -each of them represented by a Wilson line. Such differences explain the difficulty to adapt to the NLO B-JIMWLK evolution the 'collinear resummations' originally developed in the context of NLO BFKL [21][22][23][24][25], which aim at improving the convergence of the perturbative expansion for the BFKL kernel.
The collinear resummations refer to perturbative corrections, starting at NLO, which are enhanced by large, single or double, transverse logarithms. Without a proper resummation, which, strictly speaking, goes beyond the order-by-order expansion of the BFKL kernel, these large logarithms deprive the NLO BFKL formalism of its predictive power.
There is no reason to expect this lack-of-convergence problem to be attenuated by the non-linear terms in the B-JIMWLK equations: indeed, the 'collinear' corrections arise from regions in phase-space where the scattering is weak and the non-linear effects are negligible. This was anticipated in a semi-analytic study [26] and later on confirmed by the numerical observation that adding a unitarity constraint (in the form of a 'saturation boundary') to the NLO BFKL equation does not help improving the stability of the solution [27]. Very recently, while our work was being completed, this has been corroborated by a numerical study [28] of the NLO BK equation [9]: the numerical solution turns out to be unstable (the scattering amplitude decreases with increasing energy and can even turn negative) for the physically interesting initial conditions. As also shown in Ref. [28], this instability can be traced back to a large NLO correction to the BFKL kernel enhanced by a double transverse logarithm. This kind of correction, which can be associated with the choice of the reference scale in the energy logarithm, is well understood at BFKL level, where it is successfully resummed to all orders by the schemes proposed in Refs. [22][23][24][25]. It is our main objective in this paper to propose a similar resummation at the level of the BK equation.
More precisely, our goals are twofold: first, we would like to unambiguously identify the origin of the double-collinear logarithms in Feynman graphs to all orders and devise a method for their resummation; second, we would like to reformulate this resummation as a change in the kernel of the BK equation, which is energy independent. Unlike the corresponding method in the context of NLO BFKL, where the resummation is generally implemented in double Mellin space 2 [22,23], our resummation will be directly implemented in transverse coordinate space, in order to be consistent with the non-linear structure of the BK equation.
Concerning the first objective above, our main finding is that the double-collinear logs arise due to a reduction in the longitudinal phase-space for the high-energy evolution, as introduced by the condition that successive gluon emissions be strictly ordered in lifetime. The interplay between this 'kinematical constraint' and the double transverse logarithms has already been recognized in the literature [22,[29][30][31] (see [32] for a recent discussion and more references), but we are not aware of any systematic derivation of this prescription from Feynman graphs. To emphasize that this is indeed non-trivial, we notice that double collinear logs are also generated by diagrams with anti-time ordering, but they mutually cancel when all such graphs are summed together (see the discussion in Section 3 below). This observation helps understanding the peculiar way how the double transverse logs arise in the context of the NLO BK calculation in [9]. The main outcome of this diagrammatic analysis is Eq. (17), which governs the evolution in the double-logarithmic approximation (DLA): it resums to all orders the perturbative corrections in which each power of the coupling is accompanied by a double logarithm (either energy-collinear, or double collinear).
Eq. (17) however is non-local in 'rapidity' (the logarithm of the longitudinal momentum, which is our evolution variable), so it does not fully match our goals for a collinearly-improved evolution equation. 3 To cope with that, in Section 4 we demonstrate that the non-local Eq. (17) can be reformulated in a local form, modulo an analytic continuation and a reshuffling of the perturbative expansion. The new, local, Eq. (30) involves an 'improved' kernel and (for consistency) a modified initial condition, which both resum double-collinear logs to all orders.
It is then straightforward to extend this resummation to the BFKL and BK equations and thus obtain the collinearly-improved BK equation (32), which is our main result in this paper. It is 2 Note however some similarity between our strategy and that proposed in [25], where the ω-shift in Mellin space [22,23] has been approximately reformulated as an improvement of the BFKL kernel in transverse momentum space. 3 Collinearly-improved versions of the BK equation which are non-local in rapidity have been proposed too in the literature [32,33], but they suffer from some shortcomings, concerning either the systematics of the resummation (for the approach in [33]), or its feasibility in practice (for [32]). furthermore possible to promote this result to full NLL accuracy, by adding the remaining NLO BK corrections from Ref. [9]. Notice however that the NLO terms include single transverse logarithms, which may require additional resummations, as was already the case in the context of NLO BFKL [22][23][24]. Finally, in Section 5 we present the first numerical studies of the resummed BK equation (32). These studies clearly demonstrate the role of the resummation in both stabilizing and significantly slowing down the evolution: the saturation exponent extracted from the numerical solution is smaller by, roughly, a factor of two than in the absence of the resummation.

The double-logarithmic limit of the BFKL equation
In order to fix the notations and for comparison with the more refined results that we shall later obtain, it is instructive to recall the derivation of the 'naive DLA', by which we mean the version of this approximation which neglects the time-ordering of successive gluon emissions, from the leading-order (LO) BFKL equation [18][19][20]. The LO BFKL equation resums the perturbative corrections in which each power of the QCD coupling ᾱ s ≡ α s N c /π , assumed to be fixed and small, is accompanied by the energy logarithm Y ≡ ln(s/Q 2 0 ) (the 'rapidity'), with s the center-of-mass energy squared and Q 0 the characteristic transverse scale of the target. 4 In this leading-log approximation (LLA), valid when ᾱ s Y 1, it is consistent to treat the scattering and the evolution in the eikonal approximation. The LO BFKL equation can then be written as the linearized version of the BK equation [1,8], i.e. as an equation for the high-energy evolution of the scattering amplitude T x y (Y ) of a quark-antiquark dipole, with a quark leg at transverse coordinate x and an antiquark leg at transverse coordinate y, which undergoes weak scattering off a generic target (a nucleus, or a 'shockwave'): (1) This equation involves the 'dipole' version of the BFKL kernel, which describes the emission of a soft gluon with transverse coordinate z by either the quark or the antiquark leg of the dipole, followed by its reabsorption (see Fig. 1). In the limit of a large number of colors N c → ∞, the positive quantity (ᾱ s /2π )M x yz d 2 z can be interpreted [34] as the differential probability for the splitting of the original color dipole (x, y) into a pair of dipoles (x, z) and (z, y). The first two terms within the square brackets, T xz and T z y , are the 'real' terms describing the scattering of the daughter dipoles, whereas the last one, −T x y , is the 'virtual' term expressing the reduction in the probability for the parent dipole to survive at the time of scattering.
Eq. (1) is valid so long as the scattering is weak, T 1, for all the dipoles. For a dense target, such as a large nucleus, this is indeed the case provided all dipoles look small on the scale set by the target saturation momentum 5 Q 0 : (x − y) 2 Q 2 0 1, etc. In this regime, the integration over z in the r.h.s. of Eq. (1) becomes 4 This scale Q 0 is assumed to be hard enough for perturbation theory to apply: QCD . E.g., if the target is a small dipole, then 1/Q 0 is the size of that dipole. If the target is a large nucleus described by the McLerran-Venugopalan (MV) model [35], then Q 0 is the target saturation momentum in a frame where the projectile carries most of the total energy. 5 In the full BK equation, which includes unitarity corrections, this condition is is the saturation momentum of the target, which obeys Q s (0) = Q 0 and increases with Y . logarithmic when the daughter dipoles are much larger than the original one, i.e. for |x− z| |z − y| r ≡ |x − y|. For such configurations, the fast decreases of the dipole kernel, M x yz r 2 /(x − z) 4 , is partially compensated by the rapid increase of the scattering amplitudes for the daughter dipoles: T xz 2 . This last property, namely that the dipole scattering amplitude grows, roughly, like the area of the dipole in the transverse plane, is indeed satisfied by the initial condition at low energy Y 0 and is also preserved by the high energy evolution in the (doublelogarithmic) approximation discussed below.
The same arguments also imply that, in the region where the z-integration is logarithmic, the 'virtual' term T x y ∝ r 2 is much smaller than the 'real' ones and can be neglected. Its only effect will be to cut off the logarithmic phase-space at small dipole sizes To write down the double-logarithmic version of the BFKL equation, it is convenient to factor out the strong r 2 -dependence of the amplitude and write T Focusing on the logarithmically-enhanced contributions only, we can then average over azimuthal angles and impact parameters, in order to the (naive) double logarithmic accuracy, is most conveniently written in integral form, as where the upper limit 1/Q 2 0 in the integration over z 2 occurs because the proportionality law T (Y , r) ∝ r 2 holds only for r 1/Q 0 .
(E.g. for a large nucleus with saturation momentum Q 0 one has T (Y , r) 1 for r 1/Q 0 .) The integral equation above can be solved via iterations: As obvious from the previous discussion, the 'naive DLA' equation (3) resums terms of the type (ᾱ s Y ρ) n to all orders. This is a common limit of the BFKL and DGLAP evolution, which however is not very useful in practice, since valid only in a very narrow regime (namely, when ᾱ s Y ρ 1, but such that ᾱ s Y 1 and ᾱ s ρ 1). In the next section, we shall devise a more general double-logarithmic approximation, which resums all the perturbative corrections where ᾱ s is enhanced by exactly two logarithmseither Y ρ or ρ 2 .

From Feynman graphs to the DLA equation
To construct a more general DLA formalism, we shall perform a diagrammatic calculation, using light-cone (time-ordered) perturbation theory, of two successive steps in the high energy evolution of the dipole scattering amplitude. That is, we shall consider 2-gluon graphs in which one of the emitted gluons is much softer (in the sense of having a much smaller longitudinal momentum k + ) than the other one. We use conventions where the dipole projectile is an energetic right mover and we work in the light-cone gauge A + = 0. The use of 'old-fashioned' time-ordered perturbation theory is particularly convenient for our purposes, in that it allows for an economical and physically transparent classification and evaluation of the relevant Feynman graphs in the approximations of interest.
A rather compact way to organize the calculation has been devised in Ref. [36]: successive gluon emissions by the projectile, which are time-ordered and also strongly ordered in k + , are generated by repeatedly applying an evolution 'Hamiltonian' on the S-matrix for the high-energy scattering between the projectile and the target color field A − . This S-matrix is built from Wilson lines (one for each parton in the projectile), whose number keeps increasing in the course of the evolution, due to additional gluon emissions. For the original dipole, one has T x y = 1 − S x y with where t a is the SU(N c )-generator in the fundamental representation and P stands for path ordering: with increasing x + , color matrices are ordered from right to left. The Hamiltonian acts via functional differentiation w.r.t. A − and describes the emission of a soft gluon out of any of the preexisting Wilson lines, followed by its reabsorption (by either the same or a different Wilson line). After all the emissions have been produced by acting with the Hamiltonian, one must integrate over all the emission times and average over the background field A − to construct the physical scattering amplitude. This last step though (the functional average over A − ) is irrelevant for the present purposes and will be left unspecified in what follows. The first two steps in this evolution generate 2-gluon graphs like those shown in Fig. 2, where the gluon with longitudinal momentum p + is emitted first and it is harder than the gluon k + (p + > k + ). The topologies (in the sense of time-orderings) selected in Fig. 2 are quite special, in that they contribute already to LLA 6 : in Fig. 2a, the hard gluon (p + ) is real and is emitted before the soft one (k + ), but reabsorbed after it; in Fig. 2b, the hard gluon is virtual and it is both emitted and reabsorbed prior to the emission of the soft gluon, which is real. Beyond LLA, other time orderings become important as well and will be later considered (see Fig. 3).
We shall first evaluate the 2-real-gluon graph in Fig. 2a. After integrating over all emission times, within the ranges −∞ < t 1 < τ 1 < 0 and 0 < τ 2 < t 2 < ∞, one finds the following contribution In the integrals over k + and p + , the upper limit q + is the longitudinal momentum of the quark and antiquark in the original dipole, while the lower limit q + where τ p ≡ 2p + /p 2 = 1/p − is the lifetime of the hard gluon fluctuation, as determined by the uncertainty principle, and similarly for τ k . The integral over p + is logarithmic provided p + dominates both energy denominators, that is, so long as 8 Hence, to leading logarithmic accuracy for the longitudinal logarithm, one can replace τ p /(τ p + τ k ) (τ p − τ k ).
In the BFKL regime, one assumes that there is no strong hierarchy between the transverse momenta, |k| ∼ |p|, so the condition τ p > τ k is automatically satisfied when p + > k + . In that case, one can freely integrate over transverse momenta in expressions like 7 To keep expressions simple, we use the large-N c limit at intermediate steps, but some of the final results, notably the DLA equation (17), are valid for any N c . 8 For the purposes of power counting, one can use |k| ∼ |k| and |p| ∼ |p −k|; indeed, the difference between e.g. k and k is due to the scattering off the target, which is a comparatively small effect in the high transverse momenta (or small dipole sizes) regime of interest.
Eq. (6), to generate the Weizsäcker-Williams propagators of the soft gluons, according to After also summing over all possible connections for the two emitted gluons, one builds the relevant product of dipole kernels (i.e., M x yu M u yz for the sequence of emissions illustrated in Fig. 2).
However, this is strictly correct only so long as the transverse phase-space is by itself not logarithmic, meaning so long as Y ρ, where ρ ≡ ln(Q 2 /Q 2 0 ) measures the logarithmic separation in transverse scales between the original dipole, with size r ≡ 1/Q , and the target, with size 1/Q 0 . In the end, the transverse integrations in Eq. (6) are restricted to this range, e.g. Q 2 0 p 2 Q 2 (see below). For sufficiently large values of ρ, one opens the phasespace for a logarithmic integration over p 2 , which favors relatively large values |p| |k|. In this regime, the theta-function (τ p − τ k ) = (p + −k + (p 2 /k 2 )) becomes relevant and its effect is to reduce the longitudinal phase-space, roughly from Y to Y − ρ.
To the accuracy of interest, i.e. to correctly keep both the corrections of orders ᾱ s Y and ᾱ s ρ 2 generated when integrating out the hard gluon p + , the constraint τ p > τ k can be enforced directly in coordinate space, like p +ū2 > k +z2 . Here, we have anticipated that the corrections of the form ᾱ s ρ 2 come from emissions which are strongly ordered in transverse sizes, such that the daughter dipoles are much larger than the parent one. In this regime, |z − x| |z − y| |z − u| |u − x| |u − y| r = |x − y| , (9) and ū refers to any of the sizes, |u −x| or |u − y|, of the first pair of daughter dipoles, while z similarly refers to the daughter dipoles produced by the second splitting. After performing the momentum integrals in Eq. (6), summing over all the possible connections for both emitted gluons, and adding the other splitting sequence (where the gluon at z is emitted from the dipole (x, u)), one finds the following result from the 32 time-ordered graphs with two 'real' gluons (at large N c ): where ū = max (|u − x|, |u − y|) and z = max (|z − x|, |z − y|, |z − u|). Except for the theta-function enforcing time-ordering, this is recognized as the effect of two consecutive steps in the LO BFKL evolution.
To this result, one must add contributions coming from virtual graphs, evaluated to the same accuracy. The 'real-virtual' graphs in which the harder gluon (p + ) is virtual, whereas the softer one (k + ) is real, are the only ones that matter for the subsequent discussion of DLA. Consider first the 32 such graphs whose topologies (i.e. time-orderings) exist already at LLA, namely those where the two gluons have no overlap in time with each other (an example is shown in Fig. 2b). They give − ᾱ s In the BFKL context, this contribution is used to regulate the shortdistance singularities of Eq. (10) as u → x and u → y at a scale set by the original dipole size: ū r. In the present context, it plays a similar role (as anticipated in Eq. (9)), except for the fact that only the time-ordered piece of (11) is needed for that purpose. That is, albeit the virtual graphs included in Eq. (11) do not naturally involve any time ordering, it is nevertheless useful to distinguish between the respective time-ordered (TO) and anti-time-ordered (ATO) contributions, by inserting 1 = (τ p − τ k ) + (τ k − τ p ) in the integrand of Eq. (11). (Here and from now on, τ p = p +ū2 and τ k = k +z2 .) Then the TO piece must be combined with the 2-real-gluon contribution in Eq. (10), which is itself time-ordered, whereas the ATO piece is to be considered together with other virtual-real graphs, which are naturally ATO and will be discussed below.
From now on, we shall limit ourselves to the strict double- For subsequent discussions, it is important to stress that, to DLA, it is only the last emitted gluon (the one with the largest transverse size z) which contributes to scattering. Then the integrals over p + and ū are both logarithmic, as anticipated, and can be evaluated as where the logarithmic variables Y = ln(q + /k + ) and ρ = ln(z 2 /r 2 ) refer to the phase-space available to the hard gluon p + . Note that we have implicitly assumed above that Y > ρ, so that the integral over p + has indeed support for any ū ≥ r. This can be recognized as the condition for the lifetime τ k = k +z2 of the soft gluon fluctuation be (much) smaller than the 'lifetime' τ q = q + r 2 of the original dipole (the duration of the quantum process which has produced that dipole, e.g. the fluctuation of the virtual photon in DIS). To summarize, by integrating out the intermediate gluon p + , one has produced, besides the expected LLA contribution ᾱ s Y ρ, also a contribution ᾱ s ρ 2 , which can be interpreted as a NLO correction to the BFKL kernel for the emission of the soft gluon k + . This correction matches the respective piece (that enhanced by a double transverse logarithm) of the full NLO result in Ref. [9]. The last remark might suggest that the remaining 2-gluon graphs, that have not been considered so far and which correspond to other time orderings, do not contribute to order ᾱ s ρ 2 . But this is not quite true: contributions of this order arise from all the diagrams which are anti-time-ordered (ATO), in the sense that the lifetime of the hard gluon is shorter than that of the soft one (to DLA, at least). Topologically, the class includes two types of diagrams: (i) real-virtual graphs where the hard gluon is virtual and overlaps in time with the soft gluon which is real (some examples are the graphs 1a, 1b, 2a, 3a, 3b, 4a, and 4b in Fig. 3); (ii) real-real graphs where the hard gluon is emitted after, and absorbed before, the soft one (see graph 2b in Fig. 3). To these genuinely ATO diagrams, one must add the ATO pieces of the virtual-real graphs without overlap in time (see graphs 1c, 1d, 3c, and 3d in Fig. 3, which represent the ATO part of graphs like that in Fig. 2b, left over from the earlier calculations), to cancel UV divergences and introduce an effective short-distance cutoff equal to r (cf. the discussion after Eq. (11)).
When evaluating graphs of the type (i) and (ii) above mentioned, one finds that the time integrations over the overlapping region produce a factor like where the theta-function approximation in the r.h.s. holds in the double-logarithmic region. This theta-function cuts off the rapidity phase-space at the scale ρ (with ρ < Y ) and thus produces a contribution ∝ ρ 2 , as anticipated: It turns out however that all the terms of order ᾱ s ρ 2 generated by these ATO graphs exactly cancel each other. These cancellations can be understood as either the cancellation of 'infrared' (large z 2 ) logarithms between self-energy and vertex corrections, or as real vs. virtual cancellations for hard gluons whose scattering is not measured at DLA (cf. the discussion after Eq. (12)). For instance, the combinations of graphs 1a and 1b in Fig. 3, or 3a and 3b, belong to the first category, whereas 2a and 2b belong to the second. The sum of 1a and 1b leaves an uncompensated UV divergence, which is regulated at the scale r after also adding the ATO pieces of 1c and 1d (and similarly for 3a, 3b, 3c and 3d). Finally, graphs 4a and 4b mutually cancel because the color current associated with the hard gluon (and responsible for the emission of the soft gluon) has a different sign in 4a as compared to 4b. Interestingly, the order-ᾱ s ρ 2 contribution of the real-real graph 2b is such that it would exactly compensate the respective contribution of all the time-ordered real-real graphs previously discussed. This explains why, when the calculation is organized in such a way that the real-real graphs are all grouped together, as in Ref. [9] (which used the standard Feynman rules in momentum space), the net effect of order ᾱ s ρ 2 is rather seen to arise from the sum of the virtual-real graphs. The above cancellation pattern for the ATO graphs naturally generalizes to higher orders, i.e. to graphs involving an arbitrary number of strongly-ordered (in p + ) gluon emissions. To understand this, notice that graphs like those shown in Fig. 3 act as building blocks for generating via successive iterations all the diagrams that count at DLA. It is always possible to identify such subgraphs within the higher-order diagrams and then group them together in order to cancel their double-logarithmic contributions.
This leads us to the main conclusion in this section, namely the fact that the net contributions to DLA come fully from graphs which, within time-ordered perturbation theory, have exactly the same topology as the graphs contributing to LLA, but with the additional constraint that the successive gluon emissions must be strictly ordered in lifetimes. This implies that the perturbative corrections enhanced by the double logarithms Y ρ or ρ 2 can be resummed to all orders by solving a modified version of the DLA equation (3), which includes the time-ordering condition: (16) In what follows, we shall mostly use logarithmic variables, with the target scales q + 0 and Q 2 0 being the reference scales: where a step function (Y − ρ + ρ 1 ) is implicitly assumed within the integrand, to ensure that Y 1 > 0. This integral equation determines the function A(Y , ρ) for all positive values of Y and ρ, but the most interesting physical regime lies at Y > ρ.

Resummed kernel for DLA, BFKL, and BK evolutions
As compared to its 'naive' version in Eq. (3), the DLA equation with time ordering (17) is non-local in rapidity, as it can be best appreciated by rewriting it as a differential equation for the Y -evolution. This non-locality complicates the practical applications and, more importantly for our present purposes, it makes it quite tricky to extend this equation beyond DLA accuracy (albeit similar non-local versions of the BFKL and BK equations have been proposed in the literature; see Ref. [32] for a recent discussion). In what follows, we shall construct an alternative version of this equation, which is equivalent to (17) in the most interesting physical regime at Y > ρ and which is local in rapidity. The generalization of this local equation to BFKL and BK will then be almost straightforward, up to NLL accuracy.
At the mathematical level, it is more convenient to first solve the following problem: with f (0, ρ) = δ(ρ). (18) Given its solution f (Y , ρ), one can immediately construct the solution to Eq. (17) for arbitrary initial conditions according to (19) It is straightforward to solve Eq. (18) via iterations, to find where the -function arises from the longitudinal phase-space integration. For the purpose of the physical interpretation, it is useful to keep in mind that f (Y , ρ) is essentially the unintegrated gluon distribution in the dipole. The presence of the -function in the solution reflects the fact that, in order to emit a soft gluon, its lifetime τ k = k + /k 2 ⊥ must be smaller than the coherence time τ q = q + /Q 2 of the original dipole. Summing the above series we arrive at the explicit form where I 1 is the modified Bessel function. Neglecting for the moment the -function one can show that the function above admits an integral representation in the complex plane; namely, one has  (23) Here, the "characteristic function" is determined by where all the poles at γ = 1 visible in the r.h.s. are merely an artifact of expanding χ DLA (γ ) in series of ᾱ s . The resummed answer is clearly finite at γ = 1. The Jacobian J (γ ) induced by the change of variables is related to the characteristic function and reads The existence of a Mellin representation together with the exponentiation in Y (as manifest in the integrand in Eq. (23)) demonstrate that the function f (Y , ρ) obeys an evolution equation which is local in Y . Namely, Eq. (23) is tantamount to the following integral equatioñ with the kernel K DLA (ρ) defined as the inverse Mellin transform of χ DLA (γ ), that is, 9 in Eq. (21). The importance of this construction is that it can be immediately generalized to the evolution of the dipole amplitude, which can be thus rewritten as a local equation in Y . First, we define the analytic continuation of A(Y , ρ) according to (cf. Eq. (19)) This new function coincides with the physical amplitude A(Y , ρ) for Y > ρ. For with an initial condition Ã (0, ρ) which follows from Eqs. (29) and (28). For illustration, consider two interesting initial conditions, namely A(0, ρ) = 1, which has the advantage of simplicity, and A(0, ρ) = ρ, which is the limit of the McLerran-Venugopalan (MV) 9 Such a kernel (albeit written in momentum space) has already appeared in a work focusing on the collinear improvement of the NLO BFKL equation [25]. In that context, this kernel was obtained as an approximation to the equation for the BFKL eigenvalue with ω-shift [22]. model for dipole-nucleus scattering in the single scattering approximation [35]. One easily finds where we have temporarily used the notation ρ = 2 ᾱ s ρ 2 and where H α is the Struve function.
Eq. (30) is the sought-after local version of the DLA equation for the dipole amplitude: for Y > ρ, its solution coincides, by construction, with the respective physical amplitude, i.e. with the solution to the non-local equation (17). Notice that this rewriting of the DLA evolution in local form is tantamount to a complete reshuffling of the perturbation series: both the kernel in Eq. (30) and the initial condition in Eq. (31) resum double-collinear terms of the type (ᾱ s ρ 2 ) n for any n. For instance, the very first iteration of this equation generates all the terms linear in ᾱ s Y , i.e. the terms of the type ᾱ s Y ρ(ᾱ s ρ 2 ) n with n ≥ 0, that would be produced by iterating the original equation (17) to all orders. Remarkably, even though both the kernel and the initial condition exhibit oscillations as functions of ρ, their combined effect within equations like (30) or (26) yields a solution which is positive definite in the physical region Y > ρ, order by order in ᾱ s (e.g., this produces the perturbative solution (20) for f (Y , ρ)).
As we now explain, it is rather straightforward to promote this local DLA equation into a more complete equation, which includes the right BFKL and BK physics to NLL accuracy. To that aim, and starting with Eq. (30), we shall make backwards the steps leading from the LO BFKL equation (1) to the 'naive' DLA equation (3), that is: (i) We use the full expression for the dipole scattering amplitude, and more precisely its analytic continuation T (Y , ρ) ≡ e −ρÃ (Y , ρ) (which coincides with the physical amplitude for Y > ρ); (ii) We return to the use of transverse coordinates in our notations, meaning that we replace ρ = ln(1/r 2 Q 2 0 ), ρ − ρ 1 = ln(z 2 /r 2 ), T (Y , ρ) =T x y (Y ), and 2T (Y , z 2 ) →T xz (Y ) + T z y (Y ); (iii) We restore the full dipole kernel by replacing (r 2 /z 4 )dz 2 → M x yz d 2 z/π ; (iv) We reintroduce the virtual term and at the same time remove the infrared and ultraviolet cutoffs on the integral over z, since they are not needed anymore; (v) We replace the argument of the (additional) kernel K DLA , that is, ρ −ρ 1 = ln(z 2 /r 2 ), according to ln(z 2 /r 2 ) → L xzr L yzr , The last step above is the only one which is new as compared to the original discussion in Section 2 and will be thoroughly justified in a moment. We are thus led to the following equation 10 Some caution is required here, since there is a small region of integration where the product appearing under the square root can be negative; in this case, it is enough to let L xzr L yzr → |L xzr L yzr | and J 1 → I 1 (cf. Eq. (27) with ρ 2 < 0). So far, the initial condition at Y = 0 has been specified only in the weak scattering regime where T 1 and the precise normalization of T (Y = 0, r) was unessential (since the respective evolution was linear). For the purposes of the non-linear equation (32), however, we need to fix this normalization. To that aim, it is convenient to use the MV model [35], which amounts to exponen- K DLA (ρ) 1 −ᾱ s ρ 2 /2, it precisely matches the double-logarithmic term contained in the full NLO BK result [9]. The last feature makes it straightforward to formally extend Eq. (32) to full NLL accuracy: it is sufficient to add to its r.h.s. all the NLO BK corrections computed in [9], except for the double-log term that has already been included in the kernel. Let us finally remind that, strictly speaking, the solution to Eq. (32) can be trusted only for sufficiently small values of ρ Y (which in turn requires Y to be large enough, ᾱ s Y 1, in order to significantly evolve away from the 'unphysical' initial condition). In practice though, we expect the BFKL evolution encoded in Eq. (32) to eventually wash out the oscillations introduced by the initial condition at large ρ and thus progressively built a physical tail including at ρ > Y . This will be checked via numerical calculations in the next section.

Numerical tests
In this section we present a brief selection of first numerical studies which illustrate some subtle issues previously discussed (like the interplay between local and non-local evolution equations) and also some physical consequences of the resummation.  Ã (Y , ρ)). The respective solutions are supposed to coincide only at ρ < Y , where they both represent the actual physical result. This is indeed confirmed by the numerical simulations. On top of that, for ρ > Y , we see that the (physical) solution to Eq. (17) is independent of ρ and equal to A(Y , ρ) = cosh √ᾱ s Y [38], whereas its analytic continuation Ã (Y , ρ) shows (non-physical) oscillations which are inherited from the initial condition.
We now move to the collinearly-improved BK equation (32). Numerically, we solve the evolution equation following a strategy similar to the one described in the Appendix of Ref. [37]. By acting with the kernel in this equation on the power-like test function r 2γ , one can numerically extract the characteristic function ᾱ s χ (γ ) (the would-be Mellin transform of the resummed kernel). In the right-hand plot of Fig. 4, we compare the function ᾱ s χ (γ ) thus obtained for the particular value ᾱ s = 0.25 (black triangles) with the respective LO (BFKL) approximation ᾱ s χ 0 (γ ) (red squares) and with a 'NLO' approximation, 12ᾱ s χ NLO (γ ), obtained by keeping only the ᾱ s term in the expansion of K DLA , that is, This can also be understood by inspection of the DLA approximation χ DLA (γ ) in Eq. (24). At NLO accuracy, χ DLA (γ ) exhibits a cubic pole at γ = 1, the second term in the r.h.s. of Eq. (24), with a negative residue which makes the function χ NLO (γ ) unstable in the collinear limit γ → 1 (in particular, there is no saddle point on the real axis). By contrast, the all-order resummation ensures a smooth behavior near γ = 1, as already noticed after Eq. (24). For ᾱ s = 0.25, the function χ (γ ) is seen to be almost flat for γ 0.5. 11 In both cases, we found that a simple discretization of the integral and a Euler method to solve the rapidity evolution was sufficient to reach good numerical accuracy. 12 In this section, by 'NLO' we refer to the inclusion of the NLO corrections which are enhanced by double collinear logarithms, finite NLO corrections being neglected. A crude estimate of the saturation line 13 based on the DLA result in Eq. (21) yields [38] ρ s (Y ) ≡ ln which is significantly smaller than the respective LO result (no resummation) λ BFKL 4.88ᾱ s [35]. This suggests that the reduction of the longitudinal phase-space coming from time-ordering and giving rise to collinear double logs leads to a considerable reduction in the speed of the evolution. This expectation is indeed confirmed by the numerical solutions to Eq. (32). In Fig. 5, we show the results for ᾱ s = 0.25 and for an initial condition of the MV type, with A(0, ρ) = 1 (and hence Ã (0, ρ) as given in the first line of Eq. (31)). As before, the results with all-order resummation (cf. Fig. 5c) are compared to the respective predictions of LO BFKL (cf. Fig. 5a) and to the 'NLO' results obtained by using K NLO (ρ) = 1 −ᾱ s ρ 2 /2 (cf. Fig. 5b). The latter are highly unstable and physically meaningless -the evolution rapidly leads to a negative scattering amplitude -as it could have been anticipated in view of the pathological behavior of the corresponding characteristic function χ NLO (γ ) in Fig. 4. Similar instabilities have been recently observed [28] in numerical simulations of the full NLO BK equation and they have been traced back to the large double-logarithmic terms ∼ᾱ s ρ 2 in the NLO kernel, in agreement with our present findings. By contrast, the evolution with the fully 13 We recall the saturation line ρ s (Y ) is defined by the condition that T (Y , ρ) ∼ 1 when ρ = ρ s (Y ).
resummed kernel, shown in Fig. 5c, is perfectly smooth. We also see in Fig. 5c that the non-physical oscillations at ρ > Y introduced by resummation in the initial condition tend to disappear at larger rapidities. Finally, by comparing the LO results in Fig. 5a to the resummed ones in Fig. 5c, one clearly sees the anticipated reduction in the evolution speed. To more precisely characterize this reduction, we have numerically computed the target saturation momentum Q 2 s (Y ) for both the LO BFKL kernel and the fully resummed kernel, with results shown in Fig. 6 (for ᾱ s = 0.25 once again). Clearly, the growth of the saturation scale with Y is considerably reduced by the resummation: for sufficiently large Y , the saturation exponent λ s ≡ dρ s /dY approaches a value which is smaller by, roughly, a factor of 2 for the resummed kernel as compared to LO one. Remarkably, the asymptotic value which is thus obtained in the presence of resummation, namely λ s 0.55, agrees quite well with the respective DLA estimate in Eq. (33). We leave more detailed studies to a subsequent publication [38].