Collinearly improved JIMWLK evolution in Langevin form

The high-energy evolution of Wilson line operators, which at leading order is described by the Balitsky-JIMWLK equations, receives large radiative corrections en-hanced by single and double collinear logarithms at next-to-leading order and beyond. We propose a method for resumming such logarithmic corrections to all orders, at the level of the Langevin formulation of the JIMWLK equation. The ensuing, collinearly-improved Langevin equation features generalized Wilson line operators, which depend not only upon rapidity (the logarithm of the longitudinal momentum), but also upon the transverse size of the color neutral projectile to which the Wilson lines belong. This additional scale dependence is built up during the evolution, via the condition that the successive emissions of soft gluons be ordered in time. The presence of this transverse scale in the Langevin equation furthermore allows for the resummation of the one-loop running coupling corrections.


Introduction
It is by now well established that the Wilson lines are the proper degrees of freedom for describing scattering in QCD at sufficiently high energies. For instance, in the collision between a dilute projectile and a dense target -such as in deep inelastic scattering at small Bjorken x, or in proton-nucleus collisions at the LHC -, the forward scattering amplitude can be computed by associating a Wilson line to each of the partons composing the projectile. Similarly, to compute the (partonic) cross-section for particle production in a dilute-dense collision, one must associate Wilson lines to all the partons that are produced in the final state -separately in the amplitude and in the complex-conjugate amplitude. A Wilson line is a path-ordered phase built with the color field of the target which describes the change in the wavefunction of a parton (notably, its color precession) due to its multiple scattering off the target, as computed within the eikonal approximation. To obtain the physical amplitudes or cross-sections, one must also average the (gauge-invariant) product of Wilson lines describing the scattering 'in a given event' (i.e. for a given configuration of the color fields in the target) over all such configurations. The effective theory for the Color Glass Condensate (CGC) provides an appropriate weight function -a functional probability distribution -to that aim [1][2][3].
This overall picture is subjected to quantum evolution. The most important quantum corrections in the high-energy regime of interest are those associated with the emission of soft gluons, i.e. gluons which carry only a small fraction x 1 of the longitudinal momentum of their parent parton. In the leading logarithmic approximation (LLA), in which one keeps only those perturbative corrections where each power of α s is accompanied by a 'rapidity' logarithm Y ≡ ln(1/x), the high-energy evolution of the Wilson lines is JHEP08(2016)083 references). Notice also that, in order to be truly useful, such an extension should also admit a stochastic formulation, say, as a functional Langevin equation for the Wilson lines, and thus lend itself well to numerical simulations.
There are two different ways to appreciate the potential difficulties with such an extension (see also section 5 below for a more detailed discussion). If one attempts to build a collinearly-improved version of the JIMWLK Hamiltonian, then one should take guidance from the local version of the BK equation constructed in refs. [53,54]. (We recall that the rapidity Y plays the role of an evolution time, hence an equation non-local in Y , like that presented in [52], cannot be derived via a Hamiltonian principle.) However, by inspection of the results in [53,54], it is quite clear that the resummed BK kernel does not have the proper structure to be obtainable from a Hamiltonian of the Fokker-Plank type (namely, it cannot be factorized as the product of two independent emissions; see section 5 for details). This strongly suggests that it should be impossible to construct a local (in Y ) Langevin equation with collinear improvement. It furthermore rises doubts about the possibility to maintain a probabilistic interpretation for the CGC effective theory beyond leading order.
If, on the other hand, one attempts to follow the strategy in [52] -that is, to construct a non-local Langevin equation which is endowed with time ordering -, then one immediately faces the following problem: the kinematical constraints expressing time-ordering involve the transverse sizes (or momenta) of the parent projectile and of the daughter gluons; yet, this transverse-scale dependence is not encoded in the usual definition of the Wilson lines operators.
In this paper, we shall identify a strategy to circumvent this last problem; that is, we shall propose a generalization of the Wilson line operator which 'knows' about the transverse size of the parent projectile. In turn, this will allow us to construct a non-local version of the Langevin equation where time-ordering is enforced via kinematical constraints. As we shall see, the additional scale dependence of the Wilson lines is naturally and automatically built up during the evolution, via the kinematical constraints. As a result, our Wilson lines will depend upon 3 independent variables: the transverse coordinate x = (x 1 , x 2 ) of the corresponding parton, the rapidity Y , and the transverse size r of the bare projectile which initiates the evolution. The dependence upon x already exists at tree-level, while the two other ones are introduced by the high-energy evolution with time-ordering. For a more complicated projectile, like a quadrupole, which involves several transverse scales -the transverse separations between the 'valence' partons -, our construction applies so long as these various scales are commensurable with each other. Their precise values are unimportant since the scale dependence which is ultimately built in the observable (the projectile S-matrix) is merely logarithmic. Indeed, the effect of the time-ordering within the Langevin equation is to resum to all orders the radiative corrections enhanced by double collinear logarithms.
This additional scale dependence of the Wilson line operators will also enable us to resum a particular class of single collinear logarithms together with the running coupling corrections. Indeed, the transverse size of the projectile sets the scale for the collinear logs and for the running of the coupling 2 (in the case where the daughter gluons have lower 2 Previous proposals for introducing a running coupling within the leading-order JIMWLK equation [14][15][16] were not fully satisfactory precisely because of this problem: the lack of information about the projectile transverse size.

JHEP08(2016)083
transverse momenta). Our strategy for these additional resummations is similar to that followed in [54] in relation with the BK equation.
Our main new results are presented in eqs. (4.4)-(4.6) and in eqs. (4.21)-(4.22). The first set of equations expresses the Langevin evolution with time-ordering and running coupling; the second set of equations encompasses the additional resummation of single collinear logarithms. It is this last set of equations that should be numerically solved for applications to the phenomenology. Clearly, the additional transverse-scale dependence complicates the numerical problem, by enlarging the associated functional space. We nevertheless believe that this problem is tractable.
The paper is organized as follows. In section 2, we succinctly review the Langevin formulation of the leading-order JIMWLK evolution, with emphasis on its physical interpretations as either target evolution, or as the evolution of the projectile. A good understanding of the physical picture is indeed important, as it will guide us when introducing time-ordering, later on. Then, in section 3, we recall the origin of the radiative corrections enhanced by double-collinear logarithms and, notably, their relation with the condition of time-ordering for successive gluon emissions in the 'hard-to-soft' evolution of the projectile. We also clarify a rather subtle point concerning the emergence of such corrections in the context of the JIMWLK equation. (This was originally devised as 'soft-to-hard' target evolution and one may think that it should be free of double collinear logs; this is however not the case for reasons to be explained at the end of section 3.) Section 4 is the main section, which contains our new results. Besides the two versions of the collinearly-improved Langevin equation already mentioned, cf. eqs. (4.4)-(4.6) and respectively eqs. (4.21)-(4.22), we also present the version of the BK equation emerging from our eqs. (4.4)-(4.6) and compare this version with the one proposed in [52]. Finally, in section 5 we summarize our results and extract from them some lessons for the JIMWLK evolution beyond leading order.

The Langevin formulation of the leading-order JIMWLK evolution
Our starting point is the Langevin formulation [11,12] of the leading order (LO) JIMWLK equation [5][6][7][8][9][10]. This involves a functional stochastic equation on the SU(N c ) group manifold which describes the high-energy evolution of the Wilson lines in the leading logarithmic approximation (LLA). The 'Wilson lines' are unitary matrices expressing the color precession of an elementary projectile (a quark or a gluon) which scatters off a (generally strong) color background field representing the gluon distribution of a dense target. It is understood that the partonic projectile belongs to a color-neutral system, such as a quark-antiquark dipole, which is built with several such partons and which remains dilute on all the resolution scales explored by the collision and by its high-energy evolution. The ultimate goal is to compute the elastic S-matrix for the scattering between that dilute projectile and the dense target, in the eikonal approximation.
The high-energy evolution proceeds via the additional emission of soft gluons (i.e. gluons which carry a small fraction of the longitudinal momentum of their parent parton) when increasing the rapidity separation Y = ln(s/Q 2 0 ) between the projectile and the target. Here, s is the center-of-mass energy squared and Q 0 is the characteristic transverse JHEP08(2016)083 scale of the target, e.g. its saturation momentum at the rapidity Y 0 = 0 at which one starts the evolution. The LLA is valid when α s Y 1 and consists in resumming the radiative corrections of the type (α s Y ) n to all orders.
As generally with stochastic equations, the Langevin formulation of the JIMWLK evolution requires a discretization of the 'evolution time', here the rapidity Y , so we shall write Y = n , with integer n ≥ 0. Also, we shall use the 'symmetric' version of the Langevin equation [16,26], where the change in a Wilson line associated with one evolution step (n → n + 1) involves both 'left' and 'right', infinitesimal, color precessions, whose physical meaning will be shortly explained. Then the Langevin equation reads (for a Wilson line in the fundamental representation, for definiteness): where x is the transverse coordinate of the quark projectile, which is not modified by the scattering in the eikonal approximation. The other notations are as follows: where the t a 's, with a = 1, . . . , N 2 c − 1, are the generators of the SU(N c ) Lie algebra in the fundamental representation, U n is a Wilson line in the adjoint representation which describes the color precession of the emitted gluon, and is the Weizsäcker-Williams kernel for the emission of a soft gluon (the propagator of the emitted gluon in the transverse plane; see also below). Finally, ν ia n,z is a Gaussian white noise, with correlator ν ia m,x ν jb n,y = 1 δ ij δ ab δ mn δ xy . (2.5) To better understand this evolution, let us first recall the generic structure of the Wilson line describing the eikonal scattering between a quark projectile and a dense target: We work in a Lorentz frame where the dilute projectile is a right mover, while the target is a left mover. Hence, the light-cone (LC) coordinate x + plays the role of a 'time' for the projectile and, respectively, a longitudinal coordinate for the target. The quark propagates along the positive LC, x 3 = x 0 , or x − = 0; its trajectory is parametrized by x + and the fixed transverse coordinate x. The target is a shockwave localized near x + = 0. The interaction consists in the eikonal coupling between the quark color current and the LC component increasing x + , matrices are ordered from right to left.) The 'background' field A − is frozen during a given collision, by Lorentz time dilation, but it can randomly vary from event to event and physical quantities are obtained after averaging over A − . This averaging must be performed at the level of the gauge-invariant product of Wilson lines which represents a color-neutral projectile, such as a dipole. Within the CGC effective theory [1][2][3], this average is computed using the 'CGC weight function', a functional probability distribution which evolves with Y according to the JIMWLK equation. Alternatively, and equivalently, the LO JIMWLK evolution of the CGC weight function can be reformulated as the stochastic equation (2.1) for Wilson lines with open color indices. As it should be clear from the above, we work in a gauge where the component A − a is non-zero. For the physical picture below, the most useful such gauge is the projectile LC gauge A + a = 0. The physical interpretation of eq. (2.1) can be given in terms of either target or projectile evolution. For more clarity, we shall briefly describe here both points of views. But the perspective of projectile evolution turns out to be more useful for our new developments in this paper.
When the rapidity increment dY = is used to boost the target, eq. (2.1) describes the evolution of the color field A − due to the emission of a soft gluon by color sources in the target. This new gluon carries a smaller fraction of the target longitudinal momentum P − , hence it is typically delocalized over larger values of |x + |. This is consistent with eq. (2.1), which shows that the field A − evolves by extending its support towards larger values of |x + |: the 'left' and 'right' color rotations in eq. (2.1) add new layers to A − , which are located at larger (and positive) and, respectively, smaller (and negative) values of x + as compared to the color field from the previous step. It furthermore evolves by developing new correlations, as introduced by the noise term. Physically, the noise term represents the color charge density of the emitted gluon. The adjoint Wilson line U n in eq. (2.3) describes the color precession of the noise (i.e. of the emitted gluon) by the color field built in the previous steps. Via this mean field effect, the correlations built by iterating eq. (2.1) are non-linear to all orders in the color field (or gluon) distribution in the target.
When eq. (2.1) is interpreted as the evolution of the projectile, the color field of the target is not evolving anymore -for any rapidity separation Y , this is still distributed according to the CGC weight function at the original rapidity Y 0 = 0 -but the S-matrix of the projectile changes due to soft gluon emissions within the projectile wavefunction. Accordingly, the unitary matrix U † n (x) should not be viewed anymore as a genuine Wilson line -if by a 'genuine Wilson line' we understand the path-ordered phase built with the field A − according to eq. (2.6) -, but rather as a more complicated scattering operator, which refers to a multi-partonic system built with the original quark at x together with the n soft gluons generated by the evolution. These gluons are strongly ordered in k + -the longitudinal momentum for a right-mover -which decreases from the projectile towards the target. This being said, we shall still refer to the unitary operator U † n (x) as a 'Wilson line', for brevity, while keeping in mind that, when n ≥ 1, this is not a bare Wilson line anymore.
Consider now an additional step in the evolution, n → n + 1, in which an even softer gluon is being emitted -either before the collision with the shockwave (the infinitesimal color precession on the right of U † n (x)), or after the collision (the color precession on the JHEP08(2016)083 left). If the emission occurs prior to the collision, then the soft gluon itself will eventually cross the shockwave and thus acquire a color precession described by U n (z). One may expect this additional gluon to be radiated by any of the n + 1 color sources generated in the previous steps, but this is not the content of eq. (2.1). Rather, this equation describes a process where the soft gluon is emitted by the quark at x alone. This is clear e.g. from the fact that eqs. (2.2)-(2.3) involve the Weizsäcker-Williams propagators from x (the emission point) to the point z where the gluon is measured by the noise term. 3 This reflects the fact that eq. (2.1) describes the backward evolution of the projectile.
At this point, it is useful to briefly remind the difference between the 'forward' and the 'backward' evolutions. In the forward evolution, the rapidity increment is viewed as a decrease in the rapidity of the softest gluons from the projectile that can be resolved by the target. This opens the phase-space for the emission of even softer gluons, from any of the preexisting color sources. In the backward evolution, on the other hand, the increment dY = is viewed as an increase in the rapidity of the original quark; this looks 'backwards', because the rapidity of the softest gluons is fixed. Due to this boost, the quark can emit an early gluon, which becomes the new 'first' gluon in the evolution -the one which is closest in rapidity to the original quark.
This 'backward' perspective of the high-energy evolution is very convenient in practice, since quite economical: at all the steps in the evolution, there is only source of radiation -the original quark. This viewpoint is underlying other related approaches, such as the dipole picture by Mueller [66] and the original constructions of the Balitsky hierarchy [4] and of the BK equation [19]. As we shall see, this viewpoint is also convenient for the inclusion of the collinear improvement and of the running coupling corrections within the Langevin approach to the JIMWLK evolution.
As explained in the Introduction, the main virtue of this Langevin formulation is to allow for explicit numerical solutions [13][14][15][16][17][18]. However, this formulation is also useful for conceptual studies, like constructing the Balitsky equations [4] for the evolution of products of Wilson line operators (the elastic amplitudes for dilute projectiles). For the purpose of deriving differential equations, it is sufficient to keep the terms that will matter up to order after performing the average over the noise. In this expansion, however, one should keep in mind that the noise itself scales like ν ∼ 1/ √ , as visible in eq. (2.5). Hence, the left and right infinitesimal rotations in eq. (2.1) must be expanded up to quadratic order in their exponents. Moreover, the quadratic terms ∼ ( ν) 2 in this expansion can be already averaged over the noise, because to the order of interest they cannot get multiplied by noise factors coming from other Wilson lines. After some simple algebra, one finds The term of O( ) in the first line has been obtained after averaging over the noise, meaning that the gluon is both emitted and reabsorbed by the original quark at x. This term contains two pieces: (i) a real term, involving the product of two Wilson lines, which describes a process where the soft gluon has been emitted before the collision and reabsorbed after it; (ii) a virtual term, involving the quark Wilson line alone, which corresponds to processes where the gluon fluctuation has no overlap in time with the scattering process. The term linear in the noise in the second line of eq. (2.7) is the exchange term, which allows the gluon emitted by the quark at x (either after, or prior to, the collision) to be reabsorbed by some other Wilson line. Eq. (2.7) makes it clear that at each step in the high-energy evolution, a new Wilson line is generated, representing a soft gluon with a generic impact parameter.
As an application of eq. (2.7), let us use it to derive the first equation in the Balitsky hierarchy [4], that for the S-matrix of a quark-antiquark color dipole. The dipole S-matrix is constructed as where Y = n , x and y are the transverse coordinates of the quark and the antiquark, and the average sign refers to both the noise average (over the noise terms ν 1 , . . . , ν n introduced by all the evolution steps) and the CGC average over the target color field A − at Y = 0 (the field in the initial Wilson line U † 0 (x)). Using eq. (2.7) together with a similar equation for U (n+1) (y), taking the color trace and performing the average over the noise ν n+1 associated with the last evolution step, one finds has been generated by summing over the 4 possible topologies for the emission and the absorption of the soft gluon, separately for 'real' and 'virtual' contributions (see figure 1).
For what follows, it is useful to notice that the dipole kernel can be factorized as the product of two independent emissions: That is, the soft gluon is first emitted, then reabsorbed, by either the quark or the antiquark. This factorized structure naturally emerges when deriving the Balitsky equations from the JIMWLK Hamiltonian, which describes two independent, soft gluon emissions from any of the preexisting partons.
JHEP08(2016)083 As expected, the 'real' term in the r.h.s. of eq. (2.9) is the S-matrix of a quarkantiquark-gluon system. By using Fierz identities, this can be rewritten as a system of 2 quark-antiquark dipoles plus a single dipole contribution of O 1/N 2 c , which precisely cancels against the respective contribution of the 'virtual' term in eq. (2.9) (thus replacing C F → N c /2 in the coefficient of the latter). After also performing the relevant averaging, one finds Y is the average S-matrix for a system of two dipoles, defined as

Time-ordering and collinear improvement
In the physical discussion in the previous section, we have successively exposed the viewpoints of target evolution and of projectile evolution, without making any distinction between their physical consequences: the LO Langevin equation (2.1) equivalently describes target evolution with decreasing k − , or projectile evolution with increasing k + . Strictly speaking, the corresponding longitudinal phase-spaces are not exactly the same, but their difference is irrelevant at LLA. However, this difference becomes important beyond LO and to study this we need to more carefully specify the kinematics. The rapidity phase-space for the target evolution is Y − ≡ ln(P − /p − 0 ), where p − 0 is the softest (in the sense of the lowest value for the longitudinal momentum k − ) gluon in the target wavefunction which matters for the scattering with the projectile. This value p − 0 is determined by the condition that the longitudinal wavelength ∆x + 1/p − 0 of this softest gluon be comparable with the lifetime ∆x + 2q + /Q 2 of the projectile. Here q + and Q 2 are the longitudinal momentum and respectively the virtuality of the projectile; e.g., for a dipole projectile, like (2.8), one has Q 2 ∼ 1/r 2 , with r = |x − y| the dipole transverse size. This condition yields 1

JHEP08(2016)083
where s = 2P − q + is the energy squared in the center-of-mass frame. Similarly, the rapidity phase-space available for the high-energy evolution of the projectile is Y + ≡ ln(q + /q + 0 ), with q + 0 determined by a condition analogous to (3.1): where we recall that Q 2 0 is the transverse resolution scale (generally, the saturation momentum) of the un-evolved target.
The difference Y + − Y − = ln(Q 2 /Q 2 0 ) ≡ ρ is positive in most interesting physical situations, since in dilute-dense collisions the projectile looks generally 'hard' (Q 2 Q 2 0 ) on the resolution scale of the un-evolved target. In the LLA, one implicitly assumes that this difference is negligible compared to the rapidity phase-space: When this happens, one can use either Y + or Y − as the 'evolution time' Y , and both target and projectile evolutions are correctly described (to LLA) by eq. (2.1), or by the LO B-JIMWLK equations. However, there are important physical situations where the rapidity shift δY = ρ cannot be neglected. First, ρ is not necessarily small compared to Y , not even in the high-energy kinematics; e.g., in the study of DIS at small Bjorken x Bj ≡ Q 2 /s, the virtuality phase-space can be large enough forᾱ s ρ ∼ O(1). Second, the perturbative corrections to the high-energy evolution which are introduced by the transverse phasespace are 'anomalously large': they start at 4 O ᾱ s ρ 2 and not at O(ᾱ s ρ), for reasons to be shortly reviewed. Accordingly, there is a region in phase-space where Y ρ, but such that bothᾱ s Y andᾱ s ρ 2 are of order one, or larger. In that region, the evolution is still driven by the rapidity, yet some of the 'higher-order corrections' enhanced are truly of O(1) and must be resummed to all orders, for consistency.
At this level, the dissymmetry between target and projectile plays an important role. In the collinear regime where ρ is positive and large, say such thatᾱ s ρ 2 ∼ 1, the projectile evolution in k + (or Y + ) receives double-logarithmic corrections ∼ (ᾱ s ρ 2 ) n to all orders, whereas the target evolution in k − (or Y − ) does not -the respective corrections involve only single collinear logarithms, of the form (ᾱ s ρ) n . This dissymmetry can be attributed to the fact that the double collinear logarithms are related to time-ordering in the highenergy evolution: the lifetime of the soft daughter gluons should remain smaller than that of their faster parents. When ρ is not too large (ᾱ s ρ 2 1), this condition is automatically satisfied, because of the strong ordering in longitudinal momenta. But in the collinear regime, this condition might be violated by the LLA evolution of the projectile in Y + , while it is still satisfied by the LLA evolution of the target in Y − .
Recall indeed that the meaning of 'time' is different for the target and respectively the projectile. For the left-moving target, the lifetime of a fluctuation is ∆x − 2k − /k 2 ⊥ . Successive emissions are strongly ordered in the longitudinal momentum k − , which decreases from the target towards the projectile. In the collinear regime, they are typically ordered in k ⊥ as well, but in the opposite direction: the transverse momentum increases when moving from the target towards the projectile ('soft-to-hard evolution'). This k ⊥ -ordering is favored JHEP08(2016)083 by the collinear phase-space: a sequence of n gluon emissions which are simultaneously or- brings in a large contribution, of order (ᾱ s Y ρ) n , which represents the dominant effect of the LLA in this particular regime. (The approximation to LLA which consists in keeping such contributions alone is known as the double-logarithmic approximation, or DLA.) Clearly, for all such typical configurations in the evolution of the target, the lifetimes are strongly ordered, 5 as anticipated -they decrease from the target towards the projectile: ∆x − 1 ∆x − 2 · · · ∆x − n . For the right-moving projectile, on the other hand, the lifetime of a soft gluon is estimated as ∆x + 2k + /k 2 ⊥ , where k ⊥ is strongly decreasing along the cascade ('hardto-soft evolution'), for the same reason as above: the strong bias from the collinear phasespace. 6 Since k + and k ⊥ are simultaneously decreasing, the lifetime of a daughter gluon can formally be larger than that of its parent parton. This occurs indeed, in large regions of the phase-space, when the projectile evolution in this collinear regime is computed according to LLA. But this situation is unphysical and does not occur when the relevant Feynman graphs are properly computed, beyond LLA: the phase-space regions where the time-ordering would be reversed are automatically suppressed, by energy denominators (see [53] for diagrammatic calculations which demonstrate this point).
The restriction to physical evolutions which respect the proper time-ordering has the effect to reduce the rapidity phase-space for the projectile evolution, roughly from Y + to Y + − ρ, and thus introduces corrections to the LLA evolution which are enhanced by double collinear logarithms -that is, corrections of the type (ᾱ s ρ 2 ) n [52,53]. To see this, consider the first step in the high-energy evolution of a dipole, namely the emission of a soft gluon with k + q + by either the quark, or the antiquark, leg of the dipole (recall figure 1). The lifetime condition implies which is a genuine constraint if and only if k 2 ⊥ Q 2 (since ∆Y > 0 in any case). This must be true, in particular, for the softest gluon in this evolution, that having k + = q + 0 and k 2 ⊥ ∼ Q 2 0 (cf. eq. (3.2)); this implies ln(q + /q + 0 ) > ln(Q 2 /Q 2 0 ), or Y ≡ Y + > ρ. (From now on, we shall denote Y + simply as Y , since we shall only study the projectile evolution.) It is intuitively clear that one step in this constrained evolution brings an effectᾱ s (Y − ρ)ρ, instead of the naive prediction (∼ᾱ s Y ρ) of the unconstrained LLA (at DLA accuracy). Two successive such steps will produce an effect of order This can be formally interpreted as the cumulated effect of performing 2 unconstrained steps in the DLA evolution (in k + ), plus an O ᾱ s ρ 2 correction to the emission kernel 5 This condition may be violated only for rare evolutions where the transverse momenta are inversely ordered: k ⊥2 k ⊥1 ; however, such evolutions are disfavored by the phase-space in the collinear regime and, moreover, they are also cut off by gluon saturation in the target. 6 This argument is not hindered by the saturation effects, which remain small so long as Q 2 k 2 ⊥ Q 2 0 .

JHEP08(2016)083
(the term linear in Y ), and a correction of O (ᾱ s ρ 2 ) 2 to the impact factor (the last term, independent of Y ). Vice-versa, one can properly resum the corrections of the type (ᾱ s ρ 2 ) n to all orders by simply enforcing the time-ordering constraint within the equations describing the LLA evolution. This has been already demonstrated at the level of the BK equation [52,53] and in what follows we shall extend this strategy to the Langevin formulation of the JIMWLK evolution and, hence, to the ensemble of the Balitsky equations at finite N c .
But before we proceed, there is an important question that we would like to clarify: namely, does the JIMWLK evolution require time-ordering in the first place ? That this is a legitimate question, can be seen as follows: the JIMWLK equation has been originally constructed as target evolution and we have just argued that, for the target evolution in k − , the proper time-ordering is automatically built-in at LLA, including in the collinear regime. The situation of the JIMWLK equation is however special: albeit this was indeed constructed as target evolution, the successive gluon emissions in this evolution were assumed to be strongly ordered in k + , and not in k − (the would-be natural evolution variable for a left-mover). This can be checked by inspection of the manipulations in 7 [8][9][10]. This special choice was convenient for the inclusion of gluon saturation: the soft gluons which are about to be emitted can multiply scatter off the strong background field generated in the previous steps of the evolution; this field is localized in x + , hence its modes carry large components k − , which will unavoidably alter the respective momentum of the newly emitted gluon. Hence, k − is not a good 'quantum number' for ordering the evolution gluons. The k + component is better suited in that sense, since this is not modified by rescattering (the background field is independent of x − , by Lorentz time dilation). These two components are related by the mass-shell condition (the evolution gluons are nearly on-shell), 2k + k − = k 2 ⊥ . Hence, in the absence of any strong ordering in transverse momenta, the target evolution with decreasing k − is equivalent to its evolution with increasing k + . However, in the collinear regime of interest here, the transverse momenta are strongly ordered (they increase simultaneously with k + ), so the time-ordering can in fact be violated. To avoid that, one needs to amend the JIMWLK evolution by explicitly enforcing time ordering. Another argument in that sense comes from the fact that the JIMWLK evolution generates the Balitsky hierarchy, for which the necessity of time-ordering is a priori clear from the perspective of projectile evolution.

JIMWLK evolution with collinear improvement
Let us consider one step in the 'backward' evolution of a dipole projectile and rewrite the time-ordering constraint (3.3) in terms of transverse sizes, rather than transverse momenta. We recall that Q 2 1/r 2 , with r = |x − y| the transverse size of the parent dipole. Also, via the uncertainty principle, one can relate the transverse momentum of the soft gluon to the transverse sizes of the daughter dipoles: 1/k ⊥ min(|x − z|, |z − y|). The difference in JHEP08(2016)083 size between the daughter dipoles is irrelevant in practice: the condition (3.3) is non-trivial only when the daughter dipoles are sufficiently large, such that |x − z| |z − y| r. Hence, for our purposes, one can replace eq. (3.3) with This constraint, together with the condition of strong ordering in longitudinal momenta, q + k + q + 0 , can be conveniently summarized as two upper limits: one on the transverse sizes |x − z| |z − y| of the (large) daughter dipoles, the other one on the corresponding rapidities. Specifically, the condition (This can also be viewed as a necessary condition for the existence of a soft emission with transverse size |x−z|, for a given rapidity separation Y between the 'valence' quark and the target.) Also, for a given transverse size |x − z| > r of the soft gluon fluctuation, the maximal value of its rapidity Y k = ln(k + /q + 0 ) that is allowed by time ordering 8 reads, clearly, where the Θ function is intended to remind that this reduction in the rapidity becomes effective only for sufficiently large daughter dipoles.

Scale-dependent Wilson lines and the time-ordered Langevin equation
It is rather straightforward to implement the kinematical constraints (4.2) and (4.3) at the level of the BFKL or BK equation: the condition (4.2) is inserted as a Θ-function multiplying the dipole kernel, whereas the rapidity shift (4.3) is used to modify the rapidity argument of the S-matrices for the daughter dipoles [52] (see also eq. (4.12) below). But when trying to apply a similar procedure to the Langevin equation (2.1), one immediately faces the following difficulty: this equation contains no information about the overall size r of the colorless projectile. (The transverse size of the soft gluon fluctuation is properly encoded, as the distance |x − z| between the gluon and the quark.) Yet, whenever using eq. (2.1) in practice, one has in mind a well-identified projectile, whose transverse size (at tree-level) is a priori known. So, there is no conceptual difficulty to further generalize 9 the notion of 'Wilson line' by associating to it a transverse size r which represents the transverse spread of the colorless projectile to which that Wilson line belongs. For instance, the two Wilson lines which enter the definition (2.8) of the color dipole are promoted to U † n (x, r) and U n (y, r), respectively, with r = |x − y|. The additional r-dependence is generated during the evolution, via the time-ordering constraints (see below). 8 In the absence of time ordering, one has Ymax = Y − dY , which is the same as Ymax = Y to the accuracy of interest. 9 Recall that, as we argued in a previous section, within the projectile evolution picture, the 'Wilson line' represents the scattering operator for the multi-partonic system generated by the high-energy evolution of a bare quark, or gluon.

JHEP08(2016)083
A general projectile can be more complicated than just a dipole -for instance, a quadrupole is built with 4 partons (say, 2 quarks and 2 antiquarks) already at tree-level, hence it can involve up to 6 different transverse sizes. In such a case, our present proposal for a generalized Wilson line applies only if all the intrinsic transverse separations are comparable to each other (say, in such a way that the logarithms ln(r 2 ij /r 2 kl ) are of order one for all the pairs (ij) and (kl) made with the valence partons). On the other hand, if the various transverse sizes are very different from each other (at tree-level, once again), then one cannot ignore the DGLAP-like evolution inside the wavefunction of the projectile, so the problem goes beyond the scope of the high-energy evolution alone.
The fact that the proper definition of an a priori local operator at high energy requires both longitudinal and transverse resolution scales is in line with what we know about the quantum evolution of the operators. These two scales control the phase-space for the evolution, namely they fix the maximal longitudinal momentum k + max = q + and the minimal size r (or maximal transverse momentum k 2 ⊥max = 1/r 2 = Q 2 ) of the quantum fluctuations permitted by our approximations.
The above considerations lead us to propose the following generalization of the Langevin equation (2.1), which incorporates the kinematical constraints in eqs. (4.2) and (4.3): where (compare to eqs. (2.2) and (2.3)) We have introduced the notations The noise correlator is the same as before, cf. eq. (2.5). The Θ-function inside the integrands in eqs. (4.5) and (4.6) expresses the upper limit (4.2) on the transverse size of the gluon fluctuations, due to time-ordering. Also, the reduction in the rapidity argument of the adjoint Wilson line, n → n − ∆ r xz , enforces the upper limit eq. (4.3) on the rapidity of the soft gluon. Because of this reduction, the evolution described by eqs. (4.4)-(4.6) is non-local in Y .
As anticipated, the dependence upon the transverse resolution scale r is generated by the evolution, via the kinematical constraints in the Langevin equation. (The initial condition at Y = 0 has no such dependence: U † 0 (x, r) = U † 0 (x). This can be selected, e.g., according to the McLerran-Venugopalan model [67].) Note also the way how this scale dependence gets updated when moving from the parent quark at x to the soft daughter gluon at z: the corresponding scale R r xz in the gluon Wilson line is the largest among r and |x − z| since this is the characteristic distance for the transverse spread of the color charge during the lifetime of the fluctuation.

JHEP08(2016)083
In writing eqs. (4.5) and (4.6), we moved the QCD coupling √ α s inside the integral over z, to anticipate the generalization to a running coupling, which is made possible too by the presence of the reference scale r. Specifically, the experience with running coupling effects in the context of the DGLAP evolution (see e.g. the textbook [68]) instructs us to re-interpret α s in these equations as More precisely, the r.h.s. of eq. (4.8) should be understood as the one-loop running coupling, At this point, we would like to comment on a previous proposal [16] for introducing a running coupling within the Langevin equation (2.1), but without any additional scale dependence in the Wilson lines. In that case, the running of the coupling was encoded via a suitable modification of the noise correlator (2.5). As a result, the QCD coupling was effectively evaluated at the scale k 2 = k 2 ⊥ , with k ⊥ the transverse momentum of the emitted gluon. Via the uncertainty principle, this scale can be related to the transverse separation between the parent quark and the daughter gluon: k ⊥ ∼ 1/|x − z|. This particular running coupling prescription is truly correct so long as the emitted gluon is sufficiently hard, such that k 2 ⊥ Q 2 , or |x − z| r. On the other hand, this prescription artificially enhances the relatively soft emissions with |x − z| r, by associating to them the large coupling 10 α s (|x − z|) instead of the physical one α s (r). This problem is remedied by our prescription in eq. (4.8), which looks natural in the context of the modified Langevin equation (4.4), but could not be implemented in the original equation (2.1), by lack of explicit information about the projectile size r.
The proposal in eqs. (4.4)-(4.8) represents our main new result in this paper. In what follows, we shall further motivate this proposal, in particular by showing that it leads to an improved version of the BK equation which properly resums the double-logarithmic corrections associated with time-ordering. Also, we shall extend this proposal to the resummation of an important class of single-logarithmic corrections -those which in perturbation theory start already at next-to-leading order. 10 In [16] it is argued that the running of αs with the transverse momentum k ⊥ of the emitted gluon should generate the correct, 'smallest dipole', prescription at the level of the BK equation (see the discussion around eq. (4.15) below). However, this cannot be true, since, whatever the kinematics, the gluon momentum k ⊥ is controlled by the transverse sizes, |x − z| and |y − z|, of the daughter dipoles (and is insensitive to the size r of the parent dipole). Specifically, the prescription in [16] amounts to selecting k ⊥ = |p + q|/2 as the argument of the coupling to be inserted in eq. (2.11) (see eq. (31) in [16]). The integrals over p and q in (2.11) are independent from each other and they are limited in the ultraviolet by the daughter dipole sizes, via the complex exponentials. (This argument strictly holds for a fixed coupling, but the insertion of a running coupling αs(k 2 ⊥ ) cannot change the convergence properties of these integrations.) In particular, for large daughter dipoles, |x − z| |y − z| r, one necessarily has k ⊥ ∼ 1/|x − z| Q = 1/r.

The associated BK equation
Let us first derive the generalization of the BK equation corresponding to eq. (4.4). To that aim, we proceed as explained in section 2, that is, we expand the 'left' and 'right' color rotations in eq. (4.4) up to quadratic order in the noise and then use (2.5) to average the quadratic terms. The analog of eq. (2.7) reads (with the shorthand notation Θ xz ≡ Θ(n − ρ r xz )) Using this together with the corresponding Langevin equation for the antiquark Wilson line U n (y, r) and performing the average over ν ia n+1 , one finds the following evolution equation for the dipole S-matrix (with Y = n and r = |x − y|): The averaging over the noise terms ν 1 , . . . , ν n from the previous evolution steps and over the initial Wilson line U 0 is kept implicit: at this stage, this averaging is not needed and eq. (4.11) can be as well understood at operatorial level. The argument of the running coupling has been suppressed to keep the notations simple, but can be easily restored from the accompanying Θ-function, i.e., √ α s Θ xz ≡ α s (min{|x − z|, r})Θ xz . Also, we have used Fierz identities to rewrite the various S-matrices in terms of fundamental Wilson lines alone. The first two terms in the r.h.s. of eq. (4.11) refer to 'real' emissions, that is, to gluons which cross the shockwave and can interact with it: the first term describes the case where the soft gluon is emitted prior to the collision by the quark, whereas the second term similarly refers to an emission by the antiquark. The two other terms in eq. (4.11) refer to 'virtual' emissions. The first one, involving the original dipole (x, y), has a familiar structure. But the second one involves new color structures: a color quadrupole, built with four Wilson lines, together with an unusual dipole whose both legs lie at z. This last term JHEP08(2016)083 describes the situation where the gluon is exchanged between the quark and the antiquark before the collision.
Eq. (4.11) looks considerably more complicated than the usual BK equation (2.12), but this is not a problem in practice, since we do not really need to solve this equationall numerical efforts should directly focus on the Langevin equation (4.4). Here, we shall use eq. (4.11) merely to illustrate the effects of the time-ordering. Besides, we shall argue that, to the accuracy of interest, this equation can be further simplified.
Namely, the complications apparent in eq. (4.11) are largely due to the differences between the kinematical constraints associated with the quark (x) and, respectively, the antiquark (y) -e.g., the fact that, in general, the Θ-functions Θ xz and Θ yz are different from each other. Such differences prevented us from reconstructing the dipole kernel according to eq. (2.10); they also explain the appearance of the new color structures in the last term in eq. (4.11). But in reality, these differences go beyond the accuracy of interest, as anticipated above eq. (4.1) and it will be further explained.
Notice first that the 'ultraviolet' structure of eq. (4.11) -i.e., the behavior of its integrand in the limit where the daughter gluon lies very close to its parent quark (|x−z| r) or antiquark (|z − y| r) -is not affected by time-ordering. Indeed, in any of these limits, the kinematical constraints become irrelevant and eq. (4.11) reduces to the original equation (2.12). (For instance, when |x − z| r, one has ρ r xz < 0 and ρ r zy 0, hence Θ xz = Θ yz = 1, ∆ r xz = ∆ r zy = 0, and R r xz = R r zy = r.) This was to be expected, since the time-ordering can only affect the gluon emissions with relatively soft transverse momenta. This is moreover important, as it ensures that the short-distance divergences of the emission kernels at |x − z| → 0 and |z − y| → 0 are harmless, since the 'real' and 'virtual' terms mutually cancel in these limits -as in the standard BK equation (2.12).
Consider now the opposite limit, where the daughter dipoles are relatively large, |x − z| |z − y| r. In that case, the time-ordering is clearly important, but the associated kinematical constraints look identical for the emissions by the quark and the antiquark, respectively. This, together with the unitarity of the Wilson lines, implies that the quadrupole piece which appears in the last term in eq. (4.11) reduces to the original dipole (x, y), whereas the dipole piece within the same term reduces to unity. Therefore, it becomes possible to combine the various terms in the r.h.s. of eq. (4.11) in such a way to reconstruct the dipole kernel according to eq. (2.10). Then the structure of the 'improved' equation becomes quite similar to that of the original equation (2.12), namely, where we identified the dipole S-matrices which enter the 'real' term according to Notice that the two legs of this dipole live at different rapidities and it is the softest leg (here, the antiquark piece of the gluon at z) which fixes the rapidity argument of the overall S-matrix. This is appropriate since the rapidity phase-space which remains available for the JHEP08(2016)083 evolution after the emission of the first gluon is not Y , but Y −∆ r xz . The 'hat' symbol on the S-matrix is intended to emphasize that this is an operator, like the Wilson lines themselves. In particular, there is no factorization assumption in eq. (4.12): this holds for generic N c .
Eq. (4.12) has been written for the situation where |x − z| |z − y| r and hence one can identify Θ xz = Θ yz and ∆ r xz = ∆ r zy = ln[(x − z) 2 /r 2 ]. But as a matter of facts, this equation can be extended to all physical regimes, as a replacement for the more complicated equation (4.11). To that aim, it suffices to replace, in eq. (4.12), Θ xz → Θ xz Θ zy = Θ(Y − ρ xyz ), ∆ r xz → ∆ xyz , and α s (r) → α s (r min ), where and r min ≡ min |x−y|, |x−z|, |y−z| . Indeed, after these replacements, eqs. (4.12) and (4.11) coincide with each other for both very small, and very large, daughter dipoles, whereas the differences between them which occur when the three dipoles are commensurable to each other (|x − z| ∼ |z − y| ∼ r) are irrelevant to the accuracy of interest (indeed, in that case, ρ r xz and ρ r zy are both of O(1) and hence negligible compared to Y ).
The generalization of eq. (4.12) that we have just described is equivalent (to the accuracy of interest, once again) to the collinearly-improved version of the BK equation proposed in [52] as a tool to resum the double-logarithmic corrections ∼ (ᾱ s ρ 2 ) n to all orders. Furthermore, as shown by the diagrammatic analysis in [53], this equation follows directly from the QCD Feynman graphs provided one uses a convenient organization of the perturbation theory -namely, the light-cone (or time-dependent, with 'time' = x + ) perturbation theory in the projectile light-cone gauge A + = 0.

Partially resumming single logarithms
In this subsection we shall present a further refinement of the Langevin equation, which refers to a partial resummation of the radiative corrections enhanced by single collinear logarithms -that is, the corrections of the type (ᾱ s ρ) n , which represent the interplay between the BFKL and the DGLAP evolutions. Namely, as shown in ref. [54], it is relatively easy to resum a particular subset of such corrections, which fully includes the respective piece at NLO (that is, the correction of orderᾱ s ρ to the BFKL kernel), together with a part of the higher order terms. Albeit incomplete, this resummation is still useful in practice, in that it allows one to keep under control all the NLO BFKL corrections which are amplified by (double or single) transverse logarithms. In other terms, the final version of the Langevin equation to be presented below is perturbatively correct up to pure O(ᾱ s ) corrections to the BFKL kernel.
Diagrammatically, the single-logarithmic corrections of interest correspond to Feynman graphs in which one 'BFKL emission' (the emission of a gluon with small longitudinal momentum fraction) is accompanied by an arbitrary number n ≥ 1 of 'DGLAP emissions' (quarks or gluons with longitudinal momentum fractions of O(1) and whose transverse momenta are strongly ordered). The transverse momenta can be either decreasing ('collinear JHEP08(2016)083 resummation'), or increasing ('anticollinear'), along the cascade. In the context of the linear, BFKL, evolution, one was able to fully resum such corrections, by working in a double Mellin representation where the DGLAP-like corrections exponentiate [46-48, 50, 51]. Here, we shall only resum a subset of these corrections which exponentiates already in transverse coordinate (or momentum) space. This partial resummation involves the first moment of a particular linear combination of the DGLAP splitting functions, namely, P gg (z) and P qg (z) are the gluon-to-gluon and, respectively, gluon-to-quark leading-order DGLAP splitting functions, in the conventions of [68] (see eqs. (2.98)-(2.99) in that textbook). With these conventions, P gg (z) has a singular piece 2N c /z at z → 0, which has been subtracted in eq. (4.16), as this would describe a BFKL emission. The particular linear combination of splitting functions which occurs in the l.h.s. of eq. (4.16) reflects the fact that the gluon and quark distribution functions are coupled under the DGLAP evolution.
For the present purposes, we have selected the eigenvalue of the 2 ×2 matrix-valued anomalous dimension 11 which controls the growth of the dipole amplitude with increasing energy (see e.g. [48] for details). Notice also that we have defined the 'anomalous dimension' A 1 to be positive definite. We shall first describe the resummation of the single-logarithmic corrections at the level of the BK equation. For a fixed coupling, the corresponding generalization of eq. (4.12) reads (as compared to eq. (4.12), we also extend the kinematical constraints as discussed around eq. (4.15)) with r 2 < ≡ min{(x−z) 2 , (y−z) 2 }. The factor within the square brackets in the r.h.s., whose exponent features A 1 , is the result of resumming the single logarithms: the positive sign in the exponent is taken when r < r < and the negative sign otherwise. Recalling that A 1 is positive, we see that this resummation is always slowing down the evolution. Eq. (4.17) differs from the collinearly-improved version of the BK equation proposed in [54] in that the effects of the time-ordering for the collinear emissions are encoded via kinematical constraints, and not via a resummation of the kernel. Yet, these two equations (eq. (4.17) above and eq. (9) in ref. [54]) are equivalent to the accuracy of interest.
The generalization of eq. (4.17) to a one-loop running coupling reads (cf. eq. (4.8)) 11 We recall that the DGLAP anomalous dimension Pij(ω) associated with a generic partonic splitting function j → i is computed as Pij(ω) =

JHEP08(2016)083
whereb = πb/N c and the positive (negative) sign in the exponent applies when r < r < (respectively, r > r < ). The running coupling version of the factor resumming the single logarithms has been obtained via the following argument (which is familiar in the context of the DGLAP evolution; see e.g. [50] for a similar discussion). Given two widely separated transverse momentum (or virtuality) scales Q 2 Q 2 0 , the resummation of the DGLAP logs [α s ln(Q 2 /Q 2 0 )] n which is relevant for our purposes at fixed coupling reads (compare to eq. (4.17)) With a one-loop running coupling, this is replaced by as shown in eq. (4.18). Within the above equation, ρ ≡ ln(Q 2 /Λ 2 ) and similarly ρ 0 ≡ ln(Q 2 0 /Λ 2 ).
We are now prepared to present the extension of the Langevin equation which resums single-logarithms as well: this has the same general structure as shown in eq. (4.4), but with the arguments of the infinitesimal, 'left' and 'right', rotations modified to (for the case of a running coupling) α R n+1 (x, r)= 1 π z √ α s Θ (n −ρ r xz ) ᾱ s (r) α s (|x−z|) where, as before, the factor of α s under the square root is understood as α s ≡ α s (min{|x − z|, r}); also, the positive sign in the exponent is taken when r < |x−z| and the negative sign when r > |x−z|.
It is quite obvious that the version of the BK equation which corresponds to this new Langevin equation can be immediately obtained by multiplying all the Weizsäcker-Williams kernels in eq. (4.11) by the new factor appearing in eqs. (4.21)-(4.22) (or the analog factor with x → y). The ensuing equation looks considerably more complicated than eq. (4.18) above, yet they are equivalent at the accuracy of the present approximations: in the collinear regime where the daughter dipoles are large, |x − z| |z − y| r, the additional factor associated with single logarithms becomes the same for emissions by the quark and respectively the antiquark, hence it appears as a factor multiplying the dipole kernel, such as in eq. (4.18). In the opposite, anticollinear, regime, where one of the daughter dipoles is very small, e.g. |x − z| r |z − y|, then both eq. (4.18) and eq. (4.11) are dominated by processes where the soft gluonis emitted and reabsorbed by the nearby quark at x; for instance, M xyz K xxz = 1/(x − z) 2 . For these processes, the Langevin equation obviously produces the same 'single-logarithm' factor as that visible in eq. (4.18), namely [ᾱ s (|x − z|)/ᾱ s (r)] A 1 /b . The mismatch between the respective factors for the other processes (which also involve the antiquark at y) is irrelevant, since those processes are anyway negligible in the anticollinear limit.

JHEP08(2016)083
used for numerical calculations. Yet, this is useful for formal studies, e.g. as an efficient tool for deriving the equations in the Balitsky hierarchy. A collinearly-improved version of this 'dipole' Hamiltonian can be easily deduced from the results in [53,54]. Namely, in those papers, one has constructed a colllinearly-improved BK equation, which is local in rapidity and where the resummation of the double collinear logarithms is performed directly at the level of the kernel. That is, the effects of the kinematical constraints appearing in eqs. (4.12) or (4.17) have been equivalently replaced by a change in the kernel, M xyz → M xyz K DLA , with K DLA resumming powers ofᾱ s ρ 2 to all orders (compare the present eq. (4.17) with eq. (9) in [54]). The correspondingly improved version of the 'dipole' JIMWLK Hamiltonian can be obtained via a similar replacement. Yet, this suffers from the drawback that we just mentioned: it does not lend itself to a Langevin formulation.
The 'strong' level refers to a Fokker-Planck-like Hamiltonian which has a factorized structure and thus allows for a Langevin formulation -such as the leading-order JIMWLK Hamiltonian (see the discussion in [12]). In our opinion, this factorized structure is not consistent with the collinear improvement, more precisely, with the constraint of timeordering. There are two ways to see that. On one hand, the improved Langevin equations that we have here obtained are non-local in rapidity and cannot be rewritten in a local form, via manipulations similar to those in refs. [53,54]. (Indeed, the strategy in [53,54] has crucially relied on the possibility to construct exact analytic solutions for the evolution equation with time ordering, in the double logarithmic approximation. This is clearly not the case for the time-ordered Langevin equations, due to their stochastic nature.) On the other hand, the collinearly-improved BK equation obtained in [53,54], which is local in Y , cannot be generated by a Hamiltonian of a Fokker-Plank type, because the new factor K DLA in the kernel cannot be factorized as the product of two emissions, in contrast to the LO dipole kernel (recall eq. (2.11)).
These considerations suggest that the probabilistic picture underlying the CGC effective theory [1][2][3], that is heavily relying on the Fokker-Planck nature of the evolution Hamiltonian, may not survive beyond leading order. This is not necessarily a surprise: a similar situation occurs in the more familiar context of the collinear factorization, where the probabilistic interpretation of the structure functions ceases to be valid beyond leading order, notably due to the scheme dependence of the higher order corrections. The fact that the largest radiative corrections to the high-energy evolution can nevertheless be accommodated at the level of 'collinearly-improved' Langevin equations demonstrates that a generalized stochastic picture can still apply beyond leading order, while at the same time providing a suitable framework for numerical calculations.